#  MFDFA: Efficient Multifractal Detrended Fluctuation Analysis in Python

Leonardo Rydin Gorjão  Forschungszentrum Jülich, Institute for Energy and
Climate Research - Systems Analysis and Technology Evaluation (IEK-STE), 52428
Jülich, Germany  Institute for Theoretical Physics, University of Cologne,
50937 Köln, Germany  Galib Hassan  Forschungszentrum Jülich, Institute for
Energy and Climate Research - Systems Analysis and Technology Evaluation (IEK-
STE), 52428 Jülich, Germany  Institute for Theoretical Physics, University of
Cologne, 50937 Köln, Germany  Jürgen Kurths  Potsdam Institute for Climate
Impact Research, 14473 Potsdam, Germany  Institute of Physics, Humboldt
University Berlin, 12489 Berlin, Germany  Centre for Analysis of Complex
Systems, World-Class Research Center “Digital biodesign and personalised
healthcare”, Sechenov First Moscow State Medical University, 119991 Moscow,
Russia  Dirk Witthaut  Forschungszentrum Jülich, Institute for Energy and
Climate Research - Systems Analysis and Technology Evaluation (IEK-STE), 52428
Jülich, Germany  Institute for Theoretical Physics, University of Cologne,
50937 Köln, Germany

###### Abstract

Multifractal detrended fluctuation analysis (MFDFA) has become a central
method to characterise the variability and uncertainty in empiric time series.
Extracting the fluctuations on different temporal scales allows quantifying
the strength and correlations in the underlying stochastic properties, their
scaling behaviour, as well as the level of fractality. Several extensions to
the fundamental method have been developed over the years, vastly enhancing
the applicability of MFDFA, e.g. empirical mode decomposition for the study of
long-range correlations and persistence. In this article we introduce an
efficient, easy-to-use python library for MFDFA, incorporating the most common
extensions and harnessing the most of multi-threaded processing for very fast
calculations.

Software: https://github.com/LRydin/MFDFA

##  I Introduction

A common tool to unveil the nature of the scaling and fractionality of a
process, natural or computer-generated, is Multifractal Detrended Fluctuation
Analysis (MFDFA). It was initially developed by Peng et al. [1, 2] as basic
Detrended Fluctutation Analysis (DFA) and later extended to study multifractal
processes by Kandelhardt et al., giving rise to MFDFA [3]. It addresses the
question of the presence of correlations in time series and can be employed to
analyse both discrete as well as continuous-time stochastic processes. Since
its initial development in the late 90’s, it has been revisited to incorporate
several other elements, e.g. empirical mode decomposition as a method for
detrending [4, 5, 6, 7], overlapping moving windows [8, 9], and a new metric
denoted extended detrended fluctuation analysis [10, 11, 12, 13]. There are
several additional features exist, designed to study correlations of two or
more time series [14, 15], lag correlations in time series [16], and Fourier-
DFA [17], amongst others. A comprehensive study of DFA and the interplay
between trends in data and correlated noise can be found in Ref. [18]. MFDFA
has found application in various fields, such as the analysis of heartbeat
rate [19], arterial pressure [10], EEG sleep data [11, 13], physiology [20],
keystroke time series from Parkinson’s disease patients [21], cosmic microwave
radiation [22, 23], seismic activity [24, 25], sunspot activity [26],
atmospheric scintillation [27], temperature variability [28], meteorology
[29], precipitation levels [30], streamflow and sediment movement [31, 32, 33,
34, 35, 36, 7], protein folding [37], finance and econophysics [38, 39, 40,
41, 42], electricity prices [43, 44], power-grid frequency [45, 46],
epidemiology [47], music [48, 49, 50], ethology [51, 52], multifractal
harmonic signals [53], and microrheology [54].

MFDFA is a numerical algorithm designed to determine the self-similarity of a
stochastic process. Putting it simply, the algorithm examines the relation
between the diffusion of the process and its propagation in time or space.
Auto-regressive and stochastic processes with different power-law scaling will
diffuse with different rates. Fluctuation Analysis (FA) provides a method to
uncover these correlations, but fails in the presence of trends in the data,
which, for example, are particularly present in weather and climate data.
Detrending the data via polynomial fittings (DFA) allows one to uncover solely
the relation between the inherent fluctuations and the time scaling of a
process, thus circumventing the impact of non-stationarity in the data.
Likewise, other methods—as empirical mode decomposition or moving average
windows—are viable options to detrend the data. Another problem is that a
process might be driven by more than one time scale, i.e., have more than one
internal period, which can be removed either with local polynomial fittings or
EMD. Moreover, a stochastic process might be of a monofractal or multifractal
nature. By studying a continuum of power variations of DFA one extends into
MFDFA, which permits the study of the fractality of the data by comparing
power variations, i.e., a multifractal spectrum.

In this software we sought to design a computationally efficient code focused
on computational speed and usability. There are currently no flexible and
available implementations of MFDFA in python. Available are some MATLAB [55]
as well as R packages [56, 57]. There is a particularly thorough introductory
guide to MFDFA in MATLAB with a source-code by Espen A. F. Ihlen [55], which
is easy to implement but numerically inefficient. With this implementation
efficiency was sought. This was achieved by making the most out of python,
reshaping the code to allow for multi-threading, especially relying on numpy’s
polynomial, which scales easily with modern computers having more processor
cores [58]. Moreover, this library contains the most commonly applied methods
alongside with DFA and MFDFA: the added feature of empirical mode
decomposition is implemented to substitute the polynomial fittings; A moving
window is included, especially valuable for shorter time series; The extended
DFA (eDFA) method is also included, adding a second metric of fractal scaling,
especially valuable for multifractal or aperiodic time series.

In the following sections we will introduce MFDFA alongside some of the
aforementioned methods incorported into the MFDFA library. We will present two
classical applications, one with a monofractal process and one with a
multifractal noise, and show how to use MFDFA to extract their characteristics
from a single one-dimensional time series. Python code is presented to
explicate the use of the MFDFA library. Subsequently we study two real-world
time series: the sunspot time series from 1818 to 2020 which accounts for the
daily recorded sunspots and the quarter-hourly electricity trading market,
which accounts for a small volume of electricity sell and purchase at 15
minute windows in Continental Europe. Lastly we address a few details of the
library and contribute a few closing remarks.

##  II Theoretical Basis

In the following we briefly summarise the theoretical basis of Multifractal
Detrended Fluctuation Analysis. Later we detail the different included
extensions and which modifications these add to the original MFDFA algorithm.

###  II.1 Multifractal Detrended Fluctuation Analysis

Multifractal Detrended Fluctuation Analysis studies the variances of the
fluctuations of a given process by considering increasing segments of a time
series.

i) Take a time series $X(t)$ (in time or space $t$) with $N$ data points,
discretised as $X_{i}$, $i=1,2,\dots,N$. Find the “detrended” profile of the
process by defining

$Y_{i}=\sum_{k=1}^{i}\left(X_{k}-\mu_{X}\right),~{}\text{for}~{}i=1,2,\dots,N,$
(1)

i.e., the cumulative sum of $X_{i}$ subtracting the mean $\mu_{X}$ of the
data.

ii) Section the data into smaller non-overlapping segments of length $s$,
obtaining therefore $N_{s}=\text{int}(N/s)$ segments. Given the total length
of the data is not always a multiple of the segment’s length $s$, discard the
last points of the data.

iii) Consider the same data, apply the same procedure, but discard now instead
the first points of the data. One has now $2N_{s}$ segments of the time
series.

iv) To each of this segments fit a polynomial $y_{v}$ of order $m$ and
calculate the variance of the difference of the data to the polynomial fit

$F(v,s)=\frac{1}{s}\sum_{i=1}^{s}[Y_{(v-1)s+i}-y_{(v-1)s+i}]^{2},$ (2)

for $v=1,2,\dots,N_{s}$, where $y_{(v-1)s+i}$ is the polynomial fitting for
the segment $Y_{(v-1)s+i}$ of length $s$, fitted via least-squares. The order
of the polynomial $y_{v}$ can be freely chosen, giving rise to the denotes
(MF)DFA1, (MF)DFA2, $\dots$, (MF)DFA$m$, dependent on the chosen degree $m$ of
the polynomial.

v) Notice now $F(v,s)$ is a function of each variance of each $v$-segment of
data and of the different $s$-length segments chosen. Define the $q$-th order
fluctuation function by averaging over the $N_{s}$ variances of the segments
of size $s$

$F_{q}(s)=\left\\{\frac{1}{N_{s}}\sum_{v=1}^{N_{s}}[F(v,s)]^{q/2}\right\\}^{1/q}.$
(3)

The fluctuation function $F_{q}(s)$ depends on two parameters: the segment
size $s$ and the $q$-th power. The fluctuation function $F_{q}(s)$ is the
function we will focus on which the MFDFA algorithm developed extracts from
the data.

Two closely related algorithms are discussed and introduced here, DFA [1] and
MFDFA [3]. DFA is a particular case of MFDFA for the choice of $q=2$. What is
presented above is the MFDFA algorithm as according to Kantelhardt et al. [3],
for which a particular choice of $q=2$ leads to the fluctuation function
$F_{2}(s)$. The DFA fluctuation function $F_{2}(s)$ can unveil solely the
monofractal spectrum of a time series. If the examined time series $X_{i}$ is
monofractal, DFA is sufficient to describe and uncover the scaling relations
in the data. If not, one must rely on MFDFA and the study of the spectrum
unveiled by varying the $q$-th power.

We will later detail two changes: i) The first involving empirical mode
decomposition (EMD) for detrending, where the local polynomial fittings are
replaced and the trends of the data are subtracted by removing select
Intrinsic Mode Functions (IMFs) obtained via empirical mode decomposition. ii)
The second change involves substituting the non-overlapping segments with
overlapping ones.

Figure 1: Multifractal Detrended Fluctuation Analysis (MFDFA) of an exemplary
$288$ data points time series $X_{t}$. Panel a) shows, from top to bottom: the
time series; first-order polynomial fit (m = 1); third-order polynomial fit (m
= 3); EMD detrending with the slowest Intrinsic Mode Function (EMD); first-
order polynomial fit (m = 1) with a moving windows with a step of $36$ data
points (window = 36). Segments with a size of $s=72$ data points. The lines
indicate the fits, either via polynomials or EMD. Panel b) displays the
changing segment size $s$ for a first-order polynomial fit (m = 1); From top
to bottom: the time series; segmentation with $s=48$; $s=72$; $s=96$; $s=144$.
Shaded areas on both panels indicated the standard deviation of each segment.

