Exact distributions of the maximum and range of random diffusivity processes

We study the extremal properties of a stochastic process x t defined by the Langevin equation ẋt=2Dtξt , in which ξ t is a Gaussian white noise with zero mean and D t is a stochastic ‘diffusivity’, defined as a functional of independent Brownian motion B t . We focus on three choices for the random diffusivity D t : cut-off Brownian motion, D t ∼ Θ(B t ), where Θ(x) is the Heaviside step function; geometric Brownian motion, D t ∼  exp(−B t ); and a superdiffusive process based on squared Brownian motion, Dt∼Bt2 . For these cases we derive exact expressions for the probability density functions of the maximal positive displacement and of the range of the process x t on the time interval t ∈ (0, T). We discuss the asymptotic behaviours of the associated probability density functions, compare these against the behaviour of the corresponding properties of standard Brownian motion with constant diffusivity (D t = D 0) and also analyse the typical behaviour of the probability density functions which is observed for a majority of realisations of the stochastic diffusivity process.


Introduction
The statistics of extreme values (EVs) of stochastic processes has been in the focus of extensive research in the mathematical (see, e.g. [1][2][3]) and physical (see, e.g. [4][5][6][7][8][9][10][11][12][13][14][15]) literature over several decades. More recently, EV properties have also received attention in the areas of mathematical finance [18,19] in which stochastic processes represent one of the main components in the modelling of the dynamics of asset prices, of computer science [20,21], as well as of the analysis of 'records' of different kinds [19,22,23]. Apart from 'simple' EV problems asking for the maximum behaviour of a variable, 'dual' EVs of the min-max and max-min families are relevant in game theory [24] or reliability engineering [25], for which a universal Gumbel limit law emerges [26,27].
Typically, one computes several types of EVs, which are either interrelated or independent of each other, and hence, provide complementary information about the process x t under consideration. Commonly considered EVs are, for instance, the persistence probability for not crossing the initial value x 0 of the process [10], the related probability that the process does not reach a given threshold or a given point in space up to time t (i.e. the 'survival' probability) [5], or the first-passage time to a given threshold or spatial location [4][5][6][7][8][9][10][13][14][15]. For one-dimensional processes, one often considers the maximal positive M or negative M displacements and the range R (also called the span or the extent) of x t (here and henceforth we assume that x 0 = 0) on a given time interval 6 : Here T represents the length ('observation time') of the time series x t under consideration. When a random process x t evolves on a one-dimensional lattice, the range R defines another important property, namely, the number of distinct visited sites up to time T [28]. We also note that complementary characteristics of extremal values of Brownian motions such as the distribution of times between minima and maxima has been evaluated recently [29,30]. Knowledge of the EV statistics is conceptually important for the understanding of various facets of the stochastic process x t and is relevant for diverse physical phenomena and also in applications in finance, sociophysics and biology, since EVs often trigger a particular response of the system. A prominent application is molecular chemical reaction kinetics, in which a diffusing molecule hits a reaction centre [31]. For instance, during gene regulation a protein needs to diffusively search a specific binding site on the cellular DNA [32]. Recent research demonstrated that for typical biochemical situations with extremely low reactant concentrations knowledge beyond mean chemical rates [31] is essential, due to the significant separation of relevant time scales even in simple geometries [33,34]. Notably, geometry-control of reaction time scales in gene regulation [35,36] is closely related to the most likely reaction time [33,34]. We also mention that the knowledge of EVs is often very beneficial for a non-perturbative analysis of complicated functionals of x t , permitting for a construction of convergent bounds and, hence, for obtaining non-trivial exact results [37][38][39][40][41][42].
Most of the available analyses pertain to the paradigmatic process of Brownian motion, or to lattice random walks. In particular, the exact probability density function P T (M) of the maximal positive displacement M T of a one-dimensional Brownian motion has already been derived exactly in the early work by Lévy [1]: denoting the diffusion coefficient as D 0 and supposing that t ∈ (0, T) one has normalised to the positive semi-axis M 0, compare figure 1. Subsequently, the probability density function P T (R) of the range of Brownian motion was obtained [2] (see also [43][44][45]) These two series representations are exact and thus equivalent. The series in the first line highlights the asymptotic behaviour in the limit R → ∞, while the one in the second line is appropriate for the analysis in the small-R (or long-T) limit. Note also that while P T (M = 0) is finite for any finite T, the probability density function P T (R) abruptly drops to zero when R → 0, see the entire shape in figure 1. Further on, more complicated multivariate joint distributions of maxima, minima and the range were evaluated [7], while correlations between maxima or between values of the range achieved on different time intervals were studied in [46][47][48][49]. A remarkable result has recently been obtained for the distribution of the time instant at which the range of Brownian motion first reaches a prescribed value [50]. Concurrently a variety of first-passage phenomena associated with Brownian motion have been analysed using exact approaches [4][5][6][7][8][9][10][11][12][13]. On top of this several accurate approximation schemes have been analysed, permitting one to consider first-passage events in rather complicated, experimentally-relevant geometries [51,52]. However, the progress in the theoretical analysis of EV statistics for more general processes, in particular, other than standard Brownian motion and especially non-Markovian processes, remains limited, and typically only the behaviour of the expected values of the EVs is known [4][5][6][7][8][9][10]13].
In this paper we derive exact compact expressions for the probability density function P T (M) of the maximal displacement and for the probability density function P T (R) of the range for three random-diffusivity stochastic processes introduced recently in [53]. In these models, the process x t evolves in a one-dimensional system according to the Langevin equation where ξ t is standard white Gaussian noise with zero mean and covariance ξ t ξ t = δ(t − t ), D 0 is a constant scale factor, and V(B t ) is a dimensionless random diffusivity defined as a functional of independent Brownian motion B t : with the diffusivity D B . Here and henceforth, angular brackets denote averaging with respect to all possible realisations of the Brownian motion B t , while the overline corresponds to averaging over realisations of the white noise process ξ t . We note that different versions of the model in (4) corresponding to different choices of the functional V(B t ) have been extensively studied in recent years within the context of diffusion in complex heterogeneous environments [54][55][56][57][58][59][60][61][62][63][64][65], dynamics of particles involved in polymerisation processes [66,67] which can be anomalous in the non-Stokesian limit [68], as well as in the mathematical finance literature (see, e.g. [69]). We also mention that stochastically varying diffusivities were identified in simulations of diffusing proteins with fluctuating shape [70], and switching between low and high mobility states was observed in simulations of protein-crowded membranes [71] and the motion of tracer particles in the cytoplasm of mammalian cells [72]. Following [53] we define the three random diffusivity models under study here as follows: is the Heaviside step function with the property Θ(x) = 1 for x 0 and zero otherwise. In this model, the process x t undergoes standard Brownian motion once B t > 0, and it pauses at its current location whenever B t is negative. Here, the mean-squared displacement x 2 t = D 0 t shows Brownian behaviour, however, the diffusion coefficient is smaller by a factor of two than the diffusion coefficient of standard Brownian motion with V = 1.
(II) In model II we choose V(B t ) = exp(−B t /a), where a is a scale parameter of dimension length. This choice for V(B t ) corresponds to geometric Brownian motion, as assumed for the time evolution of an asset price in the paradigmatic Black-Scholes model [73]. In this model, diffusion is strongly anomalous and the mean squared displacement has an exponential time dependence, Here, the process x t accelerates when B t goes away from the origin in either direction, and we are thus facing a superdiffusive behaviour as the process x t in (4) exhibits a random ballistic growth with time.
We focus here on the generalisation of expressions (2) and (3), derived for Brownian motion, to the above defined three models of random diffusivity. We thus seek the exact expressions for probability density functions of the maximum and of the range, respectively, defined as and where P T (M) and P T (R) denote the probability density functions calculated for a given realisation of B t (and thus a given realisation of diffusivity). We note that the latter realisation-dependent distributions are evidently given by expressions (2) and (3) with D 0 T replaced by the integral D 0 T 0 V(B t )dt, implying that P T (M) and P T (R) are, respectively, a Gaussian function or an infinite sum of Gaussians with random variances. As we proceed to show, the averaging over realisations of B t can be performed exactly for the three models under study and requires only the knowledge of the moment-generating function Υ(T; λ) of the random variable D 0 This function can indeed be calculated exactly for many cases (see, e.g. [3,21] and references therein) and, in particular, for the models we study here. We proceed to show that the averaged distributions, i.e. P T (M) and P T (R), exhibit a markedly different behaviour, as compared to the distributions in (2) and (3). We will compare these predictions against the estimates of the 'typical' behaviour of these distributions (see, e.g. [74,75]), where N M and N R are normalisation constants, while p is an irrelevant auxiliary parameter of inverse length that is introduced to deal with dimensionless quantities under the logarithm but cancels anyway. We will demonstrate that their functional form is supported by some atypical realisations of the process B t . The paper is organised as follows. In section 2 we briefly summarise our main results. In two subsequent sections 3 and 4, we present the details of the derivations of our main results, analyse their asymptotic behaviour and their moments, and also estimate their effective broadness by calculating the coefficients of variation of the respective distributions. Additionally, we compare our analytical predictions with the results of numerical simulations. Section 3 is devoted to the maximum, in section 4 we consider the range. Next, in section 5 for the example of the distribution of the maximum, we will discuss its 'typical' shape which should be observed for a majority of realisations of the process B t (or for small statistical samples) and demonstrate that the exact form obtained for P T (M) (and for P T (R)) defined in (6) stems from some atypical realisations of the stochastic diffusivity process. Concluding remarks are provided in section 6.

