Applying the Horizontal Visibility Graph Method to Study Irreversibility of Electromagnetic Turbulence in Non-Thermal Plasmas

One of the fundamental open questions in plasma physics is the role of non-thermal particles distributions in poorly collisional plasma environments, a system that is commonly found throughout the Universe, e.g., the solar wind and the Earth’s magnetosphere correspond to natural plasma physics laboratories in which turbulent phenomena can be studied. Our study perspective is born from the method of Horizontal Visibility Graph (HVG) that has been developed in the last years to analyze time series avoiding the tedium and the high computational cost that other methods offer. Here, we build a complex network based on directed HVG technique applied to magnetic field fluctuations time series obtained from Particle In Cell (PIC) simulations of a magnetized collisionless plasma to distinguish the degree distributions and calculate the Kullback–Leibler Divergence (KLD) as a measure of relative entropy of data sets produced by processes that are not in equilibrium. First, we analyze the connectivity probability distribution for the undirected version of HVG finding how the Kappa distribution for low values of κ tends to be an uncorrelated time series, while the Maxwell–Boltzmann distribution shows a correlated stochastic processes behavior. Subsequently, we investigate the degree of temporary irreversibility of magnetic fluctuations that are self-generated by the plasma, comparing the case of a thermal plasma (described by a Maxwell–Botzmann velocity distribution function) with non-thermal Kappa distributions. We have shown that the KLD associated to the HVG is able to distinguish the level of reversibility that is associated to the thermal equilibrium in the plasma, because the dissipative degree of the system increases as the value of κ parameter decreases and the distribution function departs from the Maxwell–Boltzmann equilibrium.


