Long-term changes in the north – south asymmetry of solar activity : a nonlinear dynamics characterization using visibility graphs

Solar activity is characterized by complex dynamics superimposed onto an almost periodic, approximately 11-year cycle. One of its main features is the presence of a marked, time-varying hemispheric asymmetry, the deeper reasons for which have not yet been completely uncovered. Traditionally, this asymmetry has been studied by considering amplitude and phase differences. Here, we use visibility graphs, a novel tool of nonlinear time series analysis, to obtain complementary information on hemispheric asymmetries in dynamical properties. Our analysis provides deep insights into the potential and limitations of this method, revealing a complex interplay between factors relating to statistical and dynamical properties, i.e., effects due to the probability distribution and the regularity of observed fluctuations. We demonstrate that temporal changes in the hemispheric predominance of the graph properties lag those directly associated with the total hemispheric sunspot areas. Our findings open a new dynamical perspective on studying the north–south sunspot asymmetry, which is to be further explored in future work.


Introduction
Starting with the pioneering works of Schuster (1898) and Yule (1927), long-term variations in solar activity have been among the most studied examples of complex variability patterns in nature for more than a century, and their past presence has been shown in many paleo-archives (Bond et al., 2001;Agnihotri et al., 2002). A variety of complementary approaches has allowed attributing the complex modulations of the almost periodic 11-year Schwabe cycle of solar activity (and the associated twice-as-long cycle of solar magnetic field reversals) to the action of nonlinear processes (Ruzmaikin et al., 1994;Li and Li, 2007;Donner, 2008) related to the complex dynamics of the solar dynamo (Charbonneau, 2010). Importantly, variations in solar activity are known to have considerable influence on Earth system dynamics at various timescales. On the one hand, grand solar minima (Voss et al., 1996) are known to have triggered global cooling events. Hence, knowledge of solar activity is important to understand past climate variations (Bard and Frank, 2006). On the other hand, there are potentially severe space weather events due to eruptions on the solar surface; the thustriggered strong magnetic storms affect the functioning of electricity and telecommunication networks. The two aforementioned examples illustrate the importance of anticipating both short-and long-term behavior of solar activity (Kurths and Ruzmaikin, 1990;Brajša et al., 2009;Petrovay, 2010).
Because of the importance of factors such as total solar irradiation, solar wind, and geomagnetic field for the living conditions on the Earth, studying the dynamical characteristics of the solar activity cycle has recently attracted considerable interest. While some aspects of the underlying nonlinear dynamics have been revealed in the past using a variety of different analysis techniques, there are a multitude of observed phenomena that have not yet been fully understood. One of these findings is on the hemispheric asymmetry of solar activity, which manifests itself in the statistical properties of a variety of activity indicators such as sunspot numbers, areas and spatial distribution, the numbers of flares and coronal mass ejections, solar radio and Xray flux, and has been recognized to vary on multi-decadal timescales (see, e.g., Newton and Milsom, 1955;Carbonell et al., 1993;Zolotova and Ponyavin, 2006;Donner and Thiel, 2007;Donner, 2008;Li et al., 2008;Li, 2008;Zolotova et al., 2009, and references therein). Notably, it is commonly believed that the observed distinct hemispheric asymmetry is an intrinsic property associated with underlying solar magnetic field dynamics, which in turn serves as the driver of solar activity responsible for particle and electromagnetic emissions directly affecting the Earth.
Traditionally, the north-south asymmetry of solar activity has been mainly defined in terms of amplitude differences between the hemispheric values of different properties (Newton and Milsom, 1955). Recently, it has been argued that phase information should be explicitly taken into account as well, noting that a mutual time shift between the activity cycles observed separately at both solar hemispheres could provide a significant contribution to the observed asymmetry (Zolotova and Ponyavin, 2006;Donner and Thiel, 2007). As a particularly powerful tool, an extraction of the main variability components related to the approximately 11-year Schwabe cycle by means of continuous wavelet transform has been found to be widely applied (Donner and Thiel, 2007;Donner, 2008;Li et al., 2008Li et al., , 2009). However, despite these methodological advances, properly quantifying the north-south asymmetry is a challenging problem itself. Specifically, the complex dynamics of the entire solar activity cycle calls for replacing traditional linear statistical approaches with methods originated in the field of nonlinear dynamics.
Among other fundamental paradigms of nonlinear time series analysis that could be helpful for obtaining additional information on the complex dynamics beyond the north-south asymmetry of activity magnitude and phase, recently developed complex network-based approaches to time series analysis provide prospective candidates for corresponding analyses. In recent years, a variety of approaches have been introduced by various authors (see Donner et al., 2011 for a review), with recurrence networks Donner et al., 2010a) and visibility graphs (Lacasa et al., 2008) as probably the two most prominent and widely applied concepts. Specifically, some variants of visibility graph analysis have been used recently to study the statistical properties of different solar activity indicators (Yu et al., 2012;Zou et al., 2014b). However, these analyses have been restricted to univariate time series so far, not allowing for an explicit study of interhemispheric asymmetry properties. In this work, we generalize the visibility graph approach to the bivariate case in order to study the long-term (i.e., multi-decadal) changes in the asymmetry of dynamical characteristics observed at both solar hemispheres. The remainder of this paper is organized as follows: in Sect. 2, we describe the data used in this work. Subsequently, we introduce our analysis framework and discuss different methodological options. Our findings are detailed in Sect. 3, providing (i) a general discussion of the impact of statistical vs. dynamical properties of the considered data on the observed asymmetries, (ii) the obtained long-term changes revealed by different properties based on visibility graphs, and (iii) a comparison of the obtained results with findings of other recent studies. Finally, the main results of our present work are briefly summarized in Sect. 4. Some additional results relating to different visibility graph properties of artificial and observational data are provided and discussed in two appendices.