Main results
In this section we summarise our main results for the probability density functions P T (M) and P T (R), (6) and (7), for the three models under study. The parameters entering these results were defined above in the description of our models.

Model I
For model I we find that where K 0 (z) is the modified Bessel function of the second kind of the zeroth order. In turn, the exact probability density function P (I) T (R) can be written in either of two equivalent forms: (i) as we proceed to show, the analysis of the short-R behaviour (i.e. the behaviour of the left tail of the probability density function P T (R)) can be realised from the following expression, where I 0 (z) is the modified Bessel function of the first kind; (ii) in turn, the behaviour of the right tail is conveniently given by an alternative series expansion, Expressions (12) and (13) are related to each other through the Poisson summation formula.

Model II
For model II we obtain the following exact expression for the probability density function of the maximum, The detailed discussion of its rather unusual asymptotic behaviour is presented in the next section. In turn, the probability density function P (II) T (R) of the range obeys the exact expression where K iz (x) is the modified Bessel function of the second type of purely imaginary order. This latter form is suitable for the analysis of the short-R behaviour (see section 4). An alternative form appropriate for the analysis of the large-R behaviour follows from (15) via the Poisson summation formula and reads

Model III
For model III the probability density function of the maximum has the exact form where Γ(z) is the gamma function. The probability density function of the range admits the exact expansion which is suitable in the short-R limit, while the right tail of the distribution can be accessed via an equivalent expansion, The rest of the paper presents the details of the derivations of our main results and a discussion of their behaviour.

