The “Singular” Behavior of the Solar Wind Scaling Features during Parker Solar Probe–BepiColombo Radial Alignment

At the end of 2020 September, the Parker Solar Probe (PSP) and BepiColombo were radially aligned: PSP was orbiting near 0.17 au and BepiColombo near 0.6 au. This geometry is of particular interest for investigating the evolution of solar wind properties at different heliocentric distances by observing the same solar wind plasma parcels. In this work, we use the magnetic field observations from both spacecraft to characterize both the topology of the magnetic field at different heliocentric distances (scalings, high-order statistics, and multifractal features) and its evolution when moving from near-Sun to far-Sun locations. We observe a breakdown of the statistical self-similar nature of the solar wind plasma with an increase in the efficiency of the nonlinear energy cascade mechanism when moving away from the Sun. We find a complex organization of large field gradients to dissipate the excess of kinetic energy across the inertial range near the Sun, whereas the topological organization of small fluctuations is still primarily responsible for the energy transfer rate at 0.6 au. These results provide, for the first time, evidence of the different roles of dissipation mechanisms near and far away from the Sun.


Introduction
The heliosphere is permeated by the continuous flow of an ionized medium known as the solar wind, which changes its dynamical, topological, and physical properties after being ejected from the Sun (Bruno & Carbone 2016). The solar wind is a typical example of a multiscale complex system, being indeed characterized by several processes and phenomena occurring on a wide range of scales, ranging from kinetic to large-scale ones (Verscharen et al. 2019). Within this large variety of scales, a lot of attention has been paid to the so-called inertial range, corresponding to the range of scales between the injection/integral scale L, at which energy is injected into the system, and the dissipative scale η, at which energy is dissipated (Bruno & Carbone 2016). This range can be described by means of the macroscopic theory of magnetohydrodynamics (MHD), in which the solar wind is described as a single fluid, assuming that spatial scales are larger than all inherent scales, such as the Debye length or the ion gyroradius, and that the inherent frequencies are lower than the ion gyrofrequency (Biskamp 2003). A striking feature of the inertial range is its scale-invariant nature, usually associated with the idea that the energy transfer mechanism across these scales can be described via the concept of a nonlinear cascade driven by nonlinear interactions between scales (e.g., Kolmogorov 1941;Carbone 1993;Frisch 1995;Biskamp 2003), i.e., the solar wind is a turbulent medium (Bruno & Carbone 2016).
During the past 50 yr, several attempts have been made to investigate and characterize the inertial range dynamics, by means of different space missions, allowing us to explore different regions of the heliosphere, from the near-Earth (e.g., Cluster, the Magnetospheric Multiscale Mission) to the near-Sun (e.g., the Parker Solar Probe, Helios), passing through the Lagrangian point L1 (e.g., Wind, Advanced Composition Explorer), and also reaching the boundaries of the heliosphere (e.g., Voyager). In particular, the Helios mission did pioneering work in the field of solar wind turbulence in the inner heliosphere (Tu et al. 1989;Marsch & Tu 1990a, 1990bHe et al. 2013). The above missions allowed us to determine that solar wind turbulence is characterized by an anisotropic cascade (Horbury et al. 2008;Chen et al. 2010;Narita 2018), whose scaling depends on different mechanisms, such as the large-scale forcing (Velli 2003), the Alfvénicity and the compressibility (Matteini et al. 2018), and so on (Matthaeus et al. 2012;Boldyrev et al. 2015). When approaching the Sun, a transition region from the steeper scaling to the shallower ones is observed near 0.4 au (Alberti et al. 2020;Chen et al. 2020), suggesting the existence of some control parameters that determine the topology of solar wind fluctuations like cross-helicity (Matthaeus et al. 1982;Tu & Marsch 1990) and/or the role of inward versus outward perturbations (Dobrowolny et al. 1980), also marked at larger radial distances (Wawrzaszek et al. 2015(Wawrzaszek et al. , 2019. For the above reasons, there is now increasing interest in investigating the radial dependence of solar wind turbulent features, especially thanks to the large numbers of spacecraft that are monitoring the solar wind at different heliospheric locations. This will significantly increase the possibility of multispacecraft observations of the interplanetary medium variability, facilitating the exploration of temporal, spatial, and radial variations of the solar wind and transients (Telloni et al. 2021). Thus, in the near future, we could advance our understanding of solar wind dynamics and structures via multispacecraft analysis.
In this manuscript, we present the results of a multispacecraft analysis in investigating the evolution of solar wind turbulence at different heliocentric distances by observing the same solar wind plasma parcels. By using magnetic field observations from the Parker Solar Probe (PSP) and the BepiColombo spacecraft (Section 3), we characterize the underlying multifractal nature of magnetic field fluctuations when observing the same plasma parcels at different heliocentric distances. After a few notes on the theoretical background of the multifractal approach to solar wind turbulence (Section 2), we introduce a recently developed method to investigate the full multifractal spectrum (i.e., for both positive and negative statistical moments q), based on Empirical Mode Decomposition (Section 4). In Section 5, we present and discuss the main results: we show the existence of a breakdown of the statistical self-similar nature of the solar wind plasma with an increase in the efficiency of the nonlinear energy cascade mechanism when moving away from the Sun. In Section 6, we conclude that the results of this study support previous evidence of the radial dependence of solar wind scaling behavior, and suggest that it can open up a novel framework for modeling magnetic field topological changes across the heliosphere.

Some Notes on Turbulence and Multifractals
The theory of MHD turbulence (e.g., Iroshnikov 1964;Kraichnan 1965) is strictly connected to the Kolmogorov theory of turbulence (Kolmogorov 1941), which is based on the assumption that small-scale statistics obey universal scaling properties, i.e., there exists a range of scales whose statistical properties are only determined by the average energy transfer rate ò. This is related to the scale-invariant nature of the MHD equations, meaning that they are formally invariant with respect to scaling transformations (Carbone 1993;Bruno & Carbone 2016). Assuming that the energy transfer rate is preserved along the cascade, the Navier-Stokes and the MHD equations possess only one scaling exponent, h = 1/3 and h = 1/4, respectively (Kolmogorov 1941;Iroshnikov 1964;Kraichnan 1965;Carbone 1993;Frisch 1995). Since the 1960s, experiments on fluids, and later, in the 1980s, in situ measurements of solar wind, have pointed toward the absence of the existence of a unique scaling exponent. This absence relates to the singular character of Navier-Stokes and MHD equations in the limit of high Reynolds numbers (Frisch 1991;Carbone 1993;Bruno & Carbone 2016;Dubrulle 2019), reflecting the intermittent character of the cascade mechanism moving energy across scales within the inertial range (Benzi et al. 1984;Carbone 1993). This led to the development of the "multifractal formalism" to take into account the energy transfer rate not being preserved along scales, but instead fluctuating both in space and time (Mandelbrot 1974(Mandelbrot , 1982Frisch 1995). With this formalism being P(ℓ) ∼ ℓ 3− D( h) , the probability of observing a rescaling exponent h for the turbulent fluctuation D  z ℓ at the scale ℓ, the qth-order structure function must be written as being a nonlinear convex function of q (Benzi et al. 1984;Frisch 1995), with τ q/m being taken into account for intermittency correction of the energy transfer rate at different scales (Frisch 1991;Carbone 1993). In this way, all properties of the cascade are encoded in the multifractal spectrum D(h), while all information on the energy transfer mechanisms is encoded into the rescaling exponent (or singularity) h. This formalism allows us to overcome the two fundamental failures of the Kolmogorov theory of turbulence: (i) its believed universal behavior (i.e., that the energy transfer rate depends on time and space); and (ii) its assumption of considering only one possible rescaling symmetry (i.e., that associated with only one value of h = 1/m; Dubrulle 2019).
The multifractal formalism thus provides a statistical and probabilistic theory of turbulence, allowing us to gain further information on the energy transfer mechanism by looking at the singular behavior of turbulent fields. In particular, the rescaling exponent associated with the maximum value of D(h) represents the most probable scaling exponent: when q < 0, the topology of small fluctuations is emphasized, being the reflection of a smooth fractal structure of the turbulent field, while when q > 0, we explore the increasingly singular points, emphasizing the topology of extreme events with large gradients in the field fluctuations (Benzi et al. 1984;Frisch 1995).

Data
Since we are interested in investigating the radial evolution of solar wind turbulence, we used measurements provided by PSP and the BepiColombo spacecraft when they were radially aligned and orbiting near 0.17 and 0.6 au, respectively. Since we cannot evaluate the Elsässer variables, due to the lack of solar wind speed measurements, we used magnetic field measurements provided by the magnetometers on board both spacecraft. The PSP magnetic field observations are taken from the FIELDS Fluxgate Magnetometer (MAG; Bale et al. 2016) and are averaged to a 1 s cadence from their native four samples per cycle cadence (see also Alberti et al. 2020;Chen et al. 2020). The data were freely retrieved from the Space Physics Data Facility Coordinated Data Analysis Web interface. 6 The BepiColombo magnetic field data were obtained from the MPO-MAG team (Glassmeier et al. 2010;Heyner et al. 2021).
For this study, we use solar wind magnetic field components in the heliocentric Radial Tangential Normal (RTN) spacecraftcentered coordinate system, where R is the Sun-spacecraft direction pointing toward the Sun, T is the tangential axis to the orbit of the spacecraft, and N completes the right-handed triad (Fränz & Harper 2002).

Multifractal Approach
During the past 40 yr, several methods have been developed to detect singularities by introducing a suitable invariant measure and by searching for scaling-law behavior of this quantity (e.g., Chhabra et al. 1989;Consolini et al. 2021), via a Legendre transformation from the partition function scaling exponents (e.g., Paladin & Vulpiani 1987), or by moving in a conjugate space, like the wavelet transform modulus maxima (e.g., Muzy et al. 1991). Here, we used a recently proposed method, based on Empirical Mode Decomposition (EMD), to provide a new perspective in searching for multifractal features (e.g., Welter & Esquef 2013).

Empirical Mode Decomposition (EMD)
We first introduce some key relevant features of EMD, since it is the starting point of our multifractal formalism. EMD is an adaptive decomposition method that allows us to preserve the nonlinear properties of time series and to deal with their nonstationary features (Huang et al. 1998). By means of the socalled sifting process (e.g., Huang et al. 1998, for more details), a time-dependent signal s(t) can be written as Each γ j (t), named as an Intrinsic Mode Function or empirical mode, is a component of the decomposition basis, and it must satisfy two properties: (i) have the same (or differing at most by one) number of extrema and zero crossings; and (ii) the average of the upper envelope, derived via cubic spline interpolation of local maxima, must be the same as the absolute value of the average of the lower envelope, derived via cubic spline interpolation of local minima (Huang et al. 1998). The set of γ j (t) forms the decomposition basis, whose orthogonality can be checked a posteriori once the decomposition has been completed, while the procedure satisfies the completeness and convergence properties by means of the nonoscillatory residue of the decomposition ρ(t), i.e., a monotonic function of time (Huang et al. 1998). By means of the Hilbert Transform, we can write each empirical mode as modulated both in amplitude a j (t) and in phase f j (t) as thus allowing us to write each component as an oscillating function of time on a typical timescale τ j , expressed as For further details, and for a more accurate description of the method, the reader is referred to Huang et al. (1998) and Huang & Wu (2008).

The EMD-based Multifractal Formalism (EMD-DAMF)
In a similar fashion to the wavelet transform modulus maxima (e.g., Muzy et al. 1991, for more details), an EMDbased method for detecting multifractal features of time series can be built up (Welter & Esquef 2013). It allows us to evaluate singularities and fractal properties by means of the following steps.
1. For each empirical mode γ j (t), derive the instantaneous amplitude a j (t) and mean timescale τ j ; 2. For each local maximum of a j (t), define a support  j ℓ , , i.e., a closed set of points around each local maximum ℓ; 3. Define the dominant amplitude coefficients 4. Define the qth-order structure function S q (τ j ): 6. Estimate the Hölder exponents h and the multifractal spectrum D(h) by means of a Legendre transform of ζ(q) This formalism allows us to compute structure functions based on the local features of the fluctuations at different scales, as described via the empirical modes. Indeed, instead of defining increments at different a priori selected scales, the procedure allows us to exploit the local properties of empirical modes to evaluate the differences/increments as the differences between each maximum of the instantaneous amplitudes derived via the EMD and the two adjacent ones (e.g., Figure 1 in Welter & Esquef 2013). In this way, increments are not evaluated over a constant scale, as for canonical structure functions, but over a scale that is defined as the difference between two consecutive maxima. Since one of the properties of the empirical mode is to have a time-dependent frequency, the difference between two consecutive maxima is not constant, but varies around the mean timescale τ j of each empirical mode. However, the derived structure functions are equivalent to the usual definitions given via increments, unless being derived in an adaptive way (Welter & Esquef 2013). Thus, by means of the pair (h, D(h)), we can introduce some fundamental measures of complexity, allowing us to characterize the underlying fractal structure of the time series, the level of intermittency, and the role of singularities ). These measures are (Macek & Wawrzaszek 2009;Macek et al. 2011;Nguyen et  For monofractals (Δh, ΔD) = (0, 0), for symmetric fractals A = 1, and for regular (singular) fractals A < 1(A > 1), according to the Kolmogorov theory of hydrodynamic turbulence, (h, D(h), A) = (1/3, 1, 1) (Kolmogorov 1941); while according to the Iroshnikov-Kraichnan picture of MHD turbulence, (h, D(h), A) = (1/4, 1, 1) (Iroshnikov 1964;Kraichnan 1965).

Results
To find whether a radial alignment between PSP and BepiColombo could occur, we first evaluate what would be the solar wind speed measured by PSP such that it would hit BepiColombo.
As shown in the sketch presented in Figure 1, PSP and BepiColombo would observe the same plasma parcel if the travel time of the plasma to cover the spatial distance Δr were the same as the orbit time for BepiColombo to cover the spatial distance ℓ. Thus where V SW represents the solar wind speed measured by PSP and V Bepi the orbital velocity of BepiColombo. Since we can access the hourly measurements of both the PSP and BepiColombo positions, we can easily solve Equation (10) to find the time series of the expected solar wind speed measurements at any time for PSP to hit BepiColombo as The intersection between our theoretical estimates via Equation (11) and the solar wind speed measurements by PSP gives us the time instant t 0 at which the plasma parcel observed by PSP would be observed by BepiColombo at the time = + D t t V r 1 0 SW . As reported in Figure 2, this occurs at 04:59 UT on 2020 September 25. Given the measured solar wind speed V SW = 240 km s −1 and Δr = 0.4775 au, and assuming a constant solar wind as a first-order approximation (Telloni et al. 2021), the corresponding travel time is τ = 3.4528 days, thus reaching BepiColombo's orbit on 2020 September 28 at 15:51 UT. It is also interesting to note that another intersection between the two curves could occur around 19:00 UT on 2020 September 24, but due to the lack of PSP measurements, it cannot be considered meaningful. To further assess our estimates, we also follow the same procedure reported in Telloni et al. (2021; see also Zank et al. 2017) by considering pairs of 2 hr length moving windows to compare the magnetic field measurements by PSP and BepiColombo. We select the time window corresponding to the radial alignment where the mutual information (MI) between the magnetic field intensity B at PSP and BepiColombo reaches its maximum value. The MI of a pair of time series (x(t j ), y(t k )) is defined as where p(x, y) is the joint probability of observing the pair of values (x,y), while p(x) and p(y) are the independent distributions. For statistically independent time series, MI = 0, while for correlated time series, MI MI th , a threshold associated with a particular statistical significance level (e.g., 95%). We use the MI instead of the cross-correlation since it allows us to consider both linear and nonlinear relations between pairs of signals, while the cross-correlation is a linear operator. As presented in Figure 3, the MI reaches its maximum for the pair (September 25 05:00, September 28 16:00), as also depicted in Figure 2. Thus, in the following, we will consider for our analysis the magnetic field measurements made by PSP on 2020 September 25 between 04:00-06:00 UT and by BepiColombo on 2020 September 28 between 15:00-17:00. Figure 4 presents the magnetic field components in the RTN reference frame during the two selected time intervals. It can be easily observed that moving from PSP's orbit to BepiColombo, the magnetic field intensity decreases from B PSP ∼ 100 nT to B Bepi ∼ 10 nT. This is consistent with the predicted radial evolution of the magnetic field according to the Parker spiral (Parker 1958 . As a first step in our analysis, and for the global picture usually employed in MHD turbulence, we present in Figure 5 the behavior of the traces of the Fourier power spectral densities (PSDs) of each magnetic field component measured at the PSP and BepiColombo orbits. It is evident that a scaling-law behavior f −β is observed for both the PSP and BepiColombo measurements, covering at least two decades of frequencies, which is well in agreement with both Iroshnikov-Kraichnan and Kolmorogov pictures of turbulence, β being close to [3/2, 5/3]. However, the PSDs only provide information on the autocorrelation function (i.e., the second-order moment) of the time series, thus not directly allowing us to investigate interesting features of turbulence, like intermittency, self-organization, or emergent behaviors. For this reason, we evaluate the qth-order structure functions S q (τ) to gain more statistical information on the distributions of the fluctuations across scales. Our diagnosis of the statistical properties of the solar wind turbulence at both orbits starts with the investigation of the second-order structure function S 2 (τ), being directly related to the PSD, thus providing information on the variance distribution across timescales (i.e., the autocorrelation function ;Wiener 1930;Khintchine 1934). Figure 6 presents the second-order structure function for each magnetic field component and for both spacecraft. As also reported in recent works (e.g., Alberti et al. 2020;Chen et al. 2020), the near-Sun solar wind is characterized by a second-order scaling exponent ζ(2) ∼ 1/2, while moving away from the Sun (i.e., r  0.4 au), we observe ζ(2) ∼ 2/3. This variation has been attributed to the decreasing Alfvénic nature of the solar wind (e.g., Alberti et al. 2020;Chen et al. 2020). Furthermore, we clearly observe larger increments for B T and B N with respect to B R near the Sun (i.e., at PSP's orbit), which seems to suggest the existence of larger gradients in the perpendicular direction of the main field, the latter being mostly aligned along the radial direction. This seems to suggest that a 2D geometry of field fluctuations is observed, while the fluctuations are almost inhibited along the radial direction, due to the presence of a strong magnetic field. Conversely, similar increments between the three magnetic field components are observed at BepiColombo's orbit (see the stars in Figure 6). Finally, for both PSP and BepiColombo, a clear scalinglaw behavior is observed over almost two decades of scales, corresponding to the MHD/inertial regime, i.e., τ ä (5500) s, which is also consistent with previous analysis (e.g., Alberti et al. 2020;Chen et al. 2020;Telloni et al. 2021). Furthermore, a signature of large-scale dynamics at BepiColombo for τ  10 3 s is observed, which could be associated with the emergence of a large-scale forcing regime. However, the latter cannot be completely characterized, due to the reduced length of the time series (i.e., 2 hr), precluding any significant estimation over this reduced range of scales (less than one decade).
As the second step of our analysis, we evaluate the highorder scaling exponents ζ(q), as in Equation (7), across the inertial range of scales τ ä (5500) s. The behavior of ζ(q), with corresponding 95% confidence, as a function of q is reported in Figure 7. First of all, a nonlinear convex behavior is observed both for the PSP and BepiColombo measurements for positive q, thus suggesting the multifractal nature of the underlying fractal support of the field fluctuations. This points toward the existence of a non-self-similar distribution of the energy across scales, in agreement with recent studies on PSP (e.g., Chhiber et al. 2021) and the wide literature at larger heliocentric distances (e.g., Wawrzaszek et al. 2015;Bruno & Carbone 2016;Wawrzaszek et al. 2019;Alberti et al. 2020, and references therein). Furthermore, it can be easily noted that for positive q, ζ(q) Bepi > ζ(q) PSP , thus suggesting the increased role  of singularities at the PSP orbit with respect to those at the BepiColombo one (see below for a more detailed discussion). Conversely, by looking at q < 0, there seems to be the opposite condition for the (T, N) plane, while the clear linear behavior of ζ(q) with q is observed for the radial component of the magnetic field observed at PSP, lying along the expected q/3 behavior of the Kolmogorov theory of turbulence (Kolmogorov 1941). Thus, neither fluid (q/3) nor MHD Figure 6. The second-order structure function S 2 (τ) as a function of the timescales τ as derived via Equation (6). The red, green, and blue circles (stars) correspond to the Radial (R), Tangential (T), and Normal (N) components at the PSP (BepiColombo) orbits, respectively. The errors on the timescales Δτ are obtained as the standard deviations of the instantaneous timescales derived for each empirical mode, while the errors on the second-order structure function are evaluated by means of error propagation rules, i.e., since S 2 (τ) ∼ τ ζ(2) , then ΔS 2 (τ) = ζ(2) τ ζ(2)−1 Δτ. The dotted and dashed-dotted lines, referring to the Iroshnikov-Kraichnan prediction of MHD turbulence τ 1/2 and to the Kolmogorov theory of fluid turbulence τ 2/3 , are shown for reference. The vertical lines mark the extension of the inertial range we used for evaluating the scaling exponents (see the text for more details). Figure 7. The high-order scaling exponents ζ(q) as a function of the moment order q for PSP (left panel) and BepiColombo (right panel). The bar corresponds to the 95% confidence level. The red, green, and blue symbols correspond to the Radial (R), Tangential (T), and Normal (N) components, respectively. The dotted line refers to the Iroshnikov-Kraichnan prediction of MHD turbulence q/4, while the dashed-dotted line corresponds to the Kolmogorov theory of fluid turbulence q/3. (q/4) self-similar scale invariance is observed for either case, since they are instead characterized by a time-and scale-dependent energy transfer rate (Mandelbrot 1974;Benzi et al. 1984).
Once the scaling exponents ζ(q) have been evaluated, using the Legendre transform, e.g., Equations (8)-(9), we can estimate the singularities h and their spectrum D(h), as shown in Figure 8. It is evident that the most probable values are near h = 1/4 for PSP and h = 1/3 for BepiColombo, thus suggesting that both cases present features of MHD and fluid turbulence, respectively. By looking at the perpendicular plane (T, N) with respect to the radial direction, we can observe that a symmetric shape of D(h) is found for the PSP measurements, while a right-hand asymmetric form characterizes the BepiColombo measurements. This suggests that the energy transfer across scales at BepiColombo's orbit and across the radial direction at PSP's orbit is operating through an avalanching nonmultiplicative process dominated by convergent and divergent singularities (organized small and large fluctuations), respectively. Conversely, the symmetric shape in the T-N plane observed at PSP's orbit seems to point toward the existence of a multiplicative process. By looking at the radial components, we can observe an interesting situation in which a right-hand asymmetric shape is again observed at BepiColombo, thus highlighting the predominant role of organized small fluctuations on large gradients at 0.6 au; conversely, a left-hand-shaped singularity spectrum is observed at PSP's orbit, thus suggesting that along the radial direction, the near-Sun energy transfer only occurs via large gradients. This latter result has not been directly highlighted before, and deserves a deeper discussion in Section 6. The above results are summarized in Table 1 by evaluating some measures of the complexity introduced in Section 4.2. Specifically, by looking at the singularity width D = h h h max min , we can easily observe a wider range of local scaling exponents h of turbulent fluctuations at BepiColombo's orbit (Δh ∼ 0.6) with respect to PSP's one (Δh  0.6). This suggests that a more efficient cascade mechanism operates at 0.6 au, being a system more able to dissipate the excess of energy by means of small fluctuations representative of a smooth fractal field. Moreover, by evaluating the fractal width D = -D D D max min , we can easily observe that the magnetic field sampled by both PSP and BepiColombo is characterized by an underlying multifractal topology (e.g., Chhiber et al. 2021). The slightly larger values found at BepiColombo seem to support the idea of a much more efficient cascade mechanism than that at PSP (e.g., Alberti et al. 2020). Finally, by evaluating the asymmetry of the singularity spectrum = - being D(h 0 ) = 1, we can confirm that a right-hand-skewed multifractal spectrum is observed at BepiColombo's orbit, being representative of a more regular dissipation field mostly driven by small fluctuations. Conversely, a symmetric topology is found at PSP's orbit in the perpendicular plane to the main field, suggesting a concurrent

Discussion and Conclusions
In this work, we have investigated the evolution of the solar wind magnetic field scaling properties between the PSP and BepiColombo orbits during their first radial alignment when they sampled the same plasma parcel. We found that the second-order structure function, closely related to the PSD behavior, shows a steeper scaling at BepiColombo's orbit, being the scaling exponent ζ(2) ∼ 2/3, than at PSP's one, being ζ(2) ∼ 1/2. This result is in agreement with previous studies (e.g., Alberti et al. 2020;Chen et al. 2020), and can be related to the increasingly Alfvénic nature of the near-Sun solar wind, reflected in a shallower spectrum (Telloni et al. 2021). By evaluating the highorder scaling exponents ζ(q), up to the statistically accessible orders (i.e.,  q N log p max 10 ( ), N p being the number of data points; Dudok de Wit 2004), we have reported the evidence for the nonlinear convex behavior of ζ(q) as a function of q, usually related to the phenomenon of intermittency (e.g Benzi et al. 1984;Carbone 1993;Frisch 1995;Bruno & Carbone 2016), both at PSP's and BepiColombo's orbits. However, by only looking at the behavior of the scaling exponents, ζ(q) does not allow us to have a complete view of the underlying fractal topology and the mechanisms behind the energy transfer across scales. For this purpose, we used the multifractal formalism introduced in the 198os to investigate the behavior of the energy transfer across the inertial range of scales in terms of the competitive role of "regular" regions, where the energy transfer is dominated by small fluctuations, and singular regions, where the energy transfer occurs through localized bursts and large gradients (Benzi et al. 1984;Frisch 1991Frisch , 1995Dubrulle 2019;Nguyen et al. 2019). With this formalism, we are also able to characterize the nature of the inertial range in terms of the distribution of the rescaling exponents h, i.e., the presence of all possible symmetries of MHD equations, at different scales, thus providing additional information on the dependence of the energy transfer rate on scales. Our results evidence that, although the most probable values are h 0 = 1/4 for PSP and h 0 = 1/3 for BepiColombo, an entire spectrum of scaling exponents is observed at both locations. However, a clear difference emerges when looking along the different directions, especially at PSP's orbit. Indeed, by looking at the perpendicular plane (T, N) with respect to the radial direction, a symmetric shape of the multifractal spectrum D(h) is found for the PSP measurements, while a right-hand asymmetric form is observed at BepiColombo's orbit. This suggests a concurrent effect of small fluctuations and large gradients at PSP, while the dominant role of organized smooth fluctuations is observed at BepiColombo's orbit. A righthand asymmetric shape is again observed at BepiColombo's orbit along the radial direction, while a left-hand-skewed singularity spectrum is observed at PSP's orbit. This latter result, which has not been highlighted before, suggests that along the radial direction the near-Sun energy transfer only occurs via large gradients. The above results, together with the obtained complexity measures shown in Table 1, evidence that (i) a more efficient cascade mechanism operates at 0.6 au, being a system more able to dissipate the excess of energy by means of small fluctuations; (ii) the magnetic field sampled by both PSP and BepiColombo is characterized by an underlying multifractal topology (e.g., Chhiber et al. 2021); and (iii) a more regular fractal field, mostly driven by organized small fluctuations, is observed at BepiColombo's orbit, while a more singular topology characterizes the near-Sun magnetic field, suggesting a concurrent effect of both singular and regular regions in the perpendicular plane and the dominant role of singular regions along the radial direction.
To conclude our work, we try to interpret our results in terms of the possible physical processes behind the observed fractal topology near and far away from the Sun. Clearly, the energy transfer rate occurs through scale-to-scale nonlinear interactions between fluctuations (eddies in the Richardson picture) at both PSP's and BepiColombo's orbits, but with different roles of large versus small gradients. This suggests that the fractal structure of the inertial range evolves during the expansion. A heterogeneous structure, mainly dominated by large gradients, especially along the radial direction, where the main field resides and inhibits small fluctuations, is observed at PSP's orbit, while a more homogeneous structure, dominated by small fluctuations with a smooth fractal topology, is observed at BepiColombo's orbit. This could be interpreted as a different efficiency in the cascade mechanism in transferring energy across scales: a complex organization of the large field gradients is needed to dissipate the excess of kinetic energy near the Sun, whereas the topological organization of small fluctuations is primarily responsible for the energy transfer rate far away. Moreover, the asymmetry of the singularity spectra at BepiColombo's orbit and for the radial component measured at PSP suggests an avalanching nonmultiplicative process dominated by convergent and divergent singularities (organized small and large fluctuations), respectively. Conversely, the symmetric shape in the T-N plane observed at PSP's orbit seems to point toward the existence of a multiplicative process. These findings seem to suggest the role of the mean field in organizing fluctuations at different heliocentric distances, acting as a sort of order parameter (Stumpo et al. 2021). Clearly, these findings cannot only be related to a single specific parameter, like the mean field, but they could also be related to different solar wind properties evolving with distance. Indeed, as shown by Chen et al. (2020), when approaching the Sun the solar wind becomes increasingly more Alfvénic, with the Alfvénicity s =d , and the ratio between the outward and inward fluctuations increases up to |z + | 2 /|z − | 2 ∼ 15, although the scaling of their power spectra is the same. Further investigations will be needed to investigate additional features via comparisons between slow and fast solar wind streams to provide a more complete overview of the physical reasons and mechanisms behind the radial evolution of solar wind turbulence. Moreover, additional studies should be conducted to investigate what is the widest angular separation between pairs of probes, such that our findings and assumptions still hold, by using different spacecraft alignments and coordinated observations, which could also specifically be planned in the future.
The data used in this study are available at the NASA Space Physics Data Facility (SPDF): https://spdf.gsfc.nasa.gov/ index.html. The FIELDS experiment on the PSP spacecraft was designed and developed under NASA contract NNN06AA01C. We acknowledge the contributions of the FIELDS team to the PSP mission. T.A. and A.M. acknowledge