Analysis of the structural complexity of Crab Nebula observed at radio frequency using a multifractal approach

The Crab Nebula is an astrophysical system that exhibits complex morphological patterns at different observing frequencies. We carry out a systematic investigation of the structural complexity of the nebula using publicly available imaging data at radio frequency. For the analysis, we use the well-known multifractal detrended fluctuation analysis in two dimensions. We find that radio data exhibit long-range correlations, as expected from the underlying physics of the supernova explosion and evolution. The correlations follow a power-law scaling with length scales. The structural complexity is found to be multifractal in nature, as evidenced by the dependence of the generalized Hurst exponent on the order of the moments of the detrended fluctuation function. By repeating the analysis on shuffled data, we further probe the origin of the multifractality in the radio imaging data. For the radio data, we find that the probability density function is close to a Gaussian form. Hence, the multifractal behavior is due to the differing nature of long-range correlations of the large and small detrended fluctuation field values. We investigate the multifractal parameters across different partitions of the radio image and find that the structures across the image are highly heterogeneous, making the Crab Nebula a structurally complex astrophysical system. Our analysis thus provides a fresh perspective on the morphology of the Crab Nebula from a complexity science viewpoint.

Natural systems are generally complex, consisting of a large number of interlinked/integrated components of heterogenous properties [1].The individual components of such complex systems have non-linear character [2,3].Their interaction can give rise to emergent behaviors at different scales.An important property of complex systems is self-organization [4], which can be captured by fractals or multifractals.The notion of multifractality originated in the study of turbulence in fluid mechanics [5,6].Classically, multifractals are analyzed using two approaches: structure-function [7] and partition function [8].While the structure-function approach is based on moments of increment distribution, the partition-function approach mainly adopts the boxcounting formalism to determine the multifractal spectrum (α, f ).Other approaches to analyzing multifrac-tals in time-series data include the wavelet transform approach, fluctuation analysis approach, and detrended fluctuation approach.We refer interested readers to the excellent review [9] for details of these approaches.The well-established multifractal formalism is used across a wide variety of fields.It is used to study long-term persistence in river precipitation (earth sciences) [10], dynamics of stock markets (econophysics) [9], anomalies in the dynamics of brain (neuroscience) [11], social media activities (social dynamics) [12] and species abundance distributions (ecology) [13].Kantelhardt et al. [14] developed a numerical method known as the multifractal detrended fluctuation analysis (MFDFA) to calculate various multifractal parameters for characterizing multifractal scaling properties and detecting long-range correlations in noisy, non-stationary time series data of complex systems.The method is based on a generalization of the detrended fluctuation analysis (DFA) initially given by Peng et al [15].Gu and Zhou [16] extended the conventional MFDFA to the two-dimensional multifractal detrended fluctuation analysis (2D-MFDFA) to analyze higher dimensional fractals and multifractals.
The multifractal theoretical framework has long been studied by theoretical physicists and mathematicians.The multifractal spectrum, f (α), also known as Hölder spectrum, associates with each positive α the Hausdorff arXiv:2206.04717v2[nlin.AO] 9 Feb 2024 dimension of the set [17].For instance, in the case of fully developed turbulence, the Hölder spectrum is determined for the velocity of the fluid [18].The formalism of fractal or multifractal analysis [19] runs parallel to mathematical frameworks, such as multifractal decomposition of Moran fractals [20], and digraph recursive fractals [21], maximal measure [22], discrete measures [23], the theory of large deviations [24], Minkowski measurability [25], self-similar functions [26,27], etc.The present paper aims to introduce and explore a fractal theory and a multifractal analysis in the context of a complex astrophysical system, namely, the Crab Nebula.Astrophysics provides a multitude of examples of complex systems owing to the different physical interactions and evolution time scales involved.Examples of complex astrophysical systems include turbulent plasma, such as stellar flares and accretion discs, supernova remnants, spatial structures in the interstellar medium, galaxies, and the large-scale structure of the Universe [28][29][30].To measure the complexity of their dynamics or morphological structures from observational data is an intriguing and active research topic.Supernova remnants are particularly interesting for questions related to complexity because of their complicated morphological patterns that evolve with time.In this paper, we focus our attention on the Crab Nebula, a well-known remnant of a supernova recorded in 1054 AD [31].It is located about 6500 light-years from the earth [32] and has been well-studied using multi-wavelength observations from radio frequencies to gamma rays (see reviews [33,34]).We choose to study the Crab Nebula from the viewpoint of complex systems because of its historical importance and availability of imaging data with good spatial resolution in several frequency bands.Several studies have been carried out to investigate the morphological variation of the Crab Nebula across different frequencies (see e.g.[35][36][37]).However, a quantitative analysis of the morphological complexity of the Crab Nebula is lacking in the literature.
The intricate morphology of the Crab Nebula that can be easily discerned by the eye from high-resolution images naturally leads to the question of its structural complexity.A complex systems approach to the Crab Nebula views the physical system as one system of interconnected components with their organization exhibiting emergent behaviors instead of studying the dynamics of the individual components.To investigate the structural complexity of the complex system of the Crab Nebula, we adopt a multi-scale multifractal approach.A multifractal approach is very well-suited to study and analyze the structural properties of the Crab Nebula since it provides a framework for detecting and identifying complex local structures and describing local singularities.In astrophysics, multifractal formalism is applied to study the large scale structure of the universe [38,39], the behavior of galaxy clustering on large scales [40], the interstellar medium [41], Gamma-ray bursts [42,43], flux correlations of pulsars [44], flux variability in the quasar 3C 273 [45], gravitational wave signals detected by LIGO [46], to name a few.The well-established onedimensional MFDFA is used to study sunspot number fluctuations [47], investigate the statistical properties of the cosmic microwave background radiation [48], and analyze the pseudorapidity distribution of ring-like events of charged mesons [49].We systematically investigate the structural complexity of the Crab Nebula observed at radio frequency using publicly available data with the 2D-MFDFA method.Our findings provide the local scaling behavior, nature of correlations, the origin of multifractality, heterogeneous properties, and organization in the structural system of the Crab Nebula.Our results deepen our understanding of the multi-scale physics operating in the Crab Nebula.This paper is organized as follows.Section II outlines the methodology used in our analysis.We briefly introduce the concept of multifractality and describe the 2D-MFDFA algorithm and the associated physical interpretations.Section III describes the physics of the Crab Nebula and the imaging data used in our analysis.In section IV, we present our main analysis and results.We conclude with a summary and discuss the implications of our results in section V.

