Fractional relaxation noises, motions and the fractional energy balance equation

We consider the statistical properties of solutions of the stochastic fractional relaxation equation and its fractionally integrated extensions that are models for the Earth's energy balance. In these equations, the highest-order derivative term is fractional, and it models the energy storage processes that are scaling over a wide range. When driven stochastically, the system is a fractional Langevin equation (FLE) that has been considered in the context of random walks where it yields highly nonstationary behaviour. An important difference with the usual applications is that we instead consider the stationary solutions of the Weyl fractional relaxation equations whose domain is −∞ to t rather than 0 to t.

An additional key difference is that, unlike the (usual) FLEs – where the highest-order term is of integer order and the fractional term represents a scaling damping – in the fractional relaxation equation, the fractional term is of the highest order. When its order is less than $\mathrm{1}/\mathrm{2}$ (this is the main empirically relevant range), the solutions are noises (generalized functions) whose high-frequency limits are fractional Gaussian noises (fGn). In order to yield physical processes, they must be smoothed, and this is conveniently done by considering their integrals. Whereas the basic processes are (stationary) fractional relaxation noises (fRn), their integrals are (nonstationary) fractional relaxation motions (fRm) that generalize both fractional Brownian motion (fBm) as well as Ornstein–Uhlenbeck processes.

Since these processes are Gaussian, their properties are determined by their second-order statistics; using Fourier and Laplace techniques, we analytically develop corresponding power series expansions for fRn and fRm and their fractionally integrated extensions needed to model energy storage processes. We show extensive analytic and numerical results on the autocorrelation functions, Haar fluctuations and spectra. We display sample realizations.

Finally, we discuss the predictability of these processes which – due to long memories – is a past value problem, not an initial value problem (that is used for example in highly skillful monthly and seasonal temperature forecasts). We develop an analytic formula for the fRn forecast skills and compare it to fGn skill. The large-scale white noise and fGn limits are attained in a slow power law manner so that when the temporal resolution of the series is small compared to the relaxation time (of the order of a few years on the Earth), fRn and its extensions can mimic a long memory process with a range of exponents wider than possible with fGn or fBm. We discuss the implications for monthly, seasonal, and annual forecasts of the Earth's temperature as well as for projecting the temperature to 2050 and 2100.

1 Introduction

Over the last decades, stochastic approaches have rapidly developed and have spread throughout the geosciences. From early beginnings in hydrology and turbulence, stochasticity has made inroads in many traditionally deterministic areas. This is notably illustrated by stochastic parameterizations of numerical weather prediction models, e.g. Buizza et al. (1999), and the “random” extensions of dynamical systems theory, e.g. Chekroun et al. (2010).

In parallel, pure stochastic approaches have developed primarily along two distinct lines. One is the classical (integer-ordered) stochastic differential equation approach based on the Itô or Stratonovich calculus that goes back to the 1950s (see the useful review by Dijkstra, 2013). The other is the scaling strand that encompasses both linear (monofractal, Mandelbrot, 1982) and nonlinear (multifractal) models (see the review by Lovejoy and Schertzer, 2013) that are based on phenomenological scaling models, notably cascade processes. These and other stochastic approaches have played important roles in nonlinear geoscience.

Up until now, the scaling and differential equation strands of stochasticity have had surprisingly little overlap. This is at least partly for technical reasons: integer-ordered stochastic differential equations have exponential Green functions that are incompatible with wide-range scaling. However, this shortcoming can – at least in principle – be easily overcome by introducing at least some derivatives of fractional order. Once the (typically) ad hoc restriction on integer orders is dropped, the Green functions are based on “generalized exponentials” that in turn are based on fractional powers (see the review by Podlubny, 1999). The integer-ordered stochastic equations that have received the most attention are thus the exceptional, non-scaling special cases. In physics they correspond to classical Langevin equations; in geophysics and climate modelling, they correspond to the linear inverse modelling (LIM) approach that goes back to Hasselmann (1976) and later elaborated notably by Penland and Magorian (1993), Penland (1996), Sardeshmukh et al. (2000), Sardeshmukh and Sura (2009) and Newman (2013). Although LIM is not the only stochastic approach to climate, in two recent representative multi-author collections (Palmer and Williams, 2010; Franzke and O'Kane, 2017), all 32 papers shared the integer-ordered assumption (a single exception being Watkins, 2017; see also Watkins et al., 2020).

Under the title “Fractal operators”, West et al. (2003) review and emphasize that, in order to yield scaling behaviours, it suffices that stochastic differential equations contain fractional derivatives. However, when it is the time derivatives of stochastic variables that are fractional – fractional Langevin equations (FLEs) – then the relevant processes are generally non-Markovian (Jumarie, 1993), so that there is no Fokker–Planck (FP) equation describing the corresponding probabilities. Even in the relatively few cases where the FLE has been studied, the fractional terms are generally models of viscous damping, so that the highest-order terms are still integer-ordered (an exception is Watkins et al., 2020, who mention “fractionally integrated FLE” of the type studied here but without investigating its properties). Integer-ordered terms have the convenient consequence of regularizing the solutions, so that they are at least root mean square continuous; in this paper the highest-order derivatives are fractional, so that when the highest-order terms are $\le \mathrm{1}/\mathrm{2}$, the solutions are “noises”, i.e. generalized functions that must be smoothed in order to represent physically meaningful quantities.

An additional obstacle is that – as with the simplest scaling stochastic model, fractional Brownian motion (fBm, Mandelbrot and Van Ness, 1968) – we expect that the solutions will not be semi-martingales and hence that the Itô calculus used for integer-ordered equations will not be applicable (see Biagini et al., 2008). This may explain the relative paucity of mathematical literature on stochastic fractional equations (see however Karczewska and Lizama, 2009). In statistical physics, starting with Mainardi and Pironi (1996), Metzler and Klafter (2000) and Lutz (2001) helped with numerics; the FLE (and a more general “Generalized Langevin Equation”, Kou and Sunney Xie, 2004; Watkins et al., 2019) has received a little more attention as a model for (nonstationary) particle diffusion (see West et al., 2003, for an introduction, or Vojta et al., 2019, for a more recent example). These technical aspects may explain why the statistics of the resulting processes are not available in the literature.

Technical difficulties may also explain the apparent paradox of continuous-time random walks (CTRWs) and other approaches to anomalous diffusion that involve fractional equations. While CTRW probabilities are governed by the deterministic fractional-ordered generalized fractional diffusion equation (e.g. Hilfer, 2000; Coffey et al., 2012), the walks themselves are based on specific particle jump models rather than (stochastic) Langevin equations. Alternatively, a (spatially) fractional-ordered Fokker–Planck equation may be derived from an integer-ordered but nonlinear Langevin equation for a diffusing particle driven by an (infinite-variance) Levy motion (Schertzer et al., 2001).

In nonlinear geoscience, it is all too common for mathematical models and techniques developed primarily for mathematical reasons to be subsequently applied to the real world. This approach – effectively starting with a solution and then looking for a problem – occasionally succeeds, yet historically the converse has generally proved more fruitful. The proposal that an understanding of the Earth's energy balance requires the fractional energy balance equation (FEBE, Lovejoy et al., 2021, announced in Lovejoy, 2019a) is an example of the latter. First, the scaling exponent of macroweather (monthly, seasonal, interannual) temperature stochastic variability was determined ($H_{\mathrm{I}}\approx -\mathrm{0.085}\pm \mathrm{0.02}$) and shown to permit skillful global temperature predictions (Lovejoy, 2015b; Lovejoy et al., 2015; Del Rio Amador and Lovejoy, 2019), and then it was extended to regional temperatures (at $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ resolution) (Del Rio Amador and Lovejoy, 2019, 2021a, b). The latter papers showed how the long-memory high-frequency approximation to the FEBE can not only make state-of-the-art multi-month temperature forecasts, but also how the corresponding simulations generate emergent properties such as realistic El Niño events.

In parallel, the multidecadal deterministic response to external (anthropogenic, deterministic) forcing was shown to also obey a scaling law but with a different exponent (Hébert, 2017; Lovejoy et al., 2017; Procyk et al., 2020, 2022; Procyk, 2021), $H_{\mathrm{F}}\approx -\mathrm{0.5}\pm \mathrm{0.2}$. It was only then realized that the order h FEBE naturally accounts for both the high- and low-frequency global temperature exponents with $h=H_{\mathrm{I}}+\mathrm{1}/\mathrm{2}$ and $H_{\mathrm{F}}=-h$, with both empirical exponents recovered with a FEBE of order $h\approx \mathrm{0.38}\pm \mathrm{0.03}$. The realization that the FEBE fit these basic empirical facts motivated the present research into its statistical properties, including its predictability.

In the EBE, energy storage is modelled by a uniform slab of material, implying that, when perturbed, the temperature exponentially relaxes to a new thermodynamic equilibrium. However, as reviewed in Lovejoy and Schertzer (2013), both conventional global circulation models and observations show that atmospheric, oceanic and surface (e.g. topographic) structures are spatially scaling. A consequence is that the temperature relaxes to equilibrium in a power law manner. This motivated earlier approaches (van Hateren, 2013; Rypdal, 2012; Hébert, 2017; Lovejoy et al., 2017) to postulate that the climate response function (CRF) itself is scaling. However, these models require either ad hoc truncations or imply infinite sensitivity to small perturbations (Rypdal, 2015; Hébert and Lovejoy, 2015).

The FEBE instead situates the scaling in the energy storage processes; this is the physical basis for the phenomenological derivation of the FEBE proposed in Lovejoy et al. (2021), and the zeroth-order term guarantees that equilibrium is reached after long enough times. The scaling of the basic physical quantities in both time and space motivates the study of the FEBE and its fractionally integrated extensions discussed below with temperature treated as a stochastic variable. The FEBE determines the Earth's global temperature when the energy storage processes are scaling and modelled by a fractional time-derivative term. Recently, analysis of the atmospheric radiation budget has shown that, at least over some regions, the internal component of the radiative forcing may itself be scaling: this justifies the consideration of the extensions to fGn forcing.

The FEBE differs from the classical energy balance equation (EBE) in several ways. Whereas the EBE is integer-ordered and describes the deterministic, exponential relaxation of the Earth's temperature to equilibrium, the FEBE is of fractional order, and because it is both deterministic and stochastic, it unites all the forcings and responses into a single model. Whereas the stochastic part represents the forcing and response to the unresolved degrees of freedom – the “internal variability” – and is treated as a zero mean Gaussian noise, the deterministic part represents the external (e.g. anthropogenic) forcing and the forced response modelled by the total external forcing. Complementary work (Procyk et al., 2020, 2022; Procyk, 2021) uses the deterministic FEBE as the basic model for the response to external forcing, but its Bayesian parameter estimation uses the stochastic FEBE to characterize the likelihood function of the residuals assumed to be the responses to stochastic internal forcing and governed by the same equation. It thus avoids the ad hoc error models involved in conventional Bayesian parameter estimation. The result is a parsimonious, FEBE projection of the Earth's temperature to 2100 that has much lower uncertainty than the classical global circulation model alternative. This is the first time that classical general circulation model climate projections have been confirmed by an independent, qualitatively different, approach.

An important but subtle EBE–FEBE difference is that, whereas the former is an initial value problem whose initial condition is the Earth's temperature at t=0, the FEBE is effectively a past value problem whose prediction skill improves with the amount of available past data, and – depending on the parameters – it can have an enormous memory (Del Rio Amador and Lovejoy, 2021b). To understand this, recall that an important aspect of fractional derivatives is that they are defined as convolutions over various domains. To date, the main one that has been applied to physical problems is the Riemann–Liouville (and the related Caputo) fractional derivative specialized to convolutions over the interval between an initial time=0 and a later time t. With one or two exceptions, this is the domain considered in Podlubny's mathematical monograph on deterministic fractional differential equations (Podlubny, 1999) as well as in the stochastic fractional physics discussed in West et al. (2003), Herrmann (2011), Atanackovic et al. (2014), and most of the papers in Hilfer (2000) (with the partial exceptions of Schiessel et al., 2000, and Nonnenmacher and Metzler, 2000). A key point of the FEBE is that it is instead based over semi-infinite domains – here from −∞ to t – often called Weyl fractional derivatives. This is the natural range to consider for the Earth's energy balance, and it is needed to obtain statistically stationary responses. Random walk problems involving fractional equations over the domain 0 to t can be dealt with using Laplace transform techniques. In comparison, the Earth's temperature balance involves statistically stationary stochastic forcings that are more conveniently dealt with using Fourier techniques.

We have mentioned that the FEBE can be derived phenomenologically where the fractional derivative of order h term represents the energy storage processes (Lovejoy et al., 2021). In this approach order h is an empirically determined parameter with h=1 corresponding to the classical (exponential) exception. Alternatively, it may derived from a more fundamental starting point, the classical heat equation – the same starting point as the classical Budyko–Sellers energy balance models (Budyko, 1969; Sellers, 1969). Recently it was shown with the help of Babenko's operator method that the special $h=\mathrm{1}/\mathrm{2}$ FEBE – the half-ordered energy balance equation (HEBE) – could be derived analytically from the classical heat equation (Lovejoy, 2021a, b).

To obtain the HEBE, it is sufficient to follow the Budyko–Sellers approach but to avoid one of their key approximations. The Earth's atmosphere and ocean are driven by local imbalances in radiative fluxes. While Budyko–Sellers models simply redirect this flux away from the Equator, the HEBE improvement (Lovejoy, 2021a, b) is to instead use the mathematically correct radiative–conductive surface boundary conditions. When this is done in the classical energy transport equation, one obtains an important $h=\mathrm{1}/\mathrm{2}$ special case of the FEBE, the half-order EBE or HEBE. The use of half-order derivatives in the heat equation is completely classical and goes back to at least Oldham (1973), Oldham and Spanier (1972), Babenko (1986), Magin et al. (2004), and Sierociuk et al. (2013). The extension to $h\ne \mathrm{1}/\mathrm{2}$ can be obtained using the same mathematical techniques by starting with the fractional generalization of the classical heat equation, the fractional heat equation. Further generalizations are also possible and will be reported elsewhere.

The choice of a Gaussian white noise forcing was made not so much for its theoretical simplicity as for its physical realism. Using scaling to divide atmospheric dynamics into dynamical ranges (Lovejoy, 2013, 2015a, 2019b), the main ones are weather, macroweather and climate. While the temperature variability in both space and time is generally highly intermittent (multifractal), there is one exception: the temporal macroweather regime (starting at the lifetime of planetary structures – roughly 10 d – up until the climate regime at much longer scales). Macroweather is the regime over which the FEBE applies, and it has exceptionally low intermittency: temporal (but not spatial) temperature anomalies are not far from Gaussian (Lovejoy, 2018). Responses to multifractal or Levy process FEBE forcings may however be of interest elsewhere.

This paper is structured as follows. In Sect. 2 we present the fractional relaxation equation, forced by a Gaussian white noise as a natural generalization of classical fractional Brownian motion, fractional Gaussian noise and Ornstein–Uhlenbeck processes (Sect. 2.1 and 2.2). When forced by Gaussian white noises, the solutions define the corresponding fractional relaxation motions (fRm) and fractional relaxation noises (fRn). We consider further extensions to the case where the equation is forced by a scaling noise fGn (Sect. 2.3, Eqs. 21 and 22). This is equivalent to considering the fractionally integrated fractional relaxation equation with white noise forcing. In Sect. 2, we first solve the equations in terms of Green's functions and then introduce powerful Fourier techniques that yield integral representations of the second-order statistics, including autocorrelations, structure functions (Eqs. 33 and 35), Haar fluctuations and spectra (with many details in Appendix A and in Appendix B, we derive the properties of the HEBE special case). In Sect. 3, we develop both short- and long-time (asymptotic) series expansions for the statistics (Eqs. 49 and 51), and we display and discuss sample fRn and fRm processes. In Sect. 4 we discuss the problem of prediction – important for macroweather forecasting – and derive expressions for the optimum predictor (Eq. 63) and its theoretical prediction skill as a function of forecast lead time (Eq. 68). In Sect. 5 we conclude.

We could note that the paper is somewhat complex due to the necessity of developing several approaches: Fourier for the main integral representations (Sect. 2), Laplace for the asymptotic expansions (Sect. 3), and real space for the predictability results (Sect. 4).

2 The fractional relaxation equation

2.1 fRn, fRm, fGn and fBm

In the introduction, we outlined physical arguments that the Earth's global energy balance could be well modelled by the fractional energy balance equation. Taking T as the globally averaged temperature, τ as the characteristic timescale for energy storage/relaxation processes, F as the (stochastic) forcing (energy flux; power per area), and s as the climate sensitivity (temperature increase per unit flux of forcing), the FEBE can be written in Langevin form as

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}^h\left({}_a{D}_t^{h}T\right)+T=sF,\end{array}$$

where the Riemann–Liouville fractional derivative symbol ${}_a{D}_t^h$ is defined as

$$\begin{array}{}\text{(2)}& {}_a{D}_t^{h}T={\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }(\mathrm{1}-h)}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\underset{a}{\overset{t}{\int }}(t-s{)}^{-h}T(s)\mathrm{d}s;\mathrm{0}<h<\mathrm{1},\end{array}$$

where Γ is the standard gamma function. Derivatives of order ν>1 can be obtained using $\mathit{\nu }=h+m$, where m is the integer part of ν, and then applying this formula to the mth ordinary derivative. The main case studied in applications (e.g. random walks) is a=0, so that Laplace transform techniques are often used (alternatively, the somewhat different Caputo fractional derivative is used). However, here we will be interested in $a=-\mathrm{\infty }$: the Weyl fractional derivative ${}_{-\mathrm{\infty }}{D}_t^h$, which is naturally handled by Fourier techniques (Sect. 2.4 and Appendices A and B), and, in this case, this distinction is unimportant.

Since Eq. (1) is linear, by taking ensemble averages, it can be decomposed into deterministic and random components with the former driven by the mean forcing external to system 〈F〉 and the latter by the fluctuating stochastic component F−〈F〉 representing the internal forcing driving the internal variability. The deterministic part has been used to project the Earth's temperature throughout the 21st century (Procyk et al., 2020, 2022); in the following we consider the simplest purely stochastic model in which 〈F〉=0 and F=γ, where γ is a Gaussian “delta-correlated” and with unit amplitude white noise:

$$\begin{array}{}\text{(3)}& 〈\mathit{\gamma }\left(v\right)〉=\mathrm{0};〈\mathit{\gamma }\left(v\right)\mathit{\gamma }\left(u\right)〉=\mathit{\delta }(u-v).\end{array}$$

In Hébert (2017), Lovejoy et al. (2017), and Hébert et al. (2021), it was argued on the basis of an empirical study of ocean–atmosphere coupling that τ_r≈2 years, while recent work indicates a value somewhat higher, ≈5 years (Procyk et al., 2022). At high frequencies, Lovejoy et al. (2015) and Del Rio Amador and Lovejoy (2019, 2021a) showed that the value h≈0.4 reproduced the Earth's temperature at scales <τ as well as for macroweather scales (longer than the weather regime scales of about 10 d) but still <τ. Procyk et al. (2020) also used the FEBE to estimate (the global) $s=[\mathrm{0.45},\mathrm{0.67}]$ $\mathrm{K}(\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{)}^{-\mathrm{1}}$ (90 % confidence interval), and the amplitude of the radiative forcing at monthly resolution was [0.89;1.42] W m⁻² (90 % confidence interval).

When $\mathrm{0}<h<\mathrm{1}$, Eq. (1) with γ(t) replaced by a deterministic forcing is a fractional generalization of the usual (h=1) relaxation equation; when $\mathrm{1}<h<\mathrm{2}$, it is the “fractional oscillation equation”, a generalization of the usual (h=2) oscillation equation (Podlubny, 1999).

To simplify the development, we use the relaxation time τ to nondimensionalize time, i.e. to replace time by $t/\mathit{\tau }$ to obtain the canonical Weyl fractional relaxation equation:

$$\begin{array}{}\text{(4)}& ({}_{-\mathrm{\infty }}{D}_t^h+\mathrm{1})U_h=\mathit{\gamma };Q_h\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}{U}_h\left(v\right)\mathrm{d}v,\end{array}$$

for the nondimensional process U_h. The dimensional solution of Eq. (1) with nondimensional γ=sF is simply $T\left(t\right)={\mathit{\tau }}^{-\mathrm{1}}$ $U_h(t/\mathit{\tau }$), so that in the nondimensional Eq. (4), the characteristic transition “relaxation” time between dominance by the high frequency (differential) and the low frequency (U_h term) is t=1. Although we give results for the full range $\mathrm{0}<h<\mathrm{2}$ – i.e. both the “relaxation” and “oscillation” ranges – for simplicity, we refer to the solution U_h(t) as “fractional relaxation noise” (fRn) and to Q_h(t) as “fractional relaxation motion” (fRm). Note that fRn is only strictly a noise when $h\le \mathrm{1}/\mathrm{2}$.

In dealing with fRn and fRm, we must be careful of various small and large t divergences. For example, Eqs. (1) and (4) are the fractional Langevin equations corresponding to generalizations of integer-ordered stochastic diffusion equations: the classical h=1 case is the Ohrenstein–Uhlenbeck process. Since γ(t) is a “generalized function” – a “noise” – it does not converge at a mathematical instant in time, and it is only strictly meaningful under an integral sign. Therefore, a standard form of Eq. (4) is obtained by integrating both sides by order h (i.e. by differentiating by −h and assuming that differentiation and integration of order h commute):

$$\begin{array}{ll}{\displaystyle }{U}_h\left(t\right)=& {\displaystyle }-{}_{-\mathrm{\infty }}{D}_t^{-h}{U}_h+{}_{-\mathrm{\infty }}{D}_t^{-h}\mathit{\gamma }\\ {\displaystyle }=& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }\left(h\right)}}\underset{-\mathrm{\infty }}{\overset{t}{\int }}(t-v{)}^{h-\mathrm{1}}{U}_h(v)\mathrm{d}v\\ \text{(5)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }\left(h\right)}}\underset{-\mathrm{\infty }}{\overset{t}{\int }}(t-v{)}^{h-\mathrm{1}}\mathit{\gamma }(v)\mathrm{d}v,\end{array}$$

