Spatial Association Network Evolution and Variance Decomposition of Economic Sustainability Development Efficiency in China

The economy’s sustainable development has become a national strategic deployment in China. Research on the difference between the economic sustainable development efficiency (ESDE) and the spatial network will assist the government with the deployment of sustainable development strategies and the achievement of the “peak carbon dioxide emissions”. This paper designs the input–output indicator system of sustainable economic development efficiency and builds an unexpected output super-EBM-Malmquist model to measure the ESDE of 30 provinces in China from 2008–2020. According to the ranking of ESDE, the 30 provinces in China are classified into four groups by applying the quartile method, and the difference in the ESDE in different regions and the temporal variation of different provinces are studied by using the Dagum Gini coefficient and Gaussian Kernel density. Moreover, the relationship between ESDE in different provinces is studied based on the revised gravity model and social network analysis method. The connections between provinces with related relations constitute the ESDE network. Results show that (1) the average ESDE in China shows an upward trend, the eastern region is in a leading position, the central and western regions are trying to catch up with the eastern region, and the development of the northeast region is lagging behind. (2) The level of ESDE in different provinces is clearly arranged from high to low, illuminating a distinct pattern. Moreover, provinces with high levels of development are much higher than provinces with low levels of development, presenting a phenomenon of polarization. (3) The regional ESDE development imbalance is prominent, and the ESDE in the eastern region is closely related, while the connection in the western region is lower. (4) Beijing–Tianjin Urban Agglomeration and the Yangtze River Delta have significant spatial spillover effects in the association network, while the northeast, northwest, southwest and central regions have significant spatial benefit relationships. These findings provide important enlightenment for promoting the sustainable and balanced development of China’s economy.


Introduction
With the growing energy consumption, energy demand and energy use to promote economic growth at the expense of environmental sustainability [1,2], sustainable economic development is facing a huge threat. As the world's second largest economy, China's energy consumption structure has long been dominated by fossil energy, its economic model is characterized by significantly high consumption and high emissions, and its rapid economic growth is accompanied by environmental pollution and ecological degradation [3].
China's economic development has been maintaining the characteristics of crude growth with high input and high consumption and low output, and the energy structure dominated by coal and high overall carbon emissions is the reality of China's national 2 of 22 situation. According to BP statistics, coal will account for 56% of total energy consumption in 2021, and China has become the world's largest emitter of carbon dioxide. Obviously, this energy structure is unsustainable and will cause irreversible damage to the ecological environment [4]. The high consumption of primary energy, especially the large-scale exploitation of fossil energy, has led to an excessive amount of carbon dioxide emissions and increasingly serious environmental impacts, which is contrary to the principles of sustainable development [5][6][7]. Sustainable economic development is becoming a major concern for governments [8], especially after a significant increase in the exploitation of natural resources that has led to serious environmental problems [9]. In the 1970s, the economist Meadows argued in The Limits to Growth that technological progress plays a crucial role in the harmonization of economic development with resources and the environment. In the context of tightening resources and environmental constraints, "carbon peaking", "carbon neutrality" and "emission reduction", rapid economic development cannot be achieved at the expense of resources and the environment. Therefore, the sustainable development of China's economy necessarily needs to be discussed within the framework of environmental protection [10], and the study of economic sustainability under resource and environmental constraints has gradually become the focus of academic attention.
The Chinese government gives great importance to sustainable economic development. How to reconcile economic development with the environment and how to enhance the level of sustainable economic development in the future are key issues for sustainable economic development [11,12]. In 2015, China proposed supply-side structural reform, which aims to restructure the economy, achieve an optimal allocation of factors, improve the quality of economic growth and promote sustainable economic development. In the report of the 19th National Congress, General Secretary Xi Jinping emphasized the development concept of "green water and golden mountains are golden mountains and silver mountains", highlighting the importance that the government attaches to the organic integration of economic development and ecological protection. In the 14th Five-Year Plan of China, the importance of sustainable economic development was re-emphasized. In the framework of environmental protection and energy constraints, it is therefore of great theoretical and practical importance to discuss the sustainable development of China's economy.
This study measures the ESDE in China with the aim of identifying gaps in the ESDE across different regions and provinces while also examining the spatial network structure of sustainable economic development to help the government implement environmental protection, conservation and economic sustainable development strategies, and to raise China's level of sustainable development. This study designs a research framework model of "efficiency measurement-differentiation analysis-network structure". Gini coefficient and Kernel density estimation are two tools for the difference measurement of an object among different individuals, and they also provide a good visualization. In particular, the Gini coefficient method is able to calculate the development gap among the eastern, western and central regions of China, and can decompose the existing gap to distinguish that the difference in the sustainable development level of the national economy mainly comes from the gap among the three major regions or the gap within the three major regions. Moreover, Kernel density estimation is also a method to analyze differences and gaps in economic sustainability and result visualization. The social network analysis method is a method that can visualize the flow relationship of sustainable economic development among different provinces. It clarifies the closeness of the level of sustainable economic development in different regions, thus providing information for governments to take further actions that are tailored to local conditions and better promote China's sustainable economic development. Furthermore, in addition to estimating the ESDE, this study applies ESDE outcomes to further measure the regional development gap, the regional development evolution pattern, and the spatial network structure.