II. METHODOLOGY
In this section, we briefly review multifractality and complexity measures, then present the 2D-MFDFA algorithm and its associated physical interpretations.

A. Complexity and multifractality
Given y(0) = 0, a random process {y(t)} is said to be a self-affine process if it obeys the following scaling relation [50] where c is the scale factor.The scaling exponent D > 0 represents the self-affine process's fractal or self-similarity dimension.A uni-fractal or uniscaling process is characterized by a single scaling law at any scale.Data generated by various complex systems exhibit fluctuations on a broad range of scales.Such fluctuations often follow a scaling relation of the type (1) over several orders of magnitude in both equilibrium and non-equilibrium situations [14].
The theory of multifractals studies a more general relation of the type where y(t) and Γ(c) are independent random functions.
If the random process {y(t)} is multifractal, then the scaling function Γ(c) satisfies [50] Γ(c where Γ 1 , Γ 2 , . . .Γ n are n independent copies of Γ for various local scales c 1 , c 2 , . . ., c n .Using eq. ( 1) and ( 2), each local scale c κ of eq. ( 3) has the local fractal dimension D κ and hence the local scaling function follows the relation The multifractal process becomes a monofractal when ).Using eq. ( 4), we can write the scaling function Γ of eq. ( 3) as Γ = c D1+D2+  We briefly review here the algorithm and physical interpretations of the two-dimensional MFDFA method, following [14,16] 1. 2D-MFDFA algorithm Consider a two-dimensional compact space S, which is a subset of flat two-dimensional space.For our purpose, S is a rectangular region that is subdivided into M × N equal-area square pixels.Let d denote the data field, a function on S. The d is thus a two-dimensional M × N array.
The algorithm for the 2D-MFDFA involves the following steps.
1. Partition S into M s × N s disjoint square subsets, with each subset containing s 2 pixels.We refer to s as the scale size.Let us call each subset a superpixel.Let each superpixel be indexed by I, J, such that I = 1, 2, . . ., M s and J = 1, 2, . . ., N s .We assume that M, N are multiples of s so that M s = M/s and N s = N/s are positive integers.We then carry out the partitioning for different values of s.For a given s, all superpixels contain the same number of pixels.
2. Let d I,J denote the s × s data array within each superpixel I, J.The cumulative sum U I,J of the data for each superpixel is defined by where i, j = 1, . . ., s are indices for the pixels within the superpixel.
3. To subtract the trend in each superpixel, we first fit U I,J (i, j) with a linear bivariate polynomial function U fit I,J (i, j), given by where A, B, C are parameters determined by the fitting process.In general, one can also use higherorder polynomials for the fitting function.The detrended fluctuation function that is defined in the steps below will differ based on the order of the fitting function.We have chosen a linear form since it is the simplest.Even though the numerical values of our derived quantities are affected by the order of the fitting polynomial, the overall conclusions of our results remain the same.We then calculate the residual ∆U I,J (i, j) = U I,J (i, j) − U fit I,J (i, j).
Then, for a given s, the detrended fluctuation function of each superpixel, which we denote by G(I, J, s), is calculated as 4. Let q be a non-zero real number.The detrended fluctuation function of order q, denoted by F q (s), is calculated by averaging over all the superpixels as For positive integer values of q, this expression reduces to the usual definition of the q th order moment of G(I, J, s).
For the case of q = 0, eq. ( 9) cannot be used since it diverges.Instead, we use the following expression for F 0 (s) Since we use the positive square root of eq. ( 8), F q (s) is always real.