(see e.g. Karczewska and Lizama, 2009). The white noise forcing in the above is statistically stationary; the solution for U_h(t) is also statistically stationary. It is tempting to obtain an equation for the motion Q_h(t) by integrating Eq. (4) from −∞ to t to obtain the fractional Langevin equation ${}_{-\mathrm{\infty }}{D}_t^h{Q}_h+Q_h=W$, where W is the Wiener process (a standard Brownian motion) satisfying dW=γ(t)dt. Unfortunately the Wiener process-integrated −∞ to t almost surely diverges, and hence we relate Q_h to U_h by an integral from 0 to t.

In the high-frequency limit, the derivative dominates, and we obtain the simpler fractional Langevin equation

$$\begin{array}{}\text{(6)}& {}_{-\mathrm{\infty }}{D}_t^h{F}_h=\mathit{\gamma };B_h\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}{F}_h\left(v\right)\mathrm{d}v,\end{array}$$

whose solution F_h is the fractional Gaussian noise process (fGn, not to be confused with the forcing) and whose integral B_h is fractional Brownian motion (fBm). We thus anticipate that F_h and B_h are the high-frequency limits of fRn and fRm.

2.2 Green's functions

Although it will turn out that Fourier techniques are very convenient for calculating the statistics, there are also advantages to classical (real-space) approaches, and in any case they are needed for studying the predictability properties (Sect. 4). We therefore start with a discussion of Green's functions that are the classical tools for solving inhomogeneous linear differential equations:

$$\begin{array}{ll}{\displaystyle }{F}_h\left(t\right)& {\displaystyle }=\underset{-\mathrm{\infty }}{\overset{t}{\int }}{G}_{\mathrm{0},h}^{\left(\text{fGn}\right)}(t-v)\mathit{\gamma }\left(v\right)\mathrm{d}v,\\ \text{(7)}& {\displaystyle }{U}_h\left(t\right)& {\displaystyle }=\underset{-\mathrm{\infty }}{\overset{t}{\int }}{G}_{\mathrm{0},h}^{\left(\text{fRn}\right)}(t-v)\mathit{\gamma }\left(v\right)\mathrm{d}v,\end{array}$$

where $G_{\mathrm{0},h}^{\left(\text{fGn}\right)}$ and $G_{\mathrm{0},h}^{\left(\text{fRn}\right)}$ are Green's functions for the differential operators corresponding respectively to ${}_{-\mathrm{\infty }}{D}_t^h$ and ${}_{-\mathrm{\infty }}{D}_t^h+\mathrm{1}$. Note that, due to causality, all Green's functions used in this paper vanish for t<0.

$G_{\mathrm{0},h}^{\left(\text{fGn}\right)}$ and $G_{\mathrm{0},h}^{\left(\text{fRn}\right)}$ are the usual “impulse” (Dirac) response Green's functions (hence the subscript “0”). For the differential operator Ξ, they satisfy

$$\begin{array}{}\text{(8)}& \mathrm{\Xi }{G}_{\mathrm{0},h}\left(t\right)=\mathit{\delta }\left(t\right).\end{array}$$

Integrating this equation, we find an equation for their integrals G_1,h, which are thus “step” (Heaviside, subscript “1”) response Green functions satisfying

$$\begin{array}{}\text{(9)}& \mathrm{\Xi }{G}_{\mathrm{1},h}\left(t\right)=\mathrm{\Theta }\left(t\right);\mathrm{\Theta }\left(t\right)=\underset{-\mathrm{\infty }}{\overset{t}{\int }}\mathit{\delta }\left(v\right)\mathrm{d}v;{\displaystyle \frac{\mathrm{d}{G}_{\mathrm{1},h}}{\mathrm{d}t}}=G_{\mathrm{0},h},\end{array}$$

where Θ is the Heaviside (step) function (=0 for t<0, =1 for t≥0). The inhomogeneous equation

$$\begin{array}{}\text{(10)}& \mathrm{\Xi }f\left(t\right)=F\left(t\right)\end{array}$$

has a solution in terms of either an impulse or a step Green function:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}f\left(t\right)& =\underset{-\mathrm{\infty }}{\overset{t}{\int }}{G}_{\mathrm{0},h}(t-v)F\left(v\right)\mathrm{d}v\\ & =\underset{-\mathrm{\infty }}{\overset{t}{\int }}{G}_{\mathrm{1},h}(t-v)F^{\prime }\left(v\right)\mathrm{d}v;F{}^{\prime }\left(v\right)={\displaystyle \frac{\mathrm{d}F}{\mathrm{d}v}},\end{array}\end{array}$$

the equivalence being established by integration by parts with the conditions $F(-\mathrm{\infty })=\mathrm{0}$ and $G_{\mathrm{1},h}\left(\mathrm{0}\right)=\mathrm{0}$. The use of the step rather than impulse response is standard in the energy balance equation literature since it gives direct information on energy balance and the approach to equilibrium (see e.g. Lovejoy et al., 2021). The step response for the noise is also the basic impulse response function for the motion.

For fGn, Green's functions are simply the kernels of the fractional integrals

$$\begin{array}{}\text{(12)}& F_h\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }\left(h\right)}}\underset{-\mathrm{\infty }}{\overset{t}{\int }}(t-v{)}^{h-\mathrm{1}}\mathit{\gamma }(v)\mathrm{d}v\end{array}$$

obtained by integrating both sides of Eq. (6) by order h. We conclude that

$$\begin{array}{}\text{(13)}& G_{\mathrm{0},h}^{\left(\text{fGn}\right)}={\displaystyle \frac{t^{h-\mathrm{1}}}{\mathrm{\Gamma }\left(h\right)}};G_{\mathrm{1},h}^{\left(\text{fGn}\right)}={\displaystyle \frac{t^h}{\mathrm{\Gamma }(h+\mathrm{1})}};-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\le h<{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

For fRn, we now recall some classical results useful in geophysical applications. First, these Green functions are often equivalently written in terms of Mittag–Leffler functions (“generalized exponentials”), E_α,β.

$$\begin{array}{}\text{(14)}& G_{\mathrm{0},h}\left(t\right)=t^{h-\mathrm{1}}{E}_{h,h}(-t^h);E_{\mathit{\alpha },\mathit{\beta }}\left(z\right)=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{\displaystyle \frac{z^n}{\mathrm{\Gamma }(\mathit{\alpha }n+\mathit{\beta })}}\end{array}$$

$$G_{\mathrm{0},h}\left(t\right)=\sum _{n=\mathrm{1}}^{\mathrm{\infty }}(-\mathrm{1}{)}^{n+\mathrm{1}}{\displaystyle \frac{t^{nh-\mathrm{1}}}{\mathrm{\Gamma }\left(nh\right)}};\mathrm{0}<h\le \mathrm{2}$$

To lighten the notation in Eq. 14 and in the following, we suppress the superscripts for fRn and fRm processes. A convenient feature of Mittag–Leffler functions is that they can easily be integrated by any positive order α:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{G}_{\mathit{\alpha },h}\left(t\right)={}_{\mathrm{0}}{D}_t^{-\mathit{\alpha }}\left(G_{\mathrm{0},h}\right(t\left)\right)=\\ {\displaystyle }& {\displaystyle }\begin{array}{cl}{t}^{h-\mathrm{1}+\mathit{\alpha }}{E}_{h,h+\mathit{\alpha }}(-t^h)=t^{\mathit{\alpha }-\mathrm{1}}\sum _{n=\mathrm{1}}^{\mathrm{\infty }}(-\mathrm{1}{)}^{n+\mathrm{1}}{\displaystyle \frac{t^{\text{nh}}}{\mathrm{\Gamma }(\mathit{\alpha }+nh)}};& t\ge \mathrm{0}\\ =\mathrm{0};& t<\mathrm{0}\end{array}\\ \text{(15)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }\ge \mathrm{0};\mathrm{0}\le h\le \mathrm{2}\end{array}$$

(Podlubny, 1999). As mentioned, the constraint t>0 is due to causality, and physical Green functions vanish for negative arguments. In the following this will simply be assumed. With α=1, we obtain the useful formulas

$$\begin{array}{}\text{(16)}& G_{\mathrm{1},h}\left(t\right)=t^h{E}_{h,h+\mathrm{1}}(-t^h),G_{\mathrm{1},h}\left(t\right)=\sum _{n=\mathrm{1}}^{\mathrm{\infty }}(-\mathrm{1}{)}^{n+\mathrm{1}}{\displaystyle \frac{t^{nh}}{\mathrm{\Gamma }(\mathrm{1}+nh)}}.\end{array}$$

With this, we see that $G_{\mathrm{0},h}^{\left(\text{fGn}\right)}$ and $G_{\mathrm{1},h}^{\left(\text{fGn}\right)}$ are simply the first terms in the power series expansions of the corresponding fRn and fRm Green functions. The solution to Eq. (4) with the white noise forcing γ(t) is therefore

$$\begin{array}{}\text{(17)}& U_{\mathrm{0},h}\left(t\right)=\underset{-\mathrm{\infty }}{\overset{t}{\int }}{G}_{\mathrm{0},h}(t-v)\mathit{\gamma }\left(v\right)\mathrm{d}v,\end{array}$$

where for this “pure” fRn process, we have added the subscript “0” for reasons discussed below. We note that, at the origin, for $\mathrm{0}<h<\mathrm{1}$, G_0,h is singular, whereas G_1,h is regular, so that it may be advantageous to use the latter (step) response function (for example in the numerical simulations in Sect. 4). These Green function responses are shown in Fig. 1. When $\mathrm{0}<h\le \mathrm{1}$, the step response is monotonic; in an energy balance model, this would correspond to relaxation to equilibrium. When $\mathrm{1}<h<\mathrm{2}$, we see that there are overshoot and oscillations around the long-term value; it is therefore (presumably) outside the physical range of an equilibrium process.

Figure 1 The impulse (a, b) and step response functions (c, d) for the fractional relaxation range (a, c: $\mathrm{0}<h<\mathrm{1}$, and red is h=1, the exponential), the black curves, bottom to top, are for $h=\mathrm{1}/\mathrm{10}$, $\mathrm{2}/\mathrm{10}$, … $\mathrm{9}/\mathrm{10}$ and the fractional oscillation range (b, d: $\mathrm{1}<h<\mathrm{2}$, red is the integer values h=1, c, d is the exponential, and a, b h=2), and the sine function, the black curves, bottom to top are for $h=\mathrm{11}/\mathrm{10}$, $\mathrm{12}/\mathrm{10}$, …, $\mathrm{19}/\mathrm{10}$.

In order to understand the relaxation process – i.e. the approach to the asymptotic value 1 in Fig. 1 for the step response G_1,h – we need the asymptotic expansion

$$\begin{array}{}\text{(18)}& G_{\mathit{\alpha },h}\left(t\right)=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{\displaystyle \frac{(-\mathrm{1}{)}^n}{\mathrm{\Gamma }(\mathit{\alpha }-nh)}}{t}^{\mathit{\alpha }-\mathrm{1}-nh};t\gg \mathrm{1}.\end{array}$$

For α=0, 1 we obtain the special cases corresponding to impulse and step responses:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& G_{\mathrm{0},h}\left(t\right)=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}(-\mathrm{1}{)}^n{\displaystyle \frac{t^{-\mathrm{1}-nh}}{\mathrm{\Gamma }(-nh)}},\\ & G_{\mathrm{1},h}\left(t\right)=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}(-\mathrm{1}{)}^n{\displaystyle \frac{t^{-nh}}{\mathrm{\Gamma }(\mathrm{1}-nh)}};t\gg \mathrm{1},\end{array}\end{array}$$

($\mathrm{0}<h<\mathrm{1}$, $\mathrm{1}<h<\mathrm{2}$; note that the n=0 terms are 0 and 1 for G_0,h and G_1,h respectively) (Podlubny, 1999), i.e. the asymptotic expansions are power laws in t^−h rather than t^h. According to this, the asymptotic approach to the step function response (bottom row in Fig. 1) is a slow, power law process. In the FEBE, this implies for example that the classical CO₂ doubling experiment would yield a power law rather than exponential approach to a new thermodynamic equilibrium. Comparing this to the EBE, i.e. the special case h=1, we have

$$\begin{array}{}\text{(20)}& G_{\mathrm{0},\mathrm{1}}\left(t\right)=e^{-t};G_{\mathrm{1},\mathrm{1}}\left(t\right)=\mathrm{1}-e^{-t},\end{array}$$

so that when h=1, the asymptotic step response is instead approached exponentially quickly. We see that when h=1, the process is a classical Ornstein–Uhlenbeck process, so that fRn can be considered a generalization of the latter. There are also analytic formulae for fRn when $h=\mathrm{1}/\mathrm{2}$ (the HEBE) is discussed in Appendix B, notably involving logarithmic corrections.

2.3 The α-order fractionally integrated fRn and fRm processes

Before proceeding to discuss the statistics of fRn and fRm processes, it is useful to make a generalization to the fractionally integrated processes:

$$\begin{array}{}\text{(21)}& U_{\mathit{\alpha },h}={}_{-\mathrm{\infty }}{D}_t^{-\mathit{\alpha }}{U}_{\mathrm{0},h}.\end{array}$$

U_α,h is the “α-order-integrated, fractional h relaxation noise”. Combined with Green's function relation $G_{\mathit{\alpha },h}={}_{-\mathrm{\infty }}{D}_t^{-\mathit{\alpha }}{G}_{\mathrm{0},h}$ (Eq. 15; recall that $G_{\mathrm{0},h}\left(t\right)=\mathrm{0}$ for t<0), we find that U_α,h and G_α,h are respectively the fractionally integrated relaxation noises and Green's functions of the fractionally integrated fractional relaxation equation:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& ({}_{-\mathrm{\infty }}{D}_t^{\mathit{\alpha }+h}+{}_{-\mathrm{\infty }}{D}_t^{\mathit{\alpha }})U_{\mathit{\alpha },h}=\mathit{\gamma },\\ & ({}_{-\mathrm{\infty }}{D}_t^{\mathit{\alpha }+h}+{}_{-\mathrm{\infty }}{D}_t^{\mathit{\alpha }})G_{\mathit{\alpha },h}=\mathit{\delta }\left(t\right).\end{array}\end{array}$$

If the highest-order derivative is constrained to be an integer (i.e. $\mathit{\alpha }+h=\mathrm{1}$ or 2), then the equation is a standard fractional Langevin equation; for example, U could be for the velocity of a particle with fractional damping and white noise forcing, although even here, the initial conditions are usually taken to be at t=0 and not $t=-\mathrm{\infty }$. Equivalently, U_α,h is the solution of the relaxation equation but with an fGn forcing:

$$\begin{array}{}\text{(23)}& ({}_{-\mathrm{\infty }}{D}_t^h+\mathrm{1})U_{\mathit{\alpha },h}={}_{-\mathrm{\infty }}{D}_t^{-\mathit{\alpha }}\mathit{\gamma }=F_{\mathit{\alpha }}\left(t\right);\mathrm{0}\le \mathit{\alpha }<\mathrm{1}/\mathrm{2}\end{array}$$

(the Weyl fractional derivatives commute). F_α is the α-order fGn process, and the restriction $\mathit{\alpha }<\mathrm{1}/\mathrm{2}$ is needed to ensure low-frequency convergence (see below).

In the Earth's radiative balance, such fractionally integrated fRn processes arise in two physically interesting situations. The first is where the forcing itself has a long memory – e.g. it is an fGn process. Whereas the memory in a pure fRn process is purely from the high-frequency storage term, in this case, the forcing (the overall radiative imbalance) also contributes to the memory, and this has important consequences for the predictability (Sect. 4). Although the solutions U_α,h are mathematically the same whether from the fractional relaxation equation with fGn forcing (Eq. 23) or the fractionally integrated fractional relaxation equation with white noise forcing (Eq. 22), only the former is directly relevant for the Earth's energy balance. This is because the energy balance involves the response from both stochastic (internal) and deterministic (external) forcing. For the latter, it is important that, following a step function forcing, at long times, the system will approach a new state of thermodynamic equilibrium. This implies that the term in the equation that dominates at low frequencies – the lowest-order term – is of order zero, so that if F in Eq. (1) is a step function, the new equilibrium temperature (anomaly) is T=sF.

The second situation where fractionally integrated fRn processes arise is for the energy storage (even in the purely white noise forcing case). The storage process is the difference between the forcing and the response:

$$\begin{array}{}\text{(24)}& S_{\mathit{\alpha },h}=F_{\mathit{\alpha }}-U_{\mathit{\alpha },h},\end{array}$$

so that

$$\begin{array}{}\text{(25)}& S_{\mathit{\alpha },h}={}_{-\mathrm{\infty }}{D}_t^h{U}_{\mathit{\alpha },h}=U_{h-\mathit{\alpha },h}.\end{array}$$

Even when the forcing is pure white noise (α=0), the storage is an h-ordered fractionally integrated process: $S_{\mathrm{0},h}=U_{h,h}$; this corresponds to the storage following an impulse forcing. The storage following a step forcing is obtained by integration order 1: $U_{\mathrm{1}+h,h}$. Similarly, Green's function for the fRn storage following an impulse forcing is G_h,h and, following a step forcing, $G_{\mathrm{1}+h,h}$ (Fig. 2). Since it turns out that most of the pure fRn (α=0) results are readily generalized to $\mathrm{0}<\mathit{\alpha }<\mathrm{1}/\mathrm{2}$, many fractionally integrated results are given below.

Figure 2 The storage Green functions for the fractional relaxation equation (α=0): (a) impulse response (G_h,h); (b) step response ($G_{\mathrm{1}+h,h}$). Black is for $h=\mathrm{1}/\mathrm{10}$, $\mathrm{2}/\mathrm{10}$, … $\mathrm{10}/\mathrm{10}$ and red for $\mathrm{11}/\mathrm{10}$, $\mathrm{12}/\mathrm{10}$, … $\mathrm{19}/\mathrm{10}$ (to identify the curves; use the fact that at large t, they are in order of increasing h – bottom to top). For small t, $G_{h,h}\propto t^{\mathrm{2}h-\mathrm{1}}$ (Eq. 15), so that for $h\le \mathrm{1}/\mathrm{2}$, the impulse response is singular at the origin. For large t, $G_{h,h}\propto t^{h-\mathrm{1}}$ (Eq. 18), so that for h<1, the total impulse response storage decreases following the impulse; for h=1 (the EBE), it tends to unity, and for h>1, it diverges.

2.4 Statistics

In the above, we discussed fGn, fRn and their order 1 integrals fBm, fRm as well as fractional generalizations, presenting a classical (real-space) approach stressing the links with fGn and fBm. We now turn to their statistics. U_α,h(t) is a mean zero stationary Gaussian process (i.e. $\langle U_{\mathit{\alpha },h}\left(t\right)\rangle =\mathrm{0}$, where “〈.〉” indicates ensemble or statistical averaging); therefore, its statistics are determined completely by its autocorrelation function R_α,h(t), which is only a function of the lag t:

$$\begin{array}{ll}{\displaystyle }{R}_{\mathit{\alpha },h}\left(t\right)& {\displaystyle }=〈U_{\mathit{\alpha },h}(t+v)U_{\mathit{\alpha },h}\left(v\right)〉\\ \text{(26)}& {\displaystyle }& {\displaystyle }=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{G}_{\mathit{\alpha },h}(t+v)G_{\mathit{\alpha },h}\left(v\right)\mathrm{d}v.\end{array}$$

The far-right equality follows from $U_{\mathit{\alpha },h}=G_{\mathit{\alpha },h}\ast \mathit{\gamma }$ and $〈\mathit{\gamma }\left(t\right)\mathit{\gamma }\left(t^{\prime }\right)〉=\mathit{\delta }(t-t^{\prime })$ (“*” indicates “convolution”). The process can only be normalized by R_α,h(0) when there is no small-scale divergence, i.e. when

$$\begin{array}{}\text{(27)}& R_{\mathit{\alpha },h}\left(\mathrm{0}\right)=\langle U_{\mathit{\alpha },h}^{\mathrm{2}}\rangle =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{G}_{\mathit{\alpha },h}(v{)}^{\mathrm{2}}\mathrm{d}v<\mathrm{\infty };\mathit{\alpha }+h>\mathrm{1}/\mathrm{2}.\end{array}$$

When $\mathit{\alpha }+h\le \mathrm{1}/\mathrm{2}$, this diverges; in order to be normalized, the process must be averaged at a finite resolution (below).

Although it is possible to follow Mandelbrot and Van Ness (1968) and derive many statistical properties in real space, a Fourier approach is not only more streamlined, but is also more powerful. The reason for the simplicity of the Fourier approach is that the Fourier transform (FT, indicated by the tilde) of the Weyl fractional derivative is symbolically

$$\begin{array}{}\text{(28)}& (i\mathit{\omega }{)}^h\stackrel{\text{FT}}{\leftrightarrow }{}_{-\mathrm{\infty }}{D}_t^h\end{array}$$

(e.g. Podlubny, 1999). This is simply the extension of the usual rule for the FT of integer-ordered derivatives. Therefore, since U_α,h and G_α,h are respectively solutions and Green's functions of the fractionally integrated fractional relaxation equation (Eq. 22), we have

$$\begin{array}{}\text{(29)}& \begin{array}{rl}\left(\right(i\mathit{\omega }{)}^{\mathit{\alpha }+h}& +\left(i\mathit{\omega }{)}^{\mathit{\alpha }}\right){\stackrel{\mathrm{̃}}{U}}_{\mathit{\alpha },h}=\stackrel{\mathrm{̃}}{\mathit{\gamma }}\\ & \stackrel{\text{FT}}{\leftrightarrow }({}_{-\mathrm{\infty }}{D}_t^{\mathit{\alpha }+h}+{}_{-\mathrm{\infty }}{D}_t^{\mathit{\alpha }})U_{\mathit{\alpha },h}=\mathit{\gamma },\\ \left(\right(i\mathit{\omega }{)}^{\mathit{\alpha }+h}& +\left(i\mathit{\omega }{)}^{\mathit{\alpha }}\right){\stackrel{\mathrm{̃}}{G}}_{\mathit{\alpha },h}=\mathrm{1}\\ & \stackrel{\text{FT}}{\leftrightarrow }({}_{-\mathrm{\infty }}{D}_t^{\mathit{\alpha }+h}+{}_{-\mathrm{\infty }}{D}_t^{\mathit{\alpha }})G_{\mathit{\alpha },h}=\mathit{\delta },\end{array}\end{array}$$