Description of the data
We use monthly data of hemispheric sunspot areas from May 1874 to March 2013 compiled by D. Hathaway, which are available from http://solarscience.msfc.nasa.gov/greenwch. shtml. The corresponding time series are shown in Fig. 1. No additional preprocessing of the data is required in the following analysis.
The traditional way of characterizing the north-south asymmetry is to calculate either the absolute area difference or the normalized area difference (Newton and Milsom, 1955;Carbonell et al., 1993) where A N,S are the cumulated sunspot areas on the northern and southern solar hemisphere, respectively. The long-term mean behavior of both characteristics is shown in Fig. 2b, c.

Visibility graphs
In recent years, complex network approaches have tackled a great variety of challenges in various fields. Among other topics, studying time series from a network perspective has received considerable interest and has led to the development of a plethora of algorithms highlighting different aspects of the complex dynamics encoded in time series data (Zhang and Small, 2006;Xu et al., 2008;Lacasa et al., 2008;Marwan et al., 2009;Donner et al., 2010aDonner et al., , 2011. In this work, we utilize the concept of visibility graphs (VGs) (Lacasa et al., 2008;. Here, individual univariate observations are interpreted as vertices of a complex network, and edges are placed between pairs of vertices that exhibit some geometric visibility property detailed below. Let us consider a univariate time series {x(t i )} T i=1 where T denotes the time series length. As for every undirected and unweighted network, the VG is completely described by the binary T × T adjacency matrix M. The adjacency matrix has nonzero entries M i,j = 1 (corresponding to an edge between observations made at two time points t i and t j ) if and only if the convexity criterion is fulfilled for all time points t k with t i < t k < t j (Lacasa et al., 2008). An illustrative example for the construction of a VG is shown in Fig. 3a. From definition (3), it directly follows that consecutive observations are always connected. Hence, the VG does not consist of mutually disjoint subgraphs, but exhibits a giant component containing all T vertices. Furthermore, unlike most other existing approaches for transforming time series into complex network representations, the VG is not affected by the choice of certain algorithmic parameters (for example, a threshold distance ε in phase space used in the definition of recurrence networks, see Donner et al., 2010a, b). While such parameters render these alternative schemes algorithmically more complex, they increase the capability to qualitatively reconstruct (given enough data) the dynamical system under study. In turn, the resulting network characteristics typically display a certain dependence on the respective parameter.
Due to its algorithmic simplicity as a parameter-free method, for the purpose of the following analyses we prefer using the VG concept rather than other types of time series networks. Another reason for this choice is that the degree distribution of the resulting network is directly linked with the fractal properties of the underlying time series Ni et al., 2009). In recent years, the first applications to different geoscientific problems have been reported, such as the interannual changes in the frequency of landfalling hurricanes (Elsner et al., 2009), long-term variability of river discharges (Tang et al., 2010), and the scaling properties of earthquakes , wind speed records , and ocean tides . In parallel with these multiple fields of applications, conceptual problems related with the VG analysis of geophysical time series have been identified and discussed (Donner and Donges, 2012). In turn, VGs have been found useful for tackling further research questions, such as probing the time-reversal symmetry of real-world records Donges et al., 2013).