Defining Sustainable Economic Development
Sustainable development is one of the most common words used today and represents an expression to stop environmental degradation, which can lead to ongoing climate change and changes in citizens' lifestyles and habits [13]. The definition, connotation and characteristics of sustainable development have been defined in the economics literature for a long time, and although there is no unified and clear theoretical framework, it has had a significant impact [14]. Stivers [15] defines the sustainability of economic growth as "an economy in equilibrium with its ecological support system". The concept of "sustainable development" was first introduced in 1987 in the report of the United Nations World Commission on Environment and Development on the Future of Humanity: Our Common Future. Pointing to the need to change the development paradigm for the benefit of present and future generations, the report defines sustainable development as development that meets the needs of the present without compromising the ability of future generations to meet their needs [16], and this definition is widely used by international organizations and relevant scholars. Pezzey [17] argues that the core meaning of sustainable development is a development path that keeps the welfare of the population from declining, and Dailami, et al. [18] believes that sustainable economic growth can be an important element of the pace of economic development; Barro [19] defines sustainable economic growth as encompassing characteristics such as gross domestic product, the political system, income distribution, health status and religious beliefs. Brock and Taylor [20] proposes that sustainable economic growth is primarily about economic growth and environmental quality improvement. In 2005, the World Summit for Social Development adopted the '3 E's' of economic, equity and environmental as the core elements of sustainable development. The UN Sustainable Development Goals (SDGs) aim to address the social, economic and environmental dimensions of development in an integrated manner between 2015 and 2030, and to shift towards a sustainable development path. Many other scholars have defined and discussed various aspects of sustainable development, including its definition, connotation and applicability, from the perspectives of ecology and environmental economics [21,22].

Factors Influencing Sustainable Economic Development
Scholars have studied the different influencing factors of economic sustainability because of its rich connotations. Wang and Huang [23] found through bibliometric research methods that research on sustainable development in the context of COVID-19 is broad in scope and covers a wide range of disciplines, but it is lacking in depth overall. While developed countries are working on the sustainability of education, developing countries are more concerned with economic sustainability.
Khan, et al. [24] used the Environmental Kuznets Curve (EKC) to explore the non-linear link between energy intensity, financial development and environmental sustainability, showing a non-linear inverted 'U' shaped relationship between financial development and environmental sustainability. Bai, et al.'s [25] study of China's economic sustainability goals is based on a panel threshold model. The study finds that investing in renewable energy reduces greenhouse gas emissions, thereby maintaining sustainable economic growth and contributing to the achievement of the SDGs. Dabbous and Tarhini [26] applied a fixed effects panel model to find that a sharing economy contributes to sustainable economic development. Hosan, et al. [27] used a panel model to analyze the links between the demographic dividend, digital innovation, energy intensity and sustainable economic growth in 30 emerging economies from 1995 to 2018, finding that the demographic dividend and digitalization contributed to sustainable economic growth, while energy intensity was negatively associated with sustainable economic growth.

Economic Sustainable Development Efficiency (ESDE)
Economic sustainability is a multidisciplinary subject of study in all the countries of the world, and it is critical to any country's economic development and environmen-tal protection to assess a country's sustainable economic development performance and outcomes [28]. In terms of the selection of models for measuring efficiency, a number of scholars have chosen traditional radial DEA models or non-radial SBM models. Studies mostly quantify GDP as the desired output and environmental pollutants as the nondesired output, and establish a system of input-output indicators. DEA (Data Envelopment Analysis) models that can handle multiple inputs and outputs are the common tools used. Qu, et al. [29], constructed an improved super-efficient radial network DEA model to assess China's regional sustainable development performance from 2009 to 2017; the study shows that there are significant differences in the level of sustainable development in different regions, with the eastern region having the highest level of sustainable development and the central and western regions being in a steady catch-up stage of development. Many scholars consider undesired output as an output variable and construct a non-radial SBM model (Slacks-Based Measure), which can handle efficiency measurements with undesired output variables. In recent years, the SBM model has been widely used in the evaluation of green economic efficiency [30,31]. Wen, et al. [32] used a three-stage SBM model to measure the sustainable development efficiency of counties in the Yellow River Basin from 2005 to 2015, and they classified the efficiency of different regions into low, medium-low, mediumhigh and high-efficiency through quadratic plots, revealing the spatial clustering pattern of sustainable development. Swain and Ranganathan (2021) [33] employed a network perspective to analyze IAEG-SDG data on SDG indicators from the UN Sustainable Development dataset for the period 2000-2017, finding that the same criteria cannot be set for sustainable development goals in different regions and that different priority development goals should be set for different SDG groups. Wang and Jia [34] applied the DEA-SBM super efficiency model to characterize sustainable economic growth, demonstrating that improving the energy consumption mix and controlling CO 2 emissions both contribute to sustainable economic growth.
However, neither the traditional radial DEA model nor the non-radial SBM model is capable of dealing with the situation where the input and output variables have both radial and non-radial characteristics, and the information on the ratio between the target and actual values of inputs or outputs is easily lost, resulting in an overestimation or underestimation of the efficiency values [35]. EBM has been proposed as a hybrid distance model to overcome the shortcomings of the DEA and SBM models and to improve the accuracy of the evaluation results [36]. However, there are not many cases where the EBM model has been applied to evaluate the ESDE.