Introduction
In a turbulent collisionless plasma (in which Coulomb collisions are neglected), movement on a kinetic scale (spatial scales of the order of the particles Larmor radius or skindepth) occurs in a chaotic manner, and it is determined by large-scale collective behavior and also localized small-scale processes. This kind of system can be commonly found throughout the Universe. The solar wind and the Earth's magnetosphere correspond to natural plasma physics laboratories, in which plasma phenomena can be studied [1]. Some non-linear phenomena include magnetic reconnection [2], collisionless shocks [3], electromagnetic turbulence [4], collisionless wave-particle interactions [5], or plasma energization and heating [6]. One of the fundamental open questions in plasma physics is the understanding of the energy equipartition between plasma and the electromagnetic turbulence, and the role of non-thermal plasma particles distributions ubiquitous in poorly collisional plasma environments.
The representation of the plasma velocity distribution function (VDF) using the wellknown Tsallis or Kappa distributions is one of the most used approaches to model nonthermal plasma systems. First proposed by Olbert [7] and Vasyliunas [8] to fit electron measurements in the magnetosphere, it is accepted that Kappa distributions are the most common state of electrons, see e.g., [9,10], and they have been observed in space in the solar wind [11,12], the Earth's magnetosphere [13,14], or other planetary environments [15]. These distributions resolve both the quasi-thermal core and the power-law high energy tails that were measured by the κ parameter, and they correspond to a generalization of the Maxwell-Boltzmann distribution, achieved when κ → ∞. Kappa distributions have been widely studied in the framework of non-equilibrium statistical mechanism as corresponding to a class of expected probability distribution function when the system exhibits non-extensive entropy [16][17][18]. Regarding plasma physics, it has been found that, in Kappa-distributed plasmas, the non-thermal shape of the distribution function plays a key role on the details of kinetic processes, such as wave-particle interactions [19,20], which mediate the collisionless relaxation of unstable plasma populations [21][22][23]. Moreover, in a plasma with finite temperature, the random motion of the charged particles composing the plasma produces a finite level of electromagnetic fluctuations, even in the absence of instabilities. These fluctuations, known as quasi-thermal noise, can be explained by a generalization of the Fluctuation-Dissipation Theorem, see e.g., [24,25] and the references therein, and they have been studied in the case of thermal and non-thermal plasma systems. Recent results have shown that the fluctuations level in plasmas, including supra-thermal particles following a Kappa distribution, is enhanced with respect to plasma systems in thermodynamic equilibrium [25][26][27].
Regardless the nature of the distribution function (thermal or non-thermal), plasmas show a self-organized critical behavior [28], allowing for the introduction of concepts from complex systems to study this criticality. Those methods are applied both in data sets and models [29][30][31][32]. Some authors have suggested that the change of fractals and multifractals indexes could be associated with dissipative events or related to the solar cycle, proposing a relation between multifractality and physical processes in plasmas, among them solar cycle, Sun-Earth system, or theoretical models of plasmas [29,33]. Another studies show a relation between intermittency fluctuations and multifractal behavior, while the fluctuations at kinetic-scales reveal a monofractal behavior [34][35][36]. Wawrzaszek et al. [33] apply a multifractal formalism to the solar wind, suggesting a relation between the intermittency and degree of multifractality. Those studies show different time series analysis in plasmas. However, not only fractals and multifractals could be useful in the study of time series, complex networks, particularly the Visibility Graph method, allow for a simple and direct time series analysis in self-organized critical phenomena, such as earthquakes [37], macroeconomic systems [38], or biological systems [39].
The method of Visibility Graph [40] has been developed in the last years; it allows us study and analyze time series avoiding the tedium and the high computational cost that other methods offer. The visibility algorithm proceeds to map a times series into a complex network under a geometric principle of visibility, in this sense, the algorithm could be considered to be a geometric transform of the time series in which this method decomposes a time series in connections between nodes that could be repeated or not, forming a particular weave that represents the time series as a geometric object. Inside the Visibility Graph algorithm, we can use a simplification of it, the Horizontal Visibility Graph (HVG), where the nodes are connected if it is possible to draw a horizontal line between two nodes.The HVG has been applied to different systems, from earthquakes [37] and plasmas [41], to chaotic processes [42]. The visibility graph is constructed under the visibility criterion, two data (t a , y a ) and (t b , y b ) in the time series look at each other if there are data (t c , y c ), with t a < t c < t b , which satisfy the condition [40,43]: In the field of space plasma physics, the VG has been applied to solar flares [44,45] and solar wind measurements [41]. In particular, Najafi et al. [45] show a complete and detailed analysis of solar flares through a combination of two methods of complex networks: a time-based complex network that is supported by the work of Abe and Suzuki [46] and the VG method that was proposed by Telesca and Lovallo [47]. They characterize solar flares that are based in the probability distribution of connectivity and clustering coefficient, finding good agreement with the results obtained in other works with seismic data sets. In addition, Suyal et al. [41] studied the irreversibility of velocity fluctuations. Through the HVG method, they calculated the Kullback-Leibler Divergence (KLD) of the fluctuations, and found that irreversibility in solar wind velocity fluctuations show a similar behavior at different distances from the Sun, and that there is a dependence of the KLD with the solar cycle. The KLD or relative entropy value, is a measure of temporary irreversibility of data sets produced by processes that are not in equilibrium, and it gives information on the production of entropy generated by the physical system, when considering a high degree of irreversibility as a chaotic and dissipative system [40]. Under this context, the recent results by Acosta et al. [48] have suggested that the use of the HVG method can provide valuable information to characterize turbulence in collisionless plasmas, and that the KLD may be used as a proxy to establish how thermal or non-thermal are the velocity distributions of a plasma, only by looking at the magnetic fluctuations and their properties.
Here, we build a complex network that is based on the HVG technique [43] applied to magnetic field fluctuations time series that were obtained from Particle In Cell (PIC) simulations of a magnetized collisionless plasma. We analyze the degree of irreversibility of magnetic fluctuations self-generated by the plasma, comparing the case of a thermal plasma (as described by a Maxwell-Botzmann VDF) with the fluctuations that were generated by non-thermal Kappa distributions. In order to understand the degree of the irreversibility as a parameter that could be related to the shape of the particles velocity distributions, we computed the KLD for different values of the κ parameter for comparative purposes and analyzed their time evolution throughout each simulation. The paper is organized, as follows. In Section 2, the methods and techniques are described, Section 3 shows the model used to build the time series, and, in Section 4, the results are presented. Finally, in Section 5, we discuss our results and present the conclusions of our study.