so that

$$\begin{array}{}\text{(30)}& \begin{array}{rl}& {\stackrel{\mathrm{̃}}{U}}_{\mathit{\alpha },h}\left(\mathit{\omega }\right)={\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}{\left(i\mathit{\omega }{)}^{\mathit{\alpha }}\right(\mathrm{1}+\left(i\mathit{\omega }{)}^h\right)}},\\ & {\stackrel{\mathrm{̃}}{G}}_{\mathit{\alpha },h}\left(\mathit{\omega }\right)={\displaystyle \frac{\mathrm{1}}{\left(i\mathit{\omega }{)}^{\mathit{\alpha }}\right(\mathrm{1}+\left(i\mathit{\omega }{)}^h\right)}},\mathrm{0}<\mathit{\alpha }<\mathrm{1},\\ & \mathrm{0}<h<\mathrm{2}.\end{array}\end{array}$$

We see that in the limit h→0, U_α,0 is an α-order fGn process (see e.g. Eq. 23).

Now we can use the fact that the white noise γ has a flat spectrum:

$$\begin{array}{}\text{(31)}& \begin{array}{rl}〈\stackrel{\mathrm{̃}}{\mathit{\gamma }}\left(\mathit{\omega }\right)\stackrel{\mathrm{̃}}{\mathit{\gamma }}\left({\mathit{\omega }}^{\prime }\right)〉& =\mathit{\delta }(\mathit{\omega }+\mathit{\omega }{}^{\prime })\langle {\left|\stackrel{\mathrm{̃}}{\mathit{\gamma }}\left(\mathit{\omega }\right)\right|}^{\mathrm{2}}\rangle =\mathrm{2}\mathit{\pi }\mathit{\delta }(\mathit{\omega }+\mathit{\omega }{}^{\prime })\\ & \stackrel{\text{FT}}{\leftrightarrow }\langle \mathit{\gamma }\left(t\right)\mathit{\gamma }\left(t{}^{\prime }\right)\rangle =\mathit{\delta }(t-t{}^{\prime }).\end{array}\end{array}$$

The modulus (vertical bars) intervenes since for any real function f(t) we have $\stackrel{\mathrm{̃}}{f}\left(\mathit{\omega }\right)={\stackrel{\mathrm{̃}}{f}}^{\ast }(-\mathit{\omega })$, where the superscript “*” indicates a complex conjugate.

Application of Eq. (31) leads to

$$\begin{array}{}\text{(32)}& \begin{array}{rl}{R}_{\mathit{\alpha },h}\left(t\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathit{\omega }t}{E}_U\left(\mathit{\omega }\right)\mathrm{d}\mathit{\omega },\\ E_U\left(\mathit{\omega }\right)& =\langle \left|{\stackrel{\mathrm{̃}}{U}}_{\mathit{\alpha },h}\right(\mathit{\omega }){|}^{\mathrm{2}}\rangle \\ & ={\displaystyle \frac{\mathrm{1}}{{\left|\mathit{\omega }\right|}^{\mathrm{2}\mathit{\alpha }}(\mathrm{1}+(-i\mathit{\omega }{)}^h\left)\right(\mathrm{1}+\left(i\mathit{\omega }{)}^h\right)}};\end{array}\end{array}$$

i.e. the spectrum E_U is the FT of the correlation function R_α,h(t) (the Wiener–Khinchin theorem). Applying this to U_α,h, we obtain

$$\begin{array}{}\text{(33)}& R_{\mathit{\alpha },h}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\cos \left(\mathit{\omega }t\right)\mathrm{d}\mathit{\omega }}{{\left|\mathit{\omega }\right|}^{\mathrm{2}\mathit{\alpha }}(\mathrm{1}+(i\mathit{\omega }{)}^h\left)\right(\mathrm{1}+(-i\mathit{\omega }{)}^h)}}.\end{array}$$

This shows that $R_{\mathit{\alpha },h}\left(t\right)=R_{\mathit{\alpha },h}(-t)$, so that below, we only consider t≥0.

Since R_α,h(0) diverges for $\mathit{\alpha }+h<\mathrm{1}/\mathrm{2}$, we consider the integral $Q_{\mathit{\alpha },\underset{\mathrm{¯}}{h}}$ of the process (the “motion”) from which we can easily compute the average. The corresponding variance V_α,h is

$$\begin{array}{}\text{(34)}& V_{\mathit{\alpha },h}\left(t\right)=\langle Q_{\mathit{\alpha },h}(t{)}^{\mathrm{2}}\rangle ,Q_{\mathit{\alpha },h}(t)=\underset{\mathrm{0}}{\overset{t}{\int }}{U}_{\mathit{\alpha },h}(v)\mathrm{d}v.\end{array}$$

In terms of ${\stackrel{\mathrm{̃}}{U}}_{\mathit{\alpha },h}\left(\mathit{\omega }\right)$,

$$\begin{array}{ll}{\displaystyle }{V}_{\mathit{\alpha },h}\left(t\right)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{(\mathrm{1}-\cos \mathit{\omega }t)}{{\mathit{\omega }}^{\mathrm{2}}}}\langle \left|{\stackrel{\mathrm{̃}}{U}}_{\mathit{\alpha },h}\right(\mathit{\omega }){|}^{\mathrm{2}}\rangle \mathrm{d}\mathit{\omega }\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{(\mathrm{1}-\cos \mathit{\omega }t)}{{\left|\mathit{\omega }\right|}^{\mathrm{2}+\mathrm{2}\mathit{\alpha }}}}{\displaystyle \frac{\mathrm{d}\mathit{\omega }}{(\mathrm{1}+(i\mathit{\omega }{)}^h\left)\right(\mathrm{1}+(-i\mathit{\omega }{)}^h)}},\\ \text{(35)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }<\mathrm{1}/\mathrm{2},\mathrm{0}<h<\mathrm{2}.\end{array}$$

We see that at low frequencies, when $\mathit{\alpha }\ge \mathrm{1}/\mathrm{2}$, the integral diverges for all t. Also note that a series expansion for V_α,h(t) in t will only have even-ordered integer power terms.

Comparing Eqs. (33) and (35), we see that R and V are linked by the simple relation

$$\begin{array}{}\text{(36)}& R_{\mathit{\alpha },h}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{V}_{\mathit{\alpha },h}\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}.\end{array}$$

Therefore, by integrating Eq. (26) (twice), we can express V_α,h in terms of G_α,h:

$$\begin{array}{ll}{\displaystyle }{V}_{\mathit{\alpha },h}\left(t\right)=& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left(G_{\mathit{\alpha }+\mathrm{1},h}\right(t+v)-G_{\mathit{\alpha }+\mathrm{1},h}(v){)}^{\mathrm{2}}\mathrm{d}v\\ \text{(37)}& {\displaystyle }& {\displaystyle }+\underset{\mathrm{0}}{\overset{t}{\int }}{G}_{\mathit{\alpha }+\mathrm{1},h}(v{)}^{\mathrm{2}}\mathrm{d}v.\end{array}$$

This can be verified by differentiation and using

$${\displaystyle \frac{\mathrm{d}{G}_{\mathit{\alpha }+\mathrm{1},h}}{\mathrm{d}t}}=G_{\mathit{\alpha },h}.$$

The basic behaviour can be understood in the Fourier domain. First, putting t=0 in Eq. (32) (i.e. “Parseval's theorem”), we have

$$\begin{array}{ll}{\displaystyle }{R}_{\mathit{\alpha },h}\left(\mathrm{0}\right)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{E}_U\left(\mathit{\omega }\right)\mathrm{d}\mathit{\omega }\\ \text{(38)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}\mathit{\omega }}{{\left|\mathit{\omega }\right|}^{\mathrm{2}\mathit{\alpha }}(\mathrm{1}+(i\mathit{\omega }{)}^h\left)\right(\mathrm{1}+(-i\mathit{\omega }{)}^h)}},\end{array}$$

so that when $\mathit{\alpha }+h<\mathrm{1}/\mathrm{2}$, R diverges at high frequencies (small t), and hence to represent a physical process (here, the Earth's temperature), the process must be averaged over a finite-resolution τ. When $\mathit{\alpha }+h>\mathrm{1}/\mathrm{2}$, R(0) is finite and can therefore be used to obtain a normalized autocorrelation function (Eq. 27).

From Eq. (32), we may also easily obtain the asymptotic high- and low-frequency behaviours of the energy spectrum:

$$\begin{array}{}\text{(39)}& \begin{array}{rl}& E_U\left(\mathit{\omega }\right)\\ & \propto \begin{array}{cc}{\mathit{\omega }}^{-\mathrm{2}(\mathit{\alpha }+h)}+O\left({\mathit{\omega }}^{-\mathrm{2}\mathit{\alpha }-\mathrm{3}h}\right),& \mathit{\omega }\gg \mathrm{1},\\ {\mathit{\omega }}^{-\mathrm{2}\mathit{\alpha }}-\mathrm{2}\cos \left({\displaystyle \frac{\mathit{\pi }h}{\mathrm{2}}}\right){\mathit{\omega }}^{h-\mathrm{2}\mathit{\alpha }}+O\left({\mathit{\omega }}^{\mathrm{2}h-\mathrm{2}\mathit{\alpha }}\right),& \mathit{\omega }\ll \mathrm{1}.\end{array}\end{array}\end{array}$$

2.5 Finite-resolution processes

When $\mathit{\alpha }+h<\mathrm{1}/\mathrm{2}$, the process does not converge at any instant t: it is a noise, a generalized function. To represent the Earth's temperature, it must therefore be averaged at a finite-resolution τ:

$$\begin{array}{}\text{(40)}& U_{\mathit{\alpha },h,\mathit{\tau }}\left(t\right)={\displaystyle \frac{Q_{\mathit{\alpha },h}\left(t\right)-Q_{\mathit{\alpha },h}(t-\mathit{\tau })}{\mathit{\tau }}}.\end{array}$$

Applying Eqs. (34) and (40), we obtain the “resolution τ” autocorrelation:

$$\begin{array}{}\text{(41)}& \begin{array}{rl}{R}_{\mathit{\alpha },h,\mathit{\tau }}\left(\mathrm{\Delta }t\right)=& 〈U_{\mathit{\alpha },h,\mathit{\tau }}\left(t\right)U_{\mathit{\alpha },h,\mathit{\tau }}(t-\mathrm{\Delta }t)〉\\ =& {\mathit{\tau }}^{-\mathrm{2}}\langle \left(Q_{\mathit{\alpha },h}\right(t)-Q_{\mathit{\alpha },h}(t-\mathit{\tau }\left)\right)\\ & \begin{array}{ll}\begin{array}{l}\cdot \left(Q_{\mathit{\alpha },h}(t-\mathrm{\Delta }t)-Q_{\mathit{\alpha },h}(t-\mathrm{\Delta }t-\mathit{\tau })\right)\rangle \\ ={\mathit{\tau }}^{-\mathrm{2}}\frac{\mathrm{1}}{\mathrm{2}}\left(V_{\mathit{\alpha },h}\right(\mathrm{\Delta }t-\mathit{\tau })+V_{\mathit{\alpha },h}(\mathrm{\Delta }t+\mathit{\tau })\end{array},& \mathrm{\Delta }t\ge \mathit{\tau },\end{array}\\ & -\mathrm{2}{V}_{\mathit{\alpha },h}\left(\mathrm{\Delta }t\right))\\ R_{\mathit{\alpha },h,\mathit{\tau }}\left(\mathrm{0}\right)=& {\mathit{\tau }}^{-\mathrm{2}}{V}_{\mathit{\alpha },h}\left(\mathit{\tau }\right).\end{array}\end{array}$$

Alternatively, measuring time in units of the resolution $\mathit{\lambda }=\mathrm{\Delta }t/\mathit{\tau }$,

$$\begin{array}{}\text{(42)}& \begin{array}{rl}{R}_{\mathit{\alpha },h,\mathit{\tau }}\left(\mathit{\lambda }\mathit{\tau }\right)=& 〈U_{\mathit{\alpha },h,\mathit{\tau }}\left(t\right)U_{\mathit{\alpha },h,\mathit{\tau }}(t-\mathit{\lambda }\mathit{\tau })〉\\ =& {\mathit{\tau }}^{-\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(V_{\mathit{\alpha },h}\right((\mathit{\lambda }-\mathrm{1})\mathit{\tau })+V_{\mathit{\alpha },h}((\mathit{\lambda }+\mathrm{1})\mathit{\tau })\\ & -\mathrm{2}{V}_{\mathit{\alpha },h}\left(\mathit{\lambda }\mathit{\tau }\right)),\mathit{\lambda }\ge \mathrm{1}.\end{array}\end{array}$$

$R_{a,h,\mathit{\tau }}$ can be conveniently written in terms of centred finite differences:

$$\begin{array}{ll}{\displaystyle }{R}_{\mathit{\alpha },h,\mathit{\tau }}\left(\mathit{\lambda }\mathit{\tau }\right)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{\Delta }}_{\mathit{\tau }}^{\mathrm{2}}{V}_{\mathit{\alpha },h}\left(\mathit{\lambda }\mathit{\tau }\right)\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_{\mathit{\alpha },h}^{\prime \prime }\left(\mathrm{\Delta }t\right),\\ \text{(43)}& {\displaystyle }{\mathrm{\Delta }}_{\mathit{\tau }}f\left(t\right)& {\displaystyle }={\displaystyle \frac{f(t+\mathit{\tau }/\mathrm{2})-f(t-\mathit{\tau }/\mathrm{2})}{\mathit{\tau }}}.\end{array}$$

The finite-difference formula is valid for Δt≥τ. For finite τ, it allows us to obtain the correlation behaviour by replacing the second difference with a second derivative, an approximation that is very good except when Δt is close to τ. Taking the limit τ→0 in Eq. (43), we obtain the second derivative formula Eq. (36).

3 Application to fBm, fGn, fRm, and fRn

3.1 fBm and fGn

The above derivations were for noises and motions derived from differential operators whose impulse and step Green functions had convergent V_α,h(t). Before applying them to fRn and fRm, we illustrate this by applying them first to fBm and fGn.

The fBm results are obtained by using the fGn step Green function (Eq. 13) in Eq. (35) with h=0 to obtain

$$\begin{array}{}\text{(44)}& \begin{array}{rl}{V}_h^{\left(\text{fBm}\right)}\left(t\right)=& \mathrm{4}{V}_{\mathit{\alpha }=h,\mathrm{0}}\left(t\right)\\ & =\left({\displaystyle \frac{\mathrm{2}\sin \left(\mathit{\pi }h\right)\mathrm{\Gamma }(-\mathrm{1}-\mathrm{2}h)}{\mathit{\pi }}}\right)t^{\mathrm{2}h+\mathrm{1}};\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\le h<{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}.\end{array}\end{array}$$

The standard normalization and parametrization are

$$\begin{array}{}\text{(45)}& \begin{array}{rl}{N}_h=K_h& ={\left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}\sin \left(\mathit{\pi }h\right)\mathrm{\Gamma }(-\mathrm{1}-\mathrm{2}h)}}\right)}^{\mathrm{1}/\mathrm{2}}\\ & ={\left(-{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}\cos \left(\mathit{\pi }H\right)\mathrm{\Gamma }(-\mathrm{2}H)}}\right)}^{\mathrm{1}/\mathrm{2}};\\ & H=h+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}};\mathrm{0}\le H<\mathrm{1}.\end{array}\end{array}$$

This normalization turns out to be convenient not only for fBm, but also for fRm, so that for the normalized process,

$$\begin{array}{}\text{(46)}& V_H^{\left(\text{fBm}\right)}\left(t\right)=t^{\mathrm{2}h+\mathrm{1}}=t^{\mathrm{2}H};\mathrm{0}\le H<\mathrm{1},\end{array}$$

where we have introduced the standard fBm parameter $H=h+\mathrm{1}/\mathrm{2}$, so that

$$\begin{array}{}\text{(47)}& \begin{array}{rl}& \langle \mathrm{\Delta }{B}_H(\mathrm{\Delta }t{)}^{\mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}=\mathrm{\Delta }{t}^H,\\ & \mathrm{\Delta }{B}_H\left(\mathrm{\Delta }t\right)=B_H\left(t\right)-B_H(t-\mathrm{\Delta }t),\end{array}\end{array}$$

and hence H is the fluctuation exponent for fBm. Note that fBm is usually defined as the Gaussian process with V_H given by Eq. (46), i.e. with this normalization (e.g. Biagini et al., 2008).

We can now calculate the correlation function relevant for the fGn statistics. With the above normalization,

$$\begin{array}{}\text{(48)}& \begin{array}{rl}& R_{h,\mathit{\tau }}^{\left(\text{fGn}\right)}\left(\mathit{\lambda }\mathit{\tau }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\tau }}^{\mathrm{2}h-\mathrm{1}}\left(\right(\mathit{\lambda }+\mathrm{1}{)}^{\mathrm{2}h+\mathrm{1}}+(\mathit{\lambda }-\mathrm{1}{)}^{\mathrm{2}h+\mathrm{1}}\\ & -\mathrm{2}{\mathit{\lambda }}^{\mathrm{2}h+\mathrm{1}}),\mathit{\lambda }\ge \mathrm{1},-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}<h<{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},\\ & R_{h,\mathit{\tau }}^{\left(\text{fGn}\right)}\left(\mathrm{0}\right)={\mathit{\tau }}^{\mathrm{2}h-\mathrm{1}},\\ & R_{H,\mathit{\tau }}^{\left(\text{fGn}\right)}\left(\mathit{\lambda }\mathit{\tau }\right)\approx h(\mathrm{2}h+\mathrm{1})(\mathit{\lambda }\mathit{\tau }{)}^{\mathrm{2}h-\mathrm{1}}\\ & =H(\mathrm{2}H-\mathrm{1})(\mathit{\lambda }\mathit{\tau }{)}^{\mathrm{2}(H-\mathrm{1})};\mathit{\lambda }\gg \mathrm{1},\end{array}\end{array}$$

the bottom approximations are valid for large-scale ratios λ. We note the difference in sign for $H>\mathrm{1}/\mathrm{2}$ (“persistence”) and for $H<\mathrm{1}/\mathrm{2}$ (“anti-persistence”). When $H=\mathrm{1}/\mathrm{2}$, the noise corresponds to standard Brownian motion, and it is uncorrelated.

3.2 fRm and fRn

3.2.1  R_α,h(t)

Since fRm and fRn are Gaussian, their properties are determined by their second-order statistics, by V_α,h(t) and R_α,h(t). These statistics are second order in G_α,h(t) and can most easily be determined using the Fourier representation of G_α,h(t) (Sect. 2.4, Appendices A and B). The development is challenging because unlike the G_α,h(t) functions that are entirely expressed in series of fractional powers of t, V_α,h(t) and R_α,h(t) involve mixed fractional and integer power expansions; the details are given in the Appendices, and here we summarize the main results.

First, for the noises, we have

$$\begin{array}{}\text{(49)}& \begin{array}{rl}& R_{\mathit{\alpha },h}\left(t\right)=\sum _{n=\mathrm{2}}^{\mathrm{\infty }}{D}_n\mathrm{\Gamma }(\mathrm{1}-hn-\mathrm{2}\mathit{\alpha })t^{-\mathrm{1}+hn+\mathrm{2}\mathit{\alpha }}\\ & +\sum _{j=\mathrm{1},\text{odd}}^{\mathrm{\infty }}{F}_j{\displaystyle \frac{t^{j-\mathrm{1}}}{\mathrm{\Gamma }\left(j\right)}},\\ & F_j=-{\displaystyle \frac{\cos \mathit{\pi }\left(\frac{h}{\mathrm{2}}+\mathit{\alpha }\right)}{h\sin \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)\sin \left(\frac{\mathit{\pi }}{h}(j-\mathrm{2}\mathit{\alpha })\right)}},\\ & D_n=(-\mathrm{1}{)}^n{\displaystyle \frac{\sin \left(\frac{n\mathit{\pi }h}{\mathrm{2}}+\mathit{\alpha }\mathit{\pi }\right)\sin \left(\frac{(n-\mathrm{1})\mathit{\pi }h}{\mathrm{2}}\right)}{\mathit{\pi }\sin \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)}}.\end{array}\end{array}$$

At small t, the lowest-order terms dominate, and the normalized autocorrelations are thus

$$\begin{array}{}\text{(50)}& \begin{array}{rl}& R_{\mathit{\alpha },h}^{\left(\text{norm}\right)}\left(t\right)=(h+\mathit{\alpha })(\mathrm{1}+\mathrm{2}(h+\mathit{\alpha }\left)\right)t^{-\mathrm{1}+\mathrm{2}(h+\mathit{\alpha })}\\ & +O\left(t^{-\mathrm{1}+\mathrm{3}h+\mathrm{2}\mathit{\alpha }}\right),\\ & \mathit{\tau }\ll t\ll \mathrm{1},\mathrm{0}<(h+\mathit{\alpha })<\mathrm{1}/\mathrm{2},\\ & R_{\mathit{\alpha },h}^{\left(\text{norm}\right)}\left(t\right)=\mathrm{1}-{\displaystyle \frac{\left|\mathrm{\Gamma }(\mathrm{1}-\mathrm{2}(h+\mathit{\alpha }\left)\right)\right|\sin \left(\mathit{\pi }\right(h+\mathrm{2}\mathit{\alpha }\left)\right)}{\mathit{\pi }{F}_{\mathrm{1}}}}\\ & \cdot t^{-\mathrm{1}+\mathrm{2}(h+\mathit{\alpha })}+O\left(t^{-\mathrm{1}+\mathrm{3}h+\mathrm{2}\mathit{\alpha }}\right),\\ & t\ll \mathrm{1},\mathrm{1}/\mathrm{2}<(h+\mathit{\alpha })<\mathrm{3}/\mathrm{2},\\ & R_{\mathit{\alpha },h}^{\left(\text{norm}\right)}\left(t\right)=\mathrm{1}+{\displaystyle \frac{t^{\mathrm{2}}}{\mathrm{2}{F}_{\mathrm{1}}}}{F}_{\mathrm{3}}+O\left(t^{-\mathrm{1}+\mathrm{2}(h+\mathit{\alpha })}\right)\mathrm{\dots },\\ & t\ll \mathrm{1},\mathrm{3}/\mathrm{2}<(h+\mathit{\alpha })<\mathrm{2}\end{array}\end{array}$$