Spatial Variation and Social Network Analysis
In research fields as efficiency evaluation, scholars are concerned not only with overall efficiency, but also with the analysis of the differences in efficiency and the spatial structure of networks [37,38]. The analysis of spatial differences is one of the key issues in geographic studies, allowing for a precise portrayal of the spatial differences in the study population, and thus the implementation of policies corresponding to regional differences [39]. Currently, a number of studies have adopted methods such as the Gini coefficient [40,41], spatial exploratory data analysis [42] and Kernel density functions [43] to analyze spatial differences under different spaces. With the advancement of China's coordinated regional development strategy and sustainable development strategy, the spatial linkages between regions are increasingly taking the form of complex networks [44]. It is hardly enough to analyze the ESDE by relying only on the variance analysis of "attribute data". Economic efficiency issues are often characterized by the evolution of spatially linked network structures, which requires the use of social network analysis [38,44].
Social network analysis is a method based on relational data to study the complexity of the network structure between unit nodes, effectively overcoming the shortcomings of "attribute data", and has been widely used in recent years in the study of complex correlation networks in the fields of regional energy and environment [45,46]. Scholars have combined variance analysis and social network analysis to analyze the dynamic evolution of spatial differences in efficiency issues. Shen, et al. [47] applied the Theil index and Moran indices to assess the characteristics of carbon emissions in the Yangtze River Delta urban agglomeration, the Chengdu-Chongqing urban agglomeration and the Guangdong-Hong Kong-Macao urban agglomeration, combined with the social network analysis method to examine the network structure characteristics of the three urban agglomerations from three perspectives: overall, individual and correlation. Based on the impact of the spatial correlation of carbon emissions on sustainable economic development, Zhang, Tai, Cheng, Zhu and Hou [46] constructed a carbon emission efficiency network and classified carbon emission efficiency into five spatial classes: low efficiency class, low efficiency class, medium efficiency class, high efficiency class and higher efficiency level and analyzed their spatial differences and dynamic evolution. However, there is less literature on the application of the Gini coefficient, Kernel density estimation and social network analysis as holistic research methods to the research topic of economic sustainable development efficiency. Although a few studies have also verified and discussed the spatial agglomeration characteristics and spillover effects of sustainable development based on geographical proximity or economic proximity, the current research on the structural characteristics of spatially linked networks of sustainable economic development efficiency and their spatial dynamic differences still needs to be deepened.
The contribution of this study lies in several areas. Firstly, this research extends the study of sustainable economic development to cover the level of efficient use. We take 30 provinces in China as the study objects and classify the sustainable development efficiency levels into very high efficiency, high efficiency, medium efficiency and low efficiency in order to reveal the regional spatial and temporal characteristics and the mechanisms of differences, and to better analyze the spatial evolution characteristics of China's sustainable economic development efficiency and help to achieve sustainable and efficient economic development in China. Secondly, this study provides an accurate estimate for the measurement of sustainable economic development efficiency by using an over-efficiency EBM model of non-expected output, which takes into account the nonexpected output of sustainable economic development. Thirdly, this study uses the Gini coefficient, Kernel density estimation and social network analysis as a framework model for the study of spatial dynamic differences to analyze the spatial structure and evolutionary characteristics of the efficiency of sustainable economic development from the perspective of spatial correlations.

Index System and Data Sources
The core connotation of sustainable economic development is to organically integrate the environment, resources and development issues. It not only pays attention to the rapid economic growth, but also needs to pay attention to the quality of economic development. We emphasize that sustainable economic development cannot degrade the quality of the environment or harm natural resources at the cost of being sacrificed. Sustainable economic development ensures maximum economic development on the premise of maintaining the quality of natural resources and maintaining environmental stability.
Most of the literature is an economic sustainable growth index system that incorporates labor, capital, technology and energy consumption into the input indicators of production factors [48] but does not take into account the input of environmental protection and social livelihood. However, the ESDE is based on the connotation of sustainable economic development, which needs to fully consider the comprehensive economic efficiency after the cost of the environment and resources. This paper draws on the existing evaluation system based on the principles of system, science and data availability, considering the constraints of resources and the environment. On this basis, this paper designs six dimensions of input indicators, including environment, resources, social livelihood, technology, capital and labor, as well as output indicators such as expected output and undesired output, to construct an evaluation index system for the ESDE.
Considering data availability and continuity, this paper selects 30 provincial-level administrative regions in China (excluding Hong Kong, Macao, Taiwan and Tibet) as the research unit, and takes 2000-2020 as the research period. The data come from the National Bureau of Statistics, "China Statistical Yearbook", "China Population and Employment Statistical Yearbook", "China Energy Statistical Yearbook", "China Environmental Statistical Yearbook", "China Science and Technology Statistical Yearbook" and provincial statistical yearbooks. The specific index system is shown in Table 1. Among them, the capital input index is represented by the stock of fixed asset investments. Because the EBM model limits the number of indicators, this paper uses the entropy method to calculate the weight of multiple secondary indicator variables. The entropy weight method is an objective weighting method for indicators, which can deeply reflect the distinguishing ability of indicators. Compared to the subjective weighting method, the outcome of entropy weighting method is more reliable. Finally, input indicators such as capital, environmental governance, social security, technology, resources and labor are obtained, with GDP as the final expected output and COD and SO 2 as undesired outputs. The specific indicators are shown in Table 1.

Super-EBM Model
The effective unit efficiency value measured by the EBM model is 1; it is difficult to further analyze the efficiency difference of the effective evaluation unit, and the efficiency of different evaluation units cannot be judged. Andersen and Petersen [49] put forward the idea of super-efficiency planning borrowed from the design idea of Tone's [50] super-SBM model, which optimized the constraints of the EBM model with no guidance and unexpected output, and the efficiency measurement result was equal to 1 unit for secondary calculation. The formula of the super-efficient EBM model is as follows [51].
The constraints on input variables are The constraints on the expected output factors are The constraints on the unexpected output factors are The objective function is The constraints of the super-EBM model are where x ij is the data matrix of input factor indicators, and y rj is the data matrix of output factor indicators. Variable γ * is the best efficiency with constant returns; s + r is the slack variable of expected output of category r; s b− p is the slack variable of the unexpected output of type p; w + r and w b− p are the weights of expected output r and unexpected output p, respectively; ε y and ε b are parameters; ϕ is the output expansion ratio; b pj is the unexpected output p of the decision-making unit j; b pk is the unexpected output p of the decisionmaking unit k; and q is the total number of unexpected output p.
Based on the efficiency index system of ESDE in China, this paper selects 6 types of input indicators, 1 item of expected output and 2 types of unexpected output in the nonguided, constant returns and unexpected output super-efficiency EBM model to measure the efficiency of China 's sustainable economic development.

EBM-Malmquist Model
The above super-EBM model can only measure the static efficiency value of China's sustainable economic development. In order to more comprehensively reflect the changes in the efficiency, the Malmquist index can be calculated to reflect the efficiency change rate of each decision-making unit. The Malmquist index was first proposed by Malmquist, and the Malmquist index uses the distance function to reflect the rate of change of each decision-making unit [52]. This paper uses the Malmquist index to reflect the changes in the ESDE of 30 provinces in China from 2008 to 2020. The mathematical expressions of the previous period t − 1 and the current period t of the Malmquist index are shown in Formula (6) [53,54]: where x t−1 and y t−1 represent the input vector and output vector of the previous period t − 1, and x t and y t represent the input vector and output vector of the current period t.
Variable D t (x t , y t ) represents the distance between the decision-making unit in period t and the production frontier in period t − 1, and the efficiency of a decision-making unit in period t is measured by constructing the production frontier with all decision-making units in period t − 1. The value range of the Malmquist index is MI t−1,t ∈ (0, +∞). Efficiency change of sustainable economic development in China depends on comparison with 1. If MI t−1,t > 1, it means that the efficiency of the current period is improved. If MI t−1,t < 1, it means that the efficiency of the current period has declined. If MI t−1,t = 1, it means that the efficiency of the current period remains unchanged.