Probability density function of the maximal displacement
We present the details of derivations of the exact expressions for the probability density function P T (M) summarised earlier in section 2. We find it expedient to base our analysis here on the exact expressions for the first-passage time density H(t|M) that was derived for all three models under study in a recent paper [65]. An alternative approach which takes advantage of the moment-generating function (8) will be used later on in section 5 and will permit us to access the typical behaviour of the probability density.
Let S T (M) denote the survival probability, i.e. the probability that the process x t , starting at the origin x 0 = 0 at t = 0, does not reach a point M > 0 within the time interval t ∈ (0, T). This probability can be expressed as where H(t|M) is the probability density function of the event that the process x t reached the point M for the first time at the time instant t. As a consequence, the desired probability density function P T (M), which defines the probability density that the maximal positive displacement of the process x t within the time interval t ∈ (0, T) is exactly equal to M obeys In the case of standard Brownian motion is the celebrated Lev´y-Smirnov distribution, and, eventually, P T (M) is given by (2).

Model I
For model I the first-passage time density obeys [65] Differentiating the latter expression with respect to M, inserting the result into (21) and integrating it over t, we find our compact expression (11). Note that the density in (22) resembles-but is not identical to-the density P (I) T (M). Moreover, due to the presence of K 0 (z), the distribution of the maximum exhibits a different asymptotic behaviour in the limits of small and large M as compared to behaviour (2) of Brownian motion. From (11) the asymptotic limit M → ∞ produces to leading order that is, the probability density decreases with M faster, due to the additional factor 1/M, than expression (2). For the opposite limit M → 0 we find that where γ ≈ 0.5772 is the Euler-Mascheroni constant. Expression (24) implies that P (I) T (M) logarithmically diverges in this limit, while expression (2) remains bounded. Overall we observe that the probability density function (11) is shifted towards smaller M values compared to the distribution (2). In particular, the expected (with respect to the distribution in (11)) value of the maximal displacement of the process x t in model I obeys i.e. it grows with T exactly at the same rate as the expected maximum of Brownian motion, T, but has a slightly smaller prefactor (4/π 3/2 ≈ 0.72 while 2/π 1/2 ≈ 1.12). The expression in (25) can be generalised to derive the moments of the maximum of the process x t for model I of arbitrary, not necessarily integer order q 0, From this expression we also derive the coefficient of variation v (I) M of the distribution (11): By definition, this property measures the relative weight of fluctuations around the mean value. Hence, for model I these fluctuations are of nearly the same order as the expected value itself, such that P (I) T (M) is effectively broad [11,12]. Note that the distribution of the maximum of a standard Brownian motion, (2), appears to be somewhat narrower; there, the coefficient of variation v (BM)