Horizontal visibility graphs
As a notable modification of the standard VG algorithm, Luque et al. (2009) proposed utilizing a simplified criterion of horizontal visibility for transforming a time series into a complex network. Specifically, they considered two observations made at times t i and t j , respectively, to be connected in a horizontal visibility graph (HVG) if and only if for all t k with t i < t k < t j .
The algorithmic difference between HVG and VG is illustrated in Fig. 3. It can be easily seen that the edge set of the HVG associated with a given time series is a subset of the edge set of the associated VG (i.e., if the horizontal visibility criterion Eq. (4) is fulfilled, then Eq. (3) also holds, but not necessarily vice versa). Specifically, VGs are invariant under affine transformations of the entire time series, whereas HVGs are not. One notable advantage of HVGs is that they provide an even higher degree of algorithmic simplicity than standard VGs, resulting in the observation that for certain simple stochastic processes, some basic graph properties can be calculated analytically ). On the other hand, the fact that HVGs typically contain a lower number of edges increases the demands regarding the time series length relatively to the standard VG when using this approach in applications such as tests for time-reversal asymmetry (Donges et al., 2013). Nevertheless, similar to normal VGs, HVGs have been successfully applied to studying time series from various fields of sciences. Within the scope of the present manuscript, we particularly highlight the recent paper by Yu et al. (2012), who studied the multifractal properties of a solar flare index in terms of HVG characteristics.

Conditional and joint degree sequence
In the following, we will introduce some basic networktheoretic quantities that we will use for the analysis of asymmetries in hemispheric solar activity. For A N (t) and A S (t), we construct two (H)VGs with adjacency matrices M N and M S , respectively. Note that the sets of vertices are the same for both graphs, with differences exclusively in the set of edges.
Based on the thus-obtained (H)VGs, we proceed as follows: 1. From the two (H)VGs, we have two sets of neighbors, N N, The degree sequences are then defined as 2. The joint degree sequence gives the number of common neighbors of a vertex corresponding to time t in both sequences.
Notably, we can define k joint (t) as the degree sequence of a joint (horizontal) visibility graph combining the visibility criteria for two distinct time series. Here, the adjacency matrix is defined by the point-wise multiplication of the individual (H)VGs' adjacency matrices. This idea is conceptually related to the concept of joint recurrence plots encoding the simultaneous recurrence of two dynamical systems in their respective phase spaces (Romano et al., 2004). As in the latter case, generalizing joint (H)VGs to K > 2 time series is straightforward, but will not be considered here given the bivariate nature of the data under study.
3. In a similar spirit to the joint degree sequence, we can quantify the number of edges associated with time t, which connects to vertices contained in N N (t) but not in N S (t), or vice versa. More specifically, ) measures the number of neighbors that belong only to N N (t) or N S (t), respectively. In what follows, k O N,S (t) will be referred to as the conditional degree sequences. By definition, Based on the latter definitions, we can proceed in a similar way as in Eq. (2) and compute the following properties: The excess degree k(t) quantifies how much "more convex" the fluctuations of A N are in comparison with A S around a given time t (i.e., how many more or less (H)VG connections the observation of A N at time t obeys in comparison with A S ). By additionally considering the relative excess degree rel k(t) normalized by the sum of the individual degrees, we obtain a measure that does not exhibit marked sensitivity with respect to the actual degrees k N,S , which may considerably vary over time according to the statistical and dynamical characteristics of the data.
In general, we suggest that the nonlinear properties k(t) and rel k(t) complement other characteristics previously used for studying the north-south asymmetry of solar activity. Notably, our approach is conceptually related with recently developed (H)VG-based tests for time series irreversibility, which compare (among others) degree distributions obtained when considering edges to past and future observations separately (Donges et al., 2013).

Time-resolved analysis
The (H)VG-based properties introduced above can be computed separately for each point in time. However, this strategy may have certain drawbacks: on the one hand, conditional and joint degree sequences exhibit only integer values, so that considering their time variation would imply dealing with discrete-valued and possibly highly fluctuating data. On the other hand, the approximately 11-year solar activity cycle will be clearly visible in the degree sequences and, hence, the joint and conditional degree sequences. Taken together, both effects can be expected to successfully undertake any interpretation of the obtained results.
For studying long-term variations of the properties of interest, we prefer considering their mean values and spread taken over running windows in time. As in Fig. 2, in all following considerations the window size will be chosen as w = 270 months, with a mutual overlap of 12 months between subsequent time windows. This specific choice of the window size covers about one full period of the solar magnetic field polarity cycle (approximately 22 years). The dependence of our results on this choice will be further discussed in Sect. 3.3. For convenience, we specify the midpoint of the moving windows as the reference for time axis labeling for all results shown in the following.