The inherent scaling properties of the data, if the data displays power-law
correlations, can now be studied in a log-log plot of $F_{q}(s)$ versus $s$,
where the scaling of the data obeys a power-law with exponent $h(q)$ as

$F_{q}(s)\sim s^{h(q)}$ (4)

where $h(q)$ is the generalised Hurst exponent or self-similarity exponent,
which will dependent on $q$ if the data is multifractal, and relates directly
to the Hurst index [59]. The generalised Hurst exponent $h(q)$ is obtained by
finding the slope of $F_{q}(s)$ curve in the log-log plots.

If the data is monofractal, the generalised Hurst exponent $h(q)=H$ is
independent of $q$ and the generalised Hurst exponent is simply the Hurst
index $H$. On the other hand, if the data is multifractal, the dependence on
$q$ can be understood by studying the multifractal scaling exponent $\tau(q)$,
given by

$\tau(q)=qh(q)-1,$ (5)

which depends on the generalised Hurst exponent $h(q)$. Similarly, one can
construct the singularity spectrum $D(\alpha)$ as the Legendre transform [60,
61, 62, 63]. If $\tau(q)$ is sufficiently smooth, the singularity strength
$\alpha$ is given by

$\alpha=\tau^{\prime}(q)=h(q)+qh^{\prime}(q),$ (6)

from which the singularity spectrum $D(\alpha)$ can be constructed as

$D(\alpha)=q\alpha-\tau(q).$ (7)

The singularity spectrum $D(\alpha)$ describes the dimension of the subset of
the time series which is characterised by the singularity strength $\alpha$
[64]. The breadth of singularity strength $\alpha$ indicates the strength of
the multifractality of the time series, centred around the most prominent
scale of the time series, i.e., $h$. The singularity spectrum $D(\alpha)$
takes the shape of an inverted parabola with a maximum at $D(\alpha=0)=D_{0}$,
known as the box-counting or Minkowski–Bouligand dimension, or sometimes
simply fractal dimension [60]. $D(\alpha=1)=D_{1}$ is known as the information
dimension and $D(\alpha=2)=D_{2}$ the correlation dimension [65]. For a
clearer discussion of these properties, see Refs. [66, 3]. An extensive and
very illustrative representation of this can be found in Ref. [55]. For a
careful analysis of the meaning and interpretation of the generalised Hurst
coefficients extracted from (MF)DFA, see Ref. [67], where a description and
clarification is given on what are persistent and anti-persistent motions,
stationary and non-stationarity time series, among other relevant details.

####  II.1.1 Empirical mode decomposition

Empirical mode decomposition (EMD) is a method with a variety of applications
in time series analysis [68]. It seeks to extract the modes of oscillation of
a time series strictly from the data. One can harness the ability of the EMD,
i.e., the Hilbert–Huang decomposition of a time series, to obtain the trend or
trends of the time series and utilise those to transform non-stationary into
stationary data. The central concept, developed by Qian, Gu, and Zhou [6], is
to substitute the detrending method employed in the traditional MFDFA, i.e.,
polynomial fittings, by removing instead particular Intrinsic Mode functions
extracted via EMD. A sketch of the method can be seen in Fig. 1.

EMD can be summarised in a few steps: a set of intrinsic mode functions (IMFs)
are extracted from the time series, obeying: 1) the number of extrema and the
number of zero crossings must maximally differ by one. 2) for any point, the
mean value of the envelope defined by the local extrema is zero. Numerical
methods—as cubic splines—are used to find the curve that best fits “between”
the local extrema of the time series. The method is applied iteratively: i)
Obtain an IMF by finding the “best” curve between the local extrema of the
time series; ii) Subtract this IMF to the time series; iii) Repeat. Apply the
process recursively to the time series until the final IMF contains solely a
residual trend of the data.

####  II.1.2 Moving windows

The overlapping moving windows included in this library is not aimed at
detrending, but instead for the analyses of rather short time series or very
large scales in longer time series [8]. In the literature several applications
of moving average windows have been proposed as methods to remove trends and
ensure stationarity, by simply removing a windowed average to the time series
[8, 42, 9]. This is not what we do here. Here, we substitute the non-
overlapping segmentation, as explain after Eq. (1), by a moving window,
replacing the two separate segmentations by a moving window of each segment
size $s$ over the time series, as proposed by Zhou and Leung [8]. This is
particularly relevant when examining short time series, where quickly the
choice of larger lags $s$ separates the data into a small number of segments,
resulting in a poor statistics for the scaling at larger lags. A sketch of the
methods can be seen in Fig. 1.

####  II.1.3 Extended Detrended Fluctuation Analysis

A new metric of similar nature as the fluctuation function $F_{q}(s)$, given
in Eq. (3), has been proposed in Ref. [10] This measure supersedes the
$q$-order powers and takes in solely the case of DFA where $q=2$. Instead of
finding the average of the variances over each choice of segments of size $s$,
it considers the difference between the extrema of the fluctuation function at
each segment $s$. Take $F(v,s)$ as given in Eq. (2) and extract the maximum
and minimum of the variances over all windows $v$ for a certain window size
$s$

$\Delta F(s)=\max_{v}[F(v,s)]-\min_{v}[F(v,s)].$ (8)

This new metric $\Delta F(s)$ is denoted Extended Fluctuation Analysis. In
general, $\Delta F(s)$ can scale as a power law with a different exponent

$\Delta F(s)\sim s^{\beta}.$ (9)

This metric takes into account aperiodicities in the data which, in some
sense, would be accounted for as a multifractal behaviour. It can unravel a
second scaling phenomenon due to local changes of a time series’ period.

##  III Examples

To exemplify the usage of MFDFA, we first take two common examples of
stochastic processes, a fractional Ornstein–Uhlenbeck process and general
process that has a symmetric Lévy $\alpha$-stable distribution, with single
parameter $\alpha$. We will show how to extract the fluctuation function
$F_{q}(s)$ and how to interpret the plots conventionally extracted to perform
the analysis. Subsequently we test the algorithm with real-world data on
sunspot time series, following Ref. [26], and later apply the algorithm to
electricity price time series from the European Power Exchange.

###  III.1 Numerically generated data

####  III.1.1 Fractional Ornstein–Uhlenbeck process

Figure 2: Multifractal Detrended Fluctuation Analysis (MFDFA) of three
exemplary sample paths of fractional Ornstein–Uhlenbeck processes, given by Eq
(10), with Hurst indices of $H=0.3,0.5$, and $0.7$. Panel a) displays the log-
log plot of the segment size $s$ versus the fluctuation function $F_{2}(s)$,
given by Eq. (3), for $q=2$. Each line has a slope of $H+1$, as expected. The
inset shows $F_{q}(s)$ for the case of $H=0.3$ and the power variations
$q=-10,-2,2,10$. These lines are all parallel indicating that the process is
monofractal, as expected. The dashed lines indicate the theoretical expected
scaling, i.e., a slope of $H+1=0.3+1$, where the $+1$ account for the increase
in regularity due to the integration. The generalised Hurst coefficients
$h(q)$, which are simply $H+1$, are obtained by extracting the slopes of the
curves (in a log-log scale). Panel b) shows the multifractal scaling exponent
$\tau(q)$, given by Eq. (4), which exhibits a linear dependency, i.e.,
$h(q)=H$, indicating again the process is monofractal. The processes were
numerically integrated with an integration step $\Delta t=0.001$ over
$N=10^{4}$ time units ($N=10^{7}$ data points). The MFDFA algorithm ran in $1$
min $29$ s $\pm$ $1.85$ s, for $100$ segments $s$ and $40$ $q$-variation
powers, with first-order polynomial fits.

To study the scaling effects in continuous stochastic processes, three
exemplary fractional Ornstein–Uhlenbeck processes are taken, defined as [69]

$dX_{t}=-\theta X_{t}dt+\sigma dB^{H}_{t},$ (10)

with a fractional Brownian motion $B^{H}_{t}$ with the covariance function

$\mathbb{E}\left[B^{H}_{t}B^{H}_{t^{\prime}}\right]=\frac{1}{2}\left(|t|^{2H}+|t^{\prime}|^{2H}-|t-t^{\prime}|^{2H}\right).$
(11)

Eq. (10) fixed mean reverting strength $\theta=1.0$ and volatility
$\sigma=0.5$, with three distinct Hurst indices of $H=0.3$, $0.5$, and $0.7$.
Note here that the classic uncorrelated Brownian motion is given by $H=0.5$. A
fractional Brownian motion has a self-similarity exponent given by the Hurst
index $H$, thus the three choices of fractional Ornstein–Uhlenbeck should
result in a scaling of $h(q)=H+1$. The $+1$ is due to the integration, which
smooths the fluctuations and thus increases the regularity of the process. We
will now numerically integrate these processes and utilise the MFDFA library
to identify the Hurst coefficients and the presence of a monofractal vs
multifractal spectrum in the time series.

Let us exemplify how to numerically generate data and utilise the MFDFA
library Load the MFDFA library alongside with the fractional Brownian noise
generator fgn included in your python console or editor.

Listing 1: Load the MFDFA library

⬇

1from MFDFA import MFDFA

2from MFDFA import fgn

To numerically integrate an Ornstein–Uhlenbeck process, given by Eq. (10), we
utilise an Euler–Maruyama scheme with a stepsize $\Delta t=0.001$ for a total
time of $t=10^{4}$ (thus we have $10^{7}$ data points). Here exemplified is
the fractional Ornstein–Uhlenbeck process with $H=0.3$.

Listing 2: Integrate Ornstein–Uhlenbeck process

⬇

3# integration time and time sampling

4t_final = 10000

5delta_t = 0.001

6N = int(t_final/delta_t)

7

8# The parameters theta and sigma

9theta = 1

10sigma = 0.5

11

12# Initialise the array X

13X = np.zeros(N)

14

15# Generate the fractional Brownian noise

16# with a Hurst coefficient of H = 0.3

17dB = (t_final ** H) * fgn(N, H = 0.3)

18

19# Integrate the process

20for i in range(1,N):

21 X[i] = X[i-1] - theta*X[i-1]*delta_t + sigma*dB[i]

To retrieve the MFDFA spectrum of the generated time series, define the set of
$q$ power variations and the lags $s$ to examine, and call the MFDFA function.

Listing 3: Applying MFDFA

⬇

22# 100 lag s points from 3 to 1000

23lag = np.logspace(0.6,3,118).astype(int)

24lag = np.unique(lag)

25

26# q power variations, removing 0 power

27q = np.linspace(-10,10,41)

28q = q[q!=0.0]

29

30lag, fluct = MFDFA(X, lag = lag, q = q)