Model II
For model II the exact expression for the probability density function of the first-passage time is given by [65] Again, differentiating (28) over M and integrating the resulting expression over t, we arrive at our result in (14). The small-M asymptotical behaviour of expression (14) obeys Rather surprisingly, this limiting behaviour is exactly the same as that of (2) for the maximum of standard Brownian motion. In contrast, the large-M asymptotic behaviour is very different from that of standard Brownian motion and follows i.e. the right tail of the distribution P (II) T (M) is that of a log-normal distribution. In view of such a 'heavy' tail one expects that higher values of M are more likely than in case of a standard Brownian motion.
The moments of the distribution P (II) T (M) of arbitrary order q 0 can be obtained by a straightforward integration of expression (14), leading us to where q n denotes the binomial coefficient. Naturally, when q is an integer the series is truncated at n = q, as can be observed directly from the expression in the first line of (31). From (31) we have, in particular, and so on. Here we used the notation τ = D B T/a 2 . We observe that in the case of model II there is no unique time scale, in contrast to model I (and also to model III below). This is a direct consequence of the fact that the right tail of P (II) T (M) decreases with M slower than an exponential function, which gives rise to the behaviour specific to the so-called strongly anomalous superdiffusion for which a growth of the moments with time is not characterised by a unique exponent (see, e.g. [76][77][78]).
The first two expressions in (32) permit us to evaluate the coefficient of variation v (II) M of the probability density function P (II) T (M) in (14): Remarkably, v (II) M diverges exponentially, v (II) M exp(τ/4)/ √ 2 as τ → ∞ (i.e. the observation time T tends to infinity). This signifies that moments of arbitrary order are not representative of the actual behaviour and knowledge of the full distribution P (II) T (M) is crucial.

Model III
Lastly, for model III the first-passage time density is [65] Differentiating this expression with respect to M and integrating over t, we arrive at our result in (17). For small M, the gamma function in (17) tends to a constant (with corrections of order O(M 2 )), and hence, one has In turn, for M → ∞, the asymptotic behaviour of P (III) T (M) is given by i.e. the right tail of the probability density function of the maximal displacement is an exponential function and hence, is also 'heavier' than the Gaussian tail of the corresponding probability distribution of the maximum of standard Brownian motion. Evidently, expression (17) also favours higher values of the maximum M than the probability density function (2). The moments of the distribution (17) obey where the dimensionless numerical amplitude f q is given by We were unable to perform the integral in the latter equation and, hence, to derive an explicit expression for f q , except for the particular case when q is an even integer, q = 2n. In this latter case, f q is given by where A n are integers forming Sloane's sequence A126156 [79]. In particular, The numerical factor f q as a function of q is depicted in figure 2. We realise that f q turns out to be a non-monotonic function of q. Lastly, we estimate numerically the value of the coefficient of variation of implying that fluctuations around the mean value of the maximum exceed the latter such that the distribution is effectively broad. Figure 3 presents the exact probability density functions P T (M) (solid curves) and their asymptotic forms (dashed and dash-dotted curves) for the three models, highlighting the ranges of validity of the small-M limit as well as the onset of the large-M asymptotic behaviours. The results are shown on both linear and log-log scales, to highlight the asymptotic behaviour as well as the respective crossovers. Note specifically the divergence of P T (M) in the limit M → 0 for model I.

Relation between the moments of the maximum and of the random diffusivity
To close this section we present a general relation between the moments of the maximum and the moments of the random variable D 0 which holds for arbitrary q 0. This relation can be proven directly by using the definition in (21) and also a general expression for the first-passage time distribution presented in our previous work [65]. Below we will merely demonstrate the validity of (42) by establishing a relation between the moments of the maximum and the moments of the range. Using the standard 'replica trick' we find the following simple expression that connects the typical behaviour of the maximum and the typical behaviour of the random variable D 0 where γ is again the Euler-Mascheroni constant. In particular, in the special case V(B t ) ≡ 1 (i.e. when the process x t in (4) is standard Brownian motion), equation (43) reproduces the known result ln M (BM)