Horizontal Visibility Graph: Mapping Time Series to Network
We use the directed version of the Horizontal Visibility Graph method that models the time series as a directed network according to a geometric criteria that considers the magnitude of each data and then evaluates its horizontal visibility with the other data in the series in the direction of time (see Figure 1 for a graphical illustration). More precisely, first, let {x i } i=1,...,n be a time series of n data. The algorithm consists of assigning each data of the series to a node. Subsequently, two nodes i and j in the graph are connected if one can draw a horizontal line in the time series joining x i and x j that does not intersect any intermediate data height. Hence, i and j are two connected nodes if the following geometrical criterion is fulfilled within the time series [49]: Afterwards, for a graph directed in the direction of the time axis, for a given node, two different degrees are distinguished. These are the in-going degree k in , related to how many nodes see a given node i, and an out-going degree k out that is the number of nodes that node i sees [49]. With this temporal direction, causality in the network is implicit in each degree, since the input degree k in is associated with the links of a node with other nodes of the past. Meanwhile, the degree of output k out is associated with the links with nodes of the future [43]. From the properties of these connections, it can be said that, if the graph remains invariant under the reversion of time, it could be stated that the process that generated the series is conservative [49]. degree k in for in-going links and k out for out-going links of each of the n = 10 nodes are detailed. Bottom, probability distribution P in relation to degree k, where n in and n out correspond to the frequency of appearance of the degrees k in and k out , respectively, defining P in and P out .
In order to further detail the dynamics between nodes, degree distributions play a fundamental role. The degree distribution of a graph describes the probability of an arbitrary node to have degree k (i.e., k links) [50], and it becomes absolutely necessary to measure the difference between P in and P out to understand the incidence of the action of time on the process that originates the series (see Figure 1 for an illustration of both degree distributions).

Kullback-Leibler Divergence: Measuring Irreversibility
We are interested in measuring the irreversibility of the time series, since it is indicative of the presence of nonlinearities in the underlying dynamics and is associated with systems driven out-of-equilibrium [51,52]. A stationary process X(t) is said to be statically reversible if, for each N, the series {X(t 1 ), ..., X(t N )} y {X(t N ), ..., X(t 1 )} have the same probability of degree distribution [53]. In other words, we speak of reversibility if, for all the data of the series developed in its natural time and its inverse time, the same degree distribution is presented in each case. A method was proposed to measure real-valued time series irreversibility, which combines two different tools: the HVG algorithm and the Kullback-Leibler Divergence [43]. The degree of irreversibility of the series is then estimated by the Kullback-Leibler Divergence, a way to evaluate the difference between P in and P out of the associated graph. The KLD is defined by Thus, if P out > 0 and P in = 0, then D = ∞ [54]. The KLD is always non-negative and it is zero if and only if P out = P in , so, if D → 0, the system has a low degree of irreversibility. We estimate the relative entropy production of the physical process that generated the data, since this measure gives lower bounds to the entropy production [43].