5.
For different values of q, the scaling relation between F q (s) and the scale size s can be expressed as where the scaling exponent h(q) is the generalized Hurst exponent.The h(q) is a measure of selfsimilarity symmetry and correlation present in the data field.We restrict our attention here to data for which F q (s) are increasing functions of s, or positive values of h(q).This means that the data of interest contains positive correlations.
Mathematically, there is a-priori no restriction on the range of q values.However, the physical information that can be obtained may saturate beyond some range of q values, such as what we will find when h(q) asymptotes to constant values as q increases.

Physical interpretations
As described above the 2D-MFDFA algorithm determines the scaling relation between F q (s) and s, from which we can calculate h(q).For a general data field, h can depend on both q and s.From the behaviour of h(q), we can infer the following: • If the data contains long-range power-law correlations, then the dependence of F q (s) on s has a power-law form, and h(q) is independent of s.
• If h(q) remains the same for varying q, then the data field has monofractal scaling property.This means that the data or physical system contains one structural element or pattern that is invariant under size scale transformations.Then only one scaling exponent, namely h describes the selfsimilar scaling behavior of the system.
• If, however, h has a dependence on q, then the data has multifractal scaling behaviour.This implies that the data contains more than one structural element, which is invariant under size scale transformations.This further shows that small and large fluctuations scale differently with s.
• For positive values of q, the dominant contributions to the sum on the right-hand side of eq. 9 will come from superpixels containing large fluctuations (or equivalently large deviations from the fitting function U fit I,J given in step 3 in the previous subsection).So, for q > 0, h(q) effectively describes the scaling behavior of large fluctuations.In contrast, for negative values of q, the dominant contribution for the averaging in eq. 9 will come from superpixels containing small values of fluctuations.Hence in the region q < 0, h(q) effectively describes the scaling behavior of small fluctuations.
• Furthermore, given some data having multifractal properties, the value of h(q) is usually larger for q < 0, compared to q > 0. This can be explained by the fact that the contributions to the averaging in eq. 9 from small fluctuations typically vary more across the scale s, compared to large fluctuations.
We next relate h(q) to other measures of complexity used in the literature.From the standard partition function-based multifractal formalism, the classical multifractal scaling exponent τ is related to h(q) as [14] where D f is the fractal dimension of the geometric support of the object.We take D f = 2 for a two-dimensional (2D) image [16].If h is independent of q, then τ is a linear function of q.Therefore, a nonlinear shape of τ (q) is an indication of the multifractality of the system.
The Hölder exponent (also known as singularity strength), denoted by α, is calculated by taking the Legendre transformation of eq. ( 12) as where the prime denotes derivative with respect to q.We get α is a measure of the shape of h(q) since it is a function of h ′ (q).Different values of α characterise different parts of the system roughly since small, intermediate or large values of the fluctuation field will contribute to different values of α.From eq. 14 we can see that α can take both positive and negative real values depending on the sign of the second term and its relative value with respect to the first term.Note that for a monofractal system, we have h ′ (q) = 0 since h is a constant, and as a consequence, α has only one value given by α = h.Therefore, a variation of α with q signifies multifractality, and the range over which it varies quantifies the strength of multifractality.
The last quantity we use is the multifractal spectrum (also called the singularity spectrum) f (α), defined by This quantity is useful in extracting information about the symmetry between small and large fluctuations.
In order to gain an intuitive understanding of the quantities we have defined so far, let us consider two toy examples of multifractal complex systems whose h(q) are given below.
Example A. Let h(q) be given by Example A Example B Example B. Let h(q) be given by These examples have been chosen so as to capture the essence of our results in section IV. Figure 1 shows h, h ′ , τ and α versus q, and f versus α for the two examples.
We have chosen D f = 2 for both examples.If h(q) is a well-behaved monotonic function of q, such as in example A, then all the other quantities are single-valued.However, if h(q) is not monotonic, such as in example B, then h ′ and α become multi-valued.As a consequence, the multifractal spectrum f (α) is not a well-behaved function.
In case, f (α) has a well-defined peak located at α = α 0 , then we can quantify its symmetry about α 0 .Let α min be the smallest value of α, and α max denote the value of α where f (α max ) = f (α min ).Then, let χ denote the skewness parameter, given by The and symmetrical if χ = 1.If χ < 1, then it shows that the scaling of the system is dominated by large fluctuations (smaller generalized Hurst exponents).However, if χ > 1, then it shows the dominance of scaling by small fluctuations (higher generalized Hurst exponents).In this section, we present a brief description of the physical processes that are relevant to understanding the Crab Nebula.Then we give details of the imaging data that we use for our analysis.