Probability density function of the range
In this section, we first present the arguments underlying the derivation of P T (R) and evaluate general expressions which highlight the short-R and large-R behaviour, i.e. the left and right tails of P T (R), respectively. We then establish a general relation between the moments of the range and the moments of the random variable D 0 t 0 V(B t ) dt, which also permits us to link the moments of the range and the moments of the maximum in the random diffusivity processes. Lastly, we will concentrate on the particular cases and evaluate the exact forms of P T (R) for the three models under study.
The probability density function of the range R of the process (4) can be evaluated by writing down the corresponding Fokker-Planck equation for the position probability density function Π(x, t) (in which the diffusion coefficient D t is a random function of time), appropriately rescaling the time variable and then solving the resulting diffusion equation subject to adsorbing boundary conditions. The steps involved in this approach are well described, e.g. in [43][44][45]. In this procedure we find that P T (R) can be conveniently represented by two alternative forms, one of which is suitable for the analysis of the small-R behaviour of the probability density function of the range, while the second one is adapted to the large-R asymptotic behaviour.
In the first case P T (R) is given by with and where Υ(T; λ) is the moment-generating function which is defined earlier in (8). We note that in virtue of (46) the knowledge of an exact form (8) of Υ(T; λ) appears to be the key ingredient for finding exact forms of P T (R) (see also [62] for the role of this function for the analysis of the first-passage time densities). In turn, the large-λ tail of Υ(T; λ) (corresponding to such realisations of B t when D 0 T 0 V(B t )dt is small) is responsible for the behaviour of Ψ T (R) in the limit when R → 0. We proceed to show that such a behaviour can be markedly different depending on how fast Υ(T; λ) vanishes when λ → ∞. In this sense, the three models under study provide representative examples of different kinds of such a behaviour: in model I the moment-generating function Υ(T; λ) vanishes as a power-law when λ → ∞ and P T (R) approaches a constant value as R → 0, while for both models II and III Υ(T; λ) ∼ exp(− √ λ) in the leading order in the limit λ → ∞, and ln P T (R) exhibits a singular behaviour of the form ln P T (R) ∼ −1/R.
In the second case, the form appropriate for the analysis of the large-R behaviour can be obtained by the Poisson summation formula Further on, using the integral identity we cast (47) into the form In case of standard Brownian motion the latter expression reduces to the series in the first line in (3). One observes that in the limit R → ∞ the integral in the latter expression is dominated by the behaviour of Υ(T; q 2 ) in the vicinity of q = 0, which corresponds to the small-λ asymptotic behaviour of the moment-generating function in (8) (and hence, to such realisations of B t for which D 0 T 0 V(B t )dt is large). However, we find 1 − Υ(T; q 2 ) = O(q 2 ) for models I and III (while for model II there are logarithmic corrections to the q 2 -dependence), meaning that P T (R) decays sufficiently fast in all three models to ensure the existence of all moments. Hence, the precise form of the large-R tails of P T (R) cannot be, in principle, deduced from the small-q expansions of Υ(T; q 2 ) and we have to perform the corresponding integrals explicitly. In doing so, we will demonstrate below that the large-R tails of P T (R) are markedly different in all three models.
Relation (49) between P T (R) and the moment-generating function Υ(T; λ) of D 0 T 0 V(B t )dt implies a simple and quite general relation between the moments of the range and the moments of the latter random variable. Indeed, multiplying both sides of (47) by R q (q 0) we find that whenever the moments of D 0 T 0 V(B t )dt exist, the following relation holds where ζ(z) is the Riemann zeta function (note that for q = 0, 1, 2, one has to take the limit as q approaches one of these integer values). For instance, we find and so on. Next, resorting again to the usual 'replica trick' we also deduce from (50) a linear relation between the averaged logarithm of the range and the averaged logarithm of D 0 T 0 V(B t )dt, which thus connects the 'typical' behaviour of these two random variables, In (53) A ≈ 1.2824 is Glaisher's constant-a mathematical constant related to the asymptotical behaviour of the Barnes G-function (double gamma-function) [80]. The latter emerges, e.g. in the normalisation of the joint distributions of eigenvalues in Gaussian ensembles of the random matrix theory and, hence, A plays an important role in the asymptotic analysis of some characteristic properties of such ensembles (see, e.g. [81]). We emphasise that (50) and (53) where a is an auxiliary length scale which was introduced in order to get dimensionless units. While the scaling of the typical range with √ D 0 T appears quite intuitive, the proportionality factor ≈ 2.1647 in (54) is rather nontrivial. In particular, its relation to the Glaisher's constant A is surprizing. Note that the expected value of the range in (51) also scales as √ D 0 T, but the proportionality factor 4/ √ π ≈ 2.2568 is somewhat larger. As a consequence, for most of realisations of trajectories of standard Brownian motion their ranges appear to be smaller than the range averaged over all realisations, which implies that some less probable, atypical realisations dominate the value of the range. This, in turn, implies that even in the case of simple standard Brownian motion the knowledge of the full probability density function of the range is vital.
We now evaluate explicit expressions for P T (R) for the three models of random diffusivity presented in section 2 and discuss their asymptotic behaviour.