Particle in Cell Simulations: Thermal and Non-Thermal Plasma Particle Distributions
To build time series of magnetic fluctuations that are produced by a collisionless plasma, we performed PIC simulations. The simulations treat positive ions and electrons as individual particles that are self-consistently accelerated by the electric and magnetic field through the charge and current densities that are collectively produced by themselves. For our study, we consider a so-called 1.5D PIC code, which resolves the movement of the particles in one dimension, but computes the three components of the velocity of each particles and, therefore, the three components of the current density. Our code has been tested and validated in several studies, see e.g., [55,56], and technical details about the used numerical schemes can be found in [26].
We simulate a magnetized plasma composed by electrons and protons with masses m e and M p , respectively, and realistic mass ratio M p /m e ∼ 1836. We assume the warm plasma as quasineutral, in which both species have number density n 0 , such that ω pe /Ω e = 5.
Additionally, ω pe = 4πn o e 2 /m e 1/2 is the plasma frequency, Ω e = (eB 0 )/(m e c) is the electron gyro-frequency, e is the elementary charge, c the speed of light, and B 0 is the background magnetic field. Our code solves the equations in a one dimensional grid with periodic boundary conditions, and the background magnetic field aligned with the spatial grid (B 0 = B 0x ). To resolve the kinetic physics of electrons, we set up a grid of length L = 256 λ e , where λ e = ω pe /c is the electron inertial length. We divide the grid in N = 2048 cells, initially with 1000 particles per species per cell, and run the simulation up to t = 1330.72/Ω e in time steps of length dt = 0.01Ω e . For each simulation, we initialize the particles velocities following an isotropic VDF f j (v), with j = e for electrons and j = p for protons, and v represents the velocity. For the case of a plasma in thermodynamic equilibrium, f j corresponds to a Maxwell-Boltzmann distribution and, in the case of a non-thermal plasma, f j is given by a Kappa distribution. Namely: Here, α j = 2k B T j /m j 1/2 is the thermal velocity of the distribution, κ j and T j are the kappa parameter and the temperature of each species, and k B is the Boltzmann constant. Additionally, Γ corresponds to the Gamma function, and note that Kappa distributions Equation (5) becomes the Maxwell-Boltzmann distribution Equation (4) in the limit κ → ∞. However, for kappa values κ 10, the Kappa and Maxwellian VDFs are relatively similar. Following all of these consideration, for our study we run and compare the results of three different simulations with different values of the electron κ e parameter. Case 1: a plasma in thermal equilibrium with electrons following a Maxwell-Boltzmann distribution given by Equation (4); case 2: non-thermal electrons following Equation (5) with κ e = 3, representing a system far from thermodynamic equilibrium; and, case 3: a plasma with κ e = 15, also non-thermal but closer to equilibrium. In addition, to isolate the effects of thermal or non-thermal electrons, for all three cases we consider protons following a Maxwellian; i.e., κ p → ∞. Finally, for all cases, we consider a plasma with temperature T j , such that the plasma beta parameter is β j = 8πn 0 k B T j /B 2 0 = 0.01 for both species; i.e., β e = β p = 0.01.
As already mentioned, even though a collisionless isotropic plasma is a system at equilibrium according to the Vlasov Equation, the plasma will develop a certain level of magnetic fluctuations that are spontaneously produced by the motion of the charged particles [24][25][26][27]. This is precisely the situation of our study for any of the three simulations (three cases) that we have performed. Figure 2 shows the average magnetic field energy density fluctuations (δB/B 0 ) 2 as a function of time, for all three cases. To build these time series, at each time step we have computed the transverse magnetic fluctuations at the plane perpendicular to B 0 and we have averaged the magnitude of the fluctuations at each grid point. Figure 2 (left) shows the fluctuations time series for κ e = 3 (blue), κ e = 15 (purple), and the κ e → ∞ or Maxwell-Boltzmann distribution (black). As expected, we can see that the level of fluctuations increases with a decreasing value of κ e , and that the behavior of the fluctuations with κ e = 15 is fairly similar to the Maxwellian case. In addition, Figure 2 (right) presents the time series of the detrended fluctuations, where we can see that the amplitude of the fluctuations also increases as κ e decreases. In the next section, the HVG method will applied to all of these time series.