Modified Gravity Model
The theory of economic gravity believes that there is an interactive relationship between regions within a certain range, and the construction of a spatial correlation matrix is a prerequisite for network analysis. Many scholars use a revised gravity model to measure the spatial correlation of research objects [38,55,56]. However, the connection between economies is different and unidirectional. When calculating the degree of connection between regional economies, it is necessary to take into account the economic scale, population and other related factors between economic regions.
The gravity model can be used to calculate the connection strength of the 30 provinces in China. The purpose of calculating the gravity matrix is to find the spatial correlation matrix of the ESDE. The gravity model is the basis of social network analysis, and the spatial correlation matrix in social network analysis can be converted from the correlation strength calculated by the gravity model. This paper draws on the approach of Zhang and Lu [57], and uses the following revised gravity model to calculate the spatial correlation strength among the 30 provinces in China.
where C ij represents the gravitational force between province i and province j, which represents the correlation strength of the ESDE; K ij represents the gravitational coefficient between province i and province j, representing the contribution of province i to the connection between province i and province j; E i and E j represent the ESDE of province i and province j; G i and G j represent the economic development level of province i and province j, measured by GDP; P i and P j represent the population of province i and province j; D ij represents the spherical distance between the capital city of province i and province j; g i and g j represent the per capita GDP of province i and province j; g i − g j represents the economic distance of province i and province j. Using the gravity model of Formula (7), the gravity matrix of the sustainable development efficiency of an inter-provincial economy can be constructed. The average value of the row values of the matrix C ij is used as the threshold value. If the gravity value is greater than the threshold, the value is 1, which means that the province in this row has an efficiency space overflow for the province in the column. Conversely, if it is lower than the threshold, the value is 0, indicating that there is no correlation effect between the two provinces. After the above conversion, the spatial correlation matrix I of ESDE can be obtained, which is actually a (0, 1)-matrix.

Social Network Analysis
Social Network Analysis (SNA) is an analysis method based on the perspective of "relationship", including sociology, economics, geography and other disciplines [58]. It is mainly used to analyze the relationship structure and attributes of social networks, providing an "interactive" perspective and a global analysis by describing the relational schema with graph theory tools and algebraic modeling techniques [57]. Most multivariate statistical methods cannot be used to analyze relational data, and social network analysis can just make up for this limitation. The advantage of the social network analysis method is that it can accurately quantify various associations and avoid the limitations of "adjacent" or "similar" traditional spatial econometric analysis methods, thus providing a basis for the construction of some middle-level theory, and the test of empirical propositions provide quantitative tools. It is even possible to build a bridge between "macro and micro" (Liu and Jia, 2019). Therefore, this is also one of the important reasons why many experts and scholars apply social network analysis to economics, management and other fields.
Based on the spatial correlation matrix I calculated by the revised gravity model, this paper constructs a spatial correlation network based on the ESDE. Social network analysis is used to analyze the overall network and individual network characteristics of China's sustainable economic development. Among them, the overall network characteristics are mainly quantified by four indicators: the network relationship coefficient, network density, network hierarchy, and network efficiency. Individual network characteristics use point-in-degree and point-out-degree to analyze the individual spillover effects and benefit relationship.
Social network analysis is one of the main contents of this paper. Based on the calculation of ESDE, the social network analysis method is used to analyze the crossregional flow relationship of ESDE in 30 provinces in China. The connection of regions with such connections will form a geospatial network, and social network methods can visualize the strength of this connection.  (4) and (5) ( Table 2). As shown in Table 2, 40% of China's provinces had a ESDE greater than 1 in 2008; the proportion reached 56.67% in 2020, indicating that the overall level of sustainable economic development in China is on the rise.

Result Analysis
As shown on the right in Figure 1, the average of China's ESDE showed an "M"-shaped fluctuation trend during 2008-2020. However, it showed an overall upward trend, with the average efficiency rising from 0.8354 in 2008 to 0.8818 in 2020. As shown on the left in Figure 1, it can be clearly seen that the ESDE of Beijing, Shanghai, Jiangsu, Zhejiang, Fujian, Guangdong, and Hainan is at a high level. The ESDE in Liaoning, Jilin, Heilongjiang, Shaanxi, Gansu, Qinghai, Ningxia and other provinces is at a low level. As shown in Table 2, 40% of China's provinces had a ESDE greater than 1 in 2008; the proportion reached 56.67% in 2020, indicating that the overall level of sustainable economic development in China is on the rise.
As shown on the right in Figure 1, the average of China's ESDE showed an "M"shaped fluctuation trend during 2008-2020. However, it showed an overall upward trend, with the average efficiency rising from 0.8354 in 2008 to 0.8818 in 2020. As shown on the left in Figure 1, it can be clearly seen that the ESDE of Beijing, Shanghai, Jiangsu, Zhejiang, Fujian, Guangdong, and Hainan is at a high level. The ESDE in Liaoning, Jilin, Heilongjiang, Shaanxi, Gansu, Qinghai, Ningxia and other provinces is at a low level.