(note that F₃<0 for $\mathrm{3}/\mathrm{2}<h+\mathit{\alpha }<\mathrm{2}$; see Appendix A). We see that at small t, the behaviour of the normalized autocorrelations depends essentially on the sum h+α; in particular, when $h+\mathit{\alpha }<\mathrm{1}/\mathrm{2}$, the process is effectively an fGn process with an effective fluctuation exponent $H=-\mathrm{1}/\mathrm{2}+(h+\mathit{\alpha })$. This is to be expected since α+h is the highest-order term in the fractionally integrated fractional relaxation equation (Eq. 22).

3.2.2  V_α,h(t)

Integrating twice

$$V_{\mathit{\alpha },h}\left(t\right)=\mathrm{2}\underset{\mathrm{0}}{\overset{t}{\int }}\left(\underset{\mathrm{0}}{\overset{v}{\int }}{R}_{\mathit{\alpha },h}\left(u\right)\mathrm{d}u\right)\mathrm{d}v,$$

we obtain

$$\begin{array}{ll}{\displaystyle }{V}_{\mathit{\alpha },h}\left(t\right)=& {\displaystyle }\mathrm{2}\sum _{n=\mathrm{2}}^{\mathrm{\infty }}{D}_n\mathrm{\Gamma }(-\mathrm{1}-hn-\mathrm{2}\mathit{\alpha })t^{\mathrm{1}+hn+\mathrm{2}\mathit{\alpha }}\\ {\displaystyle }& {\displaystyle }+\mathrm{2}\sum _{j=\mathrm{1},\text{odd}}^{\mathrm{\infty }}{F}_j{\displaystyle \frac{t^{j+\mathrm{1}}}{\mathrm{\Gamma }(j+\mathrm{2})}},\\ \text{(51)}& {\displaystyle }& {\displaystyle }\mathrm{0}<h<\mathrm{2},\mathrm{0}\le \mathit{\alpha }<\mathrm{1}/\mathrm{2}.\end{array}$$

When $\mathrm{0}<\mathit{\alpha }+h<\mathrm{1}/\mathrm{2}$, the leading (n=2) term for V_α,h is $t^{\mathrm{1}+\mathrm{2}(h+\mathit{\alpha })}$ ($\propto V_{\mathit{\alpha }+h}^{\left(\text{fBm}\right)}$), so that the fBm coefficient can be used for normalization using $R_{\mathit{\alpha },h,\mathit{\tau }}\left(\mathrm{0}\right)={\mathit{\tau }}^{-\mathrm{2}}{V}_{\mathit{\alpha },h}\left(\mathit{\tau }\right)$. When $h+\mathit{\alpha }>\mathrm{1}/\mathrm{2}$, this normalization becomes negative, so that it cannot be used; however, in this case, $R_{\mathit{\alpha },h}\left(\mathrm{0}\right)=F_{\mathrm{1}}$ and may be used for normalization instead. For an analytic expression, convergence properties including numerical results and modified expansions converge more rapidly; see Appendix A and, for the special case $h=\mathrm{1}/\mathrm{2}$, Appendix B.

For convenience, the leading terms of the normalized V_α,h are

$$\begin{array}{ll}{\displaystyle }{V}_{\mathit{\alpha },h}^{\left(\text{norm}\right)}\left(t\right)=& {\displaystyle }{t}^{\mathrm{1}+\mathrm{2}(h+\mathit{\alpha })}+O\left(t^{\mathrm{1}+\mathrm{3}h+\mathrm{2}\mathit{\alpha }}\right)+O\left(t^{\mathrm{2}}\right),\\ \text{(52)}& {\displaystyle }& {\displaystyle }\mathrm{0}<(h+\mathit{\alpha })<\mathrm{1}/\mathrm{2},{\displaystyle }{V}_{\mathit{\alpha },h}^{\left(\text{norm}\right)}\left(t\right)=& {\displaystyle }{t}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{2}\mathrm{\Gamma }(-\mathrm{1}-\mathrm{2}(h+\mathit{\alpha }\left)\right)\sin \left(\mathit{\pi }\right(h+\mathrm{2}\mathit{\alpha }\left)\right)}{\mathit{\pi }{F}_{\mathrm{1}}}}\\ {\displaystyle }& {\displaystyle }\cdot t^{\mathrm{1}+\mathrm{2}(h+\mathit{\alpha })}+O\left(t^{\mathrm{1}+\mathrm{3}h+\mathrm{2}\mathit{\alpha }}\right),\\ {\displaystyle }& {\displaystyle }\mathrm{1}/\mathrm{2}<(h+\mathit{\alpha })<\mathrm{3}/\mathrm{2},\\ {\displaystyle }{V}_{\mathit{\alpha },h}^{\left(\text{norm}\right)}\left(t\right)=& {\displaystyle }{t}^{\mathrm{2}}+{\displaystyle \frac{F_{\mathrm{3}}}{\mathrm{12}{F}_{\mathrm{1}}}}{t}^{\mathrm{4}}+O\left(t^{\mathrm{2}(h+\mathit{\alpha })+\mathrm{1}}\right),\\ {\displaystyle }& {\displaystyle }\mathrm{3}/\mathrm{2}<(h+\mathit{\alpha })<\mathrm{2}.\end{array}$$

3.2.3 Asymptotic expansions

For multidecadal global climate projections, the relaxation time has been estimated at ≈ 5 years (Procyk et al., 2020, 2022), so that we are interested in the long-time behaviour (exploited for example in Hébert et al., 2021). For this, asymptotic expansions are needed: in Appendix A we show

$$\begin{array}{ll}{\displaystyle }{R}_{\mathit{\alpha },h}\left(t\right)=& {\displaystyle }-\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{D}_{-n}\mathrm{\Gamma }(\mathrm{1}+nh-\mathrm{2}\mathit{\alpha })t^{\mathrm{2}\mathit{\alpha }-(\mathrm{1}+nh)}\\ \text{(53)}& {\displaystyle }& {\displaystyle }+P_{\mathit{\alpha },h,+}\left(t\right),t\gg \mathrm{1},\end{array}$$

where $P_{\mathit{\alpha },h,+}\left(t\right)=\mathrm{0}$ for h<1, while for $\mathrm{1}<h<\mathrm{2}$ it has exponentially damped oscillations (see Fig. 3d and Appendix A).

Figure 3 The normalized correlation functions R_0,h for fRn corresponding to the V_0,h function in Fig. 4: $\mathrm{0}<h<\mathrm{1}/\mathrm{2}$ (a), $\mathrm{1}/\mathrm{2}<h<\mathrm{1}$ (b), $\mathrm{1}<h<\mathrm{3}/\mathrm{2}$ (c), and $\mathrm{3}/\mathrm{2}<h<\mathrm{2}$ (d). In each plot, the curves correspond to h increasing from bottom to top in units of $\mathrm{1}/\mathrm{10}$ starting from $\mathrm{1}/\mathrm{20}$ (a) to $\mathrm{39}/\mathrm{20}$ (d). For $h<\mathrm{1}/\mathrm{2}$, the resolution is important since $R_{\mathrm{0},h,\mathit{\tau }}$ diverges at small τ. In (a), $R_{\mathrm{0},h,\mathit{\tau }}$ is shown with $\mathit{\tau }={\mathrm{10}}^{-\mathrm{5}}$; they were normalized to the value at resolution $\mathit{\tau }={\mathrm{10}}^{-\mathrm{5}}$, and for $h>\mathrm{1}/\mathrm{2}$, the curves are normalized with $F_{\mathrm{3}}^{-\mathrm{1}/\mathrm{2}}$. In all cases, the large t slope is $-\mathrm{1}-h$.

For pure fRn processes, a useful formula is

$$\begin{array}{ll}{\displaystyle }{R}_{\mathrm{0},h}\left(t\right)=& {\displaystyle }\sum _{n=\mathrm{1}}^{\mathrm{\infty }}(-\mathrm{1}{)}^n{\displaystyle \frac{\mathrm{1}+\cot \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)\tan \left(\frac{n\mathit{\pi }h}{\mathrm{2}}\right)}{\mathrm{2}\mathrm{\Gamma }(-nh)}}{t}^{-(\mathrm{1}+nh)}\\ \text{(54)}& {\displaystyle }& {\displaystyle }+P_{\mathrm{0},h,+}\left(t\right),t\gg \mathrm{1},\end{array}$$

or, more generally,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{R}_{\mathit{\alpha },h}\left(t\right)={\displaystyle \frac{\mathrm{\Gamma }(\mathrm{1}-\mathrm{2}\mathit{\alpha })\sin \left(\mathit{\pi }\mathit{\alpha }\right)}{\mathit{\pi }}}{t}^{\mathrm{2}\mathit{\alpha }-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\cos \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)}{\cos \left(\frac{\mathit{\pi }h}{\mathrm{2}}-\mathit{\pi }\mathit{\alpha }\right)\mathrm{\Gamma }(\mathrm{2}\mathit{\alpha }-h)}}{t}^{\mathrm{2}\mathit{\alpha }-(\mathrm{1}+h)}+\mathrm{\dots },\\ \text{(55)}& {\displaystyle }& {\displaystyle }t\gg \mathrm{1},\mathrm{0}\le h<\mathrm{2},\mathrm{0}\le \mathit{\alpha }<\mathrm{1}/\mathrm{2}.\end{array}$$

We see that when α≠0, D₀>0, so that, as expected, the leading behaviour has no h dependence: it is only due to the long-range correlations in the forcing; we obtain the fGn result: $\approx t^{\mathrm{2}\mathit{\alpha }-\mathrm{1}}$. For pure fRn processes this reduces to

$$R_{\mathrm{0},h}\left(t\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }(-h)}}{t}^{-\mathrm{1}-h}$$

(note that $\mathrm{\Gamma }(-h)<\mathrm{0}$ for $\mathrm{0}<h<\mathrm{1}$).

Integrating R_α,h twice and doubling, we obtain

$$\begin{array}{}\text{(56)}& \begin{array}{rl}{V}_{\mathit{\alpha },h}\left(t\right)=& {\displaystyle \frac{\mathrm{2}\mathrm{\Gamma }(-\mathrm{1}-\mathrm{2}\mathit{\alpha })\sin \left(\mathit{\pi }\mathit{\alpha }\right)}{\mathit{\pi }}}{t}^{\mathrm{1}+\mathrm{2}\mathit{\alpha }}+a_{\mathit{\alpha },h}t+b_{\mathit{\alpha },h}\\ & -{\displaystyle \frac{\mathrm{1}+\cos \left(\mathit{\pi }h\right)-\sin \left(\mathit{\pi }h\right)\cot \left(\mathit{\pi }\right(h-\mathrm{2}\mathit{\alpha }\left)\right)}{\mathrm{\Gamma }(\mathrm{2}-(h-\mathrm{2}\mathit{\alpha }\left)\right)}}\\ & \cdot t^{\mathrm{1}+\mathrm{2}\mathit{\alpha }-h}+\mathrm{\dots },t\gg \mathrm{1}\end{array}\end{array}$$

(the full expansion is given in Appendix A; see Fig. 4 for plots). The constants of integration a_α,h and b_α,h are not determined since the expansion is not valid at t=0; they can be determined numerically if needed. However, in the limit α→0 (the pure fRn case), the leading term is exactly t (corresponding to ordinary Brownian motion), so that an extra a_0,h is not needed (Appendix A). When α>0, the far-left (fGn) term from the forcing dominates; at large enough t, $V_{\mathit{\alpha },h}\left(t\right)\propto t^{\mathrm{2}H}$ with $H=\mathit{\alpha }+\mathrm{1}/\mathrm{2}$, and the corresponding motion is an fBm.

Figure 4 The normalized V_0,h functions for the various ranges of h for fRm. The plots from (a) to (d) are for the ranges $\mathrm{0}<h<\mathrm{1}/\mathrm{2}$, $\mathrm{1}/\mathrm{2}<h<\mathrm{1}$, $\mathrm{1}<h<\mathrm{3}/\mathrm{2}$, and $\mathrm{3}/\mathrm{2}<h<\mathrm{2}$. Within each plot, the lines are for h increasing in units of $\mathrm{1}/\mathrm{10}$ starting at a value $\mathrm{1}/\mathrm{20}$ above the plot minimum; overall, h increases in units of $\mathrm{1}/\mathrm{10}$ starting at values $\mathrm{1}/\mathrm{20}$ (a) to $\mathrm{39}/\mathrm{20}$ (d) (e.g. for a, the lines are for $h=\mathrm{1}/\mathrm{20}$, $\mathrm{3}/\mathrm{10}$, $\mathrm{5}/\mathrm{20}$, $\mathrm{7}/\mathrm{20}$, and $\mathrm{9}/\mathrm{20}$). For all h the large t behaviour is linear (slope =1, although note the oscillations for the lower right-hand plot (d) for $\mathrm{3}/\mathrm{2}<h<\mathrm{2}$). For small t, the slopes are 1+2h ($\mathrm{0}<h\le \mathrm{1}/\mathrm{2}$) and 2 ($\mathrm{1}/\mathrm{2}\le h<\mathrm{2}$).

Using the above results, we see that there are three limiting fRn$/$fRm cases that yield fGn$/$fBm processes:

$$\begin{array}{llll}{\displaystyle }& {\displaystyle }{R}_{\mathit{\alpha },\mathrm{0}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{R}_{\mathit{\alpha }}^{\left(\text{fGn}\right)}\left(t\right),& {\displaystyle }& {\displaystyle }\mathrm{0}<\mathit{\alpha }<\mathrm{1}/\mathrm{2},h=\mathrm{0},\\ {\displaystyle }& {\displaystyle }{R}_{\mathit{\alpha },h}\left(t\right)=R_{\mathit{\alpha }}^{\left(\text{fGn}\right)}\left(t\right),& {\displaystyle }& {\displaystyle }\mathrm{0}<\mathit{\alpha }<\mathrm{1}/\mathrm{2},t\gg \mathrm{1},\\ \text{(57)}& {\displaystyle }& {\displaystyle }{R}_{\mathit{\alpha },h}\left(t\right)=R_{\mathit{\alpha }+h}^{\left(\text{fGn}\right)}\left(t\right),& {\displaystyle }& {\displaystyle }\mathrm{0}<\mathit{\alpha }+h<\mathrm{1}/\mathrm{2},t\ll \mathrm{1}.\end{array}$$

3.3 Haar fluctuations

A useful statistical characterization of the processes is by the statistics of their Haar fluctuations over an interval Δt. For an interval Δt, Haar fluctuations (based on Haar wavelets) are the differences between the averages of the first and second halves of an interval. For a process U, the Haar fluctuation is

$$\begin{array}{}\text{(58)}& \begin{array}{rl}\mathrm{\Delta }U(\mathrm{\Delta }t{)}_{\text{Haar}}=& {\displaystyle \frac{\mathrm{2}}{\mathrm{\Delta }t}}\underset{t-\mathrm{\Delta }t/\mathrm{2}}{\overset{t}{\int }}U\left(v\right)\mathrm{d}v\\ & -{\displaystyle \frac{\mathrm{2}}{\mathrm{\Delta }t}}\underset{t-\mathrm{\Delta }t}{\overset{t-\mathrm{\Delta }t/\mathrm{2}}{\int }}U\left(v\right)\mathrm{d}v.\end{array}\end{array}$$

In terms of the process at resolution $\mathrm{\Delta }t/\mathrm{2}$ (i.e. averaged at this scale), $U_{\mathrm{\Delta }t/\mathrm{2}}\left(t\right)$:

$$\begin{array}{}\text{(59)}& \mathrm{\Delta }U(\mathrm{\Delta }t{)}_{\text{Haar}}={\displaystyle \frac{\mathrm{2}}{\mathrm{\Delta }t}}(U_{\mathrm{\Delta }t/\mathrm{2}}\left(t\right)-U_{\mathrm{\Delta }t/\mathrm{2}}(t-\mathrm{\Delta }t/\mathrm{2})).\end{array}$$

Therefore,

$$\begin{array}{}\text{(60)}& \langle \mathrm{\Delta }U(\mathrm{\Delta }t{)}_{\text{Haar}}^{\mathrm{2}}\rangle ={\left({\displaystyle \frac{\mathrm{2}}{\mathrm{\Delta }t}}\right)}^{\mathrm{2}}(\mathrm{4}V(\mathrm{\Delta }t/\mathrm{2})-V\left(\mathrm{\Delta }t\right)),\end{array}$$

where V(t) is the variance of the integral of U over an interval t (Eq. 34).

Using Eq. (60), we can determine the behaviour of the root mean square (rms) Haar fluctuations; terms like $V_{\mathit{\alpha },h}\left(t\right)\propto t^{\mathit{\xi }}$ contribute $\propto t^{\mathit{\xi }/\mathrm{2}-\mathrm{1}}$ to the rms Haar fluctuation ${〈\mathrm{\Delta }{U}_{\mathit{\alpha },h}(\mathrm{\Delta }t{)}_{\text{Haar}}^{\mathrm{2}}〉}^{\mathrm{1}/\mathrm{2}}$ (the exception is when ξ=2, which contributes nothing). Applying this equation to fGn parameter h, we obtain ${〈\mathrm{\Delta }{F}_h(\mathrm{\Delta }t{)}_{\text{Haar}}^{\mathrm{2}}〉}^{\mathrm{1}/\mathrm{2}}\propto \mathrm{\Delta }{t}^H$ with $H=h-\mathrm{1}/\mathrm{2}$.

Using the results above for V_α,h, we therefore obtain the leading exponents:

$$\begin{array}{}\text{(61)}& {\displaystyle }\begin{array}{l}\begin{array}{ll}\begin{array}{ll}H=h+\mathit{\alpha }-\mathrm{1}/\mathrm{2},& \mathrm{0}<h+\mathit{\alpha }<\mathrm{3}/\mathrm{2}\\ H=\mathrm{1},& \mathrm{3}/\mathrm{2}<h+\mathit{\alpha }<\mathrm{2}\end{array},& \mathrm{\Delta }t\ll \mathrm{1}\end{array}\\ \begin{array}{ll}H=\mathit{\alpha }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},& \mathrm{\Delta }t\gg \mathrm{1}.\end{array}\end{array}\end{array}$$

Figure 5 shows that the theory agrees well with the numerics.

Figure 5 The rms Haar fluctuation plots for the pure (α=0) fRn process for $\mathrm{0}<h<\mathrm{1}/\mathrm{2}$ (a), $\mathrm{1}/\mathrm{2}<h<\mathrm{1}$ (b), $\mathrm{1}<h<\mathrm{3}/\mathrm{2}$ (c), and $\mathrm{3}/\mathrm{2}<h<\mathrm{2}$ (d). The individual curves correspond to those of Figs. 3 and 4. The small Δt slopes follow the theoretical values $h-\mathrm{1}/\mathrm{2}$ up to $h=\mathrm{3}/\mathrm{2}$ (slope=1); for larger h, the small t slopes all equal 1. Also, at large t due to dominant V≈t terms, in all cases we obtain slopes $t^{-\mathrm{1}/\mathrm{2}}$.

For the range of α, h discussed here ($\mathrm{0}\le \mathit{\alpha }<\mathrm{1}/\mathrm{2}$, $\mathrm{0}\le h\le \mathrm{2}$), H spans the range $-\mathrm{1}/\mathrm{2}$ (white noise) to 1. In comparison, fGn processes have H covering the range $-\mathrm{1}<H<\mathrm{0}$ and fBm processes have $\mathrm{0}<H<\mathrm{1}$; therefore, depending on whether the process is observed at timescales below or above the relaxation timescale (Δt=1), fractionally integrated fRn processes can mimic fGn or fBm processes. If we consider the integrals – the motions – the value of H is increased by 1 (although for Haar fluctuations, it cannot exceed H=1). Overall, from an empirical viewpoint, if over some range of scales (that may only be a factor of 100 or less), it may be quite hard to distinguish the various models, especially since the transition from low- to high-frequency scaling may be very slow (see especially Appendix B for the $h=\mathrm{1}/\mathrm{2}$ case). Recent work shows that the maximum likelihood method may be the optimum parameter estimation technique (Procyk, 2021).

3.4 Sample processes

It is instructive to view some samples of fRn and fRm processes (here we consider only α=0). For simulations, both the small- and large-scale divergences must be considered. Starting with the approximate methods developed by Mandelbrot and Wallis (1969), it took some time for exact fBm and fGn simulation techniques to be developed (Hipel and McLeod, 1994; Palma, 2007). Fortunately, for fRm and fRn, the low-frequency situation is easier since the long-time memory is much smaller than for fBm and fGn. Therefore, as long as we are careful to always simulate series a few times longer than the relaxation time and then to throw away the earliest $\mathrm{2}/\mathrm{3}$ or $\mathrm{3}/\mathrm{4}$ of the simulation, the remainder will have accurate statistics. With this procedure to take care of low-frequency issues, we can therefore use the solution for fRn in the form of a convolution and use standard numerical convolution algorithms.

We must nevertheless be careful about the high frequencies since the impulse response Green functions G_0,h are singular for h<1. In order to avoid singularities, simulations of fRn are best made by first simulating the motions Q_0,h using $Q_{\mathrm{0},h}\propto G_{\mathrm{1},h}\ast \mathit{\gamma }$ and obtaining the resolution τ fRn, using $U_{\mathrm{0},h,\mathit{\tau }}\left(t\right)=\left(Q_{\mathrm{0},h}\right(t+\mathit{\tau })-Q_{\mathrm{0},h}(t\left)\right)/\mathit{\tau }$. Numerically, this allows us to use the smoother (nonsingular) G_1,h in the convolution rather than the singular G_0,h. The simulations shown in Figs. 6–9 follow this procedure, and the Haar fluctuation statistics were analysed, verifying the statistical accuracy of the simulations.

In order to clearly display the behaviours, recall that when t≫1, we showed that all the fRn converge to Gaussian white noises and the fRm to Brownian motions (albeit in a slow power law manner). At the other extreme, for t≪1, we obtain the fGn and fBm limits (when $\mathrm{0}<h<\mathrm{1}/\mathrm{2}$) and their generalizations for $\mathrm{1}/\mathrm{2}<h<\mathrm{2}$.