Signatures of non-Gaussianity
Related to the non-negativity of the studied sunspot areas, the considered time series exhibit strongly non-Gaussian probability distribution functions (PDFs, see Fig. 4a, b), which can change themselves with time due to the complex dynamics of the solar activity cycle. Specifically, the mass of the distribution is significantly concentrated on the left of the distribution (right-skewed or right-tailed PDF), which can be quantified by the corresponding (positive) skewness where µ is the mean of the respective random variable x, σ the associated standard deviation, and E[·] denotes the expectation value (practically approximated by the arithmetic sample mean in the following). In a similar way, one could also consider the kurtosis or other higher-order moments of the PDF of hemispheric sunspot areas, which also display marked signatures of non-Gaussianity (see Fig. 1a, b in Zou et al., 2014a). Since VGs capture subtle geometric properties of a time series, we can expect that not only the observed dynamics, but also the specific shape of the PDF has a distinct effect on the network's structural features. In order to estimate the magnitude of this effect, we initially study the simplified case of time series without any serial correlations, but with a prescribed nonzero skewness S (non-Gaussian white noise). For this purpose, we utilize sequences of independent gammadistributed random variables. In order to simplify the following discussion, we consider S as a single parameter describing the relevant changes in the PDF. We emphasize that mean and variance do not affect the VG due to its invariance under affine transformations of the variable of interest. Skewness is the lowest-order statistical characteristic that affects the shape of the PDF behind such transformations, and is -as the kurtosis -fully determined by only one of the two parameters of the gamma distribution. In turn, all moments of the gamma distribution can be fully parameterized in terms of the associated mean and skewness, justifying this restriction. Our results reveal that the skewness has indeed a marked effect on the resulting VG characteristics, particularly the edge density ρ (which is directly related with the network's mean degree as k = ρ(N − 1)). In particular, we find that as S increases, ρ increases in a first-order approximation exponentially (Fig. 5). This is a direct consequence of the convexity constraint being the foundation of the VG algorithm: the higher the skewness, the larger are the local maxima of the data in comparison with the "typical" values (e.g., the median) and, hence, the more other vertices become "visible" in the associated geometric construction. Note that the range of variations is generally rather small. Residual fluctuations superimposed onto the general trend become negligible if longer time series are used. Similar observations can also be made for other global characteristics of the associated VGs (see Appendix A).
The above results confirm that when transforming the values of a given time series in even a completely monotonous (yet nonlinear) way (i.e., changing the associated PDF but retaining the dynamics), the VG properties can exhibit significant changes, implying that they are not topologically invariant. Hence, consideration of the asymmetry of the PDF (e.g., in terms of S) is indeed of potential relevance when performing a running window analysis of the VG properties related to the north-south sunspot asymmetry below: temporal changes of the skewness of the data (Fig. 4c, d) necessarily result in variations of the edge density and, hence, the entire degree sequence. By making use of the normalized quantity rel k introduced in Sect. 2.4 instead of just the excess degree k, we reduce the possible bias due to this skewness effect. Alternatively, using the HVG instead of the standard VG serves the same purpose, but has the disadvantage that typically larger windows need to be used for obtaining reliable statistics (Donges et al., 2013).