Analysis of Malmquist Index
Based on the efficiency results of the non-oriented, constant returns to scale, and unexpected output super-efficiency EBM-Malmquist model in the 30 provinces in China from 2008 to 2020, the MI index is decomposed into the EC index and the TC index. According to the Malmquist index, this paper analyzes the change trend of the ESDE of the 30 provinces (Table 3). The MI index of each province in China fluctuated in different years. It can be seen from Table 3  4.2. Analysis of Regional Differences and Dynamic Evolution 4.2.1. Estimation and Decomposition of Regional Differences in ESDE In order to describe the level of regional differences and sources of regional differences in ESDE, this paper uses the Gini coefficient and its subgroup decomposition method proposed by Dagum (1997) to decompose the ESDE from 2008 to 2020. The Gini coefficient decomposition method solves the problem wherein the traditional Gini coefficient index does not satisfy the decomposability of subgroups-not only by taking the distribution of subgroup samples into account, but also by solving the problem of overlapping between sample data and the source of regional overall differences. The new Gini coefficient decomposition method can decompose the overall difference, see the source of the difference more clearly, and overcome the shortcomings of the traditional Gini coefficient and Theil coefficient (Liu et al., 2012;Li and Zhang, 2018). The formula of the Gini coefficient decomposition method is as follows: where k is the number of regions divided by China, k = 3, which are the eastern region, central region and western region, respectively. Variables j and h are different regions in k regions, where j = 1, · · · , k, h = 1, · · · , k, j = h. Variable n is the number of provinces in the research sample, n = 30. Variable y ji (y hr ) is the economic sustainable development efficiency of province i(r) in j(h) region, and y is the average value of the ESDE. The overall Gini coefficient can be decomposed into three parts: intra-regional gap contribution G w , inter-regional gap contribution G nb , and hypervariable density G t . According to Formulas (8) and (9), the Gini coefficient and decomposition results of China's ESDE are calculated, as shown in Table 4.  (8) and (9) based on the data in Table 2.
It can be seen from Table 4 and Figure 3 that the overall Gini coefficient of China's ESDE increased by 4% in 2009, showed a slight increase, and declined from 2010 to 2012, with an average annual decline rate of 9%. After that, it rose in 2013, and then fluctuated slightly. During 2016-2017, it showed a downward trend, with a rate of decline of 10.41%, and reached its minimum value in 2017, and then showed an upward trend with a growth rate of 7.35% in 2018-2020, indicating that the overall regional differences in China's ESDE fluctuate and have shown an expanding trend in recent years.  Overall and intra-regional SDED differences; (b) inter-regional SDED differences. Data source: Table 4.
As shown on the right in Figure 3, China's ESDE is uneven among regions, and the regional differences have shown a tendency to expand after 2017. The regional differences between the eastern and central regions show periodic characteristics, showing a "V"shaped fluctuation trend during 2009-2015 and a weak "W"-shaped fluctuation trend during 2015-2020. In the sample, there is a "W"-shaped fluctuation trend. After 2009, the difference between the eastern region and the western region generally showed a fluctuating trend of "W", showing a trend of "shrink-expand-shrink-expand". After 2009, the regional differences between the central and western regions generally showed a "W" fluctuation trend, which was similar to the fluctuation trend of the differences between the eastern and western regions. Figure 4 reflects the contribution of inter-regional differences, hypervariable density and intra-regional differences to the overall differences. Prior to 2010, inter-regional differences contributed the most, and hypervariable density and intra-regional differences contributed almost equally to overall differences. After 2010, the contribution of the hypervariable density is the largest, and the hypervariable density is the main source of the overall difference. The contribution rate of the intra-regional difference remains basically unchanged, and the intra-regional contribution rate tends to expand. The contribution of intra-regional differences is significantly greater than that of inter-regional differences. Therefore, the key to solving the problem of unbalanced efficiency in China's ESDE is to reduce the efficiency difference within the region. Figure 3. (a) Overall and intra-regional SDED differences; (b) inter-regional SDED differences. Data source: Table 4.
From the perspective of evolution, the differences in the ESDEs of the eastern region reached the minimum in 2009, and it showed an expanding trend during the period from 2010 to 2014, with an average annual growth rate of 13.93%, and reached the maximum in 2014. After 2015, regional differences showed a "W"-shaped fluctuation. The differences in the ESDEs of the central region showed a trend of substantial expansion during the period from 2008 to 2010 with a growth rate of 19.80%, and a decline rate of 23.53% during the period from 2010 to 2012. It reached the minimum in 2012, and then showed a fluctuating upward trend. During 2015-2016, there was a relatively large decline, with a decline rate of 14.47%. Then, the regional differences continued to widen during 2017-2020, with a growth rate of 12.59%. The trend of differences in the ESDEs of the western region during 2008-2015 is similar to the overall trend, with an average annual growth rate of 6.55% from 2013 to 2016, and fluctuates after 2016. In general, the difference between the whole area and the central region's ESDE tends to expand after 2011, while the eastern and western regions show a fluctuating trend.
As shown on the right in Figure 3, China's ESDE is uneven among regions, and the regional differences have shown a tendency to expand after 2017. The regional differences between the eastern and central regions show periodic characteristics, showing a "V"shaped fluctuation trend during 2009-2015 and a weak "W"-shaped fluctuation trend during 2015-2020. In the sample, there is a "W"-shaped fluctuation trend. After 2009, the difference between the eastern region and the western region generally showed a fluctuating trend of "W", showing a trend of "shrink-expand-shrink-expand". After 2009, the regional differences between the central and western regions generally showed a "W" fluctuation trend, which was similar to the fluctuation trend of the differences between the eastern and western regions. Figure 4 reflects the contribution of inter-regional differences, hypervariable density and intra-regional differences to the overall differences. Prior to 2010, inter-regional differences contributed the most, and hypervariable density and intra-regional differences contributed almost equally to overall differences. After 2010, the contribution of the hy-pervariable density is the largest, and the hypervariable density is the main source of the overall difference. The contribution rate of the intra-regional difference remains basically unchanged, and the intra-regional contribution rate tends to expand. The contribution of intra-regional differences is significantly greater than that of inter-regional differences. Therefore, the key to solving the problem of unbalanced efficiency in China's ESDE is to reduce the efficiency difference within the region.