When not declaring the values of the $q$ powers, $q=2$ is assumed, thus
resulting in the conventional DFA. Likewise, not declaring the order of the
polynomial fitting, a first-order polynomial is assumed, i.e., order = 1.

In Fig. 2 the MFDFA of the three processes can be seen. In panel a) the
fluctuation function $F_{2}(s)$, with $q=2$, is shown for a polynomial
detrending of first order. This is the conventional DFA. The slopes of each
curve in the log-log plot reveal the Hurst indices of each process, i.e., the
fractional Ornstein–Uhlenbeck with Hurst $H=0.3$ scales with a slope of
$1.3=0.3+1$, the other two with $H=0.5$ and $H=0.7$ have a slope of $1.5$ and
$1.7$, respectively. The inset in panel a) shows the fluctuation function
$F_{q}(s)$, with $q=-10$, $-2$, $2$, and $10$, for the fractional
Ornstein–Uhlenbeck process with $H=0.3$. Note that the slope of all power
variations is the same, i.e., the process is monofractal, as expected. The
monofractality of the process is also evident in panel b). The multifractal
scaling exponent $\tau(q)$ is shown and is purely a linear function. Likewise,
in the inset, the generalised Hurst indices $h(q)$ for the three processes for
a set of power variations $q\in[-10,10]$ is displayed. The linear shape of
$\tau(q)$ and constant value of $h(q)$ in the inset indicates, as expected,
that these three processes are monofractal. Small deviations are seen for very
negative $q$ powers ($q\lesssim 7$), which arise due to the numeric (negative)
powering operation, which highly depends on the numerical precision of the
data.

Figure 3: Multifractal Detrended Fluctuation Analysis (MFDFA) of three
exemplary symmetric Lévy $\alpha$-stable distributed processes, with
$\alpha=1.75$, $1.25$, and $0.75$. In panel a) the fluctuation function
$F_{q}(s)$ is shown as a function of the segment size $s$ on double
logarithmic scales for $\alpha=1.25$ and different values of the power,
$q=-10,-5,-2,2,5,10$. For $q>\alpha$ the curves are not parallel, indicating
the multifractal nature of the process. Panel b) displays the generalised
Hurst exponent $h(q)$, where a clear non-linear dependency on $q$ is
observable. The inset displays the multifractal scaling exponent $\tau(q)$
displaying two clear distinct behaviours for $q<0$ and $q>\alpha$. The solid
lines indicate the theoretical expected scaling for $q<0$. The three processes
were drawn from Lévy $\alpha$-stable distributions, each with $N=10^{7}$ data
points. The MFDFA algorithm ran in $1$ min $24$ s $\pm$ $2.17$ s for $100$
segments $s$, $40$ $q$ powers, and third-order polynomial fits.

####  III.1.2 Lévy-driven process

As a second example, take a collection of Lévy distributed random variables
[70]. That is, each $X_{t}$ is drawn independently from a symmetric
$\alpha$-stable distribution, such that the probability density function of
$X(t)$ vanishes as a power-law $P(x)\sim|x|^{-(\alpha-1)}$ for large $|x|$
[70]. These processes exhibit heavy tails, ill-defined variances, and
multifractal scaling. In Fig. 3 three symmetric Lévy $\alpha$-stable
distributed processes with $\alpha=1.75,1.25$, and $0.75$ are studied with
MFDFA.

The multifractal behaviour can be identified directly in panel a), where the
fluctuation function $F_{q}(s)$ for $\alpha=1.25$ is shown. The lines of
$F_{q}(s)$ are not parallel for different positive $q$ powers, showing that
the process is not mono-fractal. In fact, the process is bi-fractal, having a
separate behaviour for $q<0$ and $q>\alpha$. For positive power variations
$q>\alpha$, the generalised Hurst exponent $h(q)$ decays like $1/q$. For
values of $q<0$, the generalised Hurst exponent $h(q)=1/\alpha$. This can be
seen clearly in Fig. 3 b), where the generalised Hurst exponent $h(q)$ is
displayed. In the inset, one notices that the multifractal scaling exponent
$\tau(q)=0$, for $q>\alpha$ (always zero for $q>2$), once again showing that
none of these processes are distinguishable for positive power variations.

In general, without the aid of the multifractal spectra, which we uncovered by
studying MFDFA for a range of $q$ values, it is not possible to distinguish
between Lévy distributed processes. The particular choice of $q=2$, i.e.,
conventional DFA, obscures the fractality of these processes, as they all show
a similar scaling for $q=2$, i.e., $h(q=2)=1/2$ for all Lévy motions,
including (non-fractional) Brownian motions (where $\alpha=2$).

###  III.2 Real-world data, empirical mode decomposition, and extended DFA

In order to evaluate the efficiency of the algorithm, we test here two real-
world data sets. Firstly, using MFDFA we will evaluate the multifractality of
sunspots time series, a recurring phenomenon on the Sun’s photosphere which
can be observed with a telescope [26]. Secondly, we will analyse the German
and Autrian spot market intraday quarter-hourly electricity price extracted
from the European Power Electricity Exchange (EPEX SPOT) [71, 72]. We
illustrate the application of two advanced features of the developed python
package, moving windows and the masking of missing data points.

####  III.2.1 Sunspots

Figure 4: Multifractal Detrended Fluctuation Analysis (MFDFA) of sunspot time
series from 1818 to 2020, by the ILSO World Data Center, Royal Observatory of
Belgium, Brussels [73]. Panel a) shows the number of sunspots registered from
1818 to 2020. Panel b) displays the fluctuation function $F_{q}(s)$ as a
function of the segment size $s$ on a double-logarithmic scale for
$q=-10,-5,-2,2,5,10$, with positive $q$ values in orange and negative $q$
values in green. Panel c) displays the generalised Hurst coefficient $h(q)$
over $q$, and the inset displays the multifractal scaling exponent $\tau(q)$,
given by Eq. (4), both highlighting the multifractal spectrum of the data
($h(q)$ is not constant over $q$, $\tau(q)$ is not linear over $q$). Panel d)
displays the singularity spectrum $D(\alpha)$ over the singularity strength
$\alpha$ which shows a large breadth of $\alpha$ spanning over $[1.25,2.25]$,
indicating the strong multifractality of the data. The MFDFA algorithm ran in
$426$ ms $\pm$ $11.8$ ms, for $70$ segments $s$, $40$ $q$ powers, and first-
order polynomial fits. The missing values were neatly removed by utilising
numpy’s masked arrays, which is integrated in MFDFA and allows the user to
simply “mask” empty or corrupted data.

The sunspot numbers, also called Wolf numbers, are a rather simple measure of
solar activity by counting in a weighted manner the number of groups of
sunspots and single sunspots visible from the Earth in the solar photosphere,
i.e. it is an integrated measure over space [74, 75]. Hence, the sunspot
numbers form a time series which has a mean period of about 11 years, but is
far from being simply periodic [62]. Solar activity is the result of complex
magneto-hydrodynamic processes in the Sun characterised by a highly complex
spatio-temporal dynamics. It is of special interest to analyse the rather long
series of sunspot numbers in order to explore some relations to the underlying
spatio-temporal system.

The emergence of sunspots has a distinct statistics and a multifractal
spectrum which has been examined in Ref. [26]. This publication has become a
reference for multifractal studies as the data from the ILSO World Data
Center, Royal Observatory of Belgium, Brussels is freely available [73]. Here,
we will focus on numerical efficiency and how to deal with missing or corrupt
data. We utilised another feature integrated in MFDFA that enables an
efficient management of missing data points. In python’s numpy arrays, missing
or corrupt values in a time series can be handled with masked data, which logs
the missing data points and takes these into account while performing
averages, sums, and power operations. When calculating averages or the
variance of a segment, or when taking powers, the masked entries are not taken
into account. For the particular application with sunspot time series, which
are recorded daily since 1818, there are $3247$ missing values, over a total
of $74145$ entries, i.e., roughly $4.4\%$ of the data is missing. To go around
this, simply use

Listing 4: MFDFA and missing data

⬇

1# Read data whichever way preferred

2data = read_data(’sunspot.csv’)

3

4# Mask missing values. For this case -1

5# is a missing entry in the record

6data[data ==-1.] = np.nan

7data = np.ma.masked_invalid(data)

8

9# Run MFDFA (choose lag and q)

10lag, fluct = MFDFA(data, lag, q)

The MFDFA will extract the variances as it is possible, taking into account
the missing values in the time series. Here we highlight that MFDFA calculated
$40$ $q$-powers over $70$ segments $s$ in $426$ ms $\pm$ $11.8$ ms.

In Fig. 4 we display the fluctuation function $F_{q}(s)$ for
$q=-10,-5,-2,2,5,10$, in panel a), for $s\in[70,1000]$ These curves are not
parallel, suggesting the time series is not monofractal. In panel b) the
generalised Hurst exponent $h(q)$ is shown as function of $q$, and similarly
the multifractal scaling exponent $\tau(q)$ in inset. The generalised Hurst
exponent $h(q)$ is not constant over $q$ and consequently the multifractal
scaling exponent $\tau(q)$ is not linear, indicating clearly the time series
in multifractal. In panel c) we display the singularity spectrum $D(\alpha)$
over the singularity strength $\alpha$, as given by Eq. (7). The singularity
strength $\alpha$ spans a wide range of values, over $[1.25,2.25]$, indicating
how strongly multifractal the time series is. Here recall that a monofractal
time series, as the fractional Ornstein–Uhlenbeck previously shown in Fig. 2,
has a very narrow range of the singularity strength $\alpha$, centred around
$H$. For the case of sunspot time series we see a wide range of $\alpha$,
indicating the various scales of the phenomenon. Note as well that $h(q)$ and
$\alpha$ are always larger than $1$, indicating that this is a non-stationary
process.

####  III.2.2 German and Austrian spot market intraday quarter-hourly
electricity price time series

Figure 5: Multifractal Detrended Fluctuation Analysis (MFDFA) of the spot
market intraday quarter-hourly German and Austrian electricity price time
series from 2015 to 2019, traded at the European Power Exchange (EPEX SPOT)
[71, 72]. Panel a) displays the price in EUR/kWh from 2015 to 2019. Panel b)
displays the log-log plot of the segment size $s$ versus the fluctuation
function $F_{q}(s)$ for $q=-10,-5,-2,2,5,10$. Orange and green markers
indicate the segments larger than two days, where purple indicate segments
between 1 and 48 hours. There two scales are studied separately. Panel c)
displays the generalised Hurst coefficient $h(q)$ over $q$, and the inset
displays the multifractal scaling exponent $\tau(q)$, given by Eq. (4). Panel
d) displays the singularity spectrum $D(\alpha)$ over the singularity strength
$\alpha$. The short-time scale (1–48 hours) displays large generalised Hurst
coefficient $h(q)$ and a very large breadth of the singularity strength
$\alpha$, indicating precisely the high volatility of the market at short time
scales. In comparison, the longer time scales ($>48$ hours) are much “milder”,
and the variations of $\alpha\in[0.67,0.92]$, which indicates the process is
both stationary, int the long run, and only moderately volatile. The MFDFA
algorithm ran in $2$ min $11$ s $\pm$ $4.43$ s for $50$ segments $s$, $40$ $q$
powers, first-order polynomial fits, and the moving window.