VG analysis of the north-south asymmetry
We construct the VGs for monthly hemispheric sunspot areas, yielding the degree sequences k N (t) and k S (t), respectively. In order to magnify potentially relevant large-scale variations in the degree sequence, we utilize running windows for averaging the degree sequences over some time period (see Sect. 2.5). Note that we do not use running windows in the VG generation, which would result in potentially severe edge effects leading to a systematic downward bias in the estimated mean degrees (Donner and Donges, 2012). In Fig. 6, we display the mean features associated with the degree sequences for our sliding windows, together with the associated window-wise standard deviations. Our results reveal two transitions between periods of positive and negative mean (absolute and relative) excess degrees, which take place at about 1925-1935 (from conditions with higher VG degrees in the Northern Hemisphere than in the Southern Hemisphere to those with the opposite behavior, where the Southern Hemisphere exhibits higher degrees than the Northern Hemisphere) and 1985-1995 (vice versa). Notably, positive (negative) excess degrees imply higher mean degrees in the northern (southern) hemisphere. In order to understand this observation, we recall that there are different factors that can contribute to a higher mean degree of a VG, including higher positive skewness (Sect. 3.1) -or, more generally, stronger asymmetry of the PDF, weaker fragmentation of the VG obtained for time series with "less irregular" (e.g., more persistent) fluctuations, or a generally higher degree of time series convexity at the "local" scale (Donner and Donges, 2012).
It is most likely that in the present case, a combination of the aforementioned potential factors is responsible for the observed changes in hemispheric predominance. Regarding the shape of the underlying PDF, Fig. 4c shows that the skewness values for the sunspot areas in both hemispheres exhibit marked long-term variability. Notably, an almost fully equivalent behavior is found for the kurtosis (see Fig. 1c in Zou et al., 2014a). This observation is explained by the fact that sunspot areas are non-negative and exhibit less marked variations of the PDF during quiescence phases of the Sun, but instead distinct tail variations related to the activity maxima of each solar cycle. In such a situation, the maxima of the distribution are mainly responsible for the resulting values of skewness, kurtosis, and possibly all other higher-order moments of the PDF observed within given time intervals so that it is to be expected that the aforementioned characteristics mainly trace differences in the maximum solar cycle amplitude.
More specifically, we find that the skewness of the northern hemispheric sunspot areas attains its maximum over the considered time period at about 1920-1925, whereas the corresponding value for the southern hemisphere shows its longterm minimum at about the same time. If the underlying dynamics were the same, this would imply a tendency towards a higher density of the VG of the northern hemispheric series and, hence, a positive excess degree. In turn, between about 1970 and 1980, the skewness is much higher for the southern hemisphere, supporting a tendency towards negative excess degrees. Remarkably, the observed transitions in the sign of the excess degrees closely follow the time periods with the strongest skewness differences (Fig. 6c, d). In general, before the first transition (1925)(1926)(1927)(1928)(1929)(1930)(1931)(1932)(1933)(1934)(1935), we find generally higher skewness in the northern hemisphere, supporting positive excess degree. In turn, the sunspot areas exhibit a more asymmetric distribution in the southern hemisphere after about 1935, which is consistent with the observed negative excess degrees until about 1990. However, for recent years of the considered record, we again find positive excess degrees despite a persistently higher asymmetry of the data in the southern hemisphere. The latter observation suggests that skewness is not a unique factor determining the hemispheric predominance in terms of VG properties (i.e., the nonlinear properties captured by VGs indeed add information beyond the "linear" statistical evaluation in terms of S).
Notably, absolute and relative excess degrees exhibit qualitatively the same long-term variability. The reason we display both quantities is that the absolute excess degree can be easily interpreted in terms of interhemispheric differences, whereas the relative excess degree partially corrects for the skewness effect and allows quantitatively assessing the rel-evance of differences between the degree sequences of both hemispheres.
We should note two limiting factors that could affect the results discussed above and their possible interpretations. On the one hand, the considered asymmetries of the underlying PDFs change over relatively broad intervals in time, disregarding that the statistical properties of the hemispheric activity series could also change abruptly during such time windows. According to this, the time intervals mentioned above provide only coarse indications of the location of the speculated changes in the underlying observations. On the other hand, due to the edge effects previously discussed by Donner and Donges (2012), the results obtained for the first and last parts of the observational records have to be taken with special care since the systematic downward bias of the corresponding degrees in both VG and HVG leads to a higher variance in the obtained estimates of excess degrees.

Dependence on window width
Regarding our findings described above, we obtain qualitatively equivalent results if the window width is varied over a reasonable range. Specifically, there are no marked changes in the long-term variability of the VG-based characteristics for w being between about 180 and 400 months. In turn, windows much smaller than 180 months do not allow clear identification of transitions in the relative behavior of both solar hemispheres. In particular, for window sizes much smaller than the period of the solar cycle, the conditional and joint degree sequences exhibit marked variability at the scale of the solar magnetic cycle (not shown), which is not visible otherwise. In this case, the fluctuations at this timescale reach the order of the long-term variations seen in Fig. 6, masking the transitions in hemispheric predominance. On the other hand, for much larger window sizes w the effects of the finite time series length become more relevant: results are averaged over a considerably large part of the records, so that not only high-frequency variations, but also low-frequency variations are smoothed out.

HVG analysis of the north-south asymmetry
In order to further identify which changes in the observed VG properties are unrelated to changes in the probability distribution function of hemispheric sunspot areas, we repeated our previous analysis using the HVG algorithm replacing the classical VG. The results shown in Fig. 7 reveal some interesting facts: first of all, all degree-related quantities obey considerably lower values and weaker overall variability than for the VG. This is to be expected since the HVG is a subgraph of the VG. Another consequence of the latter fact is that since there are fewer edges in the HVG, its statistical properties are generally prone to a higher relative uncertainty level (Donges et al., 2013). However, while the absolute degree values in the HVG typically reduce by a factor of about  Figure 7. As in Fig. 6 for the corresponding HVG properties. The gray area indicates the only considerably long time interval with statistically significant hemispheric asymmetry of the mean conditional degree.
2-4 in comparison with the VG, the absolute excess degrees are by more than 1 order of magnitude smaller.
Moreover, for the HVG-based excess degree, we do not find comparably clear indications for transitions between time periods with clear hemispheric predominance as for the VG. The only notable exception is the time period between about 1925 (corresponding to the formerly identified first transition in the VG) and 1950, where the excess degree of the HVG is significantly negative (as also observed before for the VG). Specifically, the transition in the hemispheric predominance reflected by the VGs' conditional degree sequences coincides with a sharp drop in the corresponding series for the HVG at about 1925, whereas the end of the period of significantly negative excess degrees in the HVG at about 1950 accompanies the termination of the gradual downward trend of the excess degree obtained from the VGs. Taken together, we interpret these findings such that the effect of the asymmetry of the hemispheric sunspot area values mostly dominates possible variations in dynamical characteristics. However, to this end we tentatively conclude that parts of the observed long-term changes of the VG-based excess degree cannot be explained by combining the corresponding changes in skewness and HVG-based excess degree (i.e., distribution and dynamics, respectively). One possible reason for this could be complex changes in the PDF of the sunspot areas, which go beyond fluctuations in skewness yet have a significant effect on the resulting VGs' properties. . Running means and standard deviations (window width 270 months) of the degree sequence of the VG constructed from the normalized area difference series NA. Areas shaded in gray indicate the hypothesized transitions in the dynamics of interhemispheric activity differences, which are similar to those shown by the excess degree of the VGs for A N and A S (Fig. 6).