Model I
In model I we have V(B t ) = Θ(B t ) and, hence, the exponential function of T 0 Θ(B t )dt in (46) is simply the moment-generating function of the occupation time of Brownian motion on a positive half-axis in the time interval (0, T). Explicitly one has (see [53] for more details) such that The latter expression together with (45) yields our result (12). The small-R behaviour of Ψ (I) T (R) and, hence, of the probability density function P (I) T (R) can be derived directly from (56) by taking advantage of the asymptotic expansion which holds for large values of the argument x. Inserting this expansion into (56) and performing the summation over m we find By virtue of expression (45) the asymptotic small-R expansion for P (I) T (R) is obtained from (58) by merely multiplying the latter by R and differentiating the resulting expression twice with respect to R. In doing so we arrive at a rather curious conclusion that, in contrast to the behaviour of the probability density function of the range of standard Brownian motion, P (I) T (R) does not vanish in the limit R → 0 but rather approaches the non-trivial constant value Therefore, despite the fact that the process x t in model I exhibits the 'diffusive' behaviour x 2 t = D 0 t, its rather intricate character causes significant departures from the behaviour of standard Brownian motion-here, the fact that x t may pause at the origin for a random time (having a broad distribution without even the first moment) once B t goes initially to negative values, entails a finite value of the probability density at R = 0. This behaviour is also in line with the divergence of P (I) T (M) in the limit M → 0. Note that in (59) the amplitude decays as the inverse square root of T.
To construct the asymptotic large-R expansion of P (I) T (R) we turn to the representation in (49). Using (55) we then have and, hence, P (I) T (R) is cast into the form in (13). The leading, large-R behaviour is provided by the first term in the series in the latter equation, i.e. the right tail of the distribution P (I) T (R) vanishes faster than a Gaussian function due to the additional factor 1/R.

Model II In model II the dimensionless diffusivity is governed by geometric Brownian motion: V(B t ) = exp(−B t /a).
In this case, one has (see [82] and also [83,84]) Combining (62), (45) and (46), we thus arrive at our result (15). The small-R asymptotic behaviour of Ψ (II) T (R) and, hence, of P (II) T (R) (see (15)) can be conveniently accessed by taking advantage of the Kontorovich-Lebedev-type representation in (62). Using the large-x expansion we find from (46) that Ψ (II) T (R) (see (15)) admits the following form in the limit R → 0, We notice next that the leading in the R → 0 behaviour is provided by the term with m = 0. As a consequence, we arrive at the asymptotic formula Therefore, in model II the probability density function vanishes as R → 0, in contrast to model I. Note also that the essential singularity in model II is somewhat weaker, ln P (II) T (R) ∝ 1/R, as compared to the singular behaviour specific to standard Brownian motion, for which one has ln P (BM) T (R) ∝ 1/R 2 , see (3). The analysis of the large-R asymptotic behaviour of P (II) T (R) hinges on the expansion (49). Inserting (62) into (49) and evaluating the two-fold integral we find that expression (49) for model II admits the explicit form in (16). Inspecting the latter formula we notice that the dominant large-R behaviour is provided by the term with m = 1. As a consequence, we get whose form is nearly identical (apart from numerical factors) to the asymptotic result (30) describing the right tail of the probability density function P (II) T (M).

Model III
In model III the random diffusivity is given by V(B t ) = B 2 t /a 2 and the moment-generating function is (see [53,85,86]) Hence, we find that The latter expression, together with (45) and (46), result in (18). The small-R asymptotic behaviour of P (III) T (R) can be obtained directly from (69) by noticing that the series converges very rapidly and the dominant behaviour is provided by the zeroth term, i.e.
Substituting the latter asymptotic form into (45) and differentiating, we arrive at which shows that P (III) T (R) vanishes exponentially fast when R → 0. The large-R asymptotic behaviour of the probability density function of the range in model III can be determined as follows. Using again relation (68) we get which yields our result (19). Noticing finally that in the limit R → ∞ the dominant contribution to the expansion in (19) comes from the term m = 1 we thus arrive at the asymptotic formula Figure 4 illustrates the behaviour of the probability density function P T (R) and its asymptotic forms for the three models. To emphasise the crossover behaviours as well as the asymptotic forms of the probability density functions we show the results both on linear and log-log scales. Note specifically that while models II and III exhibit a suppression of the probability density functions to zero in the limit R → 0, in model I a finite value at R = 0 is reached.
In figure 5 we compare the probability density functions P T (M) and P T (R) for the three random diffusivity models. Note the different large-M asymptotic behaviours for the different models as discussed above. The analytical results are also confronted with Monte Carlo simulations. These simulations were obtained with the Euler integration scheme applied to the Langevin equation (4). For each realisation an independent Brownian motion run is generated to compute the dimensionless random diffusivity through the specific functional V(B t ). In this way the two noises are varied simultaneously and independently. Perfect agreement with analytical formulas is observed even for a moderately large sampling with 10 000 realisations.