We will examine now a 4-year long time series sampled at 15 minutes of the
spot market intraday quarter-hourly German and Austrian electricity price [71,
72], from the 1st of January 2015 to end of December 2019, traded at the
European Power Exchange (EPEX SPOT). To the extent of our knowledge no
multifractal analysis of this particular data has been performed before, but
other multifractal studies of price time series exist [44]. In Wang et al.
[44], the authors examine different scaling properties for selected periods of
low, regular, and high electricity price for some United States of America’s
electricity markets in the year 2000 and 2001. Here we propose a different
analysis, studying the data and examining a short and long time scale of the
data without separating different activity periods.

We know that the 15 minute trading electricity market amounts to a small
volume of the overall exchange electricity sold, thus this market serves only
electricity producers which can either extract or inject power from the power-
grid system in a very fast manner ($<15$ minutes) [76, 77]. This will lead us
to explore to separate scaling phenomena in the data: A short and a long
timescale of market activities. The expectation is that the very short-time
trading is highly volatile, given the necessity of the power grid in injecting
or extracting power is a fairly speedy manner. In the long run, the quarter-
hourly market is intrinsically linked to the larger hourly and daily
electricity market, which has far less variability, as most of the power is
sold in lengthier contracts, stabilising the value of the electricity price.
Thus one expects a narrower multifractality at large temporal scales.
Multifractality is nevertheless expected, as the system exhibits very large
yet seldom negative prices, as well as an occasional four of five-fold
increase of the (positive) prices, again occurring seldom and lasting very
short periods.

In order to obtain a better statistics of the shorter time scales, we will
employ MFDFA’s moving windows method previously discussed. The moving window
method requires the input of the number of steps used to “move” the windows.
For the following example, the window parameter is set to 1, thus each
overlapping window is displaced by solely $1$ data point. This substantially
increases the computational time as each averaging operation is repeated by
the number of segments $s-1$. For Fig. 5, the total calculation lasted $2$ min
$11$ s $\pm$ $4.43$ s for the windowed mode, compared with $764$ ms $\pm$
$9.32$ ms with the conventional non-overlapping windows, for $50$ $s$
segments.

Listing 5: MFDFA and moving window

⬇

1# Read data whichever way preferred

2data = read_data(’price.csv’)

3

4# Run MFDFA (choose lag and q)

5lag, fluct = MFDFA(data, lag, q,

6 extension = {’window’ : 1})

In Fig. 5 we display the MFDFA analysis of the price time series, as
previously done for the sunpots in Fig. 4. We perform a similar analysis as
above, thus we will condense the technical details and focus on the
interpretation. Previous studies point to a clear separation of the scaling
behaviour of price time series [78, 79, 43]. They separate two time scales for
periods shorter and longer than 24 hours. These studies used pricing data from
before 2004 for the Nordic grid (Nordpool). In our analysis, we similarly
separate two scales, between 1 and 48 hours and between 48 hours and 10 days.
These are indicated in purple (short timescale) and orange and green (long
time scale). We first observe that negative $q$ powers do not exist for the
short time scale. This is not unusual, many processes do not show a
multifractal spectrum with negative $q$ values. Note that this involves taking
negative powers of the average of the variances, which is not always well
defined for short $s$ segments. This also served as a threshold to assess the
change in the fractal behaviour of the time series. For the large time scales
($>48$ hours) the negative powers are well defined, and we can identify the
full singularity spectrum $D(\alpha)$, as seen in panel c). We note that the
short time scale ($<48$ hours) has a very strong multifractality (in purple).
The singularity strength $\alpha$, which we can only extract for positive $q$
values, has a very large breadth, especially with its equivalent for the large
time scale (in orange). This is well grounded on the previous arguments of
having a very volatile market at these short time scales, thus these results
are in line with what is known about this market: The high volatility and
occasional burst—into very large electricity prices or into negative
prices—generate a wide range of the singularity strength. The long-term
stability, connected with the larger intraday and day-ahead electricity
markets, makes the process far less multifractal at large temporal scales.

##  IV The MFDFA library

The Multifractal Detrended Fluctuation Analysis library MFDFA in python
presented is a standalone package based integrally on python’s numpy [58]. It
can be found in https://github.com/LRydin/MFDFA. It harnesses numpy’s
vectorised polynomial fittings, making it possible to utilise all
computational cores in a computer’s processor(s). Additionally, EMD is
included as an extra feature, which is integrated into MFDFA by simply
installing the python library PyEMD [80]. The conventional plots associated
with multifractal analysis, i.e., $F_{q}(s)$ vs. $s$, $h(q)$ and $\tau(q)$ vs.
$q$, and $D(\alpha)$ vs. $\alpha$, are available as well and require the
plotting library matplotlib [81].

The MFDFA library accepts numpy’s masked arrays, which is particularly
convenient when dealing with time series with missing data, as exemplified in
Sec. III.2 and Fig. 4.

The MFDFA library offers a considerable speed-up in comparison with the
available MATLAB version. The library is fully developed to work with multi-
threading, which shows an increase in the performance, while handling time
series larger than $10^{5}$ data points. In Fig. 6 we display the performance
of the MFDFA library for time series of fractional Ornstein–Uhlenbeck
processes given in Eq. (10) of increasing length. The MFDFA operation scales
linearly with the number of points of the generated time series. The MFDFA
algorithm runs in under $1$ second for time series having up to $10^{5}$
datapoints, with a first-order polynomial fittings, $100$ segments $s$ and
$40$ $q$ powers, and outperforms the conventional library in MATLAB by up to a
factor of $\times 10^{3}$ in computational speed.

Estimation error and significance calculation have not been included in the
library, as the focus lied on computational speed and the inclusion of several
extra features, as discussed.

Figure 6: Speed performance of MFDFA for time series with sizes varying
between $[10^{3},10^{6}]$ data points of a fractional Ornstein–Uhlenbeck as
given by Eq. (10). Included are first-order and third-order polynomial fits,
first-order fits with extended DFA, and first-order fits with a moving window
with a step size of $5$. A comparison with the distributed MATLAB code is
included [55]. Tests ran on python 3.8.2 and MATLAB R2020b. MFDFA has a
average speed-up compared with the MATLAB code, with a five-fold speed
increase for first-order polynomial fits (m=1) and a $\times 27$-fold increase
for third-order polynomials fits (m=3). Both codes were tested for $100$
segments $s$ and $40$ $q$ powers. All tasks were performed on a laptop on two
computer cores at 2.9 GHz each.

##  V Conclusion

We have presented a numerically efficient python implementation of
Multifractal Detrended Fluctuation Analysis called MFDFA. MFDFA has found
extensive application in the past two decades, yet a reliable, all-
encompassing open-source software in python does not exist to this date. In
this library we have harnessed the most of python’s flexibility with handling
matricial operations and multi-threaded polynomial fittings. In this
implementation we have included some of the more common extensions of MFDFA,
including a simple empirical mode decomposition as a mechanism to detrend the
data, a moving window to handle very short time series, and the extended
Detrended Fluctuation Analysis, which can track a different scaling mechanism
for non-stationary time series. The MFDFA library can also handle missing
values in the data with the aid of numpy’s masked time series.

We have initially turned to two classic numerically generated stochastic
processes, fractional Ornstein–Uhlenbeck processes and Lévy-distributed
motions, and uncovered their monofractal and multifractal with MFDFA.
Subsequently we have studied two real-world time series, the sunspot count
from 1818 to 2020 and the quarter-hourly electricity price time series from
2015 to 2019. For both we performed a multifractal analysis, unveiling their
scaling properties and the strength of their multifractality. We focused on
MFDFA’s speed, the ability to handled missing data, and the integrated
overlapping moving window. The analysis displayed here covered only part of
MFDFA’s integrated options, thus we leave the user to explore the other
implemented methods, as the extended DFA and EMD detrending, as these are more
specialised to particular research fields.

We hope with this contribution we open a door to fast MFDFA calculations that
can the performed on a local machine without an extensive numerical effort and
very long time runs, thus permitting in the future to analyse larger time
series.

##  VI Acknowledgements

L.R.G. kindly thanks Francisco Meirinhos for all the help with python, and
Fabian Harang, Marc Lagunas Merino, Anton Yurchenko-Tytarenko, Dennis
Schroeder, Michele Giordano, Giulia di Nunno, and Fred Espen Benth for their
support. L.R.G and D.W gratefully acknowledge support by the Helmholtz
Association, via the joint initiative Energy System 2050 - A Contribution of
the Research Field Energy and the grant Uncertainty Quantification – From Data
to Reliable Knowledge (UQ), with grant no. ZT-I-0029, the scholarship funding
from E.ON Stipendienfonds, and the STORM - Stochastics for Time-Space Risk
Models project of the Research Council of Norway (RCN) no. 274410. This work
was performed as part of the Helmholtz School for Data Science in Life, Earth
and Energy (HDS-LEE). J.K. was financed by the Ministry of Science and Higher
Education of the Russian Federation within the framework of state support for
the creation and development of World-Class Research Center “Digital biodesign
and personalised healthcare”, no. 075-15-2020-926.