VG analysis of normalized area differences
Complementarily to our previous discussion of the excess degree, we also studied the VG for the normalized area difference series NA. Taking moving averages of the resulting degree sequences using the same sliding windows as before, we obtain estimates of the long-term variations of the corresponding mean degree (and, hence, edge density). The results shown in Fig. 8 are widely consistent with those for the excess degrees calculated from the hemispheric activity data in their tendency towards different amplitudes: between about 1910 and 1925, the edge density increases markedly, remains afterwards at a relatively stable level until about 1990 and then decreases again. Again, we could interpret this finding as more asymmetric variations of NA during the intermediate time period in comparison with the beginning and end of the series and/or temporal changes in the dynamics of NA eventually relating to a larger degree of persistence. In order to properly capture the latter property, complementary analyses using measures of long-range persistence or dynamical complexity could provide relevant information. However, such analyses are beyond the scope of the present work.
In general, we note that there is no distinct value of the mean VG degree for NA that would allow identifying two different states as there is for the excess degree of the VGs for the hemispheric areas. (In Fig. 8, we indicate k NA = 6 as some possibly reasonable value based on visual inspection rather than theoretical justification.) Hence, we do not claim the unambiguous presence of transitions in the dynamics revealed by the NA-based degree sequences. However, together with the previous information provided by the excess degrees, the VG properties associated with NA provide a largely self-consistent picture.

Discussion
When comparing the observed transition periods with those revealed by absolute area (AA) and NA (Fig. 2b, c), the corresponding properties display similar transitions from a southward-dominating (before about 1905) to a northwarddominating mode (about 1905-1980) and back (after 1980) as the VG excess degrees. However, these transitions take place considerably earlier and have the reverse orientation than those revealed by the VG degree sequences. The latter is most probably due to the fact that the classical indicators of north-south asymmetry are based on absolute (mean) amplitudes of some activity-related properties, whereas the VG-based characteristics reflect both the underlying nonlinear dynamics and PDF. In this spirit, our results indicate that the generally stronger activity in the northern hemisphere over vast parts of the 20th century has been accompanied by a weaker asymmetry and lower degree of regularity of fluctuations in comparison with the southern hemisphere, both together resulting in lower VG connectivity. A possible interpretation of the shifts in the timing of the observed transitions could be that transitions in the (relative) amplitudes correspond to non-stationarities which could mask associated changes in the dynamical characteristics, which become only visible afterwards. Complementary analyses utilizing alternative complexity measures would be helpful for investigating this idea further.
In general, we note that the transitions in AA and NA reveal information about the magnitudes of activity at both solar hemispheres. For example, activity in the southern hemisphere was considerably weaker than in the northern one between about 1950 and 1975 (cf. Fig. 1), resulting in strongly positive values of both AA and NA. The corresponding effect is further enhanced by the later termination of the solar activity maxima in the southern hemisphere during the corresponding solar cycles, which is consistent with recent findings based on the wavelet phase differences between the Schwabe cycles on both hemispheres (Donner and Thiel, 2007;Donner, 2008). Notably, the latter approach could also be extended to the corresponding cross-wavelet (amplitude) coherency in order to trace temporal changes in the dynamical properties of sunspot areas at both solar hemispheres relative to each other. However, describing these changes by just a single scalar parameter (as done in this study) would require a careful choice of a corresponding specific and sensitive characteristics. In a similar spirit, complementary application of other -traditional as well as modern nonlinearmethods of bivariate time series analysis could provide further insights into the temporal organization of solar activity on various timescales, but is clearly beyond the scope of the present work.
Beyond the results on long-term variations in nonlinear dynamics characteristics provided in this study, we suggest that a combination of information on amplitude and phase relationships with additional results on dynamic complexity has great potential for opening a new view on the north-south asymmetry of solar activity, its temporal organization and dynamics and, hence, its potential causes and underlying physical mechanisms. To this end, there is no fully established theoretical understanding of the origins of this phenomenon as well as its long-term variations. The application of modern dynamical systems-based concepts (like the one proposed in this study) for analyzing observational data -as well as model outputs for the sake of model-data intercomparisoncould be a promising strategy for filling the corresponding gap in our knowledge in the near future.