Dynamic Evolution of ESDE Regional Variation Distribution Based on 3D Kernel Density
Kernel density estimation is used to fit sample data by a smoothed peak function that describes the distribution pattern of random variables using a continuous density curve; it has the properties of robustness and weak model dependence. It has been widely used in dynamic evolution analysis.
To describe the distributional features and dynamic evolution of the regional ESDE in China, we use Gaussian Kernel functions to estimate the density distribution pattern of the ESDE in China and to assess the dynamic evolution of the distribution via time dimension. Figure 5a-d illustrates the 3D Kernel density curves for China and the three major regions of eastern, central and western China, respectively.

Dynamic Evolution of ESDE Regional Variation Distribution Based on 3D Kernel Density
Kernel density estimation is used to fit sample data by a smoothed peak function that describes the distribution pattern of random variables using a continuous density curve; it has the properties of robustness and weak model dependence. It has been widely used in dynamic evolution analysis.
To describe the distributional features and dynamic evolution of the regional ESDE in China, we use Gaussian Kernel functions to estimate the density distribution pattern of the ESDE in China and to assess the dynamic evolution of the distribution via time dimension. Figure 5a-d illustrates the 3D Kernel density curves for China and the three major regions of eastern, central and western China, respectively. Figure 5a indicates that between 2008 and 2020, the peak of the national regional ESDE shifts to the right and to then left, indicating that the overall ESDE shows a trend of a reversed U shape, which is consistent with the trend of the ESDE mean curve in previous discussions. During the sampling period, the overall curve shows a bimodal peak with a low side peak, indicating a weak polarization of ESDE. Years 2016, 2017 and 2018 appear to have higher peaks, indicating a convergence of ESDE towards high levels of development in these three years. The Kernel density estimation curves for China have an apparent bimodal pattern, with the side peaks appearing on the left side, pointing to a gradient effect and a polarization of ESDE within China.
As shown in Figure 5b, the bimodal pattern of the Kernel density curve of the ESDE in the eastern region is insignificant in the period from 2008 to 2020 and soon vanishes, evolving into a single-peaked pattern, indicating a weak polarization in the eastern region. The main peak of the ESDE in the eastern region shows a trend of "right shift-left shiftright shift-left shift", indicative of a trend of "increase-decrease-increase-decrease" in the overall ESDE. The ESDE in the eastern region presents a fluctuating status and is at an overall high efficiency level. Meanwhile, the 2D density curve widens, showing a trend of widening differences between the eastern regions.
it has the properties of robustness and weak model dependence. It has been widely used in dynamic evolution analysis.
To describe the distributional features and dynamic evolution of the regional ESDE in China, we use Gaussian Kernel functions to estimate the density distribution pattern of the ESDE in China and to assess the dynamic evolution of the distribution via time dimension. Figure 5a-d illustrates the 3D Kernel density curves for China and the three major regions of eastern, central and western China, respectively.  Table 2. Figure 5c shows that the overall curve of ESDE in the central region shows a singlepeaked pattern during the period from 2008 to 2020, indicating that there is no polarization. The curve shows a wide variation in the magnitude of the peak, with a trend of "high-lowhigh-low" and a trend of "right shift-left shift-right shift-left shift" in the position of the main peak, yet the variation is greater than in the eastern region, indicating that the overall ESDE follows an evolutionary trend of "increase-decrease-increase-decrease". Within the sampling period, high peaks were observed in 2008, 2009 and 2018, indicating a clustering of ESDE towards high levels in these three years, which suggests an overall fluctuating trend in the level of sustainable economic development in the eastern region. Furthermore, the overall ESDE in the central region is significantly disparate, with a fluctuating level of development.
As shown in Figure 5d, the Kernel density curve in the western region appears to have multiple peaks, and the ESDE in the western region has a gradient effect, with multipolarity in some years. The main peak of the Kernel density curve in the western region shows a trend of "right shift-left shift-right shift-left shift", indicating that the overall ESDE in the western region has an "increasing-decreasing-increasing-decreasing" trend. The Kernel density function center and curve in 2012 show a distinct rightward shifting trend, demonstrating a substantial increase in ESDE. Meanwhile, the curves in 2017 and 2018 show high peaks and high levels of average efficiency. Overall, ESDE has improved and grown significantly in recent years in the western region.

Overall Network Characteristics Analysis
In order to study the relationship among the level of sustainable economic development of the 30 provinces in China and the closeness of the connection among regions, this paper uses the social network analysis method to visualize the flow of sustainable economic development between different regions. This is convenient for targeting corresponding policies that should be implemented in accordance with the degree of regional correlation, thereby enhancing the level of sustainable economic development of the region through the flow between different provinces.
The ESDE spatial correlation matrix is constructed based on the modified gravity model, and the spatial correlation network of China's regional economic growth in 2008, 2012, 2016 and 2018 is visualized and mapped with ArcGIS ( Figure 6).  Table 2.
Non-neighboring provinces break the traditional geographic constraints and produce complex cross-regional connections, with inter-provincial interaction and transmission existing through both proximity and jumping paths, resulting in a spatial structure characterized by polar nucleus fission and with multiple inter-provincial connections intertwined in a typical "spider web" structure. The average number of network associations between 2008 and 2020 reaches 191, and the network structure remains comparatively stable. Among the spatially associated networks, Beijing-Tianjin, the Yangtze River Delta region and Guangdong Province have the densest network associations, while the northeast, northwest and southwest regions have sparser network associations; the spatial association effect is higher in the eastern region than in the western region overall. The spatial association network remains relatively stable, with Beijing and Tianjin, the Yangtze River Delta and Guangzhou as the core circles, while the number of spatial associations in the central and western provinces is low, and the role of these provinces in the spatial association network is weak, showing an apparent "dense in the east and sparse in the west" pattern of regional differentiation. Ucinet 6.0 software was used to produce the values of four indicators: number of association network, network density, hierarchy and network efficiency (Table 5) and to draw the plots (Figure 7). Table 5. Overall structural features of regional economic growth in China.  Table 2.