## References

  * Peng  _et al._ [1994] C.-K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E. Stanley, and A. L. Goldberger, Physical Review E 49, 1685 (1994).
  * Peng  _et al._ [1995] C. Peng, S. Havlin, H. E. Stanley, and A. L. Goldberger, Chaos: An Interdisciplinary Journal of Nonlinear Science 5, 82 (1995).
  * Kantelhardt  _et al._ [2002] J. W. Kantelhardt, S. A. Zschiegner, E. Koscielny-Bunde, S. Havlin, A. Bunde, and H. Stanley, Physica A: Statistical Mechanics and its Applications 316, 87 (2002).
  * Wu and Huang [2004] Z. Wu and N. E. Huang, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 460, 1597 (2004).
  * Wu and Huang [2009] Z. Wu and N. E. Huang, Advances in Adaptive Data Analysis 01, 1 (2009).
  * Qian  _et al._ [2011] X.-Y. Qian, G.-F. Gu, and W.-X. Zhou, Physica A: Statistical Mechanics and its Applications 390, 4388 (2011).
  * Zhang  _et al._ [2019] X. Zhang, G. Zhang, L. Qiu, B. Zhang, Y. Sun, Z. Gui, and Q. Zhang, Water 11, 891 (2019).
  * Zhou and Leung [2010] Y. Zhou and Y. Leung, Journal of Statistical Mechanics: Theory and Experiment 2010, P06021 (2010).
  * Lai  _et al._ [2019] S. Lai, L. Wan, and X. Zeng, Journal of Physics: Conference Series 1345, 042086 (2019).
  * Pavlov  _et al._ [2020a] A. N. Pavlov, A. S. Abdurashitov, A. A. Koronovskii, O. N. Pavlova, O. V. Semyachkina-Glushkovskaya, and J. Kurths, Communications in Nonlinear Science and Numerical Simulation 85, 105232 (2020a).
  * Pavlov  _et al._ [2020b] A. N. Pavlov, A. I. Dubrovsky, A. A. Koronovskii, O. N. Pavlova, O. V. Semyachkina-Glushkovskaya, and J. Kurths, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 073138 (2020b).
  * Pavlov  _et al._ [2020c] A. Pavlov, A. Dubrovsky, A. Koronovskii Jr, O. Pavlova, O. Semyachkina-Glushkovskaya, and J. Kurths, Chaos, Solitons & Fractals 139, 109989 (2020c).
  * Pavlov  _et al._ [2021] A. N. Pavlov, O. N. Pavlova, O. V. Semyachkina-Glushkovskaya, and J. Kurths, The European Physical Journal Plus 136, 10 (2021).
  * Podobnik and Stanley [2008] B. Podobnik and H. E. Stanley, Physical Review Letters 100, 084102 (2008).
  * Zhou [2008] W.-X. Zhou, Physical Review E 77, 066211 (2008).
  * Alvarez-Ramirez  _et al._ [2009] J. Alvarez-Ramirez, E. Rodriguez, and J. C. Echeverria, Physical Review E 79, 057202 (2009).
  * Chianca  _et al._ [2005] C. Chianca, A. Ticona, and T. Penna, Physica A: Statistical Mechanics and its Applications 357, 447 (2005).
  * Hu  _et al._ [2001] K. Hu, P. C. Ivanov, Z. Chen, P. Carpena, and H. Eugene Stanley, Physical Review E 64, 011114 (2001).
  * Ivanov  _et al._ [1999] P. C. Ivanov, L. A. N. Amaral, A. L. Goldberger, S. Havlin, M. G. Rosenblum, Z. R. Struzik, and H. E. Stanley, Nature 399, 461 (1999).
  * Dutta  _et al._ [2013] S. Dutta, D. Ghosh, and S. Chatterjee, Frontiers in Physiology 4, 274 (2013).
  * Madanchi  _et al._ [2020] A. Madanchi, F. Taghavi-Shahri, S. M. Taghavi-Shahri, and M. R. Rahimi Tabar, The European Physical Journal B 93, 126 (2020).
  * Movahed  _et al._ [2011] M. S. Movahed, F. Ghasemi, S. Rahvar, and M. R. R. Tabar, Physical Review E 84, 021103 (2011).
  * Movahed  _et al._ [2013] M. S. Movahed, B. Javanmardi, and R. K. Sheth, Monthly Notices of the Royal Astronomical Society 434, 3597 (2013).
  * Telesca  _et al._ [2005] L. Telesca, V. Lapenna, and M. Macchiato, Physica A: Statistical Mechanics and its Applications 354, 629 (2005).
  * Shadkhoo  _et al._ [2009] S. Shadkhoo, F. Ghanbarnejad, G. Jafari, and M. R. R. Tabar, Open Physics 7, 620 (2009).
  * Sadegh Movahed  _et al._ [2006] M. Sadegh Movahed, G. R. Jafari, F. Ghasemi, S. Rahvar, and M. R. R. Tabar, Journal of Statistical Mechanics: Theory and Experiment 2006, P02003 (2006).
  * Tanna and Pathak [2014] H. J. Tanna and K. N. Pathak, Astrophysics and Space Science 350, 47 (2014).
  * Meyer and Kantz [2019] P. G. Meyer and H. Kantz, Climate Dynamics 53, 6909 (2019).
  * Pedron [2010] I. T. Pedron, Journal of Physics: Conference Series 246, 012034 (2010).
  * Tessier  _et al._ [1996] Y. Tessier, S. Lovejoy, P. Hubert, D. Schertzer, and S. Pecknold, Journal of Geophysical Research: Atmospheres 101, 26427 (1996).
  * Matsoukas  _et al._ [2000] C. Matsoukas, S. Islam, and I. Rodriguez-Iturbe, Journal of Geophysical Research: Atmospheres 105, 29165 (2000).
  * Koutsoyiannis [2003] D. Koutsoyiannis, Hydrological Sciences Journal 48, 3 (2003).
  * Kantelhardt  _et al._ [2006] J. W. Kantelhardt, E. Koscielny-Bunde, D. Rybski, P. Braun, A. Bunde, and S. Havlin, Journal of Geophysical Research: Atmospheres 111, D01106 (2006).
  * Zhang  _et al._ [2008] Q. Zhang, C.-Y. Xu, Y. D. Chen, and Z. Yu, Hydrological Processes 22, 4997 (2008).
  * Rodríguez  _et al._ [2013] R. Rodríguez, M. C. Casas, and A. Redaño, Meteorology and Atmospheric Physics 121, 181 (2013).
  * Wu  _et al._ [2018] Y. Wu, Y. He, M. Wu, C. Lu, S. Gao, and Y. Xu, Scientific Reports 8, 16553 (2018).
  * Figueirêdo  _et al._ [2010] P. Figueirêdo, E. Nogueira, M. Moret, and S. Coutinho, Physica A: Statistical Mechanics and its Applications 389, 2090 (2010).
  * Zunino  _et al._ [2008] L. Zunino, B. Tabak, A. Figliola, D. Pérez, M. Garavaglia, and O. Rosso, Physica A: Statistical Mechanics and its Applications 387, 6558 (2008).
  * Zunino  _et al._ [2009] L. Zunino, A. Figliola, B. M. Tabak, D. G. Pérez, M. Garavaglia, and O. A. Rosso, Chaos, Solitons & Fractals 41, 2331 (2009).
  * Grech and Pamuła [2013] D. Grech and G. Pamuła, Acta Physica Polonica A 123, 529 (2013).
  * Drożdż  _et al._ [2018] S. Drożdż, R. Kowalski, P. Oświȩcimka, R. Rak, and R. Gȩbarowski, Complexity 2018, 7015721 (2018).
  * Lee  _et al._ [2018] M. Lee, J. W. Song, S. Kim, and W. Chang, Physica A: Statistical Mechanics and its Applications 512, 1278 (2018).
  * Weron  _et al._ [2004a] R. Weron, I. Simonsen, and P. Wilman, in  _The Application of Econophysics_, edited by H. Takayasu (Springer Japan, Tokyo, 2004) pp. 182–191.
  * Wang  _et al._ [2013] F. Wang, G. Liao, J. Li, X. Li, and T. Zhou, Physica A: Statistical Mechanics and its Applications 392, 5723 (2013).
  * Shalalfeh  _et al._ [2016] L. Shalalfeh, P. Bogdan, and E. Jonckheere, in  _2016 IEEE International Conference on Smart Grid Communications (SmartGridComm)_ (2016) pp. 466–471.
  * Rydin Gorjão  _et al._ [2020] L. Rydin Gorjão, R. Jumar, H. Maass, V. Hagenmeyer, G. C. Yalcin, J. Kruse, M. Timme, C. Beck, D. Witthaut, and B. Schäfer, Nature Communications 11, 6362 (2020).
  * Leung  _et al._ [2011] Y. Leung, E. Ge, and Z. Yu, Annals of the Association of American Geographers 101, 1221 (2011).
  * Jafari  _et al._ [2007] G. R. Jafari, P. Pedram, and L. Hedayatifar, Journal of Statistical Mechanics: Theory and Experiment 2007, P04012 (2007).
  * Telesca and Lovallo [2011] L. Telesca and M. Lovallo, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 467, 3022 (2011).
  * Ribeiro  _et al._ [2012] H. V. Ribeiro, L. Zunino, R. S. Mendes, and E. K. Lenzi, Physica A: Statistical Mechanics and its Applications 391, 2421 (2012).
  * Alados and Huffman [2000] C. L. Alados and M. A. Huffman, Ethology 106, 105 (2000).
  * Rutherford  _et al._ [2003] K. M. Rutherford, M. J. Haskell, C. Glasbey, R. Jones, and A. B. Lawrence, Applied Animal Behaviour Science 83, 125 (2003).
  * Li  _et al._ [2019] J. Li, X. Ma, M. Zhao, and X. Cheng, Electronics 8, 209 (2019).
  * Madanchi  _et al._ [2021] A. Madanchi, J. W. Yu, M. R. R. Tabar, W. B. Lee, and S. E. E. Rahbari, Soft Matter , Accepted Publication (2021).
  * Ihlen [2012] E. Ihlen, Frontiers in Physiology 3, 141 (2012).
  * Laib  _et al._ [2018a] M. Laib, L. Telesca, and M. Kanevski, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 033108 (2018a).
  * Laib  _et al._ [2018b] M. Laib, J. Golay, L. Telesca, and M. Kanevski, Chaos, Solitons & Fractals 109, 118 (2018b).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del R’ıo, M. Wiebe, P. Peterson, P. G’erard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature 585, 357 (2020).
  * Hurst [1951] H. E. Hurst, Transactions of the American Society of Civil Engineers 116, 770 (1951).
  * Hentschel and Procaccia [1983] H. Hentschel and I. Procaccia, Physica D: Nonlinear Phenomena 8, 435 (1983).
  * Halsey  _et al._ [1986] T. C. Halsey, M. H. Jensen, L. P. Kadanoff, I. Procaccia, and B. I. Shraiman, Physical Review A 33, 1141 (1986).
  * Kurths and Herzel [1987] J. Kurths and H. Herzel, Physica D: Nonlinear Phenomena 25, 165 (1987).
  * Meneveau and Sreenivasan [1987] C. Meneveau and K. Sreenivasan, Nuclear Physics B - Proceedings Supplements 2, 49 (1987).
  * Salat  _et al._ [2017] H. Salat, R. Murcio, and E. Arcaute, Physica A: Statistical Mechanics and its Applications 473, 467 (2017).
  * Falconer [2014] K. Falconer,  _Fractal geometry: mathematical foundations and applications_ , 3rd ed. (John Wiley & Sons, Chichester, 2014).
  * Barabási and Vicsek [1991] A.-L. Barabási and T. Vicsek, Physical Review A 44, 2730 (1991).
  * Serinaldi [2010] F. Serinaldi, Physica A: Statistical Mechanics and its Applications 389, 2770 (2010).
  * Huang  _et al._ [1998] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 903 (1998).
  * Tabar [2019] M. R. R. Tabar,  _Analysis and Data-Based Reconstruction of Complex Nonlinear Dynamical Systems_ (Springer International Publishing, 2019).
  * Applebaum [2011] D. Applebaum,  _Lévy Processes and Stochastic Calculus_, 2nd ed. (Cambridge University Press, 2011).
  * EPE [2021] European Power Exchange (EPEX SPOT) (2021), Market data https://www.epexspot.com/en/market-data.
  * Fraunhofer Institute for Solar Energy System ISE [2020] Fraunhofer Institute for Solar Energy System ISE, Energy-Charts (2020), Energy-Charts. 2015–2019 https://energy-charts.info/charts/price_spot_market/chart.htm.
  * SILSO World Data Center, Royal Observatory of Belgium, Brussels [2020] SILSO World Data Center, Royal Observatory of Belgium, Brussels, The International Sunspot Number (2020), international Sunspot Number Monthly Bulletin and online catalogue: 1818–2020. http://www.sidc.be/silso/.
  * Stix [2002] M. Stix,  _The Sun_, 2nd ed. (Springer-Verlag Berlin Heidelberg, 2002).
  * Balogh  _et al._ [2015] A. Balogh, H. Hudson, K. Petrovay, and R. von Steiger, eds.,  _The Solar Activity Cycle_, 1st ed. (Springer-Verlag New York, 2015).
  * Braun and Brunner [2018] S. M. Braun and C. Brunner, Zeitschrift für Energiewirtschaft 42, 257 (2018).
  * Narajewski and Ziel [2020] M. Narajewski and F. Ziel, Journal of Commodity Markets 19, 100107 (2020).
  * Simonsen [2003] I. Simonsen, Physica A: Statistical Mechanics and its Applications 322, 597 (2003).
  * Weron  _et al._ [2004b] R. Weron, M. Bierbrauer, and S. Trück, Physica A: Statistical Mechanics and its Applications 336, 39 (2004b), proceedings of the XVIII Max Born Symposium “Statistical Physics outside Physics”.
  * Laszuk [2017] D. Laszuk, Python implementation of empirical mode decomposition algorithm, https://github.com/laszukdawid/PyEMD (2017).
  * Hunter [2007] J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).