Conclusions
We have used visibility graphs (VGs) as a new tool for uncovering long-term changes in the asymmetry between nonlinear dynamical characteristics of the variations of sunspot areas on both solar hemispheres. This new viewpoint adds complementary information to our present knowledge of the variability of the north-south asymmetry of solar activity, which has been previously restricted to differences in the amplitudes and phases of various activity indicators. Our results indicate that the hemispheric asymmetry of the considered nonlinear dynamical characteristics has gradually changed over the period of systematic observations of sunspot areas, with two shifts in the hemispheric predominance around 1930 and 1990. The two corresponding transition periods are potentially related to known transitions in the relative magnitude of hemispheric activity, which however clearly precede the identified changes in the dynamical characteristics.
We have demonstrated that the VG properties encode both dynamical characteristics and information associated with the probability distribution function (PDF) of the data under study. Considering both aspects together has been essential for unveiling the aforementioned long-term changes in the dynamics of solar activity, which could not be obtained when studying characteristics exclusively related to dynamical (excess degrees of horizontal VGs) or distributional properties (skewness). In turn, this entanglement between statistical and dynamical aspects makes the attribution of our findings to specific physical processes an even more challenging task. Systematic applications to other solar activity indicators are necessary to fully explore the potential and limitations of the proposed methodological approach. From a conceptual perspective, we stress that there are no generic restrictions to applying VGs and related methods to such data.
Regarding the methodological advances reported in this work, using joint and excess degrees of VGs to disentangle similarities and differences between the dynamical patterns exhibited by two simultaneously observed variables provides a new approach to utilizing the powerful concept of VGs in a bivariate setting. To the best of our knowledge, by now there has only been one related bivariate approach, VG similarity (Ahmadlou and Adeli, 2012), proposed in the literature, which, however, has a distinctively different field of application (the detection of generalized synchronization) than the approach presented here and is algorithmically more complex. We emphasize that studying the asymmetries in nonlinear dynamics could be of broad interest in geosciences beyond the present application to solar physics, e.g., when considering different oscillatory patterns in the Earth's climate or biosphere.
The results of this work underline the potential (but also some methodological limitations) of complex network methods to address problems of nonlinear time series analysis. Notably, the framework of VGs utilized here is only one example of such methods. It will be the subject of future work to investigate if complementary methodological approaches such as recurrence networks can provide comparable or even additional information on dynamical asymmetries between two time series. To this end, our results mainly contribute to a better understanding of the potential and limitations of the specific methodology used here as well as to a more detailed statistical characterization of sunspot data. However, linking these new results to specific processes within the solar interior remains a challenging task.

Appendix A: Skewness effect on global VG properties
In Sect. 3.1, we showed that the edge density of a VG can depend crucially on the skewness of the underlying data. This has important implications for properly interpreting the properties of VGs. Notably, VGs have been previously considered as representing exclusively the temporal correlation properties of a time series. Regarding our results, this view needs to be corrected. In this appendix, we provide further details for various global network characteristics beyond the edge density (for explicit definitions of the properties studied below and a brief discussion of their behavior for VGs, see Donner and Donges, 2012), using the same setting of a gammadistributed white noise process as in Sect. 3.1. Figure A1 demonstrates that the global clustering coefficient C of the VG increases with increasing S (Fig. A1b), whereas transitivity T , average path length L, and assortativity coefficient R decrease (Fig. A1a, c, d). While the aforementioned decrease is continuous for L over the entire considered range of skewness values (S ∈ [−2, 2]), we find some saturation of T for high positive skewness values, and of R for strongly negative skewness values.
The observed decrease in L indicates the rising importance of local maxima for data with strongly positive skewness, which introduce many "shortcuts" in the VG related to its generally marked clustered structure (see Elsner et al., 2009;Donner and Donges, 2012;Zou et al., 2014b for some illustrative examples). Notably, additional edges are introduced as the local maxima become more and more pronounced with rising S (see Fig. 5). These new edges enter existing clusters of vertices and, hence, increase the local clustering coefficients C i within these clusters, resulting in the observed increase in C. Taking both findings together, we have an increase in the level of "small-worldedness" (Watts and Strogatz, 1998) of the VG as S rises. Finally, the decrease in R points to the rising variety of degrees of mutually connected vertices for positive S.
As a particularly remarkable observation, unlike C the network transitivity T decreases for negative S and reaches some approximately stationary value for positive skewness (Fig. A1a). This finding is not a priori expected since both measures quantify closely related characteristics associated with the presence of closed paths of length 3 in the network. Notably, the global (i.e., average local) clustering coefficient C (Watts and Strogatz, 1998) is the unweighted arithmetic mean of the contributions of all vertices (irrespective of their degree), whereas the network transitivity T implicitly gives more weight to vertices with high degrees than to those with low degrees (Barrat and Weigt, 2000). It is also worth noting that in the present example C is always much larger than T . The reason for this difference is that vertices with minimum degree k i = 2 (which are only connected to their direct neighbors) most commonly correspond to local minima of the time series (i.e., x(t i ) < x(t i−1 ), x(t i+1 )) so that the local clustering coefficient must be C i = 1. Since there are many such vertices, they contribute strongly with high values to C, but not to T (Donner and Donges, 2012). In the same spirit, we suggest that the qualitative difference between the behaviors of T and C is related to the specific type of the degree distributions of the considered VGs.
We emphasize that in case of S < 0, the PDFs of our artificial time series exhibit inflection symmetry with respect to the S > 0 case. However, since VGs treat high and low time series values in essentially different ways (Zou et al., 2014b), the resulting network properties do not exhibit a similar symmetry (Figs. 5 and A1), but reveal gradual trends in all considered global network measures across S = 0.
As the numerical example discussed here does not reflect further on more realistic time series properties, a more detailed investigation of skewness effects on VGs resulting from stochastic processes also exhibiting serial correlations is desirable yet beyond the scope of this work. In general, a more systematic inspection of changes at the local network level and their reflectance in global VG characteristics appears to be a valuable research topic for future work.
Notably, the dependence of VG properties on the skewness of the underlying time series is not present anymore when studying the same properties for the HVG.