A. Physics of the Crab Nebula
The Crab Nebula belongs to the class of supernova remnants known as pulsar wind nebulae (PWN).Such a nebula inherits its morphology and spectrum from the following factors.
1.During the supernova explosion of the progenitor star, the stellar material gets ejected (known as ejecta).The ejecta sweeps up the interstellar medium, creates a heated shock front, slows down, and eventually freely expands into the interstellar regions beyond the confines of the nebula.
2. The collapse of the core of the progenitor star results in a rotating magnetized neutron star, known as the Crab pulsar.The magnetic field causes the charged ejecta particles to accelerate and emit synchrotron radiation.
3. The Crab pulsar spins at the rate of 33 ms [51].
The pulsar emits non-thermal plasma, known as pulsar wind, which gets accelerated to high velocities by the magnetic field.The plasma pushes into the surrounding ejecta in the nebula and undergoes turbulent mixing.This leads to Rayleigh-Taylor instabilities [52] that give shape to the Crab Nebula's filamentary network and wispy structures.
The Crab Nebula appears like a bubble whose shell expands radially over time as the ejected material moves outward, while the interior structure consists of an expanding network of filaments and wispy structures when viewed across a wide range of frequencies.Freely expanding ejecta is also expected to be present beyond the visible Crab.Here we focus on radio observation of the Crab Nebula.The origin of radio emission is described below.
Radio frequencies: The emission in radio frequencies is primarily the non-thermal synchrotron radiation from the shocked pulsar wind; charged particles are accelerated due to the magnetic field of the pulsar.Visually, one can see the expanding network of filaments and finescale wispy features (Figure 2).The source of most of the thermal emission from the filaments is understood to be the photo-ionization by the hard continuum from the synchrotron PWN emission [see, e.g., 53].In addition, radiative cooling behind the shock driven by the PWN into the freely expanding ejecta also gives rise to radio emission [54].

B. Observational data of the Crab Nebula
We focus on the data of the Crab Nebula observed at radio frequency by the Very Large Array (VLA) [1].We use publicly available radio imaging data of the Crab nebula at 4.76 GHz with an angular resolution of 0.7 arcseconds from the VLA image archive [2].These data have undergone full calibration and reduction.
The image we analyze is shown in Figure 2. The surface brightness in the image is given in units of Jy beam −1 .The data has undergone full calibration and reduction.Hence we can directly apply the 2D-MFDFA algorithm to this imaging data and analyze its structural complexity.

IV. ANALYSIS AND RESULTS
We apply the 2D-MFDFA algorithm to the radio imaging data of the Crab Nebula and study its multifractal properties.Using the intensity or field value at each pixel, we first calculate the fluctuation functions F q defined by eqs. 9 and 10, and from them, we obtain the generalized Hurst exponents h(q).Then, from h(q), we obtain the derived quantities -the classical scaling exponents τ (q), [1] VLA is operated by the National Radio Astronomy Observatory (NRAO), USA.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under a cooperative agreement by Associated Universities, Inc. [2] https://archive.nrao.edu/archive/archiveimage.html the Hölder exponents α, the multifractal spectra f (α) and the skewness parameter χ.
The multifractal analysis using the 2D-MFDFA algorithm was performed using the MATLAB R2018a platform.The 2D-MFDFA algorithm calculates h(q) by the method of linear fitting using the inbuilt regression statistics regstats of MATLAB.In particular, it uses the inbuilt routine tstat based on Student's t-distribution.The code also calculates the error bars for h as the standard errors from tstat.
In order to study the spatial scaling behaviors of the multifractal variables over different spatial regions of the Crab Nebula, we repeat the 2D-MFDFA analysis on smaller partitions of the image for both radio and IR data.This is carried out by dividing the images into N a ≡ 2 a equal-area portions, where a ≥ 0 is an even integer that determines the number of partitions.We refer to a as the partition scale parameter (so as not to confuse with the size scale s).We use values a = 0, 2, 4 since the number of pixels in each partition becomes too small beyond that.The spatial scaling behavior will reveal whether the same morphological rules are followed or not at different parts of the system.It will also enable us to identify the presence of heterogeneous properties at different length scales.
The regions on the top left and bottom right of the image in Figure 2 are predominantly noise.Since our calculations require a uniform size of spatial regions, we retain these regions and perform all multifractal analyses.We present the corresponding results in Section IV A and IV B. In Section IV C, we analyze the effect of the noise on the multifractal parameters by adopting a masking procedure (described therein the section) on the radio image of Figure 2. We then compare the results for the original (unmasked) and masked radio data.