# Looking for New Physics Through the Exclusive $b\to s\nu\bar{\nu}$ Modes

Olcyr Sumensari  IJCLab, Pôle Théorie (Bat. 210), CNRS/IN2P3 et Université,
Paris-Saclay, 91405 Orsay, France

###### Abstract

The Belle-II experiment has recently measured $\mathcal{B}(B^{+}\to
K^{+}\nu\bar{\nu})$, which appears to be almost $3\sigma$ larger than its
Standard Model (SM) prediction. In this talk, I will critically revisit the
status of the SM predictions for the $B\to K^{(\ast)}\nu\bar{\nu}$ decays, and
discuss the interpretation of the recent Belle-II measurement in terms of a
general Effective Field Theory, as well as concrete models of physics beyond
the SM.

##  1 Introduction

Flavor Changing Neutral Currents (FCNCs) are powerful probes of New Physics
(NP) effects since they are loop and CKM suppressed in the Standard Model
(SM). The LHC experiments and the $B$-factories have precisely measured
several observables related to the $b\to s\mu\mu$ transition in recent years
[1]. While some of these measurements display apparent deviations from the SM
predictions [2], the interpretation of these results requires a precise
understanding of hadronic uncertainties, in particular, those related to the
so-called charm-loops that could hinder the sensitivity to NP effects [3].

Decays based on the $b\to s\nu\nu$ transition are very appealing, as they are
theoretically cleaner than the ones with muons since they are unaffected by
charm loops. They are very sensitive to the effects of physics beyond the SM,
thus offering a complementary ground to test, e.g., the NP scenarios proposed
to explain the potential discrepancies in $b\to s\mu\mu$ decays. Furthermore,
they also allow us to indirectly constrain models with large couplings to
left-handed $\tau$’s due to $SU(2)_{L}$ gauge invariance, which are difficult
to constrain otherwise [4]. Recently, the Belle-II experiment observed, for
the first time [5],

$\mathcal{B}(B^{+}\to K^{+}\nu\bar{\nu})^{\mathrm{exp}}=[2.3\pm
0.5(\mathrm{stat})^{+0.5}_{-0.4}(\mathrm{syst})]\times 10^{-5}\,,$ (1)

which turns out to be $\approx 3\sigma$ above its SM prediction [6], namely
$\mathcal{B}(B^{+}\to K^{+}\nu\bar{\nu})^{\mathrm{SM}}=(4.44\pm 0.30)\times
10^{-6}$. 111The dependence of this result on the choice of CKM inputs will be
discussed below. Note, also, that contributions from $B^{+}\to\tau^{+}(\to
K^{+}\nu_{\tau})\bar{\nu}_{\tau}$ are treated as a background in the
experimental analysis and, for this reason, are removed from our SM
predictions [7]. In the future, Belle II expects to improve this measurement,
reaching an experimental precision of $\mathcal{O}(10\%)$ with
$50\leavevmode\nobreak\ \mathrm{ab}^{-1}$ of data [8].

In this talk, we will first critically revisit the current status of the SM
predictions for the $B\to K^{(\ast)}\nu\bar{\nu}$ decays, emphasizing the
theoretical improvements needed to match the future experimental sensitivity
of Belle II [6]. Then, we will use this information to interpret the recent
$B^{+}\to K^{+}\nu\bar{\nu}$ measurement in terms of a general Effective Field
Theory (EFT) as well as concrete NP scenarios [9].

##  2 EFT approach for $B\to K^{(\ast)}\nu\bar{\nu}$ decays

Decays based on the $b\to s\nu\bar{\nu}$ transition are described by the
following effective Lagrangian,

$\displaystyle\mathcal{L}_{\mathrm{eff}}^{\mathrm{b\to
s\nu\nu}}=\dfrac{4G_{F}}{\sqrt{2}}\dfrac{e^{2}}{(4\pi)^{2}}\lambda_{t}\sum_{X=L,R}\sum_{i,j}C_{X}^{\nu_{i}\nu_{j}}\,(\bar{s}\gamma_{\mu}P_{X}b)(\bar{\nu}_{i}\gamma^{\mu}(1-\gamma_{5})\nu_{j})+\mathrm{h.c.}\,,$
(2)

where $G_{F}$ is the Fermi constant and $\lambda_{t}=V_{tb}V_{ts}^{\ast}$ is
the product of the CKM-matrix entries. In the SM, only the left-handed Wilson
coefficient appears,
$\big{[}C_{L}^{\nu_{i}\nu_{j}}\big{]}_{\mathrm{SM}}\equiv\delta_{ij}\,C_{L}^{\mathrm{SM}}$,
which is precisely known [10],

$\displaystyle C_{L}^{\mathrm{SM}}=-X_{t}/\sin^{2}\theta_{W}\,,\qquad
X_{t}=1.462(17)(2)\,,$ (3)

where NLO QCD corrections have been included [11], as well as the two-loop
electroweak contributions [12]. Using $\sin^{2}\theta_{W}=0.23141(4)$ from PDG
[1], one finally arrives at $C_{L}^{\mathrm{SM}}=-6.32(7)$, where the dominant
source of uncertainty comes from the higher order QCD corrections.

The effective Lagrangian (2) can then be used to compute the $B\to
K^{(\ast)}\nu\bar{\nu}$ branching fractions [10], which will depend on two
dominant sources of uncertainty: (i) the $B\to K^{(\ast)}$ form-factors and
(ii) the CKM matrix elements, as we discuss below.

###  2.1 Form-factors

The $B\to K^{(\ast)}$ matrix elements are parameterized in full generality in
terms of form factors which need to be computed non-perturbatively:

####  $\underline{B\to K}$:

The relevant hadronic matrix-element for the $B\to K$ transition reads

$\displaystyle\langle\bar{K}(k)|\bar{s}\gamma^{\mu}b|\bar{B}(p)\rangle$
$\displaystyle=\Big{[}(p+k)^{\mu}-\dfrac{m_{B}^{2}-m_{K}^{2}}{q^{2}}q^{\mu}\Big{]}f_{+}(q^{2})+\dfrac{m_{B}^{2}-m_{K}^{2}}{q^{2}}q^{\mu}f_{0}(q^{2})\,,$
(4)

where $q^{2}=(p-k)^{2}$, and the vector ($f_{+}$) and scalar ($f_{0}$) form-
factors satisfy the relation $f_{+}(0)=f_{0}(0)$. Only $f_{+}$ contributes to
$B\to K\nu\bar{\nu}$ since neutrino masses are negligible. These form factors
have been computed by means of Lattice QCD (LQCD) by the FNAL/MILC [13] and
HPQCD collaborations [14]. Following the same procedure as FLAG [15], we have
updated the combined fit of these form factors [6], which supersedes the
previous FLAG fit [15], in which now already obsolete HPQCD results have been
used [16]. By using these results, we can determine [6]