Appendix B: VG characteristics for sunspot areas
For the sake of completeness, we complement our previous analysis with a more detailed characterization of the obtained VGs obtained for hemispheric sunspot areas in terms of a set of complex network measures. Subsequently, we provide a comparison of our results with the corresponding properties obtained for gamma-distributed white noise (see Sect. 3.1 and Appendix A) with the same skewness as the original data. Notably, the corresponding skewness values of the entire time series are given in Fig. 4 and are close to the upper limit of corresponding range for the artificial data considered earlier in this work.
For the different global network characteristics obtained from the sunspot area data, the computed values are given in Table B1. The corresponding results allow one to draw several general conclusions.
First, the VG properties do not differ much between the two time series representing sunspot areas on the northern and southern solar hemisphere, respectively. This finding implies that the dynamical characteristics captured by the considered complex network measures do not allow for clearly distinguishing between the observed dynamics, which is consistent with the general similarity in the fluctuations of both variables on short as well as long timescales. Notably, the same applies to the properties of the corresponding HVGs.
Second, Table B1 gives an idea of the effect of the asymmetry of the sunspot area data. Specifically, VG and HVG treat local maxima and minima of a time series by −1 in essentially different ways. Hence, multiplying a time series can provide distinct results when minima and maxima are arranged in essentially different ways, as it is the case for the considered solar activity time series (Zou et al., 2014b). Notably, for the negative series, the edge densities of the VGs are somewhat lower than for the original ones, along with a higher average path length. It is worth noting that this effect originates mainly from the positive skewness since we find the reverse behavior for the HVG that is not affected by the specific PDF of the data. This finding could be explained by the relatively short durations of the phases of almost complete solar inactivity, which lead to relatively sharp maxima of the negative series. The values of the global clustering coefficient are of about the same order for the original and negative series (being slightly lower for the negative series for the VG, but higher), whereas those of network transitivity are far larger than for the original data for both VG and HVG. In a similar way, we find that whereas the original time series exhibit only very weak positive assortativity (indicating thorough mixing of connections between vertices of different degree), the negative series exhibit extremely strong assortativity, i.e., highdegree vertices are preferentially linked with other highdegree vertices, etc. This could be explained by the fact that the phases of solar inactivity (i.e., maxima of the negative series) display relatively smooth fluctuations.
Third, Table B2 shows the values of VG properties for realizations of a gamma-distributed random variable with the same skewness and length as the sunspot area data under study. A comparison with Table B1 reveals that the sunspot area time series indeed exhibits higher VG edge density and lower average path length. Transitivity and global clustering coefficient are of about the same order for sunspot data and noise, with T being slightly higher and C slightly lower for the observational data. In turn, assortativity is clearly larger for the gamma-distributed white noise, which indicates that this property is more closely related to the actual dynamics than T and C.
Notably, the differences between the VG properties of the original and negative series for the gamma-distributed white noise differ from those found for the sunspot areas but are qualitatively consistent with the results in Figs. 5 and A1. The latter is to be expected since studying the negative white noise series corresponds to considering a negative skewness with the same absolute value as before. However, the explicit values of the obtained network characteristics partly differ from those obtained for higher values of T due to the extensivity of some of the measures discussed above.