Figure 6 fRn and fRm simulations (left and right columns respectively) for $h=\mathrm{1}/\mathrm{10}$, $\mathrm{3}/\mathrm{10}$, and $\mathrm{5}/\mathrm{10}$ (top to bottom sets, all with α=0), i.e. the exponent range that overlaps with fGn and fBm. There are three simulations, each of length 2¹⁹ pixels, and each uses the same random seed with the unit scale equal to 2¹⁰ pixels (i.e. a resolution of $\mathit{\tau }={\mathrm{2}}^{-\mathrm{10}}$). The entire simulation therefore covers the range of scale $\mathrm{1}/\mathrm{1024}$ to 512 units. The fRm on the right is from the running sum of the fRn on the left. Starting at the top line of each set, we show 2¹⁰ points of the original series degraded in resolution by a factor 2⁹. Since the length is t=2⁹ units long, each pixel has resolution $\mathit{\tau }=\mathrm{1}/\mathrm{2}$. The second line of each set takes the segment of the upper line lying to the right of the dashed vertical line, $\mathrm{1}/\mathrm{8}$ of its length. It therefore spans t=0 to $t={\mathrm{2}}^{\mathrm{9}}/\mathrm{8}={\mathrm{2}}^{\mathrm{6}}$, but resolution was taken as $\mathit{\tau }={\mathrm{2}}^{-\mathrm{4}}$, and hence it is still 2¹⁰ pixels long. Since each pixel has a resolution of 2⁻⁴, the unit scale is 2⁴ pixels long: this is shown in red in the second series from the top (middle set). The process of taking $\mathrm{1}/\mathrm{8}$ and blowing up by a factor of 8 continues to the third line (length t=2³, resolution $\mathit{\tau }={\mathrm{2}}^{-\mathrm{7}}$), unit scale=2⁷ pixels (shown by the red arrows in the third series), until the bottom series which spans the range t=0 to t=1, and resolution $\mathit{\tau }={\mathrm{2}}^{-\mathrm{10}}$ with unit scale 2¹⁰ pixels (the whole series displayed). Each series was rescaled in the vertical so that its range between maximum and minimum was the same. The unit relaxation scales indicated by the red arrows mark the transition from small to large scales. Since the top series in each set has a unit scale of 2 (degraded), it is nearly a white noise (a, c, e) or (ordinary) Brownian motion (b, d, f). In contrast, the bottom series is exactly of length t=1, so that it is close to the fGn and fBm limits (left and right) with the standard exponent $H=h+\mathrm{1}/\mathrm{2}$. As indicated in the text, the second series from the top in the bottom set is most realistic for monthly temperature anomalies.

Figure 6 shows three simulations, each of length 2¹⁹ pixels, with each pixel corresponding to a temporal resolution of $\mathit{\tau }={\mathrm{2}}^{-\mathrm{10}}$, so that the unit (relaxation) scale is 2¹⁰ elementary pixels. Each simulation uses the same random seed, but they have h's increasing from $h=\mathrm{1}/\mathrm{10}$ (top set) to $h=\mathrm{5}/\mathrm{10}$ (bottom set). The fRm on the right is from the running sum of the fRn on the left. Each series has been rescaled, so that the range (maximum–minimum) is the same for each. Starting at the top line of each group, we show 2¹⁰ points of the original series degraded by a factor of 2⁹. The second line shows a blow-up by a factor of 8 of the part of the upper line to the right of the dashed vertical line. The line below is a further blow-up by a factor of 8 until the bottom line shows a $\mathrm{1}/\mathrm{512}$ part of the full simulation but at full resolution. The unit scale indicating the transition from small to large is shown by the horizontal red line in the middle-right figure. At the top (degraded by a factor of 2⁹), the unit (relaxation) scale is 2 pixels, so that the top line degraded view of the simulation is nearly a white noise (left) or (ordinary) Brownian motion (right). In contrast, the bottom series is exactly of length unity, so that it is close to the fGn limit with the standard exponent $H=h+\mathrm{1}/\mathrm{2}$. Moving from bottom to top in Fig. 6, one effectively transitions from fGn to fRn (left column) and from fBm to fRm (right column).

If we take the empirical relaxation scale for the global temperature to be 2⁷ months (≈10 years, Lovejoy et al., 2017) and we use monthly-resolution temperature anomaly data, then the nondimensional resolution is 2⁻⁷, corresponding to the second series from the top (which is thus 2¹⁰ months ≈80 years long). Since $h\approx \mathrm{0.38}\pm \mathrm{0.03}$ (Procyk et al., 2022), the second series from the top in the bottom set is the most realistic, and we can make out the low-frequency undulations that are mostly present at scales $\mathrm{1}/\mathrm{8}$ of the series (or less).

Figure 7 shows realizations constructed from the same random seed but for the extended range $\mathrm{1}/\mathrm{2}<h<\mathrm{2}$ (i.e. beyond fGn). Over this range, the top (large-scale, degraded-resolution) series are close to white noises (left) and Brownian motions (right). For the bottom series, there is no equivalent fGn or fBm process and the curves become smoother, although the rescaling may hide this somewhat (see for example the $h=\mathrm{13}/\mathrm{20}$ set, the blow-up of the far-right $\mathrm{1}/\mathrm{8}$ of the second series from the top shown in the third line). For $\mathrm{1}<h<\mathrm{2}$, also note the oscillations with frequency $\mathrm{2}\mathit{\pi }/\sin (\mathit{\pi }/h)$ (Eqs. 53 and A3): this is the fractional oscillation range.

Figure 7 The same as Fig. 6 but for $h=\mathrm{7}/\mathrm{10}$, $\mathrm{13}/\mathrm{10}$ and $\mathrm{19}/\mathrm{10}$ (top to bottom). Over this range, the top (large-scale, degraded-resolution) series is close to a white noise (a, c, e) and Brownian motion (b, d, f). For the bottom series, there is no equivalent fGn or fBm process and the curves become smoother, although the rescaling may hide this somewhat (see for example the middle $h=\mathrm{13}/\mathrm{20}$ set, the blow-up of the far-right $\mathrm{1}/\mathrm{8}$ of the second series from the top shown in the third line). Also note for the bottom two sets with $\mathrm{1}<h<\mathrm{2}$ the oscillations that have frequency $\mathrm{2}\mathit{\pi }/\sin (\mathit{\pi }/h)$: this is the fractional oscillation range.

Figure 8 shows simulations similar to Fig. 5a (fRn on the left, fRm on the right), except that instead of making a large simulation and then degrading and zooming, all the simulations were of equal length (2¹⁰ points), but the relaxation scale was changed from 2¹⁵ pixels (bottom) to 2¹⁰, 2⁵ and 1 pixel (top). Again, the top is white noise (left) and Brownian motion (right), and the bottom is (nearly) fGn (left) and fBm (right); Fig. 9 shows the extensions to $\mathrm{1}/\mathrm{2}<h<\mathrm{2}$.

Figure 8 This set of simulations is similar to Fig. 6 ((a, c, e) fRn, (b, d, f) fRm) except that instead of making a large simulation and then degrading and zooming, all the simulations were of equal length (2¹⁰ points) but resolutions $\mathit{\tau }={\mathrm{2}}^{-\mathrm{15}}$, 2⁻¹⁰, 2⁻⁵, and 1 (bottom to top). The simulations therefore spanned the ranges of scale 2⁻¹⁵ to 2⁻⁵, 2⁻¹⁰ to 1, 2⁻⁵ to 2⁵, and 1 to 2¹⁰, and the same random seed was used in each so that we can see how the structures slowly change when the relaxation scale changes. The bottom fRn, $h=\mathrm{5}/\mathrm{10}$ set is the closest to that observed for the Earth's temperature, and since the relaxation scale is of the order of a few years, the second series from the top of this set (with 1 pixel = 1 month) is close to that of monthly global temperature anomaly series. In that case the relaxation scale would be 32 months and the entire series would be ${\mathrm{2}}^{\mathrm{10}}/\mathrm{12}\approx \mathrm{85}$ years long. The top series (of total length 2¹⁰ relaxation times) is (nearly) a white noise (a, c, e) and Brownian motion (b, d, f), and the bottom is (nearly) an fGn (a, c, e) and fBm (b, d, f). The total range of scales covered here (2¹⁰×2¹⁵) is larger than in Fig. 5a and allows one to more clearly distinguish the high- and low-frequency regimes.

Figure 9 The same as Fig. 8 but for larger h values; see also Fig. 7.

4 Prediction

The initial value for Weyl fractional differential equations is effectively at $t=-\mathrm{\infty }$, so that for fRn, it is not directly relevant at finite times (although the ensemble mean is assumed=0; for fRm, the initial condition $Q_{\mathit{\alpha },h}\left(\mathrm{0}\right)=\mathrm{0}$ is important). The prediction problem is thus to use past data (say, for t<0) in order to make the most skillful prediction for t>0. We are therefore dealing with a past value rather than usual initial value problem. The emphasis on past values is particularly appropriate since in the fGn limit, the memory is so large that values of the series in the distant past are important. Indeed, prediction of fGn with a finite length of past data involves placing strong (mathematically singular) weight on the most ancient data available (see Gripenberg and Norros, 1996; Del Rio Amador and Lovejoy, 2019, 2021a, b). This is quite different from standard stochastic predictions that are based on short-memory (exponential) auto-regressive or moving-average-type processes that are not much different from initial value problems.

To deal with the small-scale divergences when $\mathrm{0}<h+\mathit{\alpha }\le \mathrm{1}/\mathrm{2}$, it is necessary to predict the finite-resolution fRn: $U_{\mathit{\alpha },h,\mathit{\tau }}\left(t\right)$. Using Eq. (40) for $U_{\mathit{\alpha },h,\mathit{\tau }}\left(t\right)$, we have