$\displaystyle\mathcal{B}(B\to K$
$\displaystyle\nu\bar{\nu})^{\mathrm{SM}}/|\lambda_{t}|^{2}=\left\\{\begin{matrix}(1.33\pm
0.04)_{K_{S}}\times 10^{-3}\,,&\\\\[5.0pt] (2.87\pm 0.10)_{K^{+}}\times
10^{-3}\,,&\end{matrix}\right.$ (5)

One should keep in mind, however, that these predictions result from an
extrapolation since LQCD results are precisely obtained only for $q^{2}$
values, which is a potential source of systematic uncertainty. Therefore, it
is fundamental to devise strategies to cross-check these results. A
possibility is to measure the ratio of $\mathcal{B}(B\to K\nu\bar{\nu})$ at
low- and high-$q^{2}$ bins [6], which would be a direct test of the form-
factor $q^{2}$-shapes, since the CKM inputs, the $f_{+}$ form-factor
normalization and also potential NP effects cancel out in such ratio. 222We
assume neutrinos to be purely left-handed.

####  $\underline{B\to K^{\ast}}$:

The situation for the $B\to K^{\ast}$ transition is far more intricate than
for $B\to K$ because there are more form factors in this case [10].
Furthermore, the results of only one LQCD study at nonzero recoil have been
reported so far, with a specific lattice setup [17]. In our study [6], we
consider the combination [18] of the available LQCD values [17] with those
obtained by using the LCSR at low $q^{2}$, which then give

$\displaystyle\mathcal{B}(B\to
K^{\ast}\nu\bar{\nu})^{\mathrm{SM}}/|\lambda_{t}|^{2}=\left\\{\begin{matrix}(5.9\pm
0.8)_{K^{\ast 0}}\times 10^{-3}\,,&\\\\[5.0pt] (6.4\pm 0.9)_{K^{\ast+}}\times
10^{-3}\,,&\end{matrix}\right.$ (6)

We should point out, however, that this result is less robust than the one for
$\mathcal{B}(B\to K\nu\bar{\nu})$, due to our limited knowledge of the $B\to
K^{\ast}$ form factors.

###  2.2 CKM matrix elements

The CKM factor $\lambda_{t}=V_{tb}V_{ts}^{\ast}$ introduces by far the largest
parametric uncertainty in $\mathcal{B}(B\to K^{(\ast)}\nu\bar{\nu})$. The
usual procedure, often adopted in the literature, is to determine $|V_{cb}|$
from tree-level processes and then to evaluate $|\lambda_{t}|$ using the CKM
unitarity [10]. Unfortunately, this procedure is intrinsically ambiguous since
the CKM coupling $|V_{cb}^{\mathrm{incl}}|$ extracted from the inclusive
semileptonic decay does not coincide with $|V_{cb}^{\mathrm{excl}}|$ obtained
from the exclusive modes [19]. The latest HFLAV [20] average values of the
inclusive $|V_{cb}|$ are $|V_{cb}^{\mathrm{incl}}|_{\mathrm{kin}}=(42.2\pm
0.8)\times 10^{-3}$ or $|V_{cb}^{\mathrm{incl}}|_{\mathrm{1S}}=(42.0\pm
0.5)\times 10^{-3}$, which are larger than $|V_{cb}^{B\to D}|=(40.0\pm
1.0)\times 10^{-3}$ obtained after combining [15] the experimental results on
the exclusive $B\to Dl\nu$ decays [20] (with $l=e,\mu$) with the LQCD form
factors [21]. This discrepancy remains true if one compares the inclusive
values with the one derived from $B\to D^{\ast}l\nu$ by HFLAV [20], namely
$|V_{cb}^{B\to D^{\ast}}|=(38.5\pm 0.7)\times 10^{-3}$. 333The $B\to D^{\ast}$
have been recently computed at nonzero recoil for the first time [22]. Note,
also, that a larger $|V_{cb}^{B\to D^{\ast}}|$ value has been recently
advocated [23]. In other words, there is a discrepancy in $|\lambda_{t}|$
depending on the particular input considered [6],

$|\lambda_{t}|\times 10^{3}=\left\\{\begin{matrix}41.4\pm
0.8\,,&\quad\;\,{\small(B\to X_{c}l\bar{\nu})}\\\\[5.0pt] 39.3\pm
1.0\,,&\quad{\small(B\to Dl\bar{\nu})}\\\\[5.0pt] 37.8\pm
0.7\,,&\quad\;\,{\small(B\to D^{\ast}l\bar{\nu})}\end{matrix}\right.$ (7)

where the inclusive value is about $1\sigma$ and $2\sigma$ larger the ones
derived from $B\to Dl\bar{\nu}$ and $B\to D^{\ast}l\bar{\nu}$ decays,
respectively, which again highlights the importance of clarifying the issue of
the determination of $|V_{cb}|$.

###  2.3 Summary

By combining the form-factor and CKM determinations determined above, we
obtain our final SM predictions at $\mathcal{O}(G_{F}^{2})$,

$\displaystyle\mathcal{B}(B^{+}\to K^{+}\nu\bar{\nu})$
$\displaystyle=4.44(14)(27)\times 10^{-6}\,,$
$\displaystyle\mathcal{B}(B^{0}\to K_{S}\nu\bar{\nu})$
$\displaystyle=2.05(7)(12)\times 10^{-6}\,,$ (8)
$\displaystyle\mathcal{B}(B^{+}\to K^{\ast+}\nu\bar{\nu})$
$\displaystyle=9.8(13)(6)\times 10^{-6}\,,$ $\displaystyle\mathcal{B}(B^{0}\to
K^{\ast 0}\nu\bar{\nu})$ $\displaystyle=9.1(13)(6)\times 10^{-6}\,,$

where the first uncertainty comes from the form factors and the second one is
dominated by the uncertainty of $|\lambda_{t}|$, where we have used
$|\lambda_{t}|_{\mathrm{excl}}=(39.3\pm 1.0)\times 10^{-3}$ as a reference
value [6]. However, we stress again that by using
$|\lambda_{t}|_{\mathrm{excl}}$ we obtain branching fractions that are about
$1.5\sigma$ smaller than the values we would get by using
$|\lambda_{t}|_{\mathrm{incl}}$. Therefore, it is fundamental to resolve this
puzzle to match the future $\mathcal{O}(10\%)$ sensitivity of Belle II for
these decays [8].

##  3 EFT description

###  3.1 Low-energy EFT

The EFT contributions to $B\to K^{(\ast)}\nu\bar{\nu}$ can be expressed via
the ratio

$\displaystyle\begin{split}\dfrac{\mathcal{B}(B\to
K^{(\ast)}\nu\bar{\nu})}{\mathcal{B}(B\to
K^{(\ast)}\nu\bar{\nu})^{\mathrm{SM}}}=\sum_{ij}\dfrac{|C_{L}^{\nu_{i}\nu_{i}}+C_{R}^{\nu_{i}\nu_{i}}|^{2}}{3|C_{L}^{\mathrm{SM}}|^{2}}-\eta_{K^{(\ast)}}\sum_{i,j}\dfrac{\mathrm{Re}[C_{R}^{\nu_{i}\nu_{j}}C_{L}^{\nu_{i}\nu_{j}\,\ast}]}{3|C_{L}^{\mathrm{SM}}|^{2}}\,,\end{split}$
(9)

where $C_{L}^{\nu_{i}\nu_{j}}\equiv C_{L}^{\mathrm{SM}}\,\delta_{ij}+\delta
C_{L}^{\nu_{i}\nu_{j}}$ and $C_{R}^{\nu_{i}\nu_{j}}\equiv\delta
C_{R}^{\nu_{i}\nu_{j}}$. The last term in this expressions is absent for $B\to
K\nu\bar{\nu}$ since $\eta_{K}=0$, but it contributes to $B\to
K^{\ast}\nu\bar{\nu}$, where $\eta_{K^{\ast}}=3.33(7)$ after integrating over
the full $q^{2}$-range [6, 9]. By using the above expressions, it is
straightforward to derive the correlation between $\mathcal{B}(B\to
K\nu\bar{\nu})$ and $\mathcal{B}(B\to K^{\ast}\nu\bar{\nu})$ which is depicted
in the left panel of Fig. 1. While the left-handed scenario is strongly
constrained by the upper limits on $\mathcal{B}(B\to K^{\ast}\nu\bar{\nu})$
from Belle data [24], the right-handed scenario can successfully explain the
excess observed by Belle II, as it would deplete the $B\to
K^{\ast}\nu\bar{\nu}$ rates [6, 25]. Another important prediction of the
latter scenario is a sizable modification of the $B\to K^{\ast}\nu\bar{\nu}$
longitudinal polarization fraction, $F_{L}\equiv F_{L}(B\to
K^{\ast}\nu\bar{\nu})$, which would be considerably lower than its SM value,
$F_{L}^{\mathrm{SM}}=0.49(7)$, as depicted by $\mathcal{R}_{F_{L}}\equiv
F_{L}/F_{L}^{\mathrm{SM}}$ in the right panel of Fig. 1 [6]. Therefore, these
$B\to K^{\ast}\nu\bar{\nu}$ observables at Belle II would be useful tests of
the observed deviation in $B\to K\nu\bar{\nu}$ decays.

###  3.2 SMEFT

We now assume that the NP scale lies well above the electroweak and write the
EFT invariant under $SU(3)_{c}\times SU(2)_{L}\times U(1)$, namely the SMEFT
[26]. Of all the $d=6$ operators in

$\mathcal{L}^{(6)}_{\mathrm{SMEFT}}\supset\sum_{i}\frac{\mathcal{C}_{i}}{\Lambda^{2}}{\mathcal{O}}_{i}\,,$
(10)

we select those relevant to our study, namely,

$\displaystyle\begin{split}\big{[}\mathcal{O}_{lq}^{(1)}\big{]}_{ijkl}&=\big{(}\overline{L}_{i}\gamma^{\mu}L_{j}\big{)}\big{(}\overline{Q}_{k}\gamma_{\mu}Q_{l}\big{)}\,,\qquad\big{[}\mathcal{O}_{lq}^{(3)}\big{]}_{ijkl}=\big{(}\overline{L}_{i}\gamma^{\mu}\tau^{I}L_{j}\big{)}\big{(}\overline{Q}_{k}\tau^{I}\gamma_{\mu}Q_{l}\big{)}\,,\\\\[3.99994pt]
\big{[}\mathcal{O}_{ld}\big{]}_{ijkl}&=\big{(}\overline{L}_{i}\gamma^{\mu}L_{j}\big{)}\big{(}\overline{d}_{kR}\gamma_{\mu}d_{lR}\big{)}\,,\end{split}$
(11)

where $Q$ and $L$ denote the quark and lepton $SU(2)_{L}$ doublet,
respectively, while $u,d,e$ stand for the quark and lepton weak singlets, and
$\\{i,j,k,l\\}$ denote flavor indices. 444We take the basis with a diagonal
down-type quark Yukawa matrix, i.e., with
$Q_{i}=[(V^{\dagger}\,u)_{Li}\,,\,d_{Li}]^{T}$. We will not consider the
Higgs-current operators, since they are tightly constrained by
$B_{s}-\overline{B_{s}}$ mixing and $\mathcal{B}(B_{s}\to\mu\mu)$ [9]. The
tree-level matching between the SMEFT Lagrangian to Eq. (2) gives

$\displaystyle\begin{split}\delta
C_{L}^{\nu_{i}\nu_{j}}&=\dfrac{\pi}{\alpha_{\mathrm{em}}\lambda_{t}}\dfrac{v^{2}}{\Lambda^{2}}\left\\{\big{[}\mathcal{C}_{lq}^{(1)}\big{]}_{ij23}-\big{[}\mathcal{C}_{lq}^{(3)}\big{]}_{ij23}\right\\}\,,\qquad\delta
C_{R}^{\nu_{i}\nu_{j}}=\dfrac{\pi}{\alpha_{\mathrm{em}}\lambda_{t}}\dfrac{v^{2}}{\Lambda^{2}}\big{[}\mathcal{C}_{ld}\big{]}_{ij23}\,.\end{split}$
(12)

From Fig. 1, we find that the SMEFT coefficient that best describes current
data is $\mathcal{C}_{ld}$. By considering a single lepton flavor of these
coefficients ($i=j$), we find that Belle-II data is compatible with
$|\mathcal{C}_{ld}|/\Lambda^{2}\approx(5\leavevmode\nobreak\
\mathrm{TeV})^{-2}$.

One of the advantages of considering the SMEFT is the correlation induced by
$SU(2)_{L}$ invariance between $b\to s\nu_{\ell}\bar{\nu}_{\ell}$ and $b\to
s\ell\ell$ modes, which allow us to probe the lepton flavor structure of these
coefficients. Given the stringent constraints [1] on
$\mathcal{B}(B_{s}\to\mu\mu)$ and $R_{K^{(\ast)}}$, we concluded that a
scenario that only affects electrons and muons cannot accommodate current data
[9]. Similarly, lepton flavor violating contributions ($i\neq j$) are tightly
constrained by the upper limits on $\mathcal{B}(B_{s}\to\ell_{i}\ell_{j})$ and
$\mathcal{B}(B\to K^{(\ast)}\ell_{i}\ell_{j})$. We find that the only scenario
that can successfully accommodate experimental data has $\tau$-lepton
coefficients [9]. Notice that such EFT scenario is also compatible with LHC
constraints on mono- and di-tau production at high-$p_{T}$ [27].

Finally, we comment on the realization of the above EFT scenario in terms of
concrete models [9, 25]. The simplest possibility to generate
$[\mathcal{C}_{ld}]_{3323}$ via the tree-level exchange of a heavy
$Z^{\prime}\sim(\mathbf{1},\mathbf{1},0)$ or a scalar leptoquark
$\widetilde{R}_{2}\sim(\mathbf{3},\mathbf{2},1/6)$, where we denote the SM
quantum numbers by $(SU(3)_{c},SU(2)_{L},U(1)_{Y})$. However, the $Z^{\prime}$
scenario would be tightly constrained as it contributes to $\Delta m_{B_{s}}$
at tree level, in such a way that the couplings to $\tau$’s would have to be
very large, being in potential tension with high-$p_{T}$ constraints. Instead,
the $\widetilde{R}_{2}$ leptoquark can accommodate the large value of
$\mathcal{B}(B\to K\nu\bar{\nu})$ while remaining compatible with other
constraints, provided that its mass satisfies $m_{\widetilde{R}_{2}}\lesssim
3\leavevmode\nobreak\ \mathrm{TeV}$ to suppress its loop contribution to
$B_{s}-\overline{B_{s}}$ mixing [9].

Figure 1: The correlation between $\mathcal{B}(B\to K\nu\bar{\nu})$ with
$\mathcal{B}(B\to K^{\ast}\nu\bar{\nu})$ (left panel) and
$\mathcal{R}_{F_{L}}\equiv F_{L}/F_{L}^{\mathrm{SM}}$ (right panel) are shown
with respect to the variation of $\delta C_{L}$ or $\delta C_{R}$. [6] The
shaded gray area correspond to $1\sigma$ and $2\sigma$ of the recent Belle II
result for $\mathcal{B}(B\to K\nu\bar{\nu})$. The red point corresponds to the
SM predictions for these observables. We also show the region of
experimentally excluded $\mathcal{B}(B\to K^{\ast}\nu\bar{\nu})$ values [24]
(gray hatched area), as well as the region not accessible within the EFT
approach (purple hatched area).

###  3.3 Beyond the SMEFT

Finally, we briefly comment on alternatives to explain the excess observed by
Belle II, which exploit that $B\to K+\mathrm{inv}$ could also receive
contributions from light invisible NP particles if they are sufficiently light
to be produced on-shell. These particles could be a fermion, such as a right-
handed neutrino [28], or light scalar/vector bosons [29] that could act as a
portal to a dark sector. In both cases, our EFT approach has to be extended to
account for these light degrees of freedom. The main prediction of these
scenarios is a different $q^{2}$-distribution that would, e.g., be peaked at
$q^{2}\simeq m_{X}^{2}$ for a two-body decay $B\to KX$, where $X$ is a light
boson, which could be explored by a dedicated analysis at Belle II [29].

##  4 Summary and outlook

The $B\to K^{(\ast)}\nu\bar{\nu}$ decays are powerful probes of NP effects as
they are suppressed in the SM. The main interest in studying these decays is
that they are theoretically cleaner than similar decays with muons, as they
are not affected by the so-called charm loops. The main theoretical
uncertainties of these decays are the hadronic factor and the parametric
uncertainty stemming from the determination of the CKM input
$\lambda_{t}=V_{tb}V_{ts}^{\ast}$. The latter is particularly problematic
given the long-standing discrepancy in $|V_{cb}|$ from exclusive and inclusive
decays, creating ambiguity in predicting $\lambda_{t}$ through CKM unitarity,
which must be solved to match the expected sensitivity for these decays at
Belle II.

Recently, the Belle-II experiment observed for the first time the $B^{+}\to
K^{+}\nu\bar{\nu}$ decays, with a branching fraction that is about $3\sigma$
larger than the SM predictions. We have shown that this deviation can be
easily accommodated by an EFT with operators coupled to third-generation
leptons, which is compatible with all available low- and high-energy
constraints. Several measurements can be made to independently verify this
excess, such as the measurement of $\mathcal{B}(B^{0}\to K_{S}\nu\bar{\nu})$,
as well as the determination of the branching fraction and longitudinal-
polarization asymmetry of $B\to K^{\ast}\nu\bar{\nu}$ decays, which can all be
studied at Belle II. Finally, another possibility is that the excess is due to
$B\to KX$, where $X$ denotes light and invisible particles. The experimental
signature of such a scenario would be a different $q^{2}$-distribution
kinematical distribution with respect to the SM, which could also be probed at
Belle II.

## Acknowledgments

This project has received funding /support from the European Union’s Horizon
2020 research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 860881-HIDDeN and No 101086085-ASYMMETRY.

## References

## References

  * [1] R. L. Workman et al. [Particle Data Group], PTEP 2022, 083C01 (2022).
  * [2] M. Algueró et al., Eur. Phys. J. C 79, no.8, 714 (2019); T. Hurth et al., Phys. Rev. D 103, 095020 (2021); W. Altmannshofer et al., Eur. Phys. J. C 81, no.10, 952 (2021).
  * [3] M. Ciuchini et al., JHEP 06, 116 (2016).
  * [4] L. Allwicher et al., JHEP 03 (2024), 049.
  * [5] I. Adachi et al. [Belle-II], [arXiv:2311.14647 [hep-ex]].
  * [6] D. Bečirević et al., Eur. Phys. J. C 83, no.3, 252 (2023).
  * [7] J. F. Kamenik et al., Phys. Lett. B 680, 471-475 (2009).
  * [8] E. Kou et al. [Belle-II], PTEP 2019, no.12, 123C01 (2019).
  * [9] L. Allwicher et al., Phys. Lett. B 848, 138411 (2024).
  * [10] A. J. Buras et al., JHEP 02, 184 (2015); W. Altmannshofer et al., JHEP 04 (2009), 022.
  * [11] G. Buchalla et al., Nucl. Phys. B 400, 225-239 (1993); G. Buchalla et al., Nucl. Phys. B 548, 309-327 (1999); M. Misiak et al., Phys. Lett. B 451, 161-169 (1999).
  * [12] J. Brod et al., Phys. Rev. D 83, 034030 (2011).
  * [13] J. A. Bailey et al., Phys. Rev. D 93, no.2, 025026 (2016).
  * [14] W. G. Parrott et al. [HPQCD], Phys. Rev. D 107, no.1, 014510 (2023).
  * [15] Y. Aoki et al. [FLAG], Eur. Phys. J. C 82, no.10, 869 (2022).
  * [16] C. Bouchard et al. [HPQCD], Phys. Rev. D 88, no.5, 054509 (2013).
  * [17] R. R. Horgan et al., Phys. Rev. D 89, no.9, 094501 (2014).
  * [18] A. Bharucha et al., JHEP 08, 098 (2016).
  * [19] A. J. Buras et al., Acta Phys. Polon. B 53 no.6, 6-A1; A. J. Buras et al., Eur. Phys. J. C 82 (2022) no.7, 615.
  * [20] Y. S. Amhis et al. [HFLAV], Phys. Rev. D 107, no.5, 052008 (2023).
  * [21] J. A. Bailey et al. [MILC], Phys. Rev. D 92, no.3, 034506 (2015); H. Na et al. [HPQCD], Phys. Rev. D 92, no.5, 054510 (2015).
  * [22] A. Bazavov et al. [Fermilab Lattice and MILC], Eur. Phys. J. C 82, no.12, 1141 (2022); J. Harrison et al., [arXiv:2304.03137 [hep-lat]]; Y. Aoki et al. [JLQCD], Phys. Rev. D 109 (2024) no.7, 074503.
  * [23] G. Martinelli, S. Simula and L. Vittorio, Eur. Phys. J. C 82 (2022) no.12, 1083.
  * [24] J. Grygier et al. [Belle], Phys. Rev. D 96, no.9, 091101 (2017).
  * [25] R. Bause et al., Phys. Rev. D 109, no.1, 015006 (2024); D. Marzocca et al., [arXiv:2404.06533 [hep-ph]].
  * [26] W. Buchmuller et al., Nucl. Phys. B 268 (1986), 621-653; B. Grzadkowski et al., JHEP 10 (2010), 085.
  * [27] L. Allwicher et al., JHEP 03 (2023), 064; L. Allwicher et al., high-$p_{T}$ Comput. Phys. Commun. 289 (2023), 108749.
  * [28] T. Felkl et al., Eur. Phys. J. C 83, no.12, 1135 (2023); S. Rosauro-Alcaraz et al., [arXiv:2404.17440 [hep-ph]]; Buras et al., [arXiv:2405.06742 [hep-ph]].
  * [29] W. Altmannshofer et al., Phys. Rev. D 109, no.7, 075008 (2024); D. McKeen et al., Phys. Rev. D 109, no.7, 075006 (2024); P. D. Bolton et al., [arXiv:2403.13887 [hep-ph]]; K. Fridell et al., [arXiv:2312.12507 [hep-ph]]; X. G. He et al., [arXiv:2403.12485 [hep-ph]].