A. Complexity of Crab nebula at the radio frequency of 4.76 GHz
The Crab Nebula image at the radio frequency of 4.76 GHz (shown in Fig. 6) consists of 650 × 740 number of pixels.
Figure 3 shows the log-log plots of the fluctuation functions F q as functions of the scale size s, for different partition scales a = 0, 2, 4 (indicated against each panel) for different values of the order of statistical moments q (shown by different colors).For all the partition scales a, it is observed that F q has power law dependence on s and the slopes of F q vary for different q.This is a signature of the presence of multiple local fractal behaviors at different local domains s, across all a.The spread of the curves along F q , i.e. ∆F q for small values of s is much larger than the spread in F q at large values of s.This shows that the role of fluctuations is more evident in FIG.3: Radio data of Crab Nebula at 4.76 GHz -The fluctuation function F q is plotted versus scale size s, for different values of q, and different partitions of the radio image corresponding to a = 0, 2, 4.
smaller local patterns/clusters of the system, which provide distinct local behaviors.However, at large values of s, this spread ∆F q → 0, where the system converges to a single behavior.
The slopes of the log-log F q -s plots in Figure 3 give the values of the generalized Hurst exponent h as a function of q, as defined by eq.(11).The left column of Figure 4 shows h versus q, ranging from -40 to 40 in steps of one.The cases of a = 0, 2, 4 are shown in the panels from top to bottom.The number of plots (shown in different colors) in each panel corresponds to the number of partitions of the radio image for each a.The vertical dashed lines in the panels showing h(q) and τ (q) correspond to q ∼ 0. We observe significant dependence of h on q for all a, indicating distinct local behaviors at various partition scales a.This shows that small and large fluctuations of the emission at the radio frequency have very different scaling natures.
As described in section II B 2, for positive (negative) values of q, h describes the scaling behavior of the segments with large (small) fluctuations.From the left column of Figure 4, we find that the values of h(q) are larger for q < 0 than that of q > 0 for all partition scales a.Hence, local behaviors of the patterns/clusters at various a are found to be sensitive to small fluctuations.These small fluctuations drive the local behaviors significantly different compared to the large fluctuations in the radio data [14].The significant spread in h as a function of q, i.e. ∆h at each a shows that scaling laws of the patterns/clusters at various a are distinctly different, which shows a heterogenous distribution of the patterns/clusters in various partitions of the system.Further, since the values of h(q) for q < 0 are larger than that of q > 0, the correlated memory in the patterns/clusters/system driven by small fluctuations persists longer as compared to that driven by large fluctuations [14,50].
Visually, in the radio image of the Crab Nebula (see Figure 2) we can discern that in general large fluctuation values track the structurally dominant network of filaments and are correlated over large length scales.In contrast, the smaller fluctuation values track the small-scale wispy structures and are correlated over shorter length scales.Hence, our finding that h(q > 0) < h(q < 0) is not surprising.It encapsulates the long-range correlations of the filamentary network and the shorter-range correlations, as well as the stronger scale dependence of the wispy structures.Further, we observe scale-dependent heterogeneous patterns in the various partitions, which show different local multifractal scaling laws of the patterns.
In all the plots of h(q), we see that it asymptotes to constant values as |q| → ∞.This implies that small fluctuation values tend towards the same scaling behavior, and so do large fluctuations.Let h ± denote the asymptotic values at q → ±∞.We observe that h − varies considerably large across the different partitions for a ≥ 2, while h + shows relatively less variations such that h ± ≈ h ± (a, m), where m is the partition index.This means that the structural properties associated with small fluctuations vary considerably large amongst the different partitions of the radio image, while large fluctuations are comparable.Further, since the values of h − > h + , the persistence of the patterns/clusters in each partition driven by small fluctuations has longer memory as compared to large fluctuations.We thus expect to see more heterogeneous or dissimilar structures with distinct local scaling laws across the different partitions of the radio image significantly contributed by small fluctuations.This corroborates what we described in the previous paragraph.We also find the largest variation of h − across the different partitions for a = 4.A closer look at the values of h ± for different partitions of the radio image will be discussed in section IV C.
In the middle column of Figure 4, we plot the classical scaling exponent τ versus q for different values of a.In all the plots, we observe that τ does not depend on q linearly.Eq. ( 12) tells us that if h is independent of q, then τ should be a linear function of q.But h is dependent on q, and it exhibits a step-like transition around q ∼ 0 (see the left column of Figure 4).Therefore, we expect that τ will also make a transition around q ∼ 0. This is indeed what we see in the plots of τ , where we observe two linear regimes with a slope change around q ∼ 0 (indicated by the vertical dashed lines).This value marks the transition of the scaling nature between large and small fluctuations exhibiting multifractal properties.
Lastly, in the right column of Figure 4, we plot the multifractal spectra f (α) versus the singularity strength α for different partition scales a.Since h(q) is a monotonous function, α is single-valued.We observe that all the curves have their maxima, denoted by f max at roughly α = α 0 ∼ 2. Also, visually it is clear that all f curves are right-skewed around α 0 , showing a high degree of complexity in the structure of each partition defined by a [55].Further, the widths of α for a fixed value of f (α) of all the partitions of each a are significantly different, which shows heterogeneous structures in the partitions with different scaling behaviors [56].The values of the skewness parameters χ for each f are calculated using eq.( 18).We will discuss χ in detail in section IV C.