Results
We apply the HVG method to study the time series of magnetic fluctuations that were obtained from the PIC simulations. When considering the Maxwellian and Kappa distributions, we follow the HVG algorithm and build complex networks for three cases: Maxwellian distribution (thermal equilibrium with electrons), κ e = 3 (non-thermal electrons), κ e = 15 (non-thermal electrons, but closer to the equilibrium), while using the original and the detrended time series (see Figure 2). First, we build the complex network and calculate the in-going and out-going degree distribution for each time series, in order to characterize this distribution for each case. In this sense, the probability distribution of the degree gives information that is related to the correlations in a process, in this case, time correlations. According to the theorem for uncorrelated time series [49], the degree distribution of the horizontal visibility graph that is associated with a bi-infinite sequence of independent and identically distributed random variables extracted from a continuous probability density is P(k) ∼ exp(−γ un k) with γ un = ln(3/2) ≈ 0.405. When the results move away from this critical value γ un we are in the presence of correlations and this value is a border between correlated and chaotic stochastic processes, i.e., γ < ln(3/2) characterizes a chaotic process, whereas γ > ln(3/2) characterizes a correlated stochastic one [42].
If we now focus on the undirected HVG, k(i) = k in (i) + k out (i), we obtain the undirected degree distribution P(k). Figure 3 show an exponential distribution for the degree distribution for the three cases studied, this is understood as short-range exponentially decaying correlations, where γ corresponds to the slope of the linear fit in degree distribution semilog. The values of the slope are computed considering the tail of the distribution [47] in Figure 3, in this case from the degree k = 5 up to the largest value of k at each plot. The values of the slope are between γ = 0.531 (κ e = 3) and γ = 0.590 (Maxwell-Boltzmann distribution), in the case of trended data and between γ = 0.520 (κ e = 3) and γ = 0.579 (Maxwell-Boltzmann distribution), in the case of detrended data sets. We observe that, for each value of the slope, it is satisfied that γ > γ un . This shows us all series corresponds to correlated stochastic processes from which we can extract consistent information. Additonally, the trend does not seem to greatly affect these correlations. Second, when considering the directed HVG, we have computed the Kullback-Leibler Divergence, D from Equation (3), for each case mentioned before. The values of the divergence D are in Figure 4 as compared to standard deviation σ (vertical bars in the figure) calculated from the algorithm to the disarrayed randomly data, see e.g., [37] and references therein. In this algorithm the original data in randomly shuffled to obtain a large number of disordered copies (in this case, 1000 copies) of the original data set, and the divergence D is computed for each copy. The vertical lines presented in Figure 4 correspond to the average value of the divergence D of all copies (central value) plus and minus a standard deviation. After this procedure, if the D value is contained inside the σ bar, the time series represents a reversible process. This is because, by randomly disarraying the series and obtaining the same results, regardless of the temporal order of the data set, by definition it indicates that the information corresponds to a reversible process. On the contrary, if D is outside the vertical range that is defined by the random copies, then the value of D is statistically significant and, therefore, it is possible to state that the data set indeed represents an irreversible process. The technique used to determine whether the data represent a reversible process consists of applying the HVG algorithm to randomly disordered copies of the data, obtaining the standard deviation σ around the average divergence computed using the disordered data (black dot and vertical lines). Figure 4 shows that the dissipative degree of the system increases as the value of κ e decreases and the distribution function departs from the Maxwell-Boltzmann equilibrium. In Figure 4 (left), the processes for κ e = 15 is reversible, case close to thermal equilibrium. Meanwhile, in Figure 4 (right), all of the distributions correspond to reversible processes by reducing the background trend. This last result could be explained due to the fact that, independent of the value of κ e , all of the considered distributions are steady state solutions of the Vlasov equation. Finally, to further analyze the relationship between the κ parameter and the KLD, we compute the time evolution of D, as shown in Figure 5. Figure 5 (right) show the same behavior found above, exhibiting a decrease in the value of the divergence for the Maxwellian distribution, whereas, for κ e = 3, this value tends to increase. That is, given the initial conditions of the simulation, κ e = 15 and Maxwellian coincide in their behavior over time towards a low degree of divergence, while κ e = 3 evidently presents a behavior to the opposite extreme.