Relation between the moments of the maximum and of the range
We derive a general expression for the moments of the range for the processes in (4). To this end, we observe that for the three models studied here the probability density function P T (R) can be formally written as where P T (M) is the corresponding probability density function of the maximum. One may expect that this relation is valid in general for an arbitrary process defined in (4) but we are not in the position to prove it here. Then, multiplying both sides of the latter expression by R q we obtain, through a simple change of the integration variable, the following intricate relation between the moments of the range and the moments of the maximum, R In particular, we get and so on. While the first relation (76) is obvious, relations (77) and (78) are non-trivial results. Lastly, comparing (76) and (50), we arrive at (42), which we presented without a derivation in the previous section.
As a direct consequence of relation (75) we can write down exact closed-form expressions for the moments of the range of arbitrary (not necessarily integer) order. In turn, the latter permit us to evaluate the coefficients of variation of the distributions of the range. For model I we find which is about 30 per cent smaller than the coefficient of variation (27) of the corresponding distribution of the maximum. The coefficient of variation of the range for model II is given explicitly by and, hence, grows as ln(2)/2 exp(τ/4) in the limit τ → ∞. This growth is thus somewhat slower than for the corresponding coefficient of variation (33) of the distribution of the maximum due to the additional numerical factor √ ln(2) ≈ 0.833. Lastly, for model III the numerical value of the coefficient of variation of i.e. is again somewhat smaller than the corresponding value v (III) M for the maximum, equation (41). We thus conclude that distributions of the range in all three models under study are narrower than the corresponding probability density functions of the maxima.

Typical behaviour of the probability density function of the maximum
As yet, our discussion of the averaged versus a typical behaviour concerned only the maximum and the range themselves. The results for P T (M) and P T (R) we have presented so far correspond to a standard way of performing the averaging, i.e. when the averaging is first performed with respect to thermal histories at a fixed realisation of the stochastic diffusivity process B t and then over all possible realisations of B t . Conversely, realisation-dependent distributions P T (M) and P T (R) are evidently random functions themselves which fluctuate from one realisation of B t to the next, and it is, of course, not clear a priori to which extent their first moments, i.e. precisely P T (M) and P T (R), are representative of the actual behaviour of these properties. In principle, it may well happen that P T (M) and P T (R) are supported by some atypical, rare realisations of B t which nonetheless provide a dominant contribution to their values. If true, in order to observe our predictions for P T (M) and P T (R) one may need very large statistical samples. Note that in figure 5 we presented a convincing evidence for the predicted functional forms but the number of realisations used to perform the averaging was sufficiently high, 10 4 . That may not be the case for experimental studies for which such a large number is beyond reach.
Following the analysis of a typical kinetic behaviour in the so-called target problem with respect to fluctuations in the starting points of searchers (see reference [74] and the recent reference [75]), we concentrate on the properties defined in (9). Here, one first performs an averaging over thermal histories at a fixed realisation of B t and then averages the logarithm of the realisation-dependent probability density over all possible realisations of B t . Because the logarithm is a slowly-varying function of its argument, one expects that its averaged behaviour is rather insensitive to rare anomalous realisations and is thus representative of a typical behaviour which should be observed for a majority of trajectories B t , or seen for small statistical samples. The resulting expression is then exponentiated to produce an estimate of typical distributions. We will be concerned here only with the typical behaviour of the distribution of the maximum-the analysis of the typical distribution of the range appears to be somewhat more involved and lengthy, but we do not expect any significant new features, as compared to the behaviour of the maximum.
Recalling that for any given realisation of B t , the probability density function P T (M) is given by (2) with and hence, Integrating the latter expression over M we find that the normalisation is given by such that, eventually, we have the following estimate of the typical behaviour Therefore, we arrive at the conclusion that if the first inverse moment of the random variable D 0 The first inverse moment of D 0 T 0 V(B t )dt can be calculated straightforwardly by simply integrating Υ(T; λ) in (8) over λ from zero to infinity. In doing so, we realise that for model I this negative moment does not exist, because Υ(T; λ) decays as 1/ √ λ in the limit λ → ∞ (see (55)) and hence, the integral diverges at the upper integration limit. On the other hand, the average of 1/(D 0 T 0 V(B t )dt) over any finite statistical sample of trajectories B t is evidently finite and hence, in virtue of (85) for such samples P The first inverse moment of D 0 T 0 V(B t )dt is finite for both models II and III. For model II we have where the integral in the right-hand side is finite for any finite value of the parameter D B T/4a 2 . The latter integral cannot be performed exactly but its behaviour can be readily understood by noticing that for any z 0 we have 2 π z coth πz 2 z + 2 π .
As a consequence, we find that the first inverse moment of D 0 T 0 V(B t )dt obeys the following two-sided inequality 1 aD 0 In the limit T → ∞ these bounds become sharp and hence, define the leading behaviour of the first inverse moment exactly. For model III the first inverse moment of D 0 T 0 V(B t )dt has the simpler form Comparison of our analytical predictions for the ensemble-averaged and the typical behaviours for models II and III against the histograms obtained from Monte Carlo simulations with just 100 realisations of the stochastic process B t is presented in figure 6. Here, for P T (M)-the former predicts higher values of the probabilities while the latter underestimates them: a trend that is confirmed by our numerical observations. This is not the case for larger values of M. Perhaps somewhat surprisingly, the heavy log-normal tail of P (II) T (M) appears in a good agreement with numerics even for such a moderately small sample size for values of M as large as 10 3 . In turn, for model III there is no significant difference between P T (M) seems to be more accurate.