$$\begin{array}{}\text{(62)}& \begin{array}{rl}{U}_{\mathit{\alpha },h,\mathit{\tau }}\left(t\right)=& {\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\left[\underset{-\mathrm{\infty }}{\overset{t}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}\right(t-v\left)\mathit{\gamma }\right(v)\mathrm{d}v\\ & -\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(-v)\mathit{\gamma }\left(v\right)\mathrm{d}v]\\ & -{\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\left[\underset{-\mathrm{\infty }}{\overset{t-\mathit{\tau }}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}\right(t-\mathit{\tau }-v\left)\mathit{\gamma }\right(v)\mathrm{d}v\\ & -\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(-v)\mathit{\gamma }\left(v\right)\mathrm{d}v]\\ =& {\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\left[\underset{-\mathrm{\infty }}{\overset{t}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}\right(t-v\left)\mathit{\gamma }\right(v)\mathrm{d}v\\ & -\underset{-\mathrm{\infty }}{\overset{t-\mathit{\tau }}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(t-\mathit{\tau }-v)\mathit{\gamma }\left(v\right)\mathrm{d}v].\end{array}\end{array}$$

Now define the predictor for t≥0 (indicated by a circumflex):

$$\begin{array}{}\text{(63)}& \begin{array}{rl}\widehat{U_{\mathit{\alpha },h,\mathit{\tau }}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}& \left[\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(t-v)\mathit{\gamma }\left(v\right)\mathrm{d}v\right.\\ & \left.-\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(t-\mathit{\tau }-v)\mathit{\gamma }\left(v\right)\mathrm{d}v\right].\end{array}\end{array}$$

To show that it is indeed the optimal predictor, consider the predictor error E_τ(t):

$$\begin{array}{ll}{\displaystyle }{E}_{\mathit{\tau }}\left(t\right)=& {\displaystyle }{U}_{\mathit{\alpha },h,\mathit{\tau }}\left(t\right)-\widehat{U_{\mathit{\alpha },h,\mathit{\tau }}}\left(t\right)\\ {\displaystyle }=& {\displaystyle }{\mathit{\tau }}^{-\mathrm{1}}\left[\underset{-\mathrm{\infty }}{\overset{t}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(t-v)\mathit{\gamma }\left(v\right)\mathrm{d}v\right.\\ {\displaystyle }& {\displaystyle }\left.-\underset{-\mathrm{\infty }}{\overset{t-\mathit{\tau }}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(t-\mathit{\tau }-v)\mathit{\gamma }\left(v\right)\mathrm{d}v\right]\\ {\displaystyle }& {\displaystyle }-{\mathit{\tau }}^{-\mathrm{1}}\left[\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(t-v)\mathit{\gamma }\left(v\right)\mathrm{d}v\right.\\ {\displaystyle }& {\displaystyle }\left.-\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(t-\mathit{\tau }-v)\mathit{\gamma }\left(v\right)\mathrm{d}v\right]\\ {\displaystyle }=& {\displaystyle }{\mathit{\tau }}^{-\mathrm{1}}\left[\underset{\mathrm{0}}{\overset{t}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(t-v)\mathit{\gamma }\left(v\right)\mathrm{d}v\right.\\ \text{(64)}& {\displaystyle }& {\displaystyle }\left.-\underset{\mathrm{0}}{\overset{t-\mathit{\tau }}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(t-\mathit{\tau }-v)\mathit{\gamma }\left(v\right)\mathrm{d}v\right].\end{array}$$

Equation (64) shows that the error depends only on γ(v) for v>0, whereas the predictor (Eq. 63) only depends on γ(v) for v<0, and hence they are orthogonal:

$$\begin{array}{}\text{(65)}& \langle E_{\mathit{\tau }}\left(t\right)\widehat{U_{\mathit{\alpha },h,\mathit{\tau }}}\left(t\right)\rangle =\mathrm{0}.\end{array}$$

This is a sufficient condition for $\widehat{U_{\mathit{\alpha },h,\mathit{\tau }}}\left(t\right)$ to be the minimum square predictor, which is the optimal predictor for stationary Gaussian processes (e.g. Papoulis, 1965). The prediction error variance is

$$\begin{array}{}\text{(66)}& \begin{array}{rl}& \langle E_{\mathit{\tau }}(t{)}^{\mathrm{2}}\rangle \\ & ={\mathit{\tau }}^{-\mathrm{2}}\left[\underset{\mathrm{0}}{\overset{t-\mathit{\tau }}{\int }}\left(G_{\mathrm{1}+\mathit{\alpha },h}\right(t-v)-G_{\mathrm{1}+\mathit{\alpha },h}(t-\mathit{\tau }-v){)}^{\mathrm{2}}\mathrm{d}v\right.\\ & \left.+\underset{t-\mathit{\tau }}{\overset{t}{\int }}{G}_{\mathrm{1}+\mathit{\alpha },h}(t-v{)}^{\mathrm{2}}\mathrm{d}v\right],\end{array}\end{array}$$

or, with a change in variables,

$$\begin{array}{ll}{\displaystyle }\langle E_{\mathit{\tau }}(t{)}^{\mathrm{2}}\rangle =& {\displaystyle }{\mathit{\tau }}^{-\mathrm{2}}{V}_{\mathit{\alpha },h}\left(\mathit{\tau }\right)\\ \text{(67)}& {\displaystyle }& {\displaystyle }-{\mathit{\tau }}^{-\mathrm{2}}\left[\underset{t-\mathit{\tau }}{\overset{\mathrm{\infty }}{\int }}\left(G_{\mathrm{1}+\mathit{\alpha },h}\right(v+\mathit{\tau })-G_{\mathrm{1}+\mathit{\alpha },h}(v){)}^{\mathrm{2}}\mathrm{d}v\right],\end{array}$$

where we have used $\langle U_{\mathit{\alpha },h,\mathit{\tau }}^{\mathrm{2}}\rangle ={\mathit{\tau }}^{-\mathrm{2}}{V}_{\mathit{\alpha },h}\left(\mathit{\tau }\right)$ (the unconditional variance).

There are numerous skill indicators, but the most popular and easy-to-interpret definition of forecast skill is the minimum square skill score or MSSS (S_k,τ, see Del Rio Amador and Lovejoy, 2021a, for a discussion of this and other indicators). For this, we obtain

$$\begin{array}{ll}{\displaystyle }{S}_{k,\mathit{\tau }}\left(t\right)=& {\displaystyle }\mathrm{1}-{\displaystyle \frac{〈E_{\mathit{\tau }}(t{)}^{\mathrm{2}}〉}{〈E_{\mathit{\tau }}(\mathrm{\infty }{)}^{\mathrm{2}}〉}}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{{\int }_{t-\mathit{\tau }}^{\mathrm{\infty }}\left(G_{\mathrm{1}+\mathit{\alpha },h}\right(u+\mathit{\tau })-G_{\mathrm{1}+\mathit{\alpha },h}(u){)}^{\mathrm{2}}\mathrm{d}u}{V_{\mathit{\alpha },h}\left(\mathit{\tau }\right)}}\\ \text{(68)}& {\displaystyle }=& {\displaystyle }{\displaystyle \frac{{\int }_{t-\mathit{\tau }}^{\mathrm{\infty }}\left(G_{\mathrm{1}+\mathit{\alpha },h}\right(v+\mathit{\tau })-G_{\mathrm{1}+\mathit{\alpha },h}(v){)}^{\mathrm{2}}\mathrm{d}v}{\left[\begin{array}{c}{\int }_{\mathrm{0}}^{\mathrm{\infty }}\left(G_{\mathrm{1}+\mathit{\alpha },h}\right(v+\mathit{\tau })-G_{\mathrm{1}+\mathit{\alpha },h}(v){)}^{\mathrm{2}}\mathrm{d}v\\ +{\int }_{\mathrm{0}}^{\mathit{\tau }}{G}_{\mathrm{1}+\mathit{\alpha },h}(v{)}^{\mathrm{2}}\mathrm{d}v\end{array}\right]}}.\end{array}$$

When $h<\mathrm{1}/\mathrm{2}$ and

$$G_{\mathrm{1},h}\left(t\right)=G_{\mathrm{1},h}^{\left(\text{fGn}\right)}\left(t\right)={\displaystyle \frac{t^h}{\mathrm{\Gamma }(\mathrm{1}+h)}},$$

we obtain the fGn result:

$$\begin{array}{ll}{\displaystyle }{S}_k& {\displaystyle }={\displaystyle \frac{{\mathit{\xi }}_h\left(\mathrm{\infty }\right)-{\mathit{\xi }}_h\left(\mathit{\lambda }\right)}{{\mathit{\xi }}_h\left(\mathrm{\infty }\right)+\frac{\mathrm{1}}{\mathrm{2}h+\mathrm{1}}}},\\ \text{(69)}& {\displaystyle }{\mathit{\xi }}_h\left(\mathit{\lambda }\right)& {\displaystyle }=\underset{\mathrm{0}}{\overset{\mathit{\lambda }-\mathrm{1}}{\int }}\left(\right(v+\mathrm{1}{)}^h-v^h{)}^{\mathrm{2}}\mathrm{d}v\end{array}$$

(Lovejoy et al., 2015), where λ is the forecast horizon (lead time) measured in the number of time steps in the future (due to the fGn scaling, it is independent of the resolution τ). The MSSS gives the fraction of the variance explained by the optimum predictor; when skill=1, the forecast is perfect.

To survey the implications, let us start by showing the τ independent results for fGn, shown in Fig. 10, which is a variant on a plot published in Lovejoy et al. (2015). We see that when $h\approx \mathrm{1}/\mathrm{2}$ (H≈1), the skill is very high; indeed, in the limit $h\to \mathrm{1}/\mathrm{2}$, we have perfect skill for fGn forecasts (this would of course require an infinite amount of past data to attain).

Figure 10 The prediction skill (S_k) for pure fGn processes for forecast horizons up to λ=10 steps (10 times the resolution). This plot is nondimensional, and it is valid for time steps of any duration. From bottom to top, the curves correspond to $h=\mathrm{1}/\mathrm{20}$, $\mathrm{3}/\mathrm{10}$, … $\mathrm{9}/\mathrm{20}$ (red, top, close to the empirical h).

Now consider the fRn skill: we will start by considering the pure (α=0) fRn case where the memory comes completely from the (high-frequency) storage, anticipating that the fGn forced case (α≠0) obtains its memory and skill from both storage and forcing. In comparison with fGn, fRn has an extra parameter, the resolution of the data, τ. Figure 11 shows curves corresponding to Fig. 10 for fRn with forecast horizon integer multiples (λ) of τ, i.e. for times t=λτ in the future but with separate curves, one for each of five τ values increasing from 10⁻⁴ to 10 by factors of 10. When τ is small, the results should be close to those of fGn, i.e. with potentially high skill, and in all cases, the skill is expected to vanish quite rapidly for τ>1 since in this limit, fRn becomes an (unpredictable) white noise (although there are scaling corrections to this).

Figure 11 Panels (a, c, e) show the skill (S_k) of pure (α=0) fRn forecasts (as in Fig. 10 for fGn) for fRn skill with $h=\mathrm{1}/\mathrm{20}$, $\mathrm{5}/\mathrm{20}$, and $\mathrm{9}/\mathrm{20}$ (top to bottom sets); λ is the forecast horizon, the number of steps of resolution τ forecast into the future. Panels (b, d, f) show the ratio (r) of the fRn to the corresponding fGn skill. Here the result depends on τ; each curve is for different values increasing from 10⁻⁴ (top, black) to 10 (bottom, purple) and increasing by a factor of 10 (the red set in the bottom plots with $\mathit{\tau }={\mathrm{10}}^{-\mathrm{2}}$; $h=\mathrm{9}/\mathrm{20}$ are closest to the empirical values).

To better understand the fGn limit, it is helpful to plot the ratio of the fRn-to-fGn skill (Fig. 11, right column). We see even with quite small values $\mathit{\tau }={\mathrm{10}}^{-\mathrm{4}}$ (top, black curves) that some skill has already been lost. Figure 12 shows this more clearly: it shows 1-time-step and 10-time-step skill ratios. To put this into perspective, it is helpful to compare this using some of the parameters relevant to macroweather forecasting. According to Lovejoy et al. (2015) and Del Rio Amador and Lovejoy (2019), the relevant empirical Haar exponent is $\approx -\mathrm{0.1}$ for the global temperature, so that $h=\mathrm{1}/\mathrm{2}-\mathrm{0.1}\approx \mathrm{0.4}$. Although direct empirical estimates of the relaxation time are difficult since the responses to anthropogenic forcing begin to dominate over the internal variability after ≈10 years, Procyk et al. (2022) have used the deterministic response to estimate a global relaxation time of ≈5 years (work in progress using maximum likelihood estimates shows that at scales of hundreds of kilometres, it is quite variable, ranging from months to decades; Procyk, 2021). For monthly-resolution forecasts, the nondimensional resolution is $\mathit{\tau }\approx \mathrm{1}/\mathrm{100}$. With these values, we see (red curves) that we may have lost ≈30 % of the fGn skill for 1-month forecasts and ≈85 % for 10-month forecasts. Comparing this with Fig. 10, we see that this implies about 60 % and 10 % skill (see also the red curve in Fig. 11, bottom set).

Figure 12 The ratio of (α=0) fRn skill to fGn skill (a: 1-step horizon, b: 10-step forecast horizon) as a function of resolution τ for h increasing from (at left) bottom to top ($h=\mathrm{1}/\mathrm{20}$, $\mathrm{2}/\mathrm{20}$, $\mathrm{3}/\mathrm{20}$ … $\mathrm{9}/\mathrm{20}$); the $h=\mathrm{9}/\mathrm{20}$ curve (close to the empirical value) is the curve that starts at the upper left of each plot.

Going beyond the $\mathrm{0}<h<\mathrm{1}/\mathrm{2}$ region that overlaps fGn, Figs. 12 and 13 clearly show that the skill continues to increase with h. We already saw (Fig. 4) that the range $\mathrm{1}/\mathrm{2}<h<\mathrm{3}/\mathrm{2}$ has rms Haar fluctuations that for Δt<0 mimic fBm, and these do indeed have higher skill, approaching unity for h near 1 corresponding to a Haar exponent $\approx \mathrm{1}/\mathrm{2}$, i.e. close to an fBm with $H=\mathrm{1}/\mathrm{2}$, i.e. a regular Brownian motion. Recall that for Brownian motion, the increments are unpredictable but the process itself is predictable (persistence). In Fig. 12, we show the skill for various h's as a function of resolution τ. Figure 14 shows that for $h<\mathrm{3}/\mathrm{2}$, the skill decreases rapidly for τ>1. Figure 15 in the fractional oscillation equation regime shows that the skill oscillates.

Figure 13 The 1-step (a) and 10-step (b) pure (α=0) fRn forecast skill as a function of h for various resolutions (τ) ranging from $\mathit{\tau }={\mathrm{10}}^{-\mathrm{4}}$ (black, left of each set) through $\mathit{\tau }={\mathrm{10}}^{-\mathrm{3}}$ (brown) 10⁻² (red), 0.1 (blue), 1 (orange), and 10 (purple). In the right set τ=1 (orange), 10 (purple) lines are nearly on top of the S_k=0 line. Again, red ($\mathit{\tau }={\mathrm{10}}^{-\mathrm{2}}$) is the more empirically relevant value for monthly data. Recall that the regime $h<\mathrm{1}/\mathrm{2}$ (to the left of the vertical dashed lines) corresponds to the overlap with fGn.

Figure 14 One-step pure (α=0) fRn prediction skills as a function of resolution for hs increasing from $\mathrm{1}/\mathrm{20}$ (bottom) to $\mathrm{29}/\mathrm{20}$ (top) every $\mathrm{1}/\mathrm{10}$. Note the rapid transition to low skill (white noise) for τ>1. The curve for $h=\mathrm{9}/\mathrm{20}$ is shown in red.

Figure 15 Same as Fig. 14 except for $h=\mathrm{37}/\mathrm{20}$ and $\mathrm{39}/\mathrm{20}$ showing the 1-step skill (black) and 10-step skill (dashed). The right-hand dashed and right-hand solid lines are for $h=\mathrm{39}/\mathrm{20}$: they clearly show that the skill oscillates in this fractional oscillation equation regime. The corresponding left lines are for $h=\mathrm{37}/\mathrm{20}$.

We may now consider the skill of the fGn-forced process (α≠0) in Fig. 16. For small τ, short lags, λ (the upper left), the contours are fairly linear along lines of constant h+α, so that, as expected, the predictability is essentially that of an fGn process but with an effective exponent h+α. At the opposite extreme, (large τ, h) the lines are fairly horizontal, indicating that the skill from the storage (i.e. from h) is negligible and that all the memory (and hence skill) comes from the forcing fGn, exponent α. The in-between resolutions and lags generally have in-between slopes. As expected, the skill from the storage drops off quickly for resolutions $\approx >\mathit{\tau }$. For h≥1, there is some waviness in the contours due to the oscillatory nature of the Green functions.

Figure 16 Contour plots of the forecast skill, with h along the horizontal axis and α along the vertical axis. The plots are for increasing nondimensional resolutions: τ=0.001, 0.01, 0.1, 1, and 10 (top to bottom), with forecasts for lags λ=1, 3, and 10 (left to right) and with contour levels (legend) varying from nearly no skill (0.03) to nearly full skill (0.98).

5 Conclusions

Ever since Budyko (1969) and Sellers (1969), the energy balance between the Earth and outer space has been modelled by the energy balance equation (EBE) based on the continuum heat equation; see North and Kim (2017) for a recent review and see Ziegler and Rehfeld (2020) for a recent regional application. It is most commonly used as a model for the globally averaged temperature, where it is usually derived by applying Newton's law of cooling applied to a uniform slab of material, a “box”. The resulting EBE is a first-order relaxation equation describing the exponential relaxation of the temperature to a new equilibrium after it has been perturbed by an external forcing. Its first-order (h=1) derivative term accounts for energy storage.

The resulting model relaxes to equilibrium much too quickly, so that to increase realism, it is usual to introduce a few interacting slabs (representing for example the atmosphere and ocean mixed layer; the Intergovernmental Panel on Climate Change recommends two such components; IPCC, 2013). However, it turns out that these h=1 box models do not use the correct surface radiative–conductive boundary conditions. If one assumes heat transport by the classical heat equation and radiative–conductive boundary conditions are used instead, one instead obtains the half-order EBE, the HEBE with $h=\mathrm{1}/\mathrm{2}$ (Lovejoy, 2021a, b), which is already close to the global empirical value ($h=\mathrm{0.38}\pm \mathrm{0.03}$, Procyk et al., 2022; Del Rio Amador and Lovejoy, 2019; see also Lovejoy et al., 2015). However, this model is only valid in the macroweather regime – for timescales of weeks and longer and, due to the spatial scaling in the atmosphere, the fractional heat equation (FHE) may be a more appropriate model than the classical one. The use of the FHE can be justified by recognizing that a realistic energy transport model involves a continuous hierarchy of mechanisms. The extension to the FHE leads directly to a fractional relaxation equation that generalizes the EBE: the fractional energy balance equation (Lovejoy, 2021a, b) (FEBE). The FEBE can also be derived phenomenologically by assuming that energy storage processes are scaling (Lovejoy, 2019; Lovejoy et al., 2021).

When forced by a Gaussian white noise, the FEBE is also a generalization of fractional Gaussian noise (fGn), and its integral (fractional relaxation motion, fRm) generalizes fractional Brownian motion (fBm). More classically, it generalizes the Orenstein–Uhlenbeck process that corresponds to the h=1 special case (i.e. the standard EBE with white noise forcing). Over the parameter range $\mathrm{0}<h<\mathrm{1}/\mathrm{2}$, the high-frequency FEBE limit (fGn) has been used as the basis of monthly and seasonal temperature forecasts (Lovejoy et al., 2015; Del Rio Amador and Lovejoy, 2019, 2021a, b); at 1-month lead times, these macroweather forecasts are similar in skill to conventional numerical models, whereas for bimonthly, seasonal and annual forecasts, they are more skillful (Del Rio Amador and Lovejoy, 2021a). For multidecadal timescales the low-frequency limit has been used as the basis of climate projections through to the year 2100 (Hébert, 2017; Lovejoy et al., 2017; Hébert et al., 2021), and more recently, the full FEBE has been used directly (Procyk et al., 2020, 2022; Procyk, 2021).

It was the success of predictions and projections with different exponents but the same theoretically derived empirical underlying FEBE h≈0.4 that, over recent years, motivated the development of the FEBE (announced in Lovejoy, 2019) and the work reported here. The statistical characterizations, correlations, structure functions, Haar fluctuations and spectra as well as the predictability properties are important for these and other FEBE applications and are derived in this paper.

While the deterministic fractional relaxation equation is classical, various technical difficulties arise when it is generalized to the stochastic case: in the physics literature, it is a fractional Langevin equation (FLE) that has almost exclusively been considered a model of diffusion of particles starting at an origin. This requires t=0 initial conditions that imply that the solutions are strongly nonstationary. In comparison, the Earth's temperature fluctuations that are associated with its internal variability are statistically stationary. This can easily be modelled with initial conditions at $t=-\mathrm{\infty }$, i.e. by using Weyl fractional derivatives. In addition, in the usual FLE, the highest-order derivative is an integer, so that sample processes are rms differentiable of order at least 1 (Watkins et al., 2020, have called the FEBE a “Fractionally Integrated FLE”). In the FEBE and the fractionally integrated extensions, the highest-order derivative is readily of order $<\mathrm{1}/\mathrm{2}$, so that sample processes are generalized functions (“noises”) and must be smoothed/averaged for physical applications.

Although EBEs were originally developed to understand the deterministic temperature response to external forcing, the temperature also responds to stochastic “internal” forcing. While the Earth's system variability is generally highly non-Gaussian (multifractal, Lovejoy, 2018), the temporal macroweather regime modelled here is the quasi-Gaussian exception. This paper therefore explores the statistics of the temperature response when it is stochastically forced by Gaussian processes, both by white noise (α=0) and by a (long-memory) fractional Gaussian noise (fGn) process. The white noise special case – “pure fRn and fRm” – is the α=0 special case; the fGn-forced case extends the parameter range to $\mathrm{0}\le \mathit{\alpha }<\mathrm{1}/\mathrm{2}$. According to work in progress using satellite and reanalysis radiances, both cases appear to be empirically relevant for modelling the Earth's energy balance.

A key novelty is therefore to consider the fractional relaxation equation (a FLE) forced by white and scaling noises starting from $t=-\mathrm{\infty }$, equivalent to Weyl's “fractionally integrated fractional relaxation equation”. In addition, the highest-order terms in standard FLEs are integer-ordered: the fractional terms represent damping and are of lower order, guaranteeing that solutions are regular functions. However, the FEBE's highest-order term is fractional, and over the main empirically significant parameter range ($\mathit{\alpha }+h<\mathrm{1}/\mathrm{2}$) the processes are noises (generalized functions): in order to represent physical processes, they must be averaged. This is conveniently handled by introducing their integrals or “motions”. We proceeded to derive their fundamental statistical properties, including series expansions about the origin and infinity. These expansions are nontrivial since they mix fractional and integer-ordered terms (Appendix A). Since the FEBE is used as the basis for macroweather predictions, the theoretical predictability skill is important in applications and was also derived.

With these stationary Gaussian forcings, the solutions are a new stationary process – fRn (α=0) and its extensions to fractionally integrated fRn processes (α>0). Over the range $\mathrm{0}<\mathit{\alpha }+h<\mathrm{1}/\mathrm{2}$, we show that the small-scale limit is an fGn, and its integral – fRm – has stationary increments and generalizes fBm. Although at long enough times the fRn (α=0) tends to a Gaussian white noise and fRm to a standard Brownian motion, this long time convergence is typically very slow (when α>0, the long time behaviours are fGn and fBm processes, parameter α).

Much of the effort was in deducing the asymptotic small- and large-scale behaviours of the autocorrelation functions that determine the statistics and in verifying these with extensive numerical simulations. An interesting exception was the $h=\mathrm{1}/\mathrm{2}$ special case, which for fGn corresponds to an exactly $\mathrm{1}/f$ noise. Here, we give the exact mathematical expressions for the full correlation functions, showing that they had logarithmic dependencies at both small and large scales. The resulting HEBE has an exceptionally slow transition from small to large scales (a factor of a million or more is needed), and empirically it is quite close to the global temperature series over scales of months, decades and possibly longer.

Beyond improved monthly and seasonal temperature forecasts and multidecadal projections, the stochastic FEBE opens up several paths for future research. One of the more promising is to apply these techniques to the spatial FEBE and generalize it in various directions. This is a follow-up on the special value $h=\mathrm{1}/\mathrm{2}$ that is very close to that found empirically and that can be analytically deduced from the classical Budyko–Sellers energy transport equation by improving the mathematical treatment of the radiative boundary conditions (Lovejoy, 2021a, b). In the latter case, one obtains a partial fractional differential equation for the horizontal space–time variability of temperature anomalies over the Earth's surface, allowing regional forecasts and projections. This has already allowed improved regional projections (Procyk, 2021) and promises better monthly and seasonal forecasts.

While the FEBE has already demonstrated its ability to project future climates, these improvements will allow for the modelling of the nonlinear albedo–temperature feedbacks needed for modelling of transitions between different past climates. Finally, FEBE-based projections have shown that, in spite of improved computer power and algorithms, conventional GCM approaches may be suffering from diminishing returns; the GCMs in the latest IPCC assessment (AR6, 2021) are even more uncertain: a range of 2–5.5 K$/$CO₂ doubling (90 % confidence) those in the previous assessment (AR5, 2013, 1.5–4.5 K per doubling) while also being somewhat warmer. The FEBE had the somewhat lower but much less uncertain range 1.6–2.4 K$/$CO₂ doubling (90 % confidence). Conventional GCM approaches attempt to explicitly model as many degrees of freedom as possible, and by the year 2030 they are expected to have kilometric-scale (“cloud-resolving”) resolutions that will model structures that live for only 15 min and then average them over decades. The FEBE (with regional and other future extensions) is, in contrast, a high-level stochastic model that accounts for the collective interactions of huge numbers of degrees of freedom (Lovejoy, 2019). It is thus a promising candidate for a new generation of climate models.

Appendix A: The small- and large-scale fRn and fRm statistics

A1 R_α,h(t) as a Laplace transform

In Sect. 2.4, we derived general statistical formulae for the autocorrelation functions of motions and noises defined in terms of Green's functions of fractional operators. Since the processes are Gaussian, autocorrelations fully determine the statistics. While the autocorrelations of fBm and fGn are well known, those for fRm and fRn are new and are not so easy to deal with since they involve quadratic integrals of Mittag–Leffler functions. In this Appendix, we derive the basic power law expansions as well as large t (asymptotic) expansions, and we numerically investigate their accuracy.

It is simplest to start with the Fourier expression for the autocorrelation function for the unit white noise forcing (Eq. 33). First convert the inverse Fourier transform (Eq. 66) into a Laplace transform. For this, consider the integral over the contour C in the complex plane:

$$\begin{array}{}\text{(A1)}& I_C\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{C}{\int }{\displaystyle \frac{e^{zt}}{z^{\mathit{\alpha }}(-z{)}^{\mathit{\alpha }}\left(\mathrm{1}+z^h\right)\left(\mathrm{1}+(-z{)}^h\right)}}\mathrm{d}z.\end{array}$$

Take C to be the closed contour obtained by integrating along the imaginary axis (this part gives R_α,h(t), Eq. 33) and closing the contour along an (infinite) semicircle over the second and third quadrants. When $\mathrm{0}<h<\mathrm{1}$, there are no poles in these quadrants, but we must integrate around a branch cut on the negative real axis. When $\mathrm{1}<h<\mathrm{2}$, we must take into account two new branch cuts and two new poles in the negative real half-plane. In a polar representation z=re^iθ, the additional branch cuts are along the rays $z=r{e}^{\pm i\mathit{\pi }/h}$, r>1, circling around the poles at $z=e^{\pm i\mathit{\pi }/h}$. The additional branch cuts give no net contribution, but the residues of the poles do make a contribution ($P_{\mathit{\alpha },h}\ne \mathrm{0}$ below). We can express both cases with the formula

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{R}_{\mathit{\alpha },h}\left(t\right)=& -{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\text{Im}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{e^{-xt}\mathrm{d}x}{x^{\mathrm{2}\mathit{\alpha }}{e}^{i\mathit{\alpha }\mathit{\pi }}(\mathrm{1}+x^h)(\mathrm{1}+x^h{e}^{i\mathit{\pi }h})}}\\ & +P_{\mathit{\alpha },h,+}\left(t\right);t>\mathrm{0}\mathrm{0}\le \mathit{\alpha }<\mathrm{1}/\mathrm{2}.\end{array}\end{array}$$

“Im” indicates the imaginary part and

$$\begin{array}{ll}{\displaystyle }{P}_{\mathit{\alpha },h,\pm }\left(t\right)& {\displaystyle }=\mathrm{0},\mathrm{0}<h<\mathrm{1},\\ {\displaystyle }{P}_{\mathit{\alpha },h,\pm }\left(t\right)& {\displaystyle }=-e^{t\cos \left(\frac{\mathit{\pi }}{h}\right)}{\displaystyle \frac{\sin (\pm \frac{\mathit{\pi }}{h}(\mathrm{1}-\mathit{\alpha })+\frac{h\mathit{\pi }}{\mathrm{2}}+t\sin (\frac{\mathit{\pi }}{h}\left)\right)}{h\sin \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)}},\\ \text{(A3)}& {\displaystyle }& {\displaystyle }\mathrm{1}<h<\mathrm{2},\end{array}$$

While the integral term is monotonic, the P_α,h term oscillates with frequency $\mathit{\omega }=\mathrm{2}\mathit{\pi }/\sin (\mathit{\pi }/h)$. P_α,h accounts for the oscillations visible in Figs. 3, 4, and 7, although since when $\mathrm{1}<h<\mathrm{2}$, $\cos (\mathit{\pi }/h)<\mathrm{1}$, they decay exponentially. When h>1, this pole contribution dominates R_α,h(t) for a wide range of t values around t=1, although as we see below, eventually at large t, power law terms come to the fore.

Comments

• a.

  When α=0, h=1, and we obtain the classical Ornstein–Uhlenbeck autocorrelation: $R_{\mathrm{0},\mathrm{1}}\left(t\right)=\frac{\mathrm{1}}{\mathrm{2}}{e}^{-\left|t\right|}$.

• b.

  In the case h=0, the process reduces to an fGn process: $R_{\mathit{\alpha },\mathrm{0}}\left(t\right)=t^{-\mathrm{1}+\mathrm{2}\mathit{\alpha }}\mathrm{\Gamma }(\mathrm{1}-\mathrm{2}\mathit{\alpha })\sin \left(\mathit{\pi }\mathit{\alpha }\right)/\left(\mathrm{4}\mathit{\pi }\right)$. There is an extra factor of 4 that comes from the small h limit ${}_{-\mathrm{\infty }}{D}_t^h+\mathrm{1}\to \mathrm{2}$.

A2 Asymptotic expansions

An advantage of writing R_α,h(t) as a Laplace transform is that we can use Watson's lemma to obtain an asymptotic expansion (e.g. Bender and Orszag, 1978). The idea is that an expansion of Eq. (A2) around x=0 can be Laplace-transformed term by term to yield an asymptotic expansion for large t.

The expansion of the integrand around x=0 can be obtained from a binomial expansion (see also Eq. A10):

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{x^{\mathrm{2}\mathit{\alpha }}{e}^{i\mathit{\pi }\mathit{\alpha }}\left(\mathrm{1}+x^h\right)\left(\mathrm{1}+x^h{e}^{i\mathit{\pi }h}\right)}}\\ & ={\displaystyle \frac{e^{-i\mathit{\pi }\mathit{\alpha }}}{e^{i\mathit{\pi }h}-\mathrm{1}}}\sum _{n=\mathrm{0}}^{\mathrm{\infty }}(-\mathrm{1}{)}^n(e^{i(n+\mathrm{1})\mathit{\pi }h}-\mathrm{1})x^{-\mathrm{2}\mathit{\alpha }+nh},x<\mathrm{1}.\end{array}\end{array}$$

This leads to

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}-{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}& \text{Im}{\displaystyle \frac{\mathrm{1}}{x^{\mathrm{2}\mathit{\alpha }}{e}^{i\mathit{\alpha }\mathit{\pi }}(\mathrm{1}+x^h)(\mathrm{1}+x^h{e}^{hi\mathit{\pi }})}}\\ & =-\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{D}_{-n}{x}^{nh-\mathrm{2}\mathit{\alpha }},\end{array}\end{array}$$

$$\begin{array}{rl}{D}_n& =(-\mathrm{1}{)}^{n+\mathrm{1}}{\displaystyle \frac{\cos \left(\right(n-\frac{\mathrm{1}}{\mathrm{2}})\mathit{\pi }h+\mathit{\alpha }\mathit{\pi })-\cos (\frac{\mathit{\pi }h}{\mathrm{2}}+\mathit{\alpha }\mathit{\pi })}{\mathrm{2}\mathit{\pi }\sin \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)}}\\ & =(-\mathrm{1}{)}^n{\displaystyle \frac{\sin (\frac{n\mathit{\pi }h}{\mathrm{2}}+\mathit{\alpha }\mathit{\pi })\sin \left(\frac{(n-\mathrm{1})\mathit{\pi }h}{\mathrm{2}}\right)}{\mathit{\pi }\sin \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)}}\end{array}$$

(note that D_−n is used in the expansion here; D_n is used below).

Therefore, taking the term-by-term Laplace transform and using Watson's lemma,

$$\begin{array}{ll}{\displaystyle }{R}_{\mathit{\alpha },h}\left(t\right)=& {\displaystyle }-\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{D}_{-n}\mathrm{\Gamma }(\mathrm{1}+nh-\mathrm{2}\mathit{\alpha })t^{\mathrm{2}\mathit{\alpha }-(\mathrm{1}+nh)}\\ \text{(A6)}& {\displaystyle }& {\displaystyle }+P_{\mathit{\alpha },h,+}\left(t\right);t\gg \mathrm{1}(\mathrm{0}<\mathit{\alpha }<\mathrm{1}/\mathrm{2}).\end{array}$$

We have included the exponentially decaying residue $P_{\mathit{\alpha },h,+}$ that contributes when $\mathrm{1}<h<\mathrm{2}$. Note that although Γ diverges for all negative integer arguments, using the identity $\mathrm{\Gamma }(\mathrm{1}+hn-\mathrm{2}\mathit{\alpha })\sin \left(\right(nh-\mathrm{2}\mathit{\alpha }\left)\mathit{\pi }\right)=-\mathit{\pi }/\mathrm{\Gamma }(\mathrm{2}\mathit{\alpha }-nh)$, we see that the product $\sin \left(\right(nh-\mathrm{2}\mathit{\alpha }\left)\mathit{\pi }\right)\mathrm{\Gamma }(\mathrm{2}\mathit{\alpha }-nh)$ is finite.

The first terms are explicitly

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}{R}_{\mathit{\alpha },h}\left(t\right)=& {\displaystyle \frac{\mathrm{\Gamma }(\mathrm{1}-\mathrm{2}\mathit{\alpha })\sin \left(\mathit{\pi }\mathit{\alpha }\right)}{\mathit{\pi }}}{t}^{\mathrm{2}\mathit{\alpha }-\mathrm{1}}\\ & -{\displaystyle \frac{\cos \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)}{\cos (\frac{\mathit{\pi }h}{\mathrm{2}}-\mathit{\pi }\mathit{\alpha })\mathrm{\Gamma }(\mathrm{2}\mathit{\alpha }-h)}}{t}^{\mathrm{2}\mathit{\alpha }-(\mathrm{1}+h)}+\mathrm{\dots },\\ & t\gg \mathrm{1}.\end{array}\end{array}$$