Year
Non-neighboring provinces break the traditional geographic constraints and produce complex cross-regional connections, with inter-provincial interaction and transmission existing through both proximity and jumping paths, resulting in a spatial structure characterized by polar nucleus fission and with multiple inter-provincial connections intertwined in a typical "spider web" structure. The average number of network associations between 2008 and 2020 reaches 191, and the network structure remains comparatively stable. Among the spatially associated networks, Beijing-Tianjin, the Yangtze River Delta region and Guangdong Province have the densest network associations, while the northeast, northwest and southwest regions have sparser network associations; the spatial association effect is higher in the eastern region than in the western region overall. The spatial association network remains relatively stable, with Beijing and Tianjin, the Yangtze River Delta and Guangzhou as the core circles, while the number of spatial associations in the central and western provinces is low, and the role of these provinces in the spatial association network is weak, showing an apparent "dense in the east and sparse in the west" pattern of regional differentiation. Ucinet 6.0 software was used to produce the values of four indicators: number of association network, network density, hierarchy and network efficiency (Table 5) and to draw the plots (Figure 7).  Figure 7 and Table 5 show that the number of inter-provincial ESDE connections in China follows a "rising-fluctuating-declining" trend, network density and network efficiency remain relatively stable, network hierarchies show a "stable-declining-stable" trend, and network efficiency declines in general, indicating that the stability of the overall network structure has enhanced. Specifically, the number of network connections rose from 177 in 2008 to 199 in 2014, reaching a peak in 2014 with a stable network structure. This was followed by a fluctuating downward trend in the number of network connections between 2015 and 2019, though the decline was minor and within the scope of standard fluctuations. The number of network connections in 2020 decreased to 189 from 194 in 2019, most likely as a result of the COVID-19 pandemic outbreak in 2020 and a decrease in inter-provincial economic sustainability connections due to the increasing spatial liquidity of economic factor resources and a minor decrease in network connections under the influence of market regulation and government macro-control policies. Importantly, the network density value reflects the closeness of the connections between the regions in the network, and the maximum number of possible connections in the network is 870. During the study period, the actual number of connections and the network density value remains low and the closeness of the connected network has to be strengthened, indicating that the spatially connected exchange and spillover effects of the inter-provincial ESDE remains weak; the inter-provincial collaborative and balanced economic sustainable development has the potential to be improved, and the economic connections between regions need to be further enhanced. The network hierarchy remains stable at 0.4343 between 2008 and 2013, yet the network hierarchy declined from 0.4337 to 0.3404 in 2014, with a downward movement towards a 'W'-shaped adjustment between 2014 and 2017, indicating that the spatially connected network hierarchy is less rigid and there is a possibility of spatial spillover between different regions. Spatial spillover is potentially possible between regions, and the interaction and interdependence between regions is growing. The network hierarchy moves to a stable level of approximately 0.3452 after 2017. The network efficiency level exceeds 0.6 across the study period; however, there is an overall downward trend from 0.7315 in 2008 to 0.7069 in 2020, implying that despite the existence of a higher number of spillover and redundant paths of interaction between regions, the cost of ESDE transmission and spillover between provinces is reduced and the paths of interaction move in a more reasonable and coordinated direction. The stability of the network is growing. Moreover, the overall network structure of ESDE is stable and continues   Table 5 show that the number of inter-provincial ESDE connections in China follows a "rising-fluctuating-declining" trend, network density and network efficiency remain relatively stable, network hierarchies show a "stable-declining-stable" trend, and network efficiency declines in general, indicating that the stability of the overall network structure has enhanced. Specifically, the number of network connections rose from 177 in 2008 to 199 in 2014, reaching a peak in 2014 with a stable network structure. This was followed by a fluctuating downward trend in the number of network connections between 2015 and 2019, though the decline was minor and within the scope of standard fluctuations. The number of network connections in 2020 decreased to 189 from 194 in 2019, most likely as a result of the COVID-19 pandemic outbreak in 2020 and a decrease in interprovincial economic sustainability connections due to the increasing spatial liquidity of economic factor resources and a minor decrease in network connections under the influence of market regulation and government macro-control policies. Importantly, the network density value reflects the closeness of the connections between the regions in the network, and the maximum number of possible connections in the network is 870. During the study period, the actual number of connections and the network density value remains low and the closeness of the connected network has to be strengthened, indicating that the spatially connected exchange and spillover effects of the inter-provincial ESDE remains weak; the inter-provincial collaborative and balanced economic sustainable development has the potential to be improved, and the economic connections between regions need to be further enhanced. The network hierarchy remains stable at 0.4343 between 2008 and 2013, yet the network hierarchy declined from 0.4337 to 0.3404 in 2014, with a downward movement towards a 'W'-shaped adjustment between 2014 and 2017, indicating that the spatially connected network hierarchy is less rigid and there is a possibility of spatial spillover between different regions. Spatial spillover is potentially possible between regions, and the interaction and interdependence between regions is growing. The network hierarchy moves to a stable level of approximately 0.3452 after 2017. The network efficiency level exceeds 0.6 across the study period; however, there is an overall downward trend from 0.7315 in 2008 to 0.7069 in 2020, implying that despite the existence of a higher number of spillover and redundant paths of interaction between regions, the cost of ESDE transmission and spillover between provinces is reduced and the paths of interaction move in a more reasonable and coordinated direction. The stability of the network is growing. Moreover, the overall network structure of ESDE is stable and continues to improve, yet has the potential to be improved.

Individual Network Characteristics Analysis
To reveal the centrality characteristics of the individual networks in the spatial association network, three indicator values were calculated for each province: in-degree, out-degree and total number of associated connections (Table 6).

Year 2008
Year 2012 Year 2016 Year 2020 In-Degree

Number of Association Network
In-Degree

Number of Association Network
In-Degree