Conclusion
Deviations from standard Brownian motion have been measured in a vast range of systems, starting with Richardson's cubic law for the relative diffusion of tracers in turbulent media in 1926 [87]. Such 'anomalous diffusion' has given rise to a rich variety of statistical models accounting for various physical aspects effecting deviations from standard Brownian motion [88][89][90]. As a particular case, random diffusivity models were introduced in the context of the modelling of complex measured NMR signals [91]. The randomness of the diffusivity can be assumed to be due to an inhomogeneous particle ensemble in a homogeneous environment, or due to identical particles in a heterogeneous environment. When the diffusivity distribution is fixed in time the dynamics resulting from such random diffusivities is captured by the framework of superstatistics [92] or grey Brownian motion [93]. When particles move in quenched environments with finite patch sizes and specific jump rules interesting dynamic effects and non-Gaussian phenomena have recently been revealed [94,95]. Originally devised to reproduce the observed crossover behaviour from non-Gaussian to Gaussian displacement statistics in systems showing a Brownian scaling of the mean squared displacement [96,97], diffusive processes with stochastically evolving diffusivities were devised as an 'annealed' approach to the motion of the test particle in a heterogeneous environment. Such diffusing diffusivity models with stationary diffusivity dynamics were analysed in terms of the mean squared displacement and the displacement distribution. Despite the different formulations several core features turn out to be robust among these models [54,55,[57][58][59]98].
Here we studied a stochastic process x t driven by white Gaussian noise, whose amplitude is being modulated by a stochastically varying diffusivity D t for three different choices: (I) cut-off Brownian motion B t with D t ∼ Θ(B t ); (II) geometric Brownian motion, D t ∼ exp(−B t ); and (III) squared Brownian motion, D t ∼ B 2 t . In contrast to the above-mentioned diffusing diffusivity models, that are all Brownian, the three choices here effect non-stationary diffusivity dynamics, and the resulting diffusion exhibits both normal and anomalous diffusive scaling.
In the analysis we concentrated on the extremal properties in terms of the maximum and the range of these three random diffusivity models. We obtain analytical expressions for the probability density functions of the maximum and the range of the processes for a given observation time. Our discussion reveals both similarities and differences of the extremal properties of these models among each other as well as compared to standard Brownian motion. In particular, we unveil that model I shows significant differences from Brownian motion while the small-maximum limit of model II coincides exactly with the Brownian behaviour. We also show that the distributions of the maximum are generically broader than the distributions of the range, as evidenced by the analysis of the coefficients of variation of the corresponding distributions. Our discussion finally unveils the difference between the ensemble and the typical behaviour of the probability density functions, an important ingredient for the analysis of finite-sized data sets.
The analysis of given stochastic time series representing a set of trajectories of diffusing test particles has more recently received considerable attention. A number of statistical observables has been introduced and discussed (see, e.g. [89,90,99,100]) to allow the physical classification of recorded data. For instance, it has been shown how to use Bayesian maximum likelihood [101,102] or machine learning [103,104] to classify a measured system and extract its physical parameters. Specifically, the power spectral analysis of single, finite-length trajectories was shown to distinguish different forms of random diffusivity models [53]. A more recent twist on data analysis of stochastic processes uses large-deviation approaches, for instance, for the time averaged mean squared displacement [105,106].
While it is not surprizing that the EV behaviour encoded in the probability density functions of the maximum and the range studied here was shown to distinguish the three, quite different, random diffusivity models investigated here, we also demonstrated that the rectified Brownian motion of model I exhibits significant differences to standard Brownian motion. It should therefore be interesting to investigate whether these two distributions allow one to distinguish between the diffusing diffusivity models encoding Brownian yet non-Gaussian motion [54,55,[57][58][59], and how these measures change for projections of higher dimensional versions of these models.