We see that when α≠0, D₀>0, so that, as expected, the leading behaviour has no h dependence: it is only due to the long-range correlations in the forcing. We obtain the fGn result t^2α−1. However, for the pure fRn case, α=0 and D₀=0, so that we obtain

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}{R}_{\mathrm{0},h}\left(t\right)=& \sum _{n=\mathrm{1}}^{\mathrm{\infty }}(-\mathrm{1}{)}^n{\displaystyle \frac{\mathrm{1}+\cot \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)\tan \left(\frac{n\mathit{\pi }h}{\mathrm{2}}\right)}{\mathrm{2}\mathrm{\Gamma }(-nh)}}{t}^{-(\mathrm{1}+nh)}\\ & +P_{\mathrm{0},h,+}\left(t\right),t\gg \mathrm{1};\end{array}\end{array}$$

i.e. the leading behaviour is $t^{-(\mathrm{1}+h)}$. Note that the leading n=1 coefficient reduces to $-\mathrm{1}/\mathrm{\Gamma }(-h)$ and that for $\mathrm{0}<h<\mathrm{1}$, $\mathrm{\Gamma }(-h)<\mathrm{0}$.

For the motions (fRm), we need the expansion of V_α,h(t); this can be obtained by integrating R_α,h twice (using Eq. 36):

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{V}_{\mathit{\alpha },h}\left(t\right)=& a_{\mathit{\alpha },h}t+b_{\mathit{\alpha },h}-\mathrm{2}\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{D}_{-n}\mathrm{\Gamma }(-\mathrm{1}+nh-\mathrm{2}\mathit{\alpha })\\ & \cdot t^{\mathrm{2}\mathit{\alpha }+\mathrm{1}-nh}+\mathrm{2}{P}_{\mathit{\alpha },h,-}\left(t\right);\\ & t\gg \mathrm{1}\mathrm{0}\le \mathit{\alpha }<\mathrm{1}/\mathrm{2},\end{array}\end{array}$$

where $P_{a,h-}$ is from the poles when $\mathrm{1}<h<\mathrm{2}$. Since the asymptotic expansion is not valid for t=0, we used the indefinite integrals of R_α,h, and hence there is a linear $a_{\mathit{\alpha },h}t+b_{\mathit{\alpha },h}$ term from the constants of integration. However, when α>0, the leading term is the t^2α+1 term from the fGn forcing, and in the pure fRn case (α=0), we can take ${lim}_{\mathit{\alpha }\to \mathrm{0}}(-\mathrm{2}{D}_{\mathrm{0}}\mathrm{\Gamma }(-\mathrm{1}-\mathrm{2}\mathit{\alpha }\left)t^{\mathrm{2}\mathit{\alpha }+\mathrm{1}}\right)=t$ so that the leading term n=0 already gives the correct fRm behaviour: $V_{\mathit{\alpha },h}\left(t\right)\approx t$, so that $a_{\mathrm{0},h}=\mathrm{0}$ (b_0,h can be determined numerically).

A3 Power series expansions about the origin

For many applications, one is interested in the behaviour of R_α,h(t) for scales of months, which is typically less than the relaxation time, i.e. t<1. It is therefore important to understand the small t behaviour. We again consider the Laplace integral for the $\mathrm{0}<h<\mathrm{1}$ case. In this case, we can divide the range of integration in Eq. (A2) into two parts for $\mathrm{0}<x<\mathrm{1}$ and x>1. For the former, we use the expansion in Eq. (A4) and, for the latter,

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{x^{\mathrm{2}\mathit{\alpha }}{e}^{i\mathit{\pi }\mathit{\alpha }}(\mathrm{1}+x^h)(\mathrm{1}+x^h{e}^{i\mathit{\pi }h})}}\\ & ={\displaystyle \frac{e^{-i\mathit{\pi }\mathit{\alpha }}}{e^{i\mathit{\pi }h}-\mathrm{1}}}\sum _{n=\mathrm{1}}^{\mathrm{\infty }}(-\mathrm{1}{)}^{n+\mathrm{1}}(e^{-i(n-\mathrm{1})\mathit{\pi }h}-\mathrm{1})x^{-\mathrm{2}\mathit{\alpha }-nh};x>\mathrm{1}.\end{array}\end{array}$$

We can now integrate each term separately using

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{e}^{-xt}{x}^{nh-\mathrm{2}\mathit{\alpha }}\mathrm{d}x=\sum _{j=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{(-\mathrm{1}{)}^{j-\mathrm{1}}}{(hn-\mathrm{2}\mathit{\alpha }+j)\mathrm{\Gamma }\left(j\right)}}{t}^{j-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }\underset{\mathrm{1}}{\overset{\mathrm{\infty }}{\int }}{e}^{-xt}{x}^{-(nh+\mathrm{2}\mathit{\alpha })}\mathrm{d}x=E_{nh+\mathrm{2}\mathit{\alpha }}\left(t\right)\\ {\displaystyle }& {\displaystyle }=\mathit{\pi }{\displaystyle \frac{t^{-\mathrm{1}+hn+\mathrm{2}\mathit{\alpha }}}{\sin (\mathit{\pi }nh+\mathrm{2}\mathit{\pi }\mathit{\alpha })\mathrm{\Gamma }(hn+\mathrm{2}\mathit{\alpha })}}\\ \text{(A11)}& {\displaystyle }& {\displaystyle }+\sum _{j=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{(-\mathrm{1}{)}^{j-\mathrm{1}}}{(hn+\mathrm{2}\mathit{\alpha }-j)\mathrm{\Gamma }\left(j\right)}}{t}^{j-\mathrm{1}},\end{array}$$

where $E_{\mathit{\beta }}\left(t\right)={\int }_{\mathrm{1}}^{\mathrm{\infty }}{e}^{-xt}{x}^{-\mathit{\beta }}\mathrm{d}x$ is the exponential integral. Adding the two integrals and summing over n, we obtain

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{R}_{\mathit{\alpha },h}\left(t\right)=& \sum _{n=\mathrm{2}}^{\mathrm{\infty }}{D}_n\mathrm{\Gamma }(\mathrm{1}-hn-\mathrm{2}\mathit{\alpha })t^{-\mathrm{1}+hn+\mathrm{2}\mathit{\alpha }}\\ & +\sum _{j=\mathrm{1}}^{\mathrm{\infty }}{F}_j{\displaystyle \frac{t^{j-\mathrm{1}}}{\mathrm{\Gamma }\left(j\right)}},\end{array}\end{array}$$

$$\begin{array}{rl}{F}_j=& {\displaystyle \frac{\mathrm{1}}{\mathit{\pi }h}}\text{Im}\left[{\displaystyle \frac{e^{-i\mathit{\alpha }\mathit{\pi }}}{e^{i\mathit{\pi }h}-\mathrm{1}}}\right(e^{i\mathit{\pi }h}\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}(-\mathrm{1}{)}^n{\displaystyle \frac{e^{i\mathit{\pi }nh}}{(n+a)}}\\ & -\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}(-\mathrm{1}{)}^n{\displaystyle \frac{\mathrm{1}}{(n+a)}})],a={\displaystyle \frac{j-\mathrm{2}\mathit{\alpha }}{h}}\end{array}$$

(we have interchanged the order of summations and used D_n from Eq. A5 with n>0).

The series for the coefficient F_j can now be summed analytically. Although the sum is a special case of the Lipchitz summation and Poisson summation formulae, the easiest method is to use the Sommerfeld–Watson transformation (e.g. Mathews and Walker, 1973) that converts an infinite sum into a contour integral that is then deformed. The Sommerfeld–Watson transformation states that for an analytic function f(z) that goes to zero at least as fast as |z|⁻¹,

$$\begin{array}{}\text{(A13)}& \sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}(-\mathrm{1}{)}^{n}f(n)=-\mathit{\pi }\sum _k{\displaystyle \frac{R_k}{\sin \mathit{\pi }{z}_k}},\end{array}$$

where z_k is the location of the poles of f(z) and R_k is the residue of the corresponding pole. In the above, take

$$\begin{array}{}\text{(A14)}& f\left(z\right)={\displaystyle \frac{e^{iz\mathit{\pi }h}}{(z+a)}}.\end{array}$$

There is a single pole at $z_{\mathrm{1}}=-a$, and the residue is $R_{\mathrm{1}}=e^{-ia\mathit{\pi }h}$; therefore,

$$\begin{array}{}\text{(A15)}& e^{i\mathit{\pi }h}\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{\displaystyle \frac{(-\mathrm{1}{)}^n{e}^{in\mathit{\pi }h}}{(n+a)}}=\mathit{\pi }{\displaystyle \frac{e^{i\mathit{\pi }h(\mathrm{1}-a)}}{\sin \mathit{\pi }a}}.\end{array}$$

The second sum needed in F_j can be obtained using h=0 in the above, so that, overall,

$$\begin{array}{}\text{(A16)}& \begin{array}{rl}{F}_j=& {\displaystyle \frac{\mathrm{1}}{h\mathit{\pi }}}\text{Im}\left[{\displaystyle \frac{e^{-i\mathit{\alpha }\mathit{\pi }}}{e^{i\mathit{\pi }h}-\mathrm{1}}}\left(\mathit{\pi }{\displaystyle \frac{e^{i\mathit{\pi }h(\mathrm{1}-a)}-\mathrm{1}}{\sin \mathit{\pi }a}}\right)\right]\\ =& {\displaystyle \frac{\mathrm{1}}{h\sin \left(\mathit{\pi }\right(j-\mathrm{2}\mathit{\alpha })/h)}}\\ & \cdot \text{Im}\left[{\displaystyle \frac{e^{-i\mathit{\pi }j}{e}^{i\mathit{\pi }(h/\mathrm{2}+\mathit{\alpha })}-e^{-i\mathit{\pi }(h/\mathrm{2}+\mathit{\alpha })}}{e^{i\mathit{\pi }h/\mathrm{2}}-e^{-i\mathit{\pi }h/\mathrm{2}}}}\right].\end{array}\end{array}$$

If j is even, then the term in the square bracket is pure real, hence F_j vanishes. Otherwise

$$\begin{array}{}\text{(A17)}& F_j=-{\displaystyle \frac{\cos \mathit{\pi }(\frac{h}{\mathrm{2}}+\mathit{\alpha })}{h\sin \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)\sin \left(\frac{\mathit{\pi }}{h}\right(j-\mathrm{2}\mathit{\alpha }\left)\right)}},j=\text{odd.}\end{array}$$

Note that F₁>0 for $h+\mathit{\alpha }>\mathrm{1}/\mathrm{2}$ (with $\mathrm{0}\le \mathit{\alpha }<\mathrm{1}/\mathrm{2}$, $\mathrm{0}\le h<\mathrm{2}$), whereas for $h+\mathit{\alpha }<\mathrm{1}/\mathrm{2}$ it is quite complicated (see below).

Comments

• 1.

  These and the following formulae are for t>0; in addition, only the even integer-ordered terms are non-zero (the sum over odd j).

• 2.

  Each integer term of the expansion F_j is itself obtained as an infinite sum, so that the overall result for R_α,h(t) is effectively a doubly infinite sum. This procedure swaps the order of the summation and apparently explains the fact that, while the expansions were derived for the case $\mathrm{0}<h<\mathrm{1}$, the final expansion is valid for $\mathrm{0}\le \mathit{\alpha }<\mathrm{1}/\mathrm{2}$ and the full range $\mathrm{0}<h<\mathrm{2}$: numerically, it accurately reproduces the oscillations when h>1.

• 3.

  The fGn correlation function is given by the single n=2 term:

  $$\begin{array}{}\text{(A18)}& \begin{array}{rl}{R}_h^{\left(\text{fGn}\right)}\left(t\right)& =D_{\mathrm{2}}\mathrm{\Gamma }(\mathrm{1}-\mathrm{2}h)t^{-\mathrm{1}+\mathrm{2}h}\\ & ={\displaystyle \frac{\sin \left(h\mathit{\pi }\right)}{\mathit{\pi }}}\mathrm{\Gamma }(\mathrm{1}-\mathrm{2}h)t^{-\mathrm{1}+\mathrm{2}h}.\end{array}\end{array}$$

  It is also proportional to the correlation function of the fGn-forced h=0; fRn process: $R_h^{\left(\text{fGn}\right)}\left(t\right)=\mathrm{4}{R}_{\mathit{\alpha }=h,\mathrm{0}}\left(t\right)$.

• 4.

  When $\mathrm{0}<\mathit{\alpha }+h<\mathrm{1}/\mathrm{2}$, R is divergent at the origin; this leading term $\mathrm{\Gamma }(-\mathrm{1}-\mathrm{2}(h+\mathit{\alpha }\left)\right)\sin \left(\mathit{\pi }\right(h+\mathit{\alpha }\left)\right)t^{-\mathrm{1}+\mathrm{2}(h+\mathit{\alpha })}/\mathit{\pi }$ is only dependent on h+α corresponding to an fGn with parameter h+α. When $\mathrm{1}/\mathrm{2}<h+\mathit{\alpha }<\mathrm{2}$, it is still the leading fractional term, but the constant F₁ dominates at small t.

• 5.

  The F_j terms diverge when $(j-\mathrm{2}\mathit{\alpha })/h$ is an integer. For example, if α=0, the overall sum over all j thus diverges for all rational h. For irrational h, the convergence properties are not easy to establish, although due to the Γ functions, these series apparently converge for all t≥0, but the convergence is rather slow.

  Figure A1 shows some numerical results for α=0 showing the convergence of the 10th-order fractional 10th-order integer power approximation ($n_{\text{max}}=j_{\text{max}}=\mathrm{10}$). Since the leading (fGn) term diverges for small t, when $h\le \mathrm{1}/\mathrm{2}$, it is more useful to consider the convergence of the difference with respect to the fGn term, i.e. $R_h^{\left(\text{fGn}\right)}\left(t\right)-R_{\mathrm{0},h,a}\left(t\right)$, where the approximation $R_{\mathrm{0},h,a}\left(t\right)$ is from the sum from n=3 to 10 and odd j≤9. Figure A1 shows the logarithm of the ratio of the approximation with respect to the true value: $r={\log }_{\mathrm{10}}|\mathrm{1}-(R_h^{\left(\text{fGn}\right)}\left(t\right)-R_{\mathrm{0},h,a}\left(t\right))/(R_h^{\left(\text{fGn}\right)}\left(t\right)-R_{\mathrm{0},h}\left(t\right)\left)\right|$ (to avoid exact rationals, 10⁻⁴ was added to the h values). From the figure we see that the approximation is satisfactory except for small h. In the next section we return to this.

