Multifractal analysis of solar flare indices and their horizontal visibility graphs

. The multifractal properties of the daily solar X-ray brightness, X l and X s , during the period from 1 January 1986 to 31 December 2007 which includes two solar cycles are examined using the universal multifractal approach and multifractal detrended ﬂuctuation analysis. Then we convert these time series into networks using the horizontal visibility graph technique. Multifractal analysis of the resulting networks is performed using an algorithm proposed by us. The results from the multifractal analysis show that multifractality exists in both raw daily time series of X-ray brightness and their horizontal visibility graphs. It is also found that the empirical K(q) curves of raw time series can be ﬁtted by the universal multifractal model. The numerical results on the raw data show that the Solar Cycle 23 is weaker than the Solar Cycle 22 in multifractality. The values of h( 2 ) from multifractal detrended ﬂuctuation analysis for these time series indicate that they are stationary and persistent, and the correlations in the time series of Solar Cycle 23 are stronger than those for Solar Cycle 22. Furthermore,


Introduction
An important aim of solar-terrestrial physics is to understand the causes of geomagnetic activity in general and geomagnetic storms in particular.Since solar flares (using X-ray measurements from GOES) are coincident with many coronal mass ejections (CMEs) (see, for example, Zhang et al., 2007), they are useful for prediction of geomagnetic storms (Park et al., 2002;Yermolaev et al., 2005) due to the shorter propagation times of solar photons.Commonly the studies on solar flares (e.g.,Howard and Tappin, 2005) have focused on a relatively small number of events, ∼ 10 per year, with large magnitudes.
Fractal methods can be used to characterise the scaling properties in each time series.Multifractals are a broad generalisation of the (geometrical) fractals.They are not only more general but also fundamental (Schertzer and Lovejoy, 2011).Multifractal analysis was initially proposed to treat turbulence data and is a useful way to characterise the spatial heterogeneity of both theoretical and experimental fractal patterns (Grassberger and Procaccia, 1983;Halsy et al., 1986).It has been applied successfully in many different fields including financial modelling (e.g., Anh et al., 2000;Canessa, 2000), biological systems (e.g., Yu et al., 2001Yu et al., , 2003Yu et al., , 2004Yu et al., , 2006;;Anh et al., 2001Anh et al., , 2002;;Zhou et al., 2005), geophysical systems (e.g., Schertzer and Lovejoy, 1987;Schmitt et al., 1992;Tessier et al., 1993Tessier et al., , 1996;;Olsson, 1995;Olsson and Niemczynowicz, 1996;Harris et al., 1996;Lovejoy et al., 1996;Deidda, 2000;Lilley et al., 2006;Kantelhardt et al., 2006;Veneziano et al., 2006;Venugopal et al., 2006;Lovejoy andSchertzer, 2006, 2010a, b;Garcia-Marin et al., 2008;Serinaldi, 2010) and high energy physics (e.g., Ratti et al., 1994).Fractal and multifractal approaches have been quite successful in extracting salient features of physical processes responsible for the near-Earth magnetospheric phenomena (Lui, 2002).As solar observational techniques improve, fine small-scale structures observed on the solar surface become more pronounced.Abramenko (2005) proposed a scaling of structure functions to Z. G. Yu et al.: Multifractal analysis of solar flare indices and their HVGs analyse multifractality and found that flare-quiet regions tend to possess a lower degree of multifractality than flaring active regions do.A method to describe the multiple scaling of the measure representation of the D st time series was provided in Wanliss et al. (2005).A prediction method based on the recurrent iterated function system in fractal theory was detailed in Anh et al. (2005) together with some evaluation of its performance.A two-dimensional chaos game representation of the D st index for prediction of geomagnetic storm events was proposed in Yu et al. (2007).Yu et al. (2009) used both multifractal detrended fluctuation analysis (MF-DFA) proposed by Kantelhardt et al. (2002) and traditional multifractal analysis to study the scaling properties of D st , a p and the solar X-ray measurements.Our group used multifractal analysis and fractional stochastic differential equations to study the AE data and geomagnetic field data (Anh et al., 2007(Anh et al., , 2008;;Yu et al., 2010).
Complex networks have been studied extensively due to their relevance to many real-world systems such as the worldwide web, the internet, energy landscapes, and biological and social systems (Song et al., 2005).After analyzing a variety of real complex networks, Song et al. (2005) found that they consist of self-repeating patterns on all length scales, i.e., they have self-similar structures.In order to unfold the self-similar property of complex networks, Song et al. (2005) calculated their fractal dimension, a known useful characteristic of complex fractal sets (Mandelbrot, 1983;Falconer, 1997), and found that the box-counting method is a proper tool for further investigations of network properties.Because a concept of metric on graphs is not as straightforward as the Euclidean metric on Euclidean spaces, the computation of the fractal dimension of networks via a box-counting approach is much more complicated than the traditional boxcounting algorithm for fractal sets in Euclidean spaces.Song et al. (2007) developed a more involved algorithm to calculate the fractal dimension of complex networks.Lee and Jung (2006) found that the behaviour of complex networks is best described by a multifractal approach.Our group proposed a new box-covering algorithm to compute the generalised fractal dimensions of a network (Wang et al., 2012).
Recent works have used network techniques to investigate time series (Donner et al., 2011;and the references therein).Inspired by the concept of visibility (de Berg et al., 2008), Lacasa et al. (2008) suggested a simple computational method to convert a time series into a graph, known as a visibility graph (VG).The constructed graph inherits several properties of the series in its structure.Thereby, periodic time series convert into regular graphs, and random series into random graphs.Moreover, fractal time series convert into scale-free networks, enhancing the fact that a power-law degree distribution of its graph is related to the fractality of the time series.These findings suggest that a visibility graph may capture the dynamical fingerprints of the process that generates the time series.Elsner et al. (2009) used the visibility network to study hurricanes in the United States.Then Luque et al. (2009) proposed the horizontal visibility graphs (HVG) which are geometrically simpler and form an analytically solvable version of VG.Xie and Zhou (2011) studied the relationship between the Hurst exponent of fractional Brownian motion and the topological properties (clustering coefficient and fractal dimension) of its converted HVG.
In this paper, we examine the multifractal properties of the daily solar X-ray brightness, X l (i.e., 1-8 Å X-rays (Watts m −2 )) and X s (i.e., 0.5-4 Å X-rays (Watts m −2 )) , during the period from 1 January 1986 to 31 December 2007, including two solar cycles (cycle 1 and cycle 2 corresponding to Solar Cycles 22 and 23 respectively), using the universal multifractal approach (Schertzer and Lovejoy 1987) and MF-DFA (Kantelhardt et al. 2002).Then we convert these time series into networks using the HVG technique proposed by Luque et al., (2009).In our recent paper (Wang et al. 2012), we proposed a new algorithm to perform multifractal analysis on different types of networks.We will apply this algorithm on the resulting HVGs of X l and X s time series.The results from these multifractal analyses confirm the existence of multifractality in the time series of X l , X s and their HVGs.

Universal multifractal approach
When a cascade proceeds over a scale ratio λ = L/ l (i.e., the ratio of the largest scale of interest to the smallest scale, L being a fixed external scale, l varying from 1 to L), we denote by λ the density of the conserved energy flux.Then its statistical moments will have the following scaling behaviour as λ → ∞ (or l → 0) (Schertzer and Lovejoy, 1987): where indicates ensemble averaging.If the curve K(q) versus q is a straight line, the data set is monofractal.However, if this curve is convex, the data set is multifractal (e.g., Garcia-Marin, 2008).The scaling of the moments can be assessed by computing M q at different scales, and plotting M q against the scale ratio λ in a log-log plane, where the powerlaw relation in Eq. ( 1) becomes linear (Serinaldi, 2010).Hence the empirical K(q) functions can be estimated from the slopes of M q against the scale ratio λ in a log-log plane.
The universal multifractal model proposed by Schertzer and Lovejoy (1987) assumes that the generator of multifractals was a random variable with an exponentiated extremal Lévy distribution.Thus, the theoretical scaling exponent function K(q) for the moments q ≥ 0 of a cascade process is obtained according to (Schertzer and Lovejoy, 1987;Ratti et al., 1994;Garcia-Marin et al., 2008;Serinaldi, 2010) in which the most significant parameter α ∈ [0, 2] is the Lévy index, which indicates the degree of multifractality (i.e., the deviation from mono-fractality).The values α = 0 and α = 2 correspond to the beta model (monofractal) and log-normal model (multifractal), respectively.C 1 ∈ [0, d], with d being the dimension of the support (d=1 in our case), describes the sparseness or inhomogeneity of the mean of the process (Garcia-Marin et al., 2008).The parameter H is called the non-conservation parameter since H = 0 implies that the ensemble average statistics depend on the scale, while H = 0 is a quantitative statement of ensemble average conservation across the scales (e.g., Ratti et al., 1994;Serinaldi, 2010).
The parameters C 1 and α can be estimated by applying the double trace moment (DTM) technique (Schmitt et al., 1992;Lavallee et al., 1993).From the estimated values of α and C 1 , and taking the value of exponent β that characterises the energy spectrum of the conserved process E(ω) ≈ ω −β , with ω being the frequency, the parameter H is given by (e.g., Lavallee et al., 1993;Serinaldi, 2010): Although the above method has been widely used to estimate the parameters H , C 1 and α in geophysical research, it is complicated and the goodness of fit of the empirical K(q) functions depends on the fit for β, and sometimes the fitting of K(q) is not satisfactory (e.g., Olsson and Niemczynowicz, 1996;Garcia-Marin et al., 2008;Serinaldi, 2010).In this paper we adopt a method which is similar to that proposed in Anh et al. (2001): If we denote K T (q) the K(q) function defined by Eq. ( 2), and K d (q) the empirical K(q) function.
We estimate the parameters by solving the least-squares optimisation problem min

Multifractal detrended fluctuation analysis
The traditional multifractal analysis has been developed for the multifractal characterisation of normalised, stationary time series.This standard formalism does not give correct results for non-stationary time series which are affected by trends.Multifractal detrended fluctuation analysis (MF-DFA), which is a generalisation of the standard detrended fluctuation analysis (DFA), is based on the identification of the scaling of the qth-order moments of the time series, which may be non-stationary (Kantelhardt et al., 2002).DFA has been used to study the classification problem of protein secondary structures (Yu et al., 2006).Movahed et al. (2006) used the MF-DFA to study sunspot fluctuations.We now summarise the MF-DFA technique.Consider a time series {X 1 , X 2 , ..., X N } of length N .For an integer s ≥ 0, we divide the time series into [N/s] segments of equal length s, where [N/s] is the integer part of N/s.In each segment j , we compute the partial sums Y (i) = i k=1 X k , i = 1, 2, ..., s, fit a local trend y j (i) to Y (i) by least squares, then compute the sample variances of the residuals: (5) Note that linear, quadratic, cubic or higher order polynomials y j (i) can be used in the local trend fitting, and the DFA is accordingly called DFA1, DFA2, DFA3,...In the following we use only DFA1.
The qth-order fluctuation function is then defined as the average over all segments: Since the segments are all of the same length, the secondorder fluctuation function F 2 (s) is equivalent to the sample variance of the entire series.This is not so for the general case q = 2.We will assume that F q (s) is characterised by a power law: The scaling function h (q) is then determined by the regression of log F q (s) on log s in some range of time scale s.
For fractional Brownian motion, Movahed et al. (2006) showed that the Hurst index H 1 = h (2) − 1.Using this relationship (or H 1 = h (2) for the stationarity case) and the estimate of h (2) from the regression of log F 2 (s) on log s, an estimate of the Hurst index H 1 , and hence the extent of long memory in the time series, is obtained.For Brownian motion (with uncorrelated increments), the scaling exponent H 1 is equal to 1/2.The range 1/2 < H 1 < 1.0 indicates the presence of long memory (persistence), while the range 0 < H 1 < 1/2 indicates short memory (anti-persistence).

Visibility graph and horizontal visibility graph of a time series
A graph (or network) is a collection of nodes, which denote the elements of a system, and links or edges, which identify the relations or interactions among these elements.
Inspired by the concept of visibility (de Berg et al., 2008), Lacasa et al. (2008) suggested a simple computational method to convert a time series into a graph, known as visibility graph (VG).A visibility graph is obtained from the mapping of a time series into a network according to the following visibility criterion: Given a time series {x 1 , x 2 , ..., x N }, two arbitrary data points x t a and x t b in the time series have visibility, and consequently become two thus a connected, unweighted network could be constructed based on a time series and is called its horizontal visibility graph.As a matter of factor, for a given time series, its horizontal visibility graph is always a subgraph of its visibility graph (Luque et al., 2009).Xie and Zhou (2011) studied the relationship between the Hurst exponent of fractional Brownian motion and the topological properties (clustering coefficient and fractal dimension) of its converted HVG.

Multifractal analysis of complex networks
The most common algorithms of traditional multifractal analysis are the fixed-size box-counting algorithms (Halsy et al., 1986).For a given measure µ with support E in a metric space, we consider the partition sum q ∈ R, where the sum is evaluated over all different nonempty boxes B of a given size in a grid covering of the support E. The exponent τ (q) is defined by and the generalised fractal dimensions of the measure are defined as D(q) = τ (q)/(q − 1), for q = 1, (10) and where Z 1, = µ(B) =0 µ(B) ln µ(B).The generalised fractal dimensions are numerically estimated through a linear regression of (ln Z (q))/(q − 1) against ln for q = 1, and similarly through a linear regression of Z 1, against ln for q = 1.The D(q) corresponding to positive values of q give relevance to the regions where the measure value is large.The D(q) corresponding to negative values of q deal with the structure and the properties of the regions where the measure value is small.Our group proposed a new box-covering algorithm to compute the generalised fractal dimensions of a network (Wang et al., 2012).For a network, we denote the matrix of shortest path lengths by B = (b ij ) N×N , where b ij is the length of the shortest path between nodes i and j .Then we use B = (b ij ) N×N as input data for multifractal analysis based on our modified fixed-size box counting algorithm as follows: i.Initially, all the nodes in the network are marked as uncovered and no node has been chosen as a seed or centre of a box.
ii Set t = 1, 2, ..., T appropriately.Group the nodes into T different ordered random sequences.More specifically, in each sequence, nodes which will be chosen as seed or centre of a box are randomly arrayed.Remark: T is the number of random sequences and is also the value over which we take the average of the partition sum Z r (q).In this study, we set T = 1000 for all the networks in order to compare them.
iii.Set the size of the box in the range r ∈ [1, d], where d is the diameter of the network.
Remark: When r = 1, the nodes covered within the same box must be connected to each other directly.When r = d, the entire network could be covered in only one box no matter which node was chosen as the centre of the box.
iv.For each centre of a box, search all the neighbours within distance r and cover all nodes which are found but have not been covered yet.
v. If no newly covered nodes have been found, then this box is discarded.
vi.For the nonempty boxes B, we define their measure as µ(B) = N B /N, where N B is the number of nodes covered by the box B, and N is the number of nodes of the entire network.
vii.Repeat (iv) until all nodes are assigned to their respective boxes.
viii.When the process of box counting is finished, we calculate the partition sum as Z r (q) = µ(B) =0 [µ(B)] q for each value of r.
ix. Repeat (iii) and (iv) for all the random sequences, and take the average of the partition sums Z r (q) = ( t Z r (q))/T , and then use Z r (q) for linear regression.
Linear regression is an essential step to get the appropriate range of r ∈ [r min , r max ] and to get the generalised fractal dimensions D q .In our approach, we run the linear regression of [ln Z r (q)]/(q −1) against ln(r/d) for q = 1, and similarly the linear regression of Z 1,r against ln(r/d) for q = 1, where Z 1,r = µ(B) =0 µ(B) ln µ(B) and d is the diameter of the network.

Results and discussion
This section examines the multifractal properties of the daily solar X-ray brightness, X l and X s , during the period from 1 January 1986 to 31 December 2007 (including two solar cycles) and their horizontal visibility graphs.
The solar X-ray data are from the GOES space environment monitors on GOES 6,7,8,9,10,11,and 12. Ion chamber detectors are used to provide whole-sun Xray fluxes for the 0.5-to-3 (0.5-to-4 prior to GOES-8) and 1-to-8 wavelength bands.These bands are referred to as the X s and X l , respectively.Hourly measurements were downloaded from the National Geophysical Data Center (NGDC, http://spidr.ngdc.noaa.gov/spidr/index.jsp)and combined, using the more recent measurements to fill, whenever possible, any gaps in the earlier ones.No attempt was made to compensate for differences in calibration between the measurements or to average them.By using measurements from all the satellites, gaps in the observations are reduced significantly.During the period covered by GOES 8 (1 March 1995 to 30 June 2003) the gaps in the measurements were reduced significantly, from 1496 to 77 h for the longer wavelength, X l , observations.Most of the remaining gaps span multiple hours, even a full 24 h.Since there are usually three satellites providing observations, the remaining gaps are probably the result of geomagnetic storm effects at Earth.Although the solar measurements may be missing during a storm, any flare(s) associated with a storm is typically observed, since it occurred hours earlier.The daily time series of solar X-ray brightness, X l and X s , are shown in Fig. 1.
We divide the raw daily X l and X s data into two time series, one for each solar cycle in the data.First we perform the universal multifractal analysis on the four time series.For calculating the empirical K(q) for the time series, we use the MATLAB program "TraceMoment.m"provided by S. Lovejoy at the web site http://www.physics.mcgill.ca/∼ eliasl/.An example for obtaining the multifractal function K(q) is shown in Fig. 2. From the plots of M q against the scale ratio λ in a log-log plane, we find the best linear fit range is from λ = 7.14 (around 1 week) to λ = 365.27(around 1 yr).The empirical K(q) curves of these time series are given in Fig. 3  Fig. 2.An example for obtaining the empirical K(q) function.From the plots, we find the best linear fit range is from λ = 7.14 (around 1 week) to λ = 365.27(around 1 yr).
(the dotted lines).In order to use the universal multifractal model (i.e., Eq. 2) to fit the empirical K(q) curves, we use the function fminsearch in MATLAB to solve the optimisation problem (Eq.4) and obtain the estimates of H , C 1 and α (we set 0.5, 0.5, 0.5 as the initial values of these three parameters, respectively).The estimated values of these three parameters are given in Table 1.We also plot the theoretical K(q) curves in Fig. 3 (the continuous lines).From Fig. 3, it can be seen that the universal multifractal model fits the empirical K(q) curves very well.From Table 1, we find that all the values of α are larger than 1.0 and smaller than 2.0, indicating that the raw daily X l and X s are multifractal.We also find that the values of α for X l are larger than those for X s , The continuous lines the fitted theroretical K(q) Fig. 3.The K(q) curves of the raw daily X l and X s data (the dotted curves), and their fitted curves (continuous lines) by the universal multifractal model.and the values of α for data in cycle 1 are larger than those for cycle 2.This fact indicates that the multifractality of X l is stronger than that of X s , and the multifractality of cycle 1 is stronger than that of cycle 2.The values of H for X l are close to zero, indicating that they correspond to a conservative field.
We also perform MF-DFA on the four time series.An example for obtaining the exponent h(2) in MF-DFA is shown in Fig. 4. The numerical results on the h(q) curves are shown in Fig. 5.The values of h(2) for these time series are also given in Table 1: they are all larger than 0.5 and smaller than 1.0, indicating that these time series are stationary and persistent.The h(2) values for cycle 2 are larger than those for cycle 1, indicating that the correlations in the time series of cycle 2 are stronger than those of cycle 1.The nonlinearity of the h(q) curves in Fig. 5 also confirms that the raw daily X l and X s are multifractal.The h(q) curves of data in cycle 2 are flatter than those in cycle 1, indicating that the multifractality reflected by the h(q) curve of the time series in cycle 2 is weaker than that in cycle 1.
De Toma et al. (2004) noted that Solar Cycle 23 (cycle 2 here) is weaker than Solar Cycle 22 (cycle 1 here) in most solar activity indices, including magnetic flux; they proposed that it was a distinct, magnetically simpler variant from pre- The slope is 0.6937 ± 0.0117 The slope is 0.6834 ± 0.0201 Fig. 5.The h(q) curves of the raw daily X l and X s data.vious cycles.Recently Kossobokov et al. (2012) found that the length (13.2 yr based on flares) and maximum number of days between solar flares (466 days) in Solar Cycle 23 are longer and larger than those (9.25 yr based on flares, 157 days) in Solar Cycle 22, respectively.They also found Cycle 23 has the longer quiet period.These differences can be reflected in the X l and X s time series and will affect our estimated K(q) and h(q) curves.Our results show that the Solar Cycle 23 is weaker than Solar Cycle 22 in multifractality reflected by the K(q) and h(q) curves based on raw data.
To gain more insight into this aspect, we next convert daily X-ray data (four time series) into their visibility graphs and horizontal visibility graphs.Because there are too many edges in the visibility graphs, their diameters are relatively small (less than 8).Hence it is not meaningful to study the fractal property of these visibility graphs.The diameters of For X l cycle 2 Fig. 6.The D(q) curves of the converted horizontal visibility graphs of the X l and X s data.
the horizontal visibility graphs (HGV) are much larger and it is meaningful to study their fractal and multifractal properties.Hence we performed multifractal analysis on the HGVs using our algorithm.The estimated D q curves of the HGVs of the four time series are shown in Fig. 6.These D q curves again confirm the multifractality of the HGVs.Furthermore, from the D q curves, we can see the multifractality, which is characterised by of X s is stronger than that of X l .This assertion is different from the multifractality reflected by the α value in the universal multifractal model and the h(q) curve for the raw data (the time series point of view).Hence network analysis of the time series reflects some properties which are not shared by the same analysis on the original time series.This suggests a potentially useful method to explore geophysical data.

Remark
As claimed in our previous work (Wang et al., 2012), we considered the generalised fractal dimensions D q to determine whether the object is multifractal from the shape of D q .In our results, an anomalous behaviour is observed: the D q curves increase at the beginning.This anomalous behaviour has also been observed in Opheusden et al. (1996), Smith andLange (1998), andFernández et al. (1999).Some reasons for this behaviour have been suggested, including that the boxes contain few elements (Fernández et al., 1999), or the small scaling regime covers less than a decade so that we cannot extrapolate the box counting results for the partition function to zero box size (Opheusden et al., 1996).Hence the results from D(q) for larger q, which give relevance to the regions where the measure value is large (hubs in the network), are more convincing in our case.

Conclusions
Multifractal analysis is a useful way to characterise the spatial heterogeneity of both theoretical and experimental fractal patterns.The numerical results from the universal multifractal approach and MF-DFA on the raw daily X l and X s data show that these time series are multifractal.The MF-DFA method shows that the multifractality of the time series in cycle 2 is weaker than that in cycle 1.It is found that the empirical K(q) curves of raw time series can be fitted very well by the universal multifractal model.The estimated values of α in this model suggest that the multifractality of the X l time series is more severe than that of the X s time series.The estimated values of H in the universal multifractal model show that X l corresponds to a conservative field.The values of h(2) from MF-DFA for these time series indicate that they are stationary and persistent, and the correlations in the data of cycle 2 are stronger than those of cycle 1.
The estimated D q curves of the horizontal visibility graphs of the four time series confirm their multifractality.The multifractality of X s is stronger than that of X l , which is different from the multifractality reflected by the α value in the universal multifractal model and the h(q) curve for raw data (the time series point of view).Hence network analysis of the time series reflects some properties which are not shared by the same analysis on the original time series.This suggests a potentially useful method to explore geophysical data.

Fig. 1 .
Fig.1.The daily solar X-ray X l and X s , during the period from 1 January 1986 to 31 December 2007 which includes two solar cycles.

Fig. 4 .
Fig. 4. Examples for obtaining the exponent h(2) in MF-DFA.The linear fit range is s = 5 to 98.

www.nonlin-processes-geophys.net/19/657/2012/ Nonlin. Processes Geophys., 19, 657-665, 2012 660 Z. G. Yu et al.: Multifractal analysis of solar flare indices and their HVGs connected
Luque et al. (2009)8)ed graph, if any other data point x t c such that t a < t c < t b fulfilsx t c < x t a + (x t b − x t a ) t c − t a t b − t a ;thus a connected, unweighted network could be constructed based on a time series and is called its visibility graph.It has been shown(Lacasa et al., 2008)that time series structures are inherited in the associated graph, such that periodic, random, and fractal series map into motif-like random exponential and scale-free networks, respectively.ThenLuque et al. (2009)proposed horizontal visibility graphs (HVG) which are geometrically simpler and analytically solvable version of VG.Given a time series {x 1 , x 2 , ..., x N }, two arbitrary data points x t a and x t b in the time series have horizontal visibility, and consequently become two connected nodes in the associated graph, if any other data point x t c such that t a < t c < t b fulfils x t c < min{x t a , x t b };

Table 1 .
The estimated values of H , C 1 and α in the universal multifractal model and h(2) in the MF-DFA for the daily solar X-ray data.Here error means the minimal value in Eq. (4).