Number of Association Network
In-Degree Beijing  6  24  30  6  24  30  6  23  29  6  24  30  Tianjin  5  15  20  5  16  21  5  12  17  3  10  13  Hebei  4  3  7  5  4  9  5  4  9  5  5  10  Shanxi  5  2  7  5  3  8  6  5  11  5  2  7  Inner Mongolia  3  1  4  6  1  7  5  2  7  5  1  6  Liaoning  5  2  7  6  2  8  4  2  6  5  1  6  Jilin  6  1  7  6  1  7  6  0  6  5  0  5  Heilongjiang  6  1  7  7  1  8  7  0  7  7  0  7  Shanghai  7  27  34  8  27  35  7  27  34  7  27  34  Jiangsu  5  22  27  5  25  30  5  27  32  6 26 32 Zhejiang Table 6 shows the spillover-benefit relationship for the 30 provinces in China. Specifically, the in-degree represents the beneficial relationship and the out-degree represents the spillover relationship. The out-degree of Beijing, Tianjin, Shanghai, Jiangsu and Zhejiang is higher than the in-degree, indicating that Beijing, Tianjin and Yangtze River Delta have more spillover relationships than beneficiary relationships in the ESDE association network, forming a significant spillover effect on other regions; the out-degree of the northeast, northwest, southwest and central regions is less than the in-degree, suggesting that the these regions have fewer spillover relationships than beneficiary relationships, receiving more spillover from other regions in the association network and showing a significant beneficiary effect. The explanation for this is that the Beijing-Tianjin Urban Agglomeration and Yangtze River Delta regions are in the early stage of economic development mode transformation, focusing on environmental management and protection in the process of fast economic development. Under the pressure of national macro-control policies, especially under the promotion of the regional coordinated development strategy and the concept of sustainable development, capital and technology conducive to sustainable economic development have continued to spread to regions such as the central and western regions and the northeast, thus showing a strong spillover effect, while other regions have a markedly beneficial relationship.

Conclusions
At present, most of the literature on efficiency measurement uses human resources, capital and energy consumption as the input indicators of the EBM model. On this basis, this paper takes "social livelihood" and "environmental protection" as an input index. In order to prevent the potential bias of a single perspective on the total evaluation, six dimensions of input indicators-environment, resources, social livelihood, science and technology, capital, and labor-were designed. This study of the ESDE reflects the unity of human-centered attributes and natural environmental attributes, which fills the gap of previous studies. Based on the ESDE evaluation indicators, the EBM-Malmquist model of super-efficiency of non-expected output was constructed, which overcomes the problem of easy loss of information on the ratio between the target and actual values of inputs or outputs in the calculation and accurately measures the provincial ESDE in China. The Gini coefficient method was applied to analyze the regional differences and sources of differences in ESDE. Combined with the Gaussian density function estimation method, the density function of the distribution of regional differences and the dynamic evolution pattern are concluded. The spatial correlation matrix of ESDE is constructed based on the modified gravity model in combination with the variance analysis. The social network method was used to assess the overall ESDE network structural characteristics, individual network characteristics and spatial network spillover effects.
This study's four main findings are as follows: (1) China's overall ESDE is on an upward trend. The ESDE in the eastern regions of Beijing, Shanghai, Jiangsu, Zhejiang and Guangdong are leading at high levels, while the ESDE in the central regions of Liaoning, Jilin and Heilongjiang are lagging behind. The western regions are mostly at a medium level of ESDE, the eastern regions are in the leading position, and the central and western regions are keeping up with the East, while ESDE in the Northeast is lagging behind. (2) The overall regional differences in China's ESDE have fluctuated and have tended to widen in recent years. On a national level, the ESDE exhibits a "gradient effect" and weak polarization. The ESDE in the eastern region does not appear to be polarized; however, there is a tendency for the differences to widen, and the ESDE in the central region is evidently disparate and fluctuates, while the ESDE in the western region appears to have a gradient effect with multi-polarization in some years. (3) The regional ESDE has shown significant spatial unevenness, with Beijing-Tianjin, the Yangtze River Delta and Guangzhou as the core of the circle structure. The role of the central and western provinces in the spatial network is weak, showing an apparent pattern of geographical differentiation with high density in the east and low density in the west, and the spatial interaction and spillover effects are also weak. Furthermore, the ESDE's inter-provincial transmission and spillover costs are decreasing, the interaction paths are moving towards a more reasonable and coordinated direction, and the stability of the network is improving. Overall, the network structure of ESDE is stable and continues to progress in a positive direction, though there is still considerable improvement potential. (4) Beijing, Tianjin and the Yangtze River Delta have higher spillover relationships than beneficiary relationships in the correlation network, resulting in significant spillover effects on other regions; the node out-degrees of the northeast, northwest, southwest and central regions are less than the node in-degrees, as the spillover relationships in the above regions are lower than the beneficiary relationships; the regions receive more spillover effects from other regions in the correlation network and show significant beneficiary effects.

Recommendation
To promote the sustainable and balanced development of China's economy, three policy recommendations are proposed on the basis of our research findings.
(1) The spillover relationship is low in central and western China, and their spillover effects need to be strengthened. The western region has greater advantages in terms of promoting clean energy compared to the eastern region; the rapid development of clean energy could significantly contribute to sustainable development. Moreover, it would be beneficial to strengthen capital and technology investment in the central and western regions, focus on clean energy development policies, and use clean energy development to enhance sustainable economic spatial spillover relations, thus supporting sustainable economic development.
(2) For economically less developed regions, the level of investment in environmental protection has to be enhanced while maintaining economic growth. These regions thus will need the support and supervision of government and national policies. Meanwhile, it would be beneficial to enhance the leading role of Beijing, Tianjin and the Yangtze River Delta in ESDE to strengthen the sustainable economic development of nearby regions. (3) Strengthen the effectiveness of the policy impact of market regulation and government macro-control, accelerate the spatial movement of economic factor resources and improve the inter-provincial economic connections for sustainable development.
Reduce channels for inter-regional spillover and redundant interactions to enhance the network structure of ESDE, lower the cost of ESDE transmission and spillover between provinces, enhance the hierarchy and stability of spatially connected networks, and comprehensively enhance China's sustainable economic development.
In future research, a more refined division of input variables and output variables in the EBM-Malmquist model can be further studied so as to include most aspects of the connotation of ESDE and increase the measurement accuracy.
Author Contributions: All authors contributed to the study conception and methodology. Software, data collection and formal analysis were performed by X.F. and Y.C. All authors have read and agreed to the published version of the manuscript.