Transitions in the ergodicity of subrecoil-laser-cooled gases

With subrecoil-laser-cooled atoms one may reach nano-Kelvin temperatures while the ergodic properties of these systems do not follow usual statistical laws. Instead, due to an ingenious trapping mechanism in momentum space, power-law-distributed sojourn times are found for the cooled particles. Here, we show how this gives rise to a statistical-mechanical framework based on infinite ergodic theory, which replaces ordinary ergodic statistical physics of a thermal gas of atoms. In particular, the energy of the system exhibits a sharp discontinuous transition in its ergodic properties. Physically this is controlled by the fluorescence rate, but more profoundly it is a manifestation of a transition for any observable, from being an integrable to becoming a non-integrable observable, with respect to the infinite (non-normalised) invariant density.

Laser cooled atoms are important for fundamental and practical applications. It is well known that Lévy statistics describes the unusual properties of the cooling processes [1,2,[4][5][6][7][8]. For subrecoil laser cooling a special atomic trap in momentum space is engineered. The most efficient cooling is found when a mean trapping time, defined more precisely below, diverges [1]. The fact that the characteristic time diverges implies that the ergodic properties of these systems must differ from those of standard gases [1,9,10]. Ergodicity is a fundamental aspect of statistical mechanics that implies that the time and ensemble averages coincide. This is found when the measurement time is made long compared to the time scale of the dynamics. However, in the context of subrecoil laser cooling this time diverges, and hence no matter how long one measures, deviations from standard ergodic theory are prominent. Given that lasers replace heat baths in many cooling experiments, what are the ergodic properties of the system? In other words, what replaces the usual ergodic statistical framework? Our goal is to show how tools of infinite ergodic theory describe the statistical properties of the ensemble and corresponding time averages of the subrecoil-laser-cooled atoms.
Infinite ergodic theory was investigated by mathematicians [11][12][13] and more recently in Physics [14][15][16][17][18][19][20][21][22][23][24][25][26][27]. In generality, infinite ergodic theory deals with a peculiar non-normalised density, describing the long time limit of a system, called below the infinite invariant density. Previous works in the field of subrecoil laser cooling [1,2] foresaw this quasi-steady state. We will see how to use this tool to investigate the ensemble averages of physical observables. However, this does not give direct information on the time averages and here we will develop physical and mathematical insights on the latter. The basic question is how to relate ensemble and time averages, even when ergodicity in its standard formulation is broken. In particular we investigate the energy of the system. Since the atoms are non-interacting, in a classical thermal setting the energy of the atoms per degree of freedom would be k B T /2, as a consequence of Maxwell's velocity distribution. At variance with this, we will show, that the energy of a subrecoil-laser-cooled gas is obtained under certain conditions with a non-normalisable invariant density. A sharp transition is exposed, in the statistical properties of the energy, when the fluorescence rate, given by R(v) ∝ |v| α in the vicinity of zero velocity, is varied, more precisely when α = 3. Since experimental work demonstrates the capability of a variation of α from α = 2 to α = 4 or, in principle, α = 6 etc, the rich phase diagram of ergodic properties which we find seems to be within reach of experimental investigations. This new type of transition is related to the fact that the energy observable can switch from being an integrable observable, with respect to the infinite density, to being non-integrable.
Let ρ(v, t) be the probability density function (PDF) of the speed v > 0 of the atoms at time t. A master equation governs its evolution with typical gain and loss terms Here R(v) = 1/τ (v) is the fluorescence rate and f (v) denotes the speed PDF after the atom experiences a spontaneous jolt. The specific forms of these functions, investigated previously [1], are crucial for our analysis, but now we treat the problem in generality. It is arXiv:2104.03816v1 [cond-mat.stat-mech] 8 Apr 2021 natural to consider the long time limit of the PDF ρ(v, t). For this aim, as usual, we consider the steady-state condition ∂ρ(v, t)/∂t = 0. This gives the time-independent solution ρ * (v) = τ (v)f (v)/Z, where Z, if it exists, is a time independent constant obtained from the normalization condition. We treat the opposite non-normalizable case which is highly relevant.
Indeed, for subrecoil laser cooling, once the atom is close to a halt, the rate of change of speed is becoming small, and that drives the system to low temperatures, e.g. nano-Kelvins. More specifically, when the fluorescence time τ (v) is given by τ (v) ∼ cv −α for v → 0, then ρ * (v) becomes non-normalizable, for α > 1. Importantly, as we will shortly explain, the non-normalizable function ρ * (v) is not an object we should ignore. An analysis of the master equation (see SM) and following the footsteps of [1,2] yields For the case under study, γ < 1 or α > 1, the integration over I(v) diverges, due to the small v behavior of τ (v), and hence I(v) is called an infinite invariant density. Note that in this case, on the left-hand side of Eq. (2), we multiply the normalized density ρ(v, t) with a function Z(t) that is increasing with time, therefore it is not surprising that on the right-hand side we find in the long-time limit a nonnormalizable function. The question we tackle is what is the physical meaning of this state? And how do we use the infinite invariant density to construct a non-trivial ergodic theory for the gas?
One can argue that the infinite density, as a stand alone, cannot be the full picture, as it is nonnormalizable. Indeed, for long but finite times, the density of the velocity of the particles can be described by two regimes. As the density is evolving, it is shrinking in its width peaking on zero speeds. In the inner region, the density ρ(v, t) exhibits a scaling solution, given by [1] This describes velocities of order 1/t γ , i.e. an inner region of the packet, which for t → ∞ goes to zero (see SM). The scaling solution Eq. (3) and the infinite density Eq. (2) are not separable and both together yield a complete description of the packet. Mathematically, the two solutions match at intermediate velocities (see SM). Eq.
(3) predicts that the full width at half maximum of the velocity PDF decays like t −1/2 for α = 2 and as t −1/4 for α = 4. These theoretical predictions were indeed observed in the laboratory with Cesium [29,30], see also quantum Monte Carlo simulation and experiments with Helium in [1,10]. However, the scaling solution Eq. (3) exhibits what might seem naively as an unphysical feature. If we consider the realistic case α = 2, we find that the scaling function as a stand alone predicts that the second moment of the velocity, namely the kinetic energy, is infinite, and in this sense the system cannot be considered cold. This is due to the fat tail of the scaling function for large v. As explained below, this issue is cured, using the non-normalized solution Eq. (2). In general, we need to classify observables based on whether they are integrable or non-integrable with respect to I(v). These two classes of observables have vastly different ergodic properties, which are now investigated.
Ensemble and time averages. Let v(t) be the random speed process of an atom. We now consider generic observables O[v(t)], and we will study their time and ensemble , which is the kinetic energy when we set m/2 = 1 and the indicator function (2), we find in the long-time limit This equation shows that the non-normalized density I(v) is used in the computation of ensemble averages, in a way reminiscent of the standard averaging performed for equilibrated systems with normalized densities. The only condition we have is that the observable is integrable with respect to I(v), namely that the integral exists. Since I(v) ∼ v −α for small velocity, the kinetic energy is an integrable observable for α < 3. However it is non-integrable otherwise. The critical value α = 3 or γ = 1/3 will mark an ergodicity transition for this observable. In contrast, the indicator function is always an integrable observable, provided that v a > 0. A goal of ergodic theories is to relate the ensemble and the time averages, denoted with an over-line, According to the standard ergodic hypothesis O(t)/ O → 1 in the limit of longmeasurement-times. In our case, the time averages remain random, and we will soon investigate their fluctuations. To start, we consider an ensemble of paths and average over time and then over the ensemble (5) Considering the long-time limit and using Eq.
where we used 0 < γ < 1. This resembles the standard calculation of the time averages, performed when the invariant density is normalizable. We conclude that for integrable observables and if γ < 1 Thus we established a relation between the time average and the ensemble average. The latter is obtained using the infinite density I(v) and thus this invariant density is not merely a tool for the calculation of the ensemble average, but rather it gives also information on the time average. We see that when γ → 1 we approach standard ergodic behavior. Physically this is related to the observation made below, that the mean time between velocity updates diverges for γ < 1.
The time averages are functionals of the stochastic path v(t) and we now develop a machinery to explore their statistical properties. Here we will focus on the kinetic energy observable, due to its physical importance. To analyse the latter we recall the stochastic process under study. Initially we draw the speed v from f (v), the particle then remains in this state for a random time denotedτ . This timeτ is exponentially distributed with the life time τ (v) that depends on the velocity. The process is then renewed, namely after waiting for a timeτ we draw a new velocity from f (v) and hence a new life time etc. In simulations below, following [1], we use a uni- This holds under two conditions, namely that the energy is an integrable observable, so 1/3 < γ, and γ < 1, otherwise we are in the regime of standard ergodic theory. Interestingly, when we take γ → 1/3 the prefactor in Eq. (8) diverges and becomes v max independent. This marks the entry into the phase where the energy is non-integrable. After the calculation of the expectation of the time averaged energy with the infinite density I(v), the real challenge is to obtain the distribution of E k .
We rewrite the time average E k (t) = S(t)/t, where S(t) is the action. The key idea is to investigate the distribution of the action, and with this to infer the ergodic properties of the process. We perform this task with a new form of coupled continuous-time random walks, which in turn is a variation of the well known Lévy walk [31][32][33][34][35][36]. In the time interval of observation (0, t) we have N (t) random renewal events, and the pairs of waiting times and velocities are labeled (v i ,τ i ), where i = 1, · · · , N (t) (v 1 is the initial condition). We rewrite the action S(t) = , exhibits a wide range of physical behaviors. When γ = 1/α > 1 ergodicity holds and the PDF approaches a delta function in the long-measurement-time limit, see finite time simulations in sub-plot a with α = 0.8 (details in SM). When α = 1.25, the energy is integrable with respect to the infinite measure, and P (Υ) is non-trivial though it has a peak close to unity, since at this stage we are not too far from the ergodic phase (plot b). For α = 2 and α = 4 we find that the fluctuations of the time average, are described by a half-Gaussian distribution, presented in c. Here the energy observable is integrable in the first case while it is not in the second, further note that the Mittag-Leffler distribution is half-Gaussian when α = 2 since γ = 1/2. When α = 6 the energy observable is nonintegrable and P (Υ) diverges for Υ → 0, indicating very large deviations from usual ergodic behavior (plot d). In the limit α → ∞ we get a discontinuous behavior, a sharp cutoff at Υ = 3 (see plot e, and simulation with α = 50) . Finally, the EB parameter, subplot f , exhibits a cusp at γ = 1/α = 1/3. This value marks a transition for the ergodic behavior of the system, from an energy which is integrable with respect to the infinite measure to non-integrable. In the figure simulations match theoretical predictions without fitting.
is reminiscent of a biased random walk process. The increments s i > 0 are constrained since the measurement Here t B is the so called backward recurrence time [37,38], the time elapsing between the last update in the process and the measurement time t. Similarly s B (t) = (v N (t)+1 ) 2 t B (t) is the contribution to the action, from the last time interval in the sequence.
To advance the theory we need the joint PDF of action increments s and waiting timesτ denoted φ(s,τ ). This is obtained from when 0 ≤ s ≤ v 2 maxτ , and zero otherwise. Here the waiting timesτ and action increments s are clearly correlated. Integrating over s we find the PDF of the waiting times ψ(τ ) whose fat tail reads ψ(τ ) ∼ constτ −1−γ with const = c γ γΓ(1 + γ)/v max . As mentioned, the divergence of the mean waiting time found for 0 < γ < 1 signals special ergodic properties [9].
Let P (S, t) be the PDF of the action at time t. The first step in the analysis is to relate this PDF to Eq. (9). Using techniques borrowed from random walk theory, employing the renewal property of the process and the convolution theorem, we establish this relation using Laplace transforms. LetP (u, p) = ∞ 0 ∞ 0 dSdt exp(−uS −pt)P (S, t) be the double Laplace transform S → u and t → p of the PDF P (S, t). Then we derive a Montroll-Weiss [39] like equation which readŝ Hereφ(u, p) is the double Laplace transform of φ(s,τ ) and similarly for the pairΦ(u, p) and (11) This term stems from the contribution to the action s = v 2 t B from the last increment in the sequence, namely from the backward recurrence time t B , while exp(−R(v)t B ) is the probability of not jolting in that time interval. We now promote the use of the Montroll-Weiss like tool Eq. (10) in the context of ergodicity, following the next steps.
To investigate the ergodic properties of the process we define the dimensionless random variable Υ = E k / E k . We first focus on the case when the observable is integrable, namely 1/3 < γ < 1. In standard ergodic theories, found here if 1 < γ, the distribution of Υ will approach a delta function centred on unity. By defintion Υ = S(t)/(t E k (t) ), and here an analysis of Eq. (10) is useful, since it yields the distribution of S and then the distribution of the sought after Υ. We find the universal behavior that the PDF of Υ denoted ML(Υ) is given by the Mittag-Leffler (ML) law Here l γ,1 (.) is the one sided Lévy PDF. This law which replaces Birkhoff's ergodic theory, is a concrete manifestation of the Darling-Kac theorem [11,12]. Our physical approach was to show, how this law is related to the laser cooling parameters, i.e. to γ. This law is valid for any observable of interest, provided that it is integrable, for example we verified this numerically for the energy observable Fig. 1 but also for the indicator function. In the limit γ → 1 ML(Υ) reduces to a delta function, as expected.
We now address the case when the energy in a nonintegrable observable, namely 0 < γ < 1/3. The calculation of the mean energy cannot be made with the infinite density, since the result will diverge. Instead, here the ensemble average kinetic energy is found using the scaling solution Eq. (3). This means that the energy is now sensitive to the inner part of the velocity packet This behavior is universal in the sense that this result does not depend on v max , unlike Eq. (8). In this case, a simple time integration gives the relation between time and ensemble average E k (t) ∼ E k (t) /(1 − 2γ) which clearly differs from the generic behavior found for integrable observables Eq. (7) .
To characterize the fluctuations of the time averages it is useful to define the ergodicity breaking parameter [40] which is zero in the long-time limit, if the process is ergodic, namely for γ > 1. In the phase where energy is integrable EB = [2Γ 2 (1 + γ) − Γ(1 + 2γ)]/Γ(1 + 2γ), for 1/3 < γ < 1, which as mentioned is universal in the sense that it goes beyond the observable of interest. For example this describes also the fluctuations of the number of renewals N (t).
Returning to the case γ < 1/3, we need to find the second moment of the action, which is obtained in Laplace p space after differentiation of Eq. (10) This expression is analysed in the small p limit, and we find asymptotically for large times and for 0 < γ < 1/3 This expression clearly differs from the EB parameter found in the Darling-Kac phase 1/3 < γ < 1. When γ → 0 we find EB = 4/5. In this limit the particle maintains a constant speed for (nearly) all the observation time so E k = v 2 and hence using Eq. (3), lim γ→0 EB = v 4 − v 2 2 / v 2 2 = 4/5. The EB parameter versus 0 < γ < 1 is plotted in Fig. 1f and it exhibits a clear transition when γ = 1/3. Thus switching from the case when energy is integrable to non-integrable, manifests itself in nontrivial fluctuations of the time averages.
We have investigated semi-analytically the distribution of the time averages also in the non-integrable phase γ < 1/3 namely α > 3. Here, clearly the Mittag-Leffler law is not valid any more. In Fig. 1 we present some of the main results of this mathematically challenging domain. For example, we find lim γ→0 P (Υ) = Υ −3/2 /(2 √ 3) for Υ < 3, otherwise P (Υ) = 0. This diverges for small Υ and also has a non-analytical cutoff at Υ = 3. This result can be explained, as we did for the EB parameter, noting that the atom maintains a constant speed for practically all the duration of the measurement i.e. in this limit Υ = v 2 / v 2 and the PDF of v is the uniform f (v). For the experimentally relevant case α = 4 (γ = 1/4), we find that the distribution of Υ > 0 is half-Gaussian, see Fig. 1(c). Interestingly this case also marks a transition: the PDF of Υ diverges at the origin for any α > 4 and vanishes there for 3 < α < 4. This means that for α > 4 the most likely time average is actually found for cases where it is much smaller than the ensemble average. Finally, for α = 6, a case also considered relevant for experiments, we find the solution in terms of a Fox-H function [4], P (Υ) , which perfectly matches the simulation presented in Fig.  1(d).
We have seen how the non-normalizable infinite invariant measure I(v) replaces in some sense the Maxwellian velocity distribution for thermal gases, which is of course also invariant but perfectly normalizable. We showed how the ergodic properties of the system exhibit three phases controlled by the fluorescence rate parameter γ. The divergence of the mean trapping time at γ = 1 is a well known factor in the change of ergodic properties of a vast number of physical systems. We highlighted a second novel transition, associated with the changeover of an observable from being integrable to non-integrable, with respect to the non-normalised state I(v), found at γ = 1/3 for the energy observable. Similar behaviors can be found for other observables and as the applications of infinite ergodic theory expand, either for stochastic or deterministic systems, we expect this type of transition to be wide spread.
Acknowledgments: The support of Israel Science Foundation's grant 1898/17 is acknowledged (EB). This work was supported by the JSPS KAKENHI Grant No 240 18K03468 (TA). GR thanks Tony Albers for helpful discussions.

SUPPLEMENTARY MATERIAL The Master Equation
We briefly discuss the derivation of the non-normalized quasi steady-state I(v), Eq. (2) in the main text, and the scaling solution g(x), Eq. (3), see also [1,2]. Let ρ(v, t) be the distribution of the speed of the atom. As mentioned in the text, after a recoil event, the parent PDF of v is f (v). As is well known, this does not imply that the steady-state velocity is given by f (v) since the fluorescence rate R(v) depends on v. Note that in simulations, in the text and following [1] f (v) was uniform, however, a variation will not change dramatically the main conclusions of our work. Following [1] the transition from state v to v is given by where the first term describes the losses and the second the gains. Here τ (v) = 1/R(v) is the life time of an atom in state v. A time independent solution of this equation reads ρ * (v) = Constτ (v)f (v). If this solution is normalizable, we can determine the constant, and then we get the usual steady-state of the system. However, we consider the opposite situation, when τ (v) ∼ cv −α for v → 0 and α > 1 and, as mentioned in the text, in current experiments α = 2, 4, 6 are relevant. In this model the life time τ (v) becomes long when the velocity is small, and hence this describes trapping of atoms at small velocities, which leads eventually to efficient cooling. To solve the Master equation in the long-time limit, we invoke two types of solutions, which we later match. For long times and not too small v we have The numerator has the structure of the usual steady-state and b and ξ are obtained by matching. On the other hand for small v, of the order of t −1/α we have a scaling solution ρ(v, t) ∼ t 1/α g(vt 1/α ). Inserting this ansatz into the Master equation, we obtain an integral equation for g(x). Following this procedure and matching the two solutions one arrives at Eqs. (2,3) in the main text, which are the starting point of our work. We extensively verified these solutions with numerical simulations. Specifically one finds the exponent ξ = 1/α [1].

Figure preparation
To prepare plots, we simulate the process, using an initial condition for the speed v drawn from a uniform distribution in the interval (0, 1) so v max = 1. The waiting time is a random variable following the exponential distribution with mean τ (v). In numerical simulations, a uniform random variable x on [0, 1] is used and we transform x into y = −τ (v) ln x which gives the random time between collision updates. The speed of the atoms is fixed between these events. The next velocity and waiting time are determined in the same way as the above. Recall that τ (v) = cv −α and in simulations we use c = 1, which does not alter the main conclusions of the paper.
We now provide further details of the figure.
• Fig. 1(a). PDF P (Υ) of the normalized time averaged energy Υ = E k (t)/ E k (t) is plotted for α = 4/5. The measurement time is t = 10 9 . In the limit of infinite measurement time, this PDF will approach a delta function, since here we are considering the normal ergodic phase γ > 1 (because γ = 1/α = 5/4). The number of trajectories used was 10 4 . To obtain Υ, we calculate Y (t) = S(t)/t, where S(t) is the action defined in the text. This is performed for each trajectory and the mean is obtained using an average over the ensemble of trajectories. Y (t) divided by the mean gives the random variable Υ. Repeating many times we find the histogram for Υ.
• Fig. 1(b). Same as the previous plot, the PDF P (Υ) of the time averaged energy v 2 (t), but now α = 5/4 and hence γ = 4/5, the measurement time is t = 10 6 and the number of trajectories 10 6 . Since 1/3 < γ < 1 the theoretical PDF of Υ is the Mittag-Leffler distribution with index 4/5, see Eq. (12) in the Letter. As explained in the text, one sided Lévy stable functions l γ,1 (x) are used to plot this PDF. Specifically l γ,1 (x) is by definition the inverse Laplace transform of exp(−u γ ), and x > 0. In turn these functions are implemented in programs like Mathematica, hence it is easy to plot the theoretical prediction for P (Υ). As shown in the figure, this perfectly matches the simulations without fitting.
• Fig. 1(c). PDF P (Υ) for α = 2 and α = 4. The measurement time is fixed as t = 10 11 , the number of trajectories is 10 6 . For α = 2 note that the Mittag-Leffler PDF, defined in the text, ishalf-Gaussian. While the calculation of the PDF of Υ for α = 4, which also gives a half-Gaussian distribution, is based on the Mellin transform (see below). As in other sub-plots theory and simulations perfectly match.
• Fig. 1(d). Same as the other plots however now α = 6, t = 10 15 and the number of realizations 10 6 . The results of simulations, presented as a histogram, perfectly match the theory given in the text. As we will soon explain, the latter is expressed in terms of a Fox H-function, which gives P (Υ) = C Γ( 3 4 ) H 1,0 1,1 CΥ ( 1 3 , 2 3 ) (− 1 4 , 1) with C = 3/(4Γ(5/3)). After rewriting this expression in terms of a Meijer G-function we can plot the result using Mathematica.
• Fig. 1(e). Same as the other plots, now for simulations with α = 50. This in turn is compared with the α → ∞ limit of our theory presented in the main text.
• Fig. 1(f ). EB parameter for the energy observable v 2 (t). The process is the same as in other figures though now we vary γ, namely we present EB versus γ. The EB parameter is defined in Eq. (14) in the main text. The number of the trajectories used to produce the figure is 10 6 except for γ > 1 where it is 10 5 . The method of calculating the time average is the same as in the previous plots. The theoretical predictions, plotted in the figure, are given in the text and they read where Γ(.) is the Gamma function. Here we have three phases: energy is non-integrable with respect to the infinite measure 0 < γ < 1/3, the energy is integrable 1/3 < γ < 1 (in the Letter this was called also the Darling-Kac phase), and finally 1 < γ where standard ergodic theory holds and the EB parameter is equal zero.

DETAILS ON THE DERIVATION OF P (Υ)
We provide some details on the derivation of P (Υ) focusing on the phase where the energy observable is nonintegrable with respect to the infinite measure. The calculations are lengthy, and hence here we provide only an outline. We plan a longer publication which will provide further technical details.
The double Laplace transform of the action propagator P (S, t) is given by P (u, p) in Eq. (10) in the Letter. In the limit t → ∞, and equivalently p → 0, these quantities take corresponding scaling forms, which inserted into the Laplace transform relation yields where f α (x) and g α (y) are the scaling function in action space and in the corresponding Laplace space, respectively. Substituting x = S t 1−2/α and setting p = 1 turns Eq. (20) into the integral equation relating g α (y) and f α (x) via this convolution transform with kernel For the scaling function g α (y), with α > 3, i.e. in the non-Darling-Kac phase, we obtain from the p → 0 limit of P (u, p), Eq. (10) in the text, the following exact form The goal is to obtain from this by inversion of Eq.
While in principle along these steps the problem of calculating the limit distribution P α (Υ) is solved, a fully analytic solution is available only in two cases, α = 4 and α = ∞. For α = 4 the integrals in Eq. (23) can be evaluated by residue calculus yielding after some calculations the simple result g 4 (y) = 1 1 + y , with Mellin transform g 4 (s) = Γ(s)Γ(1 − s). Since the Mellin transform K α (s) of the integral kernel is also known in full generality, K α (s) = Γ(s)Γ(1 − (1 − 2/α)s), we get the quotient in Eq. (25), and in addition we can invert from Mellin space to obtain f 4 (x) = 1 √ π e − x 2 4 . The scaling factor C 4 = x f4 follows from the general relation between the n-th derivative g (n) α (y = 0) and the moments x n fα of f α (x) g (n) α (0) = K (n) α (0) x n fα = (−1) n Γ(1 + n(1 − 2 α )) x n fα , which follows directly from Eq. (21). For α = 4 we get C 4 = x f4 = 2 √ π , so that we obtain for P 4 (Υ) according to Eq. (26) a half Gaussian distribution as exact result, which is displayed together with the results from numerical simulations in Fig. 1 c. We can proceed similarly for α → ∞, because all integrals and Mellin transforms are known exactly also in this case yielding and eventually which is also plotted in Fig. 1e. For the experimentally also relevant case α = 6 we can still get from Eq. (23) by residue calculus a fully analytic expression for g 6 (y), but the result is very lengthy involving 3rd roots etc., and it does not simplify as in the case α = 4. Therefore one cannot analytically calculate its Mellin transform. This led us to find an approximate scaling function g ≈ 6 (y), which deviates only little from the true function g 6 (y), but for which the Mellin transform is known analytically. The general ideas is to match in g ≈ 6 (y) the true small y-behavior and the large y-asymptotics of g 6 (y), which we know analytically from an analysis of Eq. (23). The simple form shares with the exact function g 6 (y) the identical first and second derivative at y = 0, and the exponent of the asymptotic behavior for y → ∞. The relative deviation of g ≈ 6 (y) from g 6 (y) is strictly less than 4.2%. This function g ≈ 6 (y) can be Mellin transformed [4], and leads via Eq. (25) and subsequent rescaling to an analytical expression for P ≈ 6 (Υ), which can be expressed in terms of a Fox-H function as .
As mentioned, after rewriting this expression in terms of a Meijer G-function [ [4], p.629], we can plot the result with Mathematica, as shown in Fig. 1 d, with no visible deviations from the result obtained by numerical simulations of the process.