Discussion and Conclusions
In this study, we have modeled a turbulent plasma as a complex network, applying the method that is known as Horizontal Visibility Graph to study the reversibility on magnetic fluctuations. We have developed algorithms to build HVGs from magnetic field fluctuations time series obtained from PIC simulations of collisionless magnetized plasmas. We have analyzed three cases for the time series: a time series of a plasma far from the thermodynamic equilibrium (κ e = 3), a time series closer to the thermodynamic equilibrium (κ e = 15), and a Maxwell-Boltzmann distribution, representing a plasma in a thermal equilibrium. For these three time series, we have computed the degree distribution of the connectivity, which gives information that is associated with the time correlations in the distribution and the KLD, which provides information that is related to the reversibility of the time series.
In the case of the degree probability distribution, we have found an exponential behavior for all cases analyzed, i.e., a short-range correlations for all time series (Kappa and Maxwell-Boltzmann distributions). Our results show that the decaying critical exponent γ is the largest for the Maxwellian-Boltzmann distribution, and it decreases with a decreasing Kappa value. Moreover, for κ e = 3, the critical exponent is closer to the limit value γ un = ln(3/2) proposed by Lacasa and Toral [42] in which the time series becomes uncorrelated, being chaotic for smaller values (γ < γ un ). These results suggest a lower time correlation for κ e = 3 than the Maxwell-Boltzmann distribution, which is consistent with the fact that, in collisionless plasmas out of thermodynamic equilibrium long-range interactions dominate [18]. As already mentioned, in all cases, our simulations correspond to isotropic plasmas that are steady-state solutions of the Vlasov equation. Thus, the electromagnetic fluctuations correspond to spontaneous emissions of a system that is composed by discrete charged particles in random motion. Consequently, the fluctuations provide information regarding the smallest scales where fast short-range interaction dominate. In the case of a plasma system, these scales are strongly related with the Debye length λ D .
Inside the Debye sphere (a sphere of radius λ D ), particles interact individually and, outside the Debye sphere, long-range collective interactions dominate. This is directly related to the correlations between the particles that produce the magnetic fluctuations, which depend on the shape of the velocity distribution function. In the case of Kappa distributions, the Debye length of the plasma is a decreasing function of κ that collapses to zero for κ = 3/2 [57]. Therefore, in plasmas described by a Kappa VDF, the short-range correlations are less effective, since the Debye length is smaller. Outside the Debye sphere, the thermal energy dominates the potential energy, and the correlations are practically dissolved [11]. In contrast, since the Debye length is greater, in a Maxwellian plasma the short-range correlations dominate, as they decay faster, both temporally and spatially. Regarding our results, this is reflected in the gamma value that seems to behave as a increasing function of κ e .
In the case of the KLD, for both the original and detrended time series, we have obtained low values of the divergence D for all cases, which is consistent with plasmas in steady state according to the Vlasov equation. However, the method gas shown to be sensitive enough to distinguish higher values of irreversibility for the Kappa distribution than the Maxwell-Boltzmann case. The irreversibility that is associated to the Kappa distributions is related to the non-extensive nature of these distributions [58], showing an increase in the value of the KLD for decreasing values of κ. The increase in the value of the KLD indicates a larger value of the entropy in the system. For Kappa distributions following the dynamics of a non-collisional plasma, particles lose individuality and interact collectively increasing entropy [18]. On the other hand, the Maxwell-Boltzmann distribution shows low values for the KLD, being consistent with the Gibbs-Boltzmann entropy. The Maxwell-Boltzmann distribution is related to low values of entropy, in contrast to non-thermal Kappa distributions where it is possible to find a higher (non-extensive) entropy, associated to electromagnetic long-range interactions that dominate the dynamics in the plasma.
In summary, when only considering the limited information provided by the time series, our results seem to indicate a robust relation between the shape of the VDF (as given by the Debye length and its dependence on κ) and the nature of the correlations dominating the magnetic field fluctuations time series represented by γ [42]. The connectivity probability distribution shows how the Kappa distribution for low values of κ tends to be an uncorrelated time series, while the Maxwell-Boltzmann distribution shows a stochastic time series behavior. Furthermore, we can see that the KLD that is associated to the HVG is able to distinguish the level of reversibility in time series that were obtained from PIC simulations, and this reversibility seems to be associated to the thermal equilibrium in the plasma. Our results suggest a high sensitivity of the HVG algorithm and a relationship between KLD, κ, and the entropy of the system. The technique that is applied here has allowed us to address the role of non-thermal particles distributions in poorly collisional plasma environments. We expect all of these features to provide a framework in which complex networks analysis may be used as a relevant tool to characterize turbulent plasma systems, and also as a proxy to identify the nature of electron populations in space plasmas at locations where direct in-situ measurements of particle fluxes are not available.