B. Origin of the multifractal nature
We now probe the physical origin of the multifractal nature that is obtained for the radio data.Multifractality can arise due to two different physical reasons-(1) distribution-induced multifractality caused by a heavytailed probability density function (PDF) of the data, and (2) correlation-induced multifractality caused by differing nature of correlations for small and large fluctuations, such as linear and nonlinear correlations [9,14].One can analyze the origin of multifractality as follows.Shuffling the data destroys all the long-range correlations while the PDF remains unchanged.On the other hand, multifractality due to a heavy-tailed PDF cannot be removed by shuffling [57].If the multifractality is solely due to correlations, the shuffled data will show monofractality [57].If the data contains both heavy-tailed PDF and linear/nonlinear correlations, then the shuffled data will show multifractality smaller than that of the original data [58].First, in order to isolate the effect of correlation, we randomly shuffle the data to destroy all spatial correlations.We will use superscripts 'orig' and 'shuf' to indicate quantities obtained from the original and shuffled data, respectively.Suppose F orig q and F shuf q are the fluctuation functions corresponding to the original and shuffled data with their respective generalized Hurst ex-ponents h orig and h shuf .Then the ratio where h corr = h orig −h shuf .We can analyse the origin of multifractality in the data from h corr as: (i) if h corr = 0 such that h orig = h shuf , then the cause of multifractality is the fatness of the PDF, and (ii) if h corr = h corr (q) ̸ = 0 and h shuf = h shuf (q), then the origin of multifractality is due to both q-dependent long-range correlations coming from short and long-range fluctuations as well as a broad PDF in the data, and (iii) if h shuf = 2, then the multifractality in the data is due to only long-range correlations [14,59].
In the leftmost panel of Figure 5, we show the PDF for the radio data.We find that it is close to Gaussian.We calculate the multifractal quantities h(q), τ (q) and f (α) for the shuffled data.The second from the left panel of Figure 5 shows h(q) for the two cases of original and shuffled data.We indeed find h shuf (q) has weak qdependence with h shuf (q) → constant → 2 for the shuffled data.Hence, h corr → h orig (q) − 2. Similarly, we plot τ (q) and f (α) for the two cases of original and shuffled data, respectively, in the third from left and right most panels of Figure 5.We previously found that the scaling exponents τ (q) at all the considered partition scales a = 0, 2, 4 are nonlinear functions of q (see middle panels of Figure 4).This non-linearity shows the presence of intrinsic multifractality [60].When we shuffled the original field values and repeated the same multifractal analysis, we found that the τ (q) curve now becomes linear (second from the right panel in Figure 5).This shows that the shuffling process destroys all the intrinsic non-linear correlations.
Since the PDF of the radio data is close to Gaussian, we simulated an uncorrelated Gaussian random field to have a handle on what to expect of the multifractal parameters for uncorrelated Gaussian data (see Appendix A).By repeating the same analysis we performed with Crab Nebula data, we have calculated h(q), τ (q) and f (α) for one realization of an uncorrelated Gaussian random field.We find that for this uncorrelated field, τ (q) curves show monofractal behavior.From this result in appendix A, we expect the shuffled radio data to behave close to a monofractal and hence to find h(q) to be almost constant with value 2. This is in agreement with what we observe in the panels of column 2 in Figure 5, where we find τ (q) curves to be linear thereby indicating monofractal behavior.We show the breaking of intrinsic non-linear correlations responsible for intrinsic multifractality by the shuffling process.Therefore we conclude that the multifractality of the Crab nebula radio data originates from long-range non-linear correlations of large and small fluctuations in the data.FIG.4: Multifractal parameters for the radio image of the Crab nebula.Column 1: Generalized Hurst exponents h versus q for a = 0, 2, 4. Column 2: The classical scaling exponents τ versus q for a = 0, 2, 4. Column 3: Multifractal spectra f versus the singular exponents α for a = 0, 2, 4. Different colours represent N a = 2 a partitions for each value of a.The vertical dashed lines in the panels for h and τ represent the value q ∼ 0, where the slopes make a transition.
FIG. 5: Left: PDF of the radio data.The x-axis is in units of Jy beam −1 .The other three panels show h(q) (second from left), τ (q) (third from left) and f (α) (right) of the radio data for the original and shuffled cases.The black dashed horizontal line shows h = 2 for reference.