• 6.

  For $\mathit{\alpha }+h>\mathrm{1}/\mathrm{2}$, when t=0, the only nonzero term is from the constant F₁: $R_{\mathit{\alpha },h}\left(\mathrm{0}\right)=F_{\mathrm{1}}$. This gives the normalization constant. Comparing with Eq. (27), we therefore have

  $$\begin{array}{}\text{(A19)}& \begin{array}{rl}{R}_{\mathit{\alpha },h}\left(\mathrm{0}\right)=& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{G}_{\mathit{\alpha },h}(u{)}^{\mathrm{2}}\mathrm{d}u=F_{\mathrm{1}}\\ =& -{\displaystyle \frac{\cos \mathit{\pi }(\frac{h}{\mathrm{2}}+\mathit{\alpha })}{h\sin \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)\sin \left(\frac{\mathit{\pi }}{h}\right(\mathrm{1}-\mathrm{2}\mathit{\alpha }\left)\right)}},\\ & \mathit{\alpha }+h>\mathrm{1}/\mathrm{2},\begin{array}{l}\mathrm{0}\le \mathit{\alpha }<\mathrm{1}/\mathrm{2},\\ \mathrm{1}/\mathrm{2}<h<\mathrm{2}.\end{array}\end{array}\end{array}$$

  Similarly, when $\mathit{\alpha }+h>\mathrm{3}/\mathrm{2}$, for the quadratic the squared integral of $G_{\mathit{\alpha },h}^{\prime }$ is finite, and it gives the coefficient of the t² term, so that

  $$\begin{array}{}\text{(A20)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{G}_{\mathit{\alpha },h}^{\prime }(s{)}^{\mathrm{2}}\mathrm{d}s& =-{\displaystyle \frac{F_{\mathrm{3}}}{\mathrm{\Gamma }\left(\mathrm{3}\right)}}\\ & ={\displaystyle \frac{\cos \left(\frac{\mathit{\pi }}{\mathrm{2}}\right(h+\mathrm{2}\mathit{\alpha }\left)\right)}{\mathrm{2}h\sin \left(\frac{\mathit{\pi }h}{\mathrm{2}}\right)\sin \left(\frac{\mathit{\pi }}{h}\right(\mathrm{3}-\mathrm{2}\mathit{\alpha }\left)\right)}},\\ & h+\mathit{\alpha }>{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}.\end{array}\end{array}$$
• 7.

  The expression for V_α,h(t) can be obtained by integrating twice (Eq. 36).

• 8.

  In the special cases $h=\mathrm{1}/m$, with m a positive integer, F_j is independent of j, and the integer-powered series can be summed, yielding a result proportional to cosh t. However, this large t divergence is cancelled out by the fractional term, and the result is finite (this partial cancellation is discussed in the next subsection). The special important case $h=\mathrm{1}/\mathrm{2}$ is dealt with in Appendix B.

Figure A1 This shows the logarithm of the relative error in the $R_{\mathrm{0},h}^{(\mathrm{10},\mathrm{10})}\left(t\right)$ approximation (i.e. with 10 fractional terms and 10 integer-ordered terms) with respect to the deviation from the fGn R_0,h(t) $r={\log }_{\mathrm{10}}|\mathrm{1}-(R_h^{\text{fGn}}\left(t\right)-R_{\mathrm{0},h}^{(\mathrm{10},\mathrm{10})}\left(t\right))/(R_h^{\text{fGn}}\left(t\right)-R_{\mathrm{0},h}\left(t\right)\left)\right|$. The lines are for $h=\mathrm{2}/\mathrm{10}$, $\mathrm{4}/\mathrm{10}$, …, $\mathrm{16}/\mathrm{10}$, $\mathrm{18}/\mathrm{10}$ (excluding the exponential case h=1), from left to right (note that convergence is only for irrational h, and therefore an extra 10⁻⁴ was added to each h). For the low h values the convergence is particularly slow.

A4 A convenient approximation

The expansion for R_α,h is the sum of a fractional-ordered and an integer-ordered series. Partial sums appear to converge (Fig. A1), albeit slowly. For simplicity, we consider the case of primary interest, a pure fRn process (α=0). Examination of partial sums shows that the integer-ordered and fractional-ordered terms tend to cancel, the difficulty being due to the coefficient of the integer-ordered terms $j\approx hn+\mathrm{2}\mathit{\alpha }$ that comes from the exponential integral and that can be large when $j\approx hn+\mathrm{2}\mathit{\alpha }$. This suggests an alternative way of expressing the series:

$$\begin{array}{ll}{\displaystyle }{R}_{\mathrm{0},h}\left(t\right)& {\displaystyle }=\sum _{n=\mathrm{2}}^{\mathrm{\infty }}{D}_n{E}_{nh}\left(t\right)+\sum _{j=\mathrm{1}}^{\mathrm{\infty }}{C}_j{\displaystyle \frac{(-\mathrm{1}{)}^{j-\mathrm{1}}}{\mathrm{\Gamma }\left(j\right)}}{t}^{j-\mathrm{1}},\\ \text{(A21)}& {\displaystyle }{C}_j& {\displaystyle }=\sum _{n=\mathrm{2}}^{\mathrm{\infty }}{\displaystyle \frac{D_n}{(hn+j)}},\end{array}$$

where D_n is given by Eq. (A5) and the n sums start at n=2 since D₁=0. C_j can be expressed as

$$\begin{array}{}\text{(A22)}& \begin{array}{rl}{C}_j=& -{\displaystyle \frac{i{e}^{-ih\mathit{\pi }}}{\mathrm{2}\mathit{\pi }h(e^{ih\mathit{\pi }}-\mathrm{1})}}(-(e^{ih\mathit{\pi }+e^{\mathrm{2}ih\mathit{\pi }}}\left)\mathrm{\Phi }\right(-\mathrm{1},\mathrm{1},\mathrm{1}+{\displaystyle \frac{j}{h}})\\ & +\mathrm{\Phi }(e^{ih\mathit{\pi }},\mathrm{1},\mathrm{1}+{\displaystyle \frac{j}{h}})+e^{\mathrm{3}ih\mathit{\pi }}\mathrm{\Phi }(e^{-ih\mathit{\pi }},\mathrm{1},\mathrm{1}+{\displaystyle \frac{j}{h}})),\end{array}\end{array}$$

where Φ is the Hurwitz–Lerch phi function $\mathrm{\Phi }(z,s,a)=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{z}^n(n+a{)}^{-s}$.

We can also expand the exponential integral:

$$\begin{array}{}\text{(A23)}& E_{nh}\left(t\right)=\mathit{\pi }{\displaystyle \frac{t^{-\mathrm{1}+hn}}{\sin \left(\mathit{\pi }nh\right)\mathrm{\Gamma }\left(hn\right)}}+\sum _{j=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{(-\mathrm{1}{)}^{j-\mathrm{1}}}{(hn-j)\mathrm{\Gamma }\left(j\right)}}{t}^{j-\mathrm{1}}.\end{array}$$

For the j_max and n_max partial sums, we have

$$\begin{array}{}\text{(A24)}& \begin{array}{rl}{R}_{\mathrm{0},h}^{(n_{max},j_{max})}\left(t\right)=& \sum _{n=\mathrm{2}}^{n_{max}}{D}_n\mathrm{\Gamma }(\mathrm{1}-nh)t^{-\mathrm{1}+hn}\\ & +\sum _{j=\mathrm{1}}^{j_{max}}{F}_{j,n_{max}}{\displaystyle \frac{(-\mathrm{1}{)}^{j-\mathrm{1}}}{\mathrm{\Gamma }\left(j\right)}}{t}^{j-\mathrm{1}},\\ & F_{j,n_{max}}=C_j+\sum _{n=\mathrm{2}}^{n_{max}}{\displaystyle \frac{D_n}{hn-j}}.\end{array}\end{array}$$

Now define the (j_max, n_max) approximation by

$$\begin{array}{}\text{(A25)}& R_{\mathrm{0},h,n_{max},j_{max}}\left(t\right)={\displaystyle \frac{R_{\mathrm{0},h}^{(n_{max}+\mathrm{1},j_{max})}\left(t\right)+R_{\mathrm{0},h}^{(n_{max},j_{max})}\left(t\right)}{\mathrm{2}}}.\end{array}$$

This has the effect of adding half the next-highest n term and is more accurate; overall, j_max and n_max may now be taken to be much smaller than in the previous approximation. For example, putting n_max=2, j_max=1, we get with the partial sum

$$\begin{array}{}\text{(A26)}& R_{\mathrm{0},h,\mathrm{2},\mathrm{1}}\left(t\right)=R_h^{\left(\text{fGn}\right)}\left(t\right)+{\displaystyle \frac{D_{\mathrm{3}}}{\mathrm{2}}}\mathrm{\Gamma }(\mathrm{1}-\mathrm{3}h)t^{-\mathrm{1}+\mathrm{3}h}+F_{\mathrm{1}},\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{F}_{\mathrm{1}}=C_{\mathrm{1}}+{\displaystyle \frac{D_{\mathrm{2}}}{\mathrm{2}h-\mathrm{1}}}+{\displaystyle \frac{D_{\mathrm{3}}}{\mathrm{2}(\mathrm{3}h-\mathrm{1})}},\\ \text{(A27)}& {\displaystyle }& {\displaystyle }{D}_{\mathrm{2}}={\displaystyle \frac{\sin \left(\mathit{\pi }h\right)}{\mathit{\pi }}},D_{\mathrm{3}}=-{\displaystyle \frac{\sin \left(\mathit{\pi }h\right)(\mathrm{1}+\mathrm{2}\cos (\mathit{\pi }h\left)\right)}{\mathit{\pi }}}.\end{array}$$

Figure A2 The solid line is the constant term F₁, the long dashes are the coefficients $\frac{D_{\mathrm{3}}}{\mathrm{2}}\mathrm{\Gamma }(\mathrm{1}-\mathrm{3}h)$ of the fractional power, and the short dashes are the coefficients of the linear term: $F_{\mathrm{2}}=C_{\mathrm{2}}+\frac{D_{\mathrm{2}}}{\mathrm{2}h-\mathrm{2}}+\frac{D_{\mathrm{3}}}{\mathrm{2}(\mathrm{3}h-\mathrm{2})}$. We can see that the contribution of the linear term (used in the $R_{\mathrm{0},h,\mathrm{2},\mathrm{2}}\left(t\right)$ approximation) for h<1 and t<1 is fairly small, whereas for $\mathrm{1}<h<\mathrm{2}$, it is larger and the $R_{\mathrm{0},h,\mathrm{2},\mathrm{2}}\left(t\right)$ approximation is significantly better than the $R_{\mathrm{0},h,\mathrm{2},\mathrm{1}}\left(t\right)$ approximation (see Fig. A3).

To understand the behaviour, Fig. A2 shows the behaviour of the coefficient of the $t^{-\mathrm{1}+\mathrm{3}h}$ term $\frac{D_{\mathrm{3}}}{\mathrm{2}}\mathrm{\Gamma }(\mathrm{1}-\mathrm{3}h)$, the constant term F₁ and the coefficient of the next integer (linear in t) term

$$F_{\mathrm{2}}=C_{\mathrm{2}}+{\displaystyle \frac{D_{\mathrm{2}}}{\mathrm{2}h-\mathrm{2}}}+{\displaystyle \frac{D_{\mathrm{3}}}{\mathrm{2}(\mathrm{3}h-\mathrm{2})}}.$$

Up until the end of the fGn region ($h=\mathrm{1}/\mathrm{2}$), the $t^{-\mathrm{1}+\mathrm{3}h}$ and F₁ terms have opposite signs and tend to cancel. In addition, we see that for $t\approx <\mathrm{1}$ and h<1, they dominate over the (omitted) linear term. Figure A3 shows that the $R_{\mathrm{0},h,\mathrm{2},\mathrm{1}}$ approximation is surprisingly good for h<1 and is still not so bad for $\mathrm{1}<h<\mathrm{2}$. This approximation is thus useful for monthly-resolution macroweather temperature fields that have relaxation times of years or longer and where h is mostly over the range $\mathrm{0}<h<\mathrm{1}/\mathrm{2}$ but over some tropical ocean regions can increase to as much as h≈1.2 (Del Rio Amador and Lovejoy, 2021a). Figure A3 shows that the (2,1) approximation is reasonably accurate for $t\approx <\mathrm{1}$, especially for h<1.

Figure A3 This shows the logarithm of the relative error in the (2,1) approximation with respect to the deviation from the fGn R_h(t) ($r={\log }_{\mathrm{10}}|\mathrm{1}-(R_h^{\text{fGn}}\left(t\right)-R_{\mathrm{0},h,\mathrm{2},\mathrm{1}}\left(t\right))/(R_h^{\text{fGn}}\left(t\right)-R_{\mathrm{0},h}\left(t\right)\left)\right|$). For h<1, t<0 is of the order ≈30 %, whereas for h>1, it of the order 100 %. The h=1 (exponential) curve is not shown, although when t<0, the error is of order 60 %.

Appendix B: The $h=\mathrm{1}/\mathrm{2}$ special case

When α=0, $h=\mathrm{1}/\mathrm{2}$, and the high-frequency fGn limit is an exact “$\mathrm{1}/f$ noise” (spectrum ω⁻¹); it has both high- and low-frequency divergences. The high-frequency divergence can be tamed by averaging, but not the low-frequency divergence, so that fGn is only defined for $h<\mathrm{1}/\mathrm{2}$. However, for fRn, the low frequencies are convergent over the whole range $\mathrm{0}<h<\mathrm{2}$, and for $h=\mathrm{1}/\mathrm{2}$ we find that the correlation function has a logarithmic dependence at both small and large scales. This is associated with particularly slow transitions from high- to low-frequency behaviours. The critical value $h=\mathrm{1}/\mathrm{2}$ corresponds to the HEBE (Lovejoy, 2021a, b), where it was shown that the value $h=\mathrm{1}/\mathrm{2}$ could be derived analytically from the classical Budyko–Sellers energy balance equation. Therefore, $R_{\mathit{\alpha },\mathrm{1}/\mathrm{2}}\left(t\right)$ and $V_{\mathit{\alpha },\mathrm{1}/\mathrm{2}}\left(t\right)$ characterize the statistics of the temperature response of the classical heat equation response to an fGn-order α forcing.

It is possible to obtain exact analytic expressions for $R_{\mathit{\alpha },\mathrm{1}/\mathrm{2}}\left(t\right)$, $V_{\mathit{\alpha },\mathrm{1}/\mathrm{2}}\left(t\right)$ and the Haar fluctuations; we develop these in this Appendix; for some early results, see Mainardi and Pironi (1996).

The starting point is the Laplace expression (Eq. A2) with $h=\mathrm{1}/\mathrm{2}$:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{R}_{\mathit{\alpha },h}\left(t\right)=& -{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\text{Im}{e}^{-i\mathit{\alpha }\mathit{\pi }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{e^{-xt}\mathrm{d}x}{x^{\mathrm{2}\mathit{\alpha }}(\mathrm{1}+x^{\mathrm{1}/\mathrm{2}})(\mathrm{1}+i{x}^{\mathrm{1}/\mathrm{2}})}}\\ =& -{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }\sqrt{\mathrm{2}}}}\text{Im}{e}^{-i\mathit{\pi }\mathit{\alpha }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{x}^{-\mathrm{2}\mathit{\alpha }}\\ & \cdot \left({\displaystyle \frac{e^{i\mathit{\pi }/\mathrm{4}}}{\mathrm{1}+x^{\mathrm{1}/\mathrm{2}}}}+{\displaystyle \frac{e^{-i\mathit{\pi }/\mathrm{4}}}{\mathrm{1}+x}}-{\displaystyle \frac{e^{i\mathit{\pi }/\mathrm{4}}{x}^{\mathrm{1}/\mathrm{2}}}{\mathrm{1}+x}}\right)e^{-xt}\mathrm{d}x.\end{array}\end{array}$$

We require the following Laplace transforms:

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}{L}_{\mathrm{1}}\left(t\right)=& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{e^{-xt}}{x^{\mathrm{2}\mathit{\alpha }}(\mathrm{1}+x^{\mathrm{1}/\mathrm{2}})}}\mathrm{d}t\\ =& e^{-t-\mathrm{2}i\mathit{\pi }\mathit{\alpha }}\left(\mathrm{\Gamma }\right(\mathrm{1}-\mathrm{2}\mathit{\alpha }\left)\mathrm{\Gamma }\right(\mathrm{2}\mathit{\alpha },-t)-i\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}-\mathrm{2}\mathit{\alpha }\right)\\ & \cdot \mathrm{\Gamma }\left(\mathrm{2}\mathit{\alpha }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},-t\right)),\\ L_{\mathrm{2}}\left(t\right)=& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{e^{-xt}}{x^{\mathrm{2}\mathit{\alpha }}(\mathrm{1}+x)}}\mathrm{d}t=e^t\mathrm{\Gamma }(\mathrm{1}-\mathrm{2}\mathit{\alpha })\mathrm{\Gamma }(\mathrm{2}\mathit{\alpha },t),\\ L_{\mathrm{3}}\left(t\right)=& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{e^{-xt}{x}^{\mathrm{1}/\mathrm{2}}}{x^{\mathrm{2}\mathit{\alpha }}(\mathrm{1}+x)}}\mathrm{d}t=e^t\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}-\mathrm{2}\mathit{\alpha }\right)\mathrm{\Gamma }\left(\mathrm{2}\mathit{\alpha }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},t\right),\end{array}\end{array}$$

where we have introduced the incomplete gamma function: $\mathrm{\Gamma }(a,z)={\int }_z^{\mathrm{\infty }}{u}^{a-\mathrm{1}}{e}^{-u}\mathrm{d}u$ (with a branch cut in the complex plane from −∞ to 0). The general result is thus

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}{R}_{\mathit{\alpha },\mathrm{1}/\mathrm{2}}\left(t\right)=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\left(\sin \mathit{\pi }\mathit{\alpha }\right(L_{\mathrm{1}}\left(t\right)+L_{\mathrm{2}}\left(t\right)-L_{\mathrm{3}}\left(t\right))\\ & +\cos \mathit{\pi }\mathit{\alpha }(-L_{\mathrm{1}}(t)+L_{\mathrm{2}}(t)+L_{\mathrm{3}}(t\left)\right)).\end{array}\end{array}$$

Figure B1 shows plots $R_{\mathit{\alpha },\mathrm{1}/\mathrm{2}}\left(t\right)$ over 8 orders of magnitude in t, indicating the generally very slow convergence to the asymptotic behaviour (shown as straight lines on the right).

Figure B1 $R_{\mathit{\alpha },\mathrm{1}/\mathrm{2}}$ for α increasing from 0 (pure fRn) to $\mathrm{8}/\mathrm{10}$ in steps of $\mathrm{1}/\mathrm{10}$ (on the right: bottom to top). The α=0 curve has a logarithmic divergence at small t (far left). Recall from the section that at large t, $R_{\mathrm{0},\mathrm{1}/\mathrm{2}}\approx t^{-\mathrm{3}/\mathrm{2}}$ and that for α>0, $R_{\mathit{\alpha },\mathrm{1}/\mathrm{2}}\approx t^{\mathrm{2}\mathit{\alpha }-\mathrm{1}}$; for α=0, $\mathrm{1}/\mathrm{5}$, and $\mathrm{2}/\mathrm{5}$ the theoretical asymptotes of the leading terms are indicated for reference.

Figure B1 also shows the singular small t behaviour of the pure fRn case (α=0). In this limit both L₁ and L₂ are singular – they both yield logarithmic small-scale divergences. Pure fRn is of special interest and yields the somewhat simpler result:

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}{R}_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(t\right)=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(e^{-t}\text{erf}i\sqrt{t}-e^t\text{erf}c\sqrt{t})\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\left(e^{t}Ei\right(-t)+e^{-t}Ei(t\left)\right),\\ Ei\left(z\right)=& -\underset{-z}{\overset{\mathrm{\infty }}{\int }}{e}^{-u}{\displaystyle \frac{\mathrm{d}u}{u}},\text{erf}i\left(z\right)=-i\left(\text{erf}\right(iz\left)\right),\\ \text{erf}\left(z\right)=& {\displaystyle \frac{\mathrm{2}}{\sqrt{\mathit{\pi }}}}\underset{\mathrm{0}}{\overset{z}{\int }}{e}^{-s^{\mathrm{2}}}\mathrm{d}s.\end{array}\end{array}$$

We can use these results to obtain small and large t expansions:

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}{R}_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(t\right)=& -\left({\displaystyle \frac{\mathrm{2}{\mathit{\gamma }}_E+\mathit{\pi }+\mathrm{2}\log t}{\mathrm{2}\mathit{\pi }}}\right)+{\displaystyle \frac{\mathrm{2}\sqrt{t}}{\sqrt{\mathit{\pi }}}}-{\displaystyle \frac{t}{\mathrm{2}}}\\ & -\left({\displaystyle \frac{\mathrm{3}+\mathrm{2}{\mathit{\gamma }}_E+\mathit{\pi }+\mathrm{2}\log t}{\mathrm{4}\mathit{\pi }}}\right)t^{\mathrm{2}}\\ & +O\left(t^{\mathrm{3}/\mathrm{2}}\right),t\ll \mathrm{1},\end{array}\end{array}$$

$$\begin{array}{rl}& R_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\sqrt{\mathit{\pi }}}}{t}^{-\mathrm{3}/\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}{t}^{-\mathrm{2}}+{\displaystyle \frac{\mathrm{15}}{\mathrm{8}\sqrt{\mathit{\pi }}}}{t}^{-\mathrm{7}/\mathrm{2}}+O\left(t^{-\mathrm{4}}\right),\\ & t\gg \mathrm{1},\end{array}$$

where γ_E is Euler's constant=0.57… (the asymptotic formula can be obtained as a special case of Eq. A8 in Appendix A, but not the logarithmic small-scale divergence).

To obtain the corresponding results for $V_{\mathrm{0},\mathrm{1}/\mathrm{2}}$, use

$$V_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(t\right)=\mathrm{2}\underset{\mathrm{0}}{\overset{t}{\int }}\left(\underset{\mathrm{0}}{\overset{v}{\int }}{R}_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(u\right)\mathrm{d}u\right)\mathrm{d}v.$$

The exact $V_{\mathrm{0},\mathrm{1}/\mathrm{2}}$ is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{V}_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(t\right)=G_{\mathrm{3},\mathrm{4}}^{\mathrm{2},\mathrm{2}}\left[t\left|\begin{array}{ccc}\mathrm{2},& \mathrm{2},& \mathrm{5}/\mathrm{2}\\ \mathrm{2},& \mathrm{2},& \mathrm{0},& \mathrm{5}/\mathrm{2}\end{array}\right.\right]\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{e^t}{\mathit{\pi }}}\left(\text{Shi}\right(t)-\text{Chi}(t\left)\right)+\left(e^{-t}\text{erf}i\right(\sqrt{t})-e^t\text{erf}(\sqrt{t}\left)\right)\\ \text{(B6)}& {\displaystyle }& {\displaystyle }+t\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_E-\mathrm{1}}{\mathit{\pi }}}\right)-\mathrm{4}\sqrt{{\displaystyle \frac{t}{\mathit{\pi }}}}+{\displaystyle \frac{(\mathrm{1}+t)\log t}{\mathit{\pi }}}+\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_E}{\mathit{\pi }}},\end{array}$$

where $G_{\mathrm{3},\mathrm{4}}^{\mathrm{2},\mathrm{2}}$ is the MeijrG function, Chi is the CoshIntegral function and Shi is the SinhIntegral function. The expansions are

$$\begin{array}{ll}{\displaystyle }{V}_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(t\right)=& {\displaystyle }-{\displaystyle \frac{t^{\mathrm{2}}\log t}{\mathit{\pi }}}+{\displaystyle \frac{\mathrm{191}-\mathrm{156}{\mathit{\gamma }}_E-\mathrm{78}\mathit{\pi }}{\mathrm{144}\mathit{\pi }}}+{\displaystyle \frac{\mathrm{16}}{\mathrm{15}\sqrt{\mathit{\pi }}}}{t}^{\mathrm{5}/\mathrm{2}}\\ \text{(B7)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{t^{\mathrm{3}}}{\mathrm{6}}}-{\displaystyle \frac{t^{\mathrm{4}}\log t}{\mathrm{12}\mathit{\pi }}}+O\left(t^{\mathrm{3}/\mathrm{2}}\right),t\ll \mathrm{1},{\displaystyle }{V}_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(t\right)=& {\displaystyle }t+{\displaystyle \frac{\mathit{\pi }+\mathrm{2}{\mathit{\gamma }}_E}{\mathit{\pi }}}+{\displaystyle \frac{\mathrm{2}\log t}{\mathit{\pi }}}-{\displaystyle \frac{\mathrm{4}}{\sqrt{\mathit{\pi }}}}{t}^{\mathrm{1}/\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathit{\pi }}}}{t}^{-\mathrm{1}/\mathrm{2}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}{t}^{-\mathrm{2}}+{\displaystyle \frac{\mathrm{15}}{\mathrm{4}\sqrt{\mathit{\pi }}}}{t}^{-\mathrm{3}/\mathrm{2}}+O\left(t^{-\mathrm{4}}\right),t\gg \mathrm{1}.\end{array}$$

We can also work out the variance of the Haar fluctuations:

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}& 〈\mathrm{\Delta }{U}_{\mathrm{0},\mathrm{1}/\mathrm{2}}^{\mathrm{2}}(\mathrm{\Delta }t{)}_{\text{Haar}}〉\\ & ={\displaystyle \frac{\mathrm{\Delta }{t}^{\mathrm{2}}\log \mathrm{\Delta }t}{\mathrm{4}\mathit{\pi }}}+{\displaystyle \frac{\mathrm{6}\mathit{\pi }+\mathrm{12}{\mathit{\gamma }}_E-\log \mathrm{16}+\mathrm{960}\log \mathrm{2}}{\mathrm{240}\mathit{\pi }}}\\ & +{\displaystyle \frac{\mathrm{512}(\sqrt{\mathrm{2}}-\mathrm{2})}{\mathrm{240}\sqrt{\mathit{\pi }}}}\mathrm{\Delta }{t}^{\mathrm{1}/\mathrm{2}}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{3}}}+O\left(\mathrm{\Delta }{t}^{\mathrm{3}/\mathrm{2}}\right),\mathrm{\Delta }t\ll \mathrm{1},\end{array}\end{array}$$

$$\begin{array}{rl}& 〈\mathrm{\Delta }{U}_{\mathrm{0},\mathrm{1}/\mathrm{2}}^{\mathrm{2}}(\mathrm{\Delta }t{)}_{\text{Haar}}〉\\ & =\mathrm{4}\mathrm{\Delta }{t}^{-\mathrm{1}}-{\displaystyle \frac{\mathrm{32}\sqrt{\mathrm{2}}}{\sqrt{\mathit{\pi }}}}\mathrm{\Delta }{t}^{-\mathrm{3}/\mathrm{2}}+{\displaystyle \frac{\mathrm{3}\mathrm{\Delta }{t}^{-\mathrm{2}}\log \mathrm{\Delta }t}{\mathit{\pi }}}+O\left(\mathrm{\Delta }{t}^{-\mathrm{2}}\right),\\ & \mathrm{\Delta }t\gg \mathrm{1}.\end{array}$$

Figure B2 shows numerical results for α=0 and $h=\mathrm{1}/\mathrm{2}$. The transition between small and large t behaviour is extremely slow; the nine orders of magnitude depicted in the figure are barely enough. The extreme low ($R_{\mathrm{1}/\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$ (dashed) asymptotes on the left to a 0 slope (a square root logarithmic limit, Eq. B8) and to a $-\mathrm{3}/\mathrm{4}$ slope on the right. The rms Haar fluctuation (black) changes slope from H=0 to $-\mathrm{1}/\mathrm{2}$ (left to right). Figure B2 also shows the logarithmic derivative of the rms Haar (black) compared to a regression estimate over 2 orders of magnitude in scale (dashed; a factor of 10 smaller and 10 larger than the indicated scale was used; this represents a possibly empirically accessible range). This figure underlines the gradualness of the transition from H=0 to $H=-\mathrm{1}/\mathrm{2}$. If empirical data were available only over a factor of 100 in scale, depending on where this scale was with respect to the relaxation timescale (unity in the plot), the rms Haar fluctuations could have any slope in the range 0 to $-\mathrm{1}/\mathrm{2}$, with only small deviations.

Figure B2 The logarithmic derivative of the rms Haar fluctuations of $U_{\mathrm{0},\mathrm{1}/\mathrm{2}}$ (solid) in Fig. B1 compared to a regression estimate over 2 orders of magnitude in scale (dashed; factors of 10 smaller and 10 larger than the indicated scale were used). This plot underlines the gradualness of the transition from slopes 0 to −0.5 corresponding to apparent H=0 to $H=-\mathrm{1}/\mathrm{2}$ scaling. Over a range of 100 or so in scale there is approximate scaling but with exponents that depend on the range of scales covered by the data. If data were available only over a factor of 100 in scale, the rms Haar fluctuations could have any slope in the fGn range 0 to $-\mathrm{1}/\mathrm{2}$, with only small deviations.