C. Heterogeneity in the spatial structures
We now study the heterogeneity in the structures of the Crab Nebula observed at the radio frequency of 4.76GHz by examining integrated quantities obtained from h, τ and χ as functions of the scale parameter.
To address the effect of noise present in the top left and bottom right regions of the image in Figure 2, we repeat the calculations of the quantities h, τ and χ after masking these regions.We construct an elliptical mask of appropriate dimensions [82] and assign zero values to the regions outside the ellipse while retaining the original field values inside the ellipse (see Figure 6).The partitions in the image are numbered for the cases of a = 2 (orange) and a = 4 (yellow).In the latter case, the partitions are numbered in a nested manner so that se-quential groups of four are mapped hierarchically to the corresponding larger partition of a = 2.All calculations are then repeated for the masked image.
We first examine the asymptotic values h± for the different partitions indexed by i a (for a given partition scale a).The plots of h± versus i a are shown in Figure 7 for a = 2 (left panel) and a = 4 (right panel).Red stars indicate values for the unmasked original image, while blue indicates values for the masked images.For a = 2, there is an increase in h for i a = 1 and a slight increase in h for i a = 3 after masking.Similarly, for a = 4, there is a change in h for i a = 1, 2, 4 and a slight change in i a = 10, 11, 12 after masking.These changes are because the mask artificially introduces two distinct regions in each affected region -one that belongs to the Crab nebula and another that has predominantly zero values.
From the bottom panels of Figure 7, we can see that for a = 2, the radio image (both original and masked) has larger values of h − and smaller values of h + .This captures the wider variation of scaling behavior between large and small fluctuation values of the radio data, already discussed earlier.This is corroborated by visual inspection of the image.If we closely look at the smaller partition scale structures (a = 2) of the Crab nebula observed at radio (top panel of Figure 7), we can easily see that the radio image has more heterogenous structures across the partitions in the form of filaments and wisps.Moreover, the filamentary network of the radio image has longer-range correlations.This can be explained by the fact that the materials that emit synchrotron radiation at radio frequencies, which are the charged particles of the supernova ejecta and pulsar wind, are distributed over larger spatial scales.At a = 4, h − for the radio data varies randomly.This implies that the structures across the different parts of the Crab Nebula are highly heterogeneous when observed at radio frequency.These highly heterogeneous structures at various length scales make the Crab Nebula a structurally complex physical system.
Next, we focus on the values of h in the regime around q = 0, where it exhibits a transition.Let the average of h(q) over q in a suitable regime −q c < q < q c , with q c is a positive variable, be given by h ≡ In the limit q c → ∞, h will be dominated by the asymptotic values, and we simply get h ∼ (h − + h + )/2.In this limit, we lose the information of the variation or shape of h(q) in the intermediate q values.So we focus on finite q c and choose q c = 40 as the value of q, where h becomes effectively constant for each partition scale a.This is akin to weighted averaging, where the weight kernel is the top hat function centered at q = 0. Mathematically there is no restriction on taking large |q|.Since q takes integer values here, we need the discretized version of eq.19.Let i q = 1, . . ., N q index the values of q, where N q = 81 denotes the total number of q values considered.Let i a = 1, . . ., N a index each partition of the image for a given a.Then for each i a , the discrete form of eq.19 is For τ , there are no asymptotically constant values.Hence, there is no natural cut-off value of q c to define the range −q c < q < q c over which to average τ .We will choose the same q c that is defined for h.So we have τia Let χ ia denote the skewness parameter χ in each partition i a .Let Xia denote either of the three quantities hia , τia or χ ia .Then we average them over all the partitions for each a, as Xia .
Figure 8 shows ⟨h⟩ (top panel), ⟨τ ⟩ (middle panel) and ⟨χ⟩ (bottom panel) for the radio image as functions of partition scales a.The error bars shown are the standard errors.We observe that ⟨h⟩ varies with a for the radio data.This result correlates with the value of (h − +h + )/2 seen in Figure 7, with the larger contributions coming from h − .This implies that the radio image has stronger small-scale scaling behaviors, showing the presence of fine-scale structures with well-defined local scaling laws in the radio data.
Further, let h 0 denote the largest value of ⟨h⟩ chosen from the different a's.Let us normalize ⟨h⟩ by h 0 and denote If ⟨h⟩ norm ∈ (0.5, 1), then this implies a strong positive correlation in the system [14], with long memory of persistence in the states of the system [10,50].From the top panel of Figure 8, for the radio data, h 0 = 2.28, which is the value of ⟨h⟩ at a = 0. Using this, we find that ⟨h⟩ norm lies in the range [0.92, 1] for all a for the radio data.This value is greater than 0.5 and close to unity.This implies strong positive correlations [14] in the Crab Nebula structural system.Further, the decrease in ⟨h⟩ as a increases show that finer structures in the smaller partitions of the system lose their correlations and local structural memories, which is expected.In the asymptotic limit lim where γ is the power-law exponent.This shows a hierarchical organization of the heterogenous structures in the Crab Nebula structural system [61][62][63].This behavior is similar to Carpenter's law as predicted in astrophysical data analyses [63][64][65].
The middle panel of Figure 8 shows that the average classical scaling exponent ⟨τ ⟩ also varies with a, which again implies heterogeneity in the structure of the underlying system.From the bottom panel of Figure 8, ⟨χ⟩ > 2 for all a, which implies the richness of complexity of the Crab nebula at radio frequency.This shows the dominance of scaling of small fluctuations or the presence of fine structures in the Crab Nebula.The random distribution of the ⟨χ⟩ parameter also implies a heterogeneous property of the structure visible in the Crab Nebula.This can be accounted for by the non-equilibrium dynamics of the underlying system.

V. CONCLUSIONS
The Crab Nebula is one of the most studied astrophysical systems across the entire electromagnetic spec-trum from radio wavelengths to ultrahigh-energy γ-rays [66,67].In astrophysics, several recent studies confirm the complexity of the Crab Nebula's morphology across the different bands of the spectrum [68][69][70][71][72][73].In this paper, we offer a fresh perspective on the complex morphology of Crab Nebula using multifractal analysis from the vantage point of complex systems science.
We have investigated the structural complexity of the Crab Nebula at the radio frequency of 4.76 GHz using publicly available imaging data.In particular, we have analyzed the data using the 2D-MFDFA approach.We have investigated the structural system at various length scales for local scaling behaviors.Our results confirm that the structure of the Crab Nebula has a rich complexity characterized by a multifractal nature.We have studied the origin of multifractality in the radio data.We have found that the multifractal property of the radio data arises from the difference in the nature of longrange correlations (and consequent scaling with the scale parameter s) of large and small field values in the data.
Our results show a strong positive correlation in the radio frequency.The long-range correlations of large and small values in the data can be traced back to the materials that emit at this frequency.The variations of the multifractal parameters, namely the generalized Hurst exponents h calculated from the fluctuation functions F q , the classical scaling exponents τ and the skewness parameter χ across varying length scales show heterogeneity in the structure across different parts of the visible Crab Nebula.This can be accounted to the Crab Nebula being a non-equilibrium, highly dynamic structural system.
What we have done in this paper is essentially quantified how the physics of the evolution of the Crab Nebula as a supernova remnant manifests in the language of complex systems.The 2D-MFDFA approach we have used in the study is found to be a good method to characterize the complexity of an astrophysical system with complex morphological features.Our multifractal results inform us that the rich fine structures in the Crab Nebula ruled by local scaling behaviors show self-organization at different length scales of the structural system.Longrange correlations manifest as self-similar fractals and are identified as signatures of self-organized criticality [74,75].Our results thus deepen our understanding of the multi-scale physics that is operating in the Crab Nebula.An example of a complex system with multifractal characteristics similar to the Crab Nebula is the human brain.Multifractality captures the heterogeneous and multiscale interaction rules in the brain networks [76]; for instance, the community structures at different structural levels in the human brain connectome [77] ruled by scaling behaviors that are signatures of self-organization [78].This paper represents a proof of concept of the insights that can be obtained from viewing the Crab Nebula as a complex system.However, it has some limitations.First, our analysis considers the Crab as a two-dimensional object and ignores its three-dimensional structure.It will be interesting to study the full 3D structure using reconstructions such as those given by [79].Secondly, the Crab Nebula is a rapidly evolving system on astrophysical time scales.The data that we analyze here are observations taken at specific times, and we do not address the question of the time evolution of the properties related to complexity.In order to extract the full potential of the method, our work can be extended in the following directions.Given the wealth of data available for this astrophysical object, an analysis across the full range of frequencies will be valuable and will be carried out in the future.It opens up a slew of questions that will enhance our understanding of supernova remnants, and we plan to carry out a systematic investigation in the near future.
The methodology used in this paper can also be used to understand other astrophysical complex systems such as the planetary nebula NGC 1514, whose complex kinematics with multiple structures was studied [80].As pointed out in the review paper [81], pulsar wind nebulae, f which the Crab nebula is a prototype, show a variety of properties and morphologies at different evolutionary ages.It will be interesting to extend the present analysis to study the multifractal properties of a wide variety of pulsar wind nebulae and look for their universal properties.knowledges funding from the Science and Engineering Research Board of the Department of Science and Technology, Government of India, under the MATRICS scheme, bearing project reference no MTR/2021/000766.PK acknowledges funding from the Department of Atomic Energy, Government of India, under the project 12-R&D-TFR-5.02-0700.