CO2 emissions are first aggravated and then alleviated with economic growth in China: a new multidimensional EKC analysis

CO2 emissions have become a topical issue worldwide, but few studies have considered the spatial effect of income on carbon emissions and explored the relationship between CO2 emissions and income by establishing direct, indirect, and total environmental Kuznets curves (EKCs). Using an annual panel dataset collected over the 1997–2017 period in China, this study first analyzed the spatiotemporal evolutionary process of CO2 emissions and subsequently developed direct, indirect, and total EKC-based spatial Durbin model (SDM) and partial derivative approach. These results indicate that, first, CO2 emissions have characteristic positive spatial autocorrelation, with gravity centers that have shifted westward. Second, the direct EKC forms a line, while the total EKC resembles a lying-S shape as well as the total EKC, which indicates that compared to local economic growth, neighboring growth plays a very different role in impacting local CO2 emissions. Furthermore, neighboring economic growth seems to have stronger impacts on local emissions, and the turning point of the total EKC comes much earlier than that of the conventional EKC due to the spillover effects of economic growth. Finally, the growth of the population, as well as the rise of energy intensity, can stimulate CO2 emissions in both local and neighboring regions. Industrialization seems to have a nonsignificant impact on emission changes due to the offsetting effects of the positive direct and negative indirect impacts of the share of secondary industry. Improvements in local urbanization may lead to an increase in emissions, while neighboring improvements may have stronger restricting effects; thus, urbanization improvement is beneficial to emissions reduction. This study provides more scientific information from both local and neighboring perspectives, which may differ from conventional results but still be beneficial for emissions reduction policy-makers to introduce corresponding measures.


Introduction
Issues regarding climate change have attracted increasing attention since the Kyoto Protocol came to fruition on February 16, 2005 (Zhang andSun 2016). According to the Intergovernmental Panel on Climate Change (IPCC) reports, global warming is mainly caused by greenhouse gas emissions, in which carbon dioxide (CO 2 ) is the most prevalent (Griggs and Noguer 2006;Karl 2004). In fact, the energyrelated CO 2 emissions of China account for approximately 30% of the global CO 2 emissions in recent years (Shan et al. 2018) due to its greatest amount of energy consumption in the world. As one of the Paris Agreement members, China plays an important role in reducing global greenhouse emissions and mitigating climate change. In the 14th five-year plan, China plans to achieve the peak of carbon emissions by 2030 and encourage some appropriate regions to reach the peak in advance.
Since the reform and opening up, China's economy has achieved rapid growth, and in 2010 GDP successfully surpassed Japan, becoming the world's second largest economy. In order to promote economic development, China has promoted economic growth through investment and the establishment of new factories (Zhao et al. 2022), which has led to a large release of carbon emissions. Based on the diversity of economic development in different regions of China, adjusting economic structure and energy structure can promote the reduction of carbon emissions, but this behavior may slow the pace of economic growth . Therefore, the contradiction between China's economic development and carbon emission reduction is deepening (Zhao et al. 2022). In order to achieve CO 2 emission reduction, as well as sustainable economic development, the Chinese government must obtain more detailed information on the CO 2 emission process and the corresponding economic growth with consideration for spatial dependence.
In this case, it is of great significance to explore the relationship between economic growth and carbon emissions. According to the related literature, the environmental Kuznets curve (EKC) is an effective analysis tool to study the relationship between CO 2 emissions and economic growth (Auci and Trovato 2018;Luzzati et al. 2018). Using different econometric approaches, many scholars have investigated the relationship between CO 2 emissions and economic growth with corresponding datasets and obtained useful conclusions (Lin 2009;Xu et al. 2012;Zou et al. 2014). However, CO 2 emissions have obvious spatial characteristics (Ding et al. 2019;Li et al. 2019a, b); thus, spatial econometric models should be used to explore the so-called inverted-U curve. Indeed, many scholars have studied the relationship between CO 2 emissions and economic growth (Chen and Lee 2020;Fong et al. 2020;You and Lv 2018), which is beneficial to understanding the real EKC from different perspectives. However, direct, indirect, and total EKCs have not been defined with the partial derivative approach in previous studies. This study attempts to study the relationship between CO 2 emissions and economic growth from direct, indirect, and total EKC perspectives.
The novelty of this paper arises from the following three aspects. First, using a spatial panel data model and the principle of partial derivatives, the direct and indirect (spillover) effects as well as the total effect of all explanatory variables are measured, which reveals the spatial interaction effects among all variables in the model. Second, this study introduces the cubic term of Ln(GDP); we analyze the interaction process between economic growth and carbon emissions in more detail, which is conducive to improving and deepening the existing literature conclusions. Finally, based on the direct, indirect, and total effects of Ln(GDP), (Ln(GDP)) 2 , and (Ln(GDP)) 3 , this paper develops direct, indirect, and total EKCs for CO 2 emissions, which can reveal the relationship between CO 2 emissions and economic growth from local, neighboring, and entire perspectives.
The rest of the paper is organized as follows. The "Literature review" section reviews the related literature. The "Data and methodologies" section introduces the data descriptions and adopted methods. The "Empirical analyses and results" section shows the empirical results and analyses. The "Conclusion" section provides some conclusions and suggestions resulting from this paper.

Literature review
The environmental Kuznets curve (EKC) hypothesis stems from the initial work of Kuznets (1955) and was defined by Panayotou (1993). The concerns of the EKC hypothesis have always been discussed in academic circles due to its importance in pollution reduction policy formulation (Sarkodie and Strezov 2018). The EKC hypothesis posits that at earlier economic development stages, environmental quality will decrease as economic growth occurs, while after a certain turning point of income per capita, environmental quality will continually improve (Du et al. 2018;Sinha and Shahbaz 2018). Usama et al. (2016) examined the EKC with the vector error correction model and confirmed that the formulation of the EKC hypothesis was supported by the value of renewable energy consumption. Ilhan and Usama (2015) investigated whether better governance and corruption control helped form the inverted U-shaped EKC of CO 2 emissions in Cambodia during the period of 1996-2012 based on the generalized method of moments and the two-stage least squares and stated that no evidence was provided to confirm the existence of the EKC in this country. Nutakor et al. (2020) examined the causal relationship between economic growth and CO 2 in Rwanda for 1960-2014 data with the vector autoregression (VAR) model and found bidirectional causality between the real GDP and the corresponding CO 2 emissions. In addition, many other scholars have used such methods to analyze this issue, such as the autoregressive distributed lag (ARDL) approach (Ahmad et al. 2017;Alshehry and Belloumi 2017;Balaguer and Cantavella 2016;Sinha and Shahbaz 2018), cointegration analysis and vector error correction model (VECM) (Ben Jebli and Ben Youssef 2015;Azam and Khan 2016;Zoundi 2016), fully modified ordinary least squares, and pairwise Granger causality method (Alper and Onur 2016). Using various approaches, these studies have tested the EKC hypothesis and achieved many useful conclusions for policymakers, but these nonspatial methods could not identify the spatial interaction effect among different variables.
The first law of geography states that all things are related, and nearer things are more related than distant things (Tobler 1993). Both CO 2 emissions and economic growth have obvious characteristics of spatial dependence (Jorge et al. 2020;Kang et al. 2016). Considering the spatial dependence of household CO 2 emissions, Wang et al. (2018) used a geographically weighted regression (GWR) model to examine the spatial effects of urbanization, energy intensity, energy structure, and income on household CO 2 emissions. Similarly, using the dataset of CO 2 emissions in five North African countries during 1990-2014, Mahmood et al. (2020) found significant spatial dependence in the CO 2 emission distribution. From global and local perspectives, Liu et al. (2021) revealed the spatial aggregation features of CO 2 emissions at the city level and visualized them through mapping. In addition, Meng et al. (2017) noted that the spillover (indirect) effect played an important role in driving regional CO 2 emissions in China, and these results indicated that CO 2 emissions had significant spatial dependence.
Given the notable spatial dependence of CO 2 emissions, spatial econometric models, which can calculate the spatial interaction effects of variables in a model, should be adopted to test the existence of the EKC (Li et al. 2019a, b;Yang et al. 2019). In recent years, some scholars have indeed discussed the relationship between CO 2 emissions and economic growth with spatial econometric methods. For example, considering the spatial interaction effect, Balado-Naves et al. (2018) employed a panel dataset composed of 173 countries over the 1990-2014 period to examine the EKC for CO 2 emissions and found that most countries supported the standard EKC and that some areas appeared to have an inverted-U relationship between neighboring per capita income and national per capita emissions. Munir et al. (2020) found that some previous results about the EKC for CO 2 emissions were biased due to their exclusion of crosssectional dependence, and Munir et al. applied a new panel test of Granger noncausal relationships to obtain much different conclusions (Munir et al. 2020). Additionally, using a spatial panel model to avoid the coefficient estimation error, Kang et al. (2016) stated that the shape of the relationship between CO 2 emissions and economic growth followed an inverted-N trajectory and confirmed that spatial spillover effects significantly affected the EKC shape.
Previous studies on the EKC for CO 2 emissions with spatial econometric models have supported many beneficial explorations and obtained scientific results by considering the spatial interaction effects (Wang and Ye 2017). However, these papers do not pay enough attention to calculating the direct, indirect, and total effects of real GDP per capita on CO 2 emissions per capita. Moreover, most previous studies still obtained the EKC only from the estimated coefficients of spatial econometric models; for instance, Balado-Naves et al. (2018) estimated the spatial EKC with the coefficients of the spatially lagged GDP and (GDP) 2 based on the spatial lag of X model (SLX) and the spatial Durbin error model (SDEM). For SLX and SDEM, the direct and indirect effects of an explanatory variable are equal to the coefficients of this variable and its spatially lagged variable, while for other models, such as SDM and spatial autoregressive (SAR) models, the direct and indirect effects are not equal to corresponding coefficients (Elhorst 2014). Thus, we should develop the general direct, indirect, and total EKCs based on the direct, indirect, and total effects of GDP, (GDP) 2 , and (GDP) 3 , but not their coefficients, and these three types of effects can be calculated with a partial derivative approach (Elhorst 2014).
It is enough to obtain the EKC from the estimated coefficients of non-spatial econometric models, because the coefficients can reflect the effect of independent variables on dependent variable. Though the estimated coefficients of spatial econometric models cannot always represent the effects of corresponding independent variables on the dependent variable, if we explore the effects of independent variable on the dependent variable with some type of spatial econometric models, this may lead to erroneous conclusions, and a partial derivative interpretation of the impact from changes to the variables of different model specifications represents a more valid basis for testing this hypothesis (Lesage and Pace 2009). In their book, LeSage and Pace present specific process of proof on it.
Furthermore, Elhorst (2014) also proved that only for the OLS (non-spatial model) and SEM (spatial error model) the coefficient can reflect the effect of X on Y. In the SLX (spatial lagged of X) model SDEM (spatial Durbin error model), the coefficients of X are equal to the direct effect of X, and the coefficients of WX are equal to the indirect effect of X, while in other spatial models, such as SAR, SAC, SDM, and GNS, the coefficients of X (WX) are not equal to the direct (indirect) effects. Direct, indirect, and total effect of independent variables can be obtained from the partial differential method.
In summary, we only can obtain the conventional EKC (direct effect EKC) from the coefficient of independent variables with non-spatial model, while we can explore the direct, indirect, and total effect EKCs with spatial model other than SEM. In order to get these three types of EKCs, the direct, indirect, and total effects of independent variables should be obtained by spatial econometric models with the partial differential method. The direct effects, not the coefficient, of independent variables can accurately reflect the effects of independent variables on dependent one, and the indirect effects of independent variables, not the coefficient of WX, can accurately reflect the effects from spatial adjacency region. Therefore, we can always obtain the direct, indirect, and total effect EKCs from the direct, indirect, and total effects which can be calculated with the partial differential method, not the independent variables' coefficients.
Under this circumstance, considering the importance of the spatial interaction effects of GDP and other influential factors on CO 2 emissions, this paper first calculates the direct, indirect, and total effects of GDP on CO 2 emissions based on econometric models and the partial derivative approach and then develops the direct, indirect, and total EKCs. This study is a critical supplement to the existing literature related to the EKC hypothesis because the conclusions of the paper can provide more information to policymakers focusing on CO 2 emission reduction from local, neighboring, and nationwide perspectives.

Dependent variables
Carbon dioxide emission is the dependent variable of this paper; it is the key step to test EKC hypothesis. Emissions from previous studies (Wang et al. 2011;Pao et al. 2012) are usually estimated, not directly measured. Wang et al. (2011) estimated carbon dioxide emissions data for 28 provinces based on energy consumption classification for each fuel category, using historical data. Pao et al. (2012) and others used the nonlinear gray Bernoulli model (NGBM) to make multistep predictions on the three indicators of carbon emissions, energy consumption, and real output, and predicted carbon emissions, energy consumption, and real output from 2009 to 2020. However, in many cases, emission estimates are relatively uncertain (Shan et al. 2016;Marland 2012;Gurney et al. 2005). Shan et al. (2016) pointed out that this uncertainty may come from basic energy statistics, accounting scope, fuel carbon content, and other potential sources (Marland 2008; Andres et al. 2014). These uncertainties have led to extensive estimates of CO 2 emissions by different world energy research institutions (Shan et al. 2016).
The uncertainty of China's carbon dioxide emission estimation mainly includes two sources (Shan et al. 2016). The first is the uncertainty of energy statistics. Existing related research collected energy consumption data released by the National Bureau of Statistics of China (Zheng et al. 2015;Lindner and Guan 2014), and comparing the total energy consumption of China's 30 provinces with China's national consumption, found a 20% gap between the total energy consumption of China's 30 provinces and China's national consumption. Meanwhile, Guan et al. (2012), based on two publicly available official energy data sets for 2010, found a gap of 1.4 billion tons between China's 30 provinces' CO 2 emissions and China's national CO 2 emissions. Shan et al. (2016) summarized the existing research and found that the reason for this gap may be due to the different statistical standards (Pan et al. 2013) and units (Shan et al. 2016) applied in different provinces and the whole country. The second source of uncertainty is the difference in estimated emission factors. By reviewing 2368 articles on China's carbon emissions from 2004 to 2014, Shan et al. (2016) found that most of the previous studies collected emission factors from the Intergovernmental Panel on Climate Change or the National Development and Reform Commission of China (NDRC), and less than 1% of the studies used emission factors based on experiments and field measurements (Chen et al. 2009(Chen et al. , 2013Yu et al. 2014). Studies have shown that emission factors from different sources can differ by up to 40% (Liu et al. 2015;Shan et al. 2016).
Based on this, this paper draws on the practice of Shan et al. (2016) to measure China's carbon dioxide emissions using "apparent energy consumption" and updated emission factors (Liu et al. 2015). Using this method to measure carbon dioxide emissions, the new CO 2 emission inventory helps to reduce the uncertainty of China's CO 2 emissions and presents obvious emission-socioeconomic characteristics across regions. Among them, "apparent energy consumption" refers to the estimated energy consumption based on the mass balance of energy (Andres et al. 2014;Shan et al. 2016); the updated emission factor refers to the emission factor based on 602 coal samples from China's largest coal area, which is more in line with China's actual conditions than the IPCC emission factor (Mi et al. 2017;Liu et al. 2015). See Shan et al. (2016) for more detailed operation process.
CO 2 emission data can be obtained from the China Carbon Emissions Account and Data Set (CEADS) website (http:// www. ceads. net/ (Mi et al. 2017)). Because differences in population size lead to changes in carbon emissions, most studies use total carbon emissions to represent carbon emissions ). This paper uses per capita carbon emissions to measure carbon emissions, which can avoid the impact of population on carbon emissions and is more reliable. Table 1 lists the specific factors affecting environmental variables in selected previous studies. According to Table 1, the influencing factors of carbon emissions mainly include economy, population size, urbanization level, industrial structure, technology and international trade factors (Sarkodie and Strezov 2018).

Independent variables
The core explanatory variable of this paper is the level of economic development. Based on previous studies, per capita GDP is used to characterize . Compared with GDP, per capita GDP can eliminate the impact of population size on the level of economic development. In order to control the impact of other factors other than income on carbon emissions, the population, urbanization rate, the proportion of the secondary industry, per capita GDP energy consumption, and foreign trade are taken as control variables based on the existing research (shown in Table 2). The total population is represented by the total population at the end of the year ; the urbanization level was characterized by the ratio of urban population to total population (Charfeddine 2017;Talbi 2017); and the proportion of the secondary industry is characterized by the ratio of the output value of the secondary industry to the total output value of the first, second, and third industries . Per capita GDP energy consumption as a proxy variable for the technical effect, using the ratio of total energy consumption to GDP and total population product (Sinha and Shahbaz 2018); the level of international trade is characterized by unit of gross domestic product trade (Sinha and Shahbaz 2018).

Data preprocessing and description
Data from 30 Chinese provinces for 1997 to 2017 were collected. Due to a lack of data available in some regions in China, including Tibet, Taiwan, Hong Kong, and Macao, we removed these regions from our empirical analyses. Moreover, the data for our variables are all presented in a natural logarithmic form to eliminate heteroscedasticity. Table 2 shows the descriptive statistics of all variables used in the study. Figure 1 shows the real GDP, energy consumption, and CO 2 emission data collected in China. The real GDP grew exponentially from 1997 to 2017, while the fossil energy consumption and CO 2 emissions grew with the economic growth of China. During the period of 1997-2000, CO 2 emissions and energy consumption changed gradually and then increased sharply until 2012. It should be noted that the energy consumption of China first decreased in 2013, and after 2013, both CO 2 emissions and energy consumption gently increased, mainly due to the implementation of the CO 2 emissions reduction strategy.  (2018) Ecological footprint Per capita GDP/GDP 2 , trade openness (ratio of foreign trade to GDP), human capital, biocapacity Talbi (2017) CO 2 Energy consumption of fuel, intensity energy of road transport, real per capita GDP, urbanization, fuel rate Sinha and Shahbaz (2018) CO 2 Per capita GDP/GDP 2 , per capita renewable energy generation, per capita energy consumption, volume of foreign trade, total factor productivity (TFP) Charfeddine (2017) Ecological footprint GDP, Energy use, Urbanization, Fertility, Life expectancy Olale et al. (2018) CO 2 Per capita GDP/GDP 2 , real trade openness (ratio of foreign trade to GDP)   Table 3 shows the results of 4 types of unit root tests for all variables. From the results, we can find that only the variable of Ln(IND) does not pass the Im-Pesaran-Shin test, while it passes the other 3 types of unit root tests, so we can conclude that all data in logarithmic form are stable. Therefore, using these stable data, we can perform empirical analyses with spatial econometric models.

Standard deviation elliptic-gravity center model
A standard deviational ellipse (SDE) can present the overall characteristics of the spatial distribution of geographical elements (Du et al. 2019;Lefever 1926), such as the CO 2 emissions in China. SDE has four main basic parameters: mean center, major axis, minor axis, and azimuth. The mean center is the gravity center (GC) of all points of CO 2 emission data, which reveals their relative positions, while the movement of the GC can show the dynamic evolutionary trajectories of CO 2 emissions over time. The major and minor axes denote their directions and scopes, and the azimuth represents the emissions' main trend. The formulas for these parameters of the SDE-GC model are as follows: where (X, Y) represents the coordinate of the GC for CO 2 emissions; (x i , y i ) is the geographical coordinate of province i in China; (x * i , y * i ) denotes the coordinate deviation between province i and the GC; i , as the weight, denotes the per capita CO 2 emissions in province i ; x , y denotes the standard deviations along the X axis and Y axis, respectively; and denotes the elliptic azimuth defined by rotating clockwise in the direction due north to the long axis of the ellipse.

Global Moran's index (GMI)
The global Moran's index (GMI) can be used to measure the spatial autocorrelation of CO 2 emissions. The results of the GMI ranges from − 1 to 1. If its value is significantly positive, CO 2 emissions have a positive spatial autocorrelation, denoting that regions with high or low emissions are located close together. By contrast, if the value is significantly negative, CO 2 emissions have a negative spatial autocorrelation, denoting that highly and less polluted regions are located close together. If the value of the GMI is 0, no spatial autocorrelation exists in the emissions data, indicating that CO 2 emissions are stochastic. In addition, if the absolute value of the GMI is closer to 1, the degree of spatial autocorrelation is stronger. The GMI model can be described by Eqs. (4) and (5).
where x i and x j are the CO 2 emissions values for the ith and jth provinces, respectively, and x is the average value of emissions for all provinces. w ij is the element of W , and W is a n × n comprehensive spatial weight matrix defined based on a spatial geographic distance weight matrix and spatial economic distance weight matrix. We apply the comprehensive spatial weight matrix of the geographic and economic distance ( W ) as our spatial weight matrix due to its geographic and economic features. The matrix can be defined as follows. W g andW e are the n × n spatial geographic distance and economic distance weight matrices, respectively, while w g ij and w e ij are their elements between provinces i and j , respectively. w g ij can be defined as the reciprocal of the distance between provinces i and j , while w e ij is the reciprocal of the GDP (2016) difference between provinces i and j . Thus, w ij can be determined from the dot product of w g ij and w e ij . In addition, we also use the z score, defined as Eq. (5), to test the significance of the spatial autocorrelation.

Spatial panel model
Given the spatial dependence of CO 2 emissions, spatial econometric models other than nonspatial regression models should be proposed. We use a spatial panel model to quantitatively investigate the relationships between the dependent and independent variables (Lesage and Pace 2009;Shao et al. 2020;Zeng et al. 2019). Using spatial panel models to achieve the direct, indirect, and total EKCs, our analysis involves three steps: introducing the adopted spatial econometric model (SDM), estimation methods, and tests; measuring the direct, indirect, and total effects of Ln(GDP), Ln(GDP) 2 , and Ln(GDP) 3 ; and developing direct, indirect, and total effect EKCs. Specifically, the process is as follows. First, introducing the SDM, both spatially lagged dependent and independent variables are included in the SDM proposed by Lesage and Pace (2009). The SDM can be written as Eqs. (7)-(9).
Where Y , the dependent variable (shown as Eq. (4)), represents CO 2 emissions and is a constant parameter. W is the comprehensive spatial weight matrix, (W • Y) is the spatial interaction effect of CO 2 emissions, and is the coefficient of (W • Y) . The independent variable vector X (shown as Eq. (5)) denotes factors affected on CO 2 emissions, and is the coefficient of these factors. Similarly, (W • X) denotes the spatial interaction effects of influential factors found in neighboring/related provinces on local dependent variables ( Y ), and is the coefficient of these interactive effects. In specification (3), i is an index for the cross-sectional dimension (province) with  1997,1998, ⋯ , 2017 . i and t are the individual-and time-specific effects, respectively. is the error term, which is an independently and identically distributed error term with a zero mean and variance of 2 .
Second, estimation methods are preformed and a process is selected. The ordinary least square (OLS) cannot be used in spatial econometric models due to its biased and inconsistent estimator (Anselin 1988), so the maximum likelihood (ML) method was proposed by Anselin (1988) to estimate spatial econometric models. For the model selection procedure, we first determine whether nonspatial models (OLSs) or spatial models (spatial lag model (SLM) or spatial error model (SEM)) are more suitable (shown as Appendix Fig. 5) with Lagrange multiplier (LM) and robust LM tests. Then, we compare SDM with SLM and SEM based on the Wald and likelihood ratio (LR) test (shown as Appendix Fig. 6). Consequently, the optimal model can be determined; then, on the basis of the optimal model, we estimated the pooled, individual, time-period, and two-way fixed models and found the best model by comparison.
Third, we measured the direct, indirect, and total effects of Ln(GDP), Ln(GDP) 2 , and Ln(GDP) 3 and developed direct, indirect, and total effect EKCs. While calculating these three effects is our main focus, the direct effect is not necessarily the value of the estimated parameters ( ), and the indirect effect is not necessarily the value of the parameters ( ) due to the nonlinear features of the economic system (Elhorst 2014). As the volatility of dependent and/or independent variables can spread to neighboring or related regions through a special path, they may affect other provinces' dependent variables (Anselin 1988). On the other hand, interaction effects may increase or decrease and then return to other provinces. According to the above analyses, direct, indirect, and total effects should be estimated correctly. While these could not be calculated based on the estimated parameters, we tried to measure them using partial differential methods (Elhorst 2014;Lesage and Pace 2009). Thus, we must address two problems: measuring the three effects and testing their significance.
To calculate direct, indirect, and total effects, Eq. (7) is changed to the following: where R is a rest term containing a constant parameter, error term, and individual and/or time-period specific effects. The matrix of partial derivatives of the expected value of y with respect to the kth independent variable of X is determined by Eq. (11).
where w ij is the element of W , elements of the diagonal line in the matrix denote direct effects, and nondiagonal elements denote indirect effects. In general, we only estimate the average value of diagonal and nondiagonal elements as direct and indirect effects, respectively. Thus, the direct effect ( DE ), indirect effect ( IDE ), and total effect ( TE ) of x k are defined by Eqs. (12)-(14), respectively.
where y 1, y 2, ⋯ y n (n = 30) denotes the CO 2 emissions in the 1st, 2nd, ⋯ , and 30th provinces in China, while x 1k, x 2k, ⋯ x nk (n = 30) represents the kth dependent variable in the 1st, 2nd, ⋯ , and 30th provinces. Here, we mainly focus on the dependent variables of Ln(GDP), (Ln(GDP)) 2 , and (Ln(GDP)) 3 . While it is also important to test the significance of direct, indirect, and total effects, as there is no way to compute these values from the estimated parameters or their T values, we attempt simulation from variance-covariance matrices based on the results of maximum likelihood estimation for Eq. (7). We can simulate the parameter formation process due to their presumed normal distribution. We simulated this 1000 times to obtain these three interaction effects and their response T values with MATLAB.
Finally, another key goal is to develop direct, indirect, and total effect EKCs. Based on the direct effects of Ln(GDP), (Ln(GDP)) 2 , and (Ln(GDP)) 3 , the direct EKC can be developed, and the graph of the direct EKC can be obtained with the MATLAB program. Similarly, the indirect and total effect EKCs can also be developed based on the indirect and total effects of Ln(GDP), (Ln(GDP)) 2 , and (Ln(GDP)) 3 . Moreover, the long-run direct, indirect, and total EKCs can be achieved by expanding the sample period through which all features of the three types of EKCs can be observed clearly.

Empirical analyses and results
Based on the above data and adopted methods in this study, we can achieve the spatial dependence features of CO 2 emissions with GMI, their spatiotemporal evolutionary characteristics with the SDE-GC model, and the direct, indirect, and total EKC formations with SDM.

Spatiotemporal evolutionary process
We mapped the gravity centers and standard deviation ellipses of per capita CO 2 emissions in China from 1997 to 2017 with the SDE-GC model to analyze their dynamic spatiotemporal evolutionary process (Fig. 2).
As shown in Fig. 2, the GC geographical coordinates of per capita CO 2 emissions in China 1997 to 2017 ranged from 111.338 to 113.665 E and 35.806 to 36.214 N, from Changzhi City in 1997 to Linfen City in 2017. During the 1997-2005 period, the GC vibrated strongly, and the transfer direction was not clear; after that, it moved northwest. The long axis of the ellipse extended from 1153.158 in 1997 to 1186.588 km in 2017, while the short axis extended from 1020.422 to 1059.361 km. The total moving distance of the gravity center is 354.91 km, and the average annual shift speed is 17.75 km/a. Thus, we can obtain three results as follows. First, CO 2 emissions are mainly concentrated in eastern China because they are developed regions and most energy-intensive industries are distributed there. Second, during the 1997-2017 period, the GC of CO 2 emissions moved toward the west. In 2000, China implemented the strategy of developing the western region, and some industries, especially energy-intensive industries, moved to western China; thus, the GC of the CO 2 emissions moved west. Last, the extension of the long and short axes of the ellipses reveals that the CO 2 emissions have been increasing continuously, and the range is also expanding. Figure 3 presents the distribution of CO 2 emissions of China in 2017. High CO 2 emissions are mainly concentrated in Beijing-Tianjin-Hebei (BTH), Yangtze River D, the Pearl River Delta, and Xinjiang Province. The high emission regions in the BTH region include Hebei and Hebei's neighboring provinces, such as Liaoning, Inner Mongolia, Shanxi, Henan, and Shandong provinces, and those in Yangtze River Delta contain Jiangsu and Zhejiang provinces. Among the BTH regions, Beijing and Tianjin have lower carbon emissions than Hebei. Shanghai has lower carbon emissions than its seven provinces in eight regions of the Yangtze River Economic Belt. Pearl River Delta mainly refers to Guangdong Province. BTH, the Yangtze River Delta, and the Pearl River Delta are the most developed regions in China, and most of the CO 2 emissions are produced therein. Xinjiang province is considered a high-emission region mainly due to the transfer of energy-intensive industries from eastern China.

Spatial distribution of CO 2 emissions
The ten low CO 2 emissions regions can be divided into three categories. The first category contains municipalities directly under the central government, including Beijing, Tianjin, Shanghai, and Chongqing, whose energy-intensive industries have been transferred out. The second category contains Qinghai, Gansu, and Ningxia provinces located in Northwest China and Jilin province located in Northeast China, which are less-developed regions. The last category includes Yunnan and Hainan provinces, which have excellent ecological environments and developed tourism industries.

Spatial dependence of CO 2 emissions
Another spatial feature of CO 2 emissions is their spatial dependence, which can be measured by GMI. GMI cannot only reflect the spatial agglomeration characteristics of CO 2 emissions but it assists in the selection of the best (spatial or nonspatial) model. Table 4 provides the GMI results for CO 2 emissions in China, and all the GMI values from 1997 to 2017 are significantly positive at the 1% level, which means that the province-level CO 2 emissions present obvious positive spatial agglomeration. It can be seen that there is a significant spatial correlation between carbon emissions, but the existing literature does not take into account the spatial effect of economic growth on carbon emissions only considers the relationship between the two from the perspective of non-spatial effect.
As a result, the CO 2 emissions present obvious spatial dependence characteristics, which means that provinces with high CO 2 emissions are likely to be distributed together, and low-emission regions are generally adjacent. These findings are very similar to those of . Based on the results, we have to take spatial dependence into consideration to achieve a more scientific EKC. Based on this, this paper also discusses the impact of economic growth in the region and adjacent regions on carbon emissions in the region, and compares the size of the direct and indirect effects of space in order to deepen the existing research results.

Direct, indirect, and total EKCs
Considering the significant spatial dependence of CO 2 emissions, a spatial panel model should be adopted to test the EKC hypothesis. These methods can not only measure the spatial interaction effects but can also be used to study the relationship between CO 2 emissions and economic growth from both direct and indirect perspectives; subsequently, we can obtain the three types of EKCs, namely, direct, indirect and total effect EKCs.

Spatial econometric model selection
Based on the model selection process above, this study first compares and tests the models for the nonspatial model (OLS) and the spatial models for SLM and SEM. Subsequently, it compares SDM with SLM and SEM to obtain the optimal model. Finally, based on the optimal model, the specific effects, including the pooled, individual fixed, timeperiod fixed, and two-way fixed effects, should be compared. Table 5 displays the test results of the LM and robust LM for the OLS, SLM, and SEM. From the results of pooled effect and time-period fixed of OLS, both LM and robust LM tests for the SLM and SEM are all significant, so spatial models (SLM or SEM) are better than nonspatial model (OLS). From the results of the individual fixed effects of OLS, both LM and robust LM for SEM are significant, so SEM is better than OLS, while from the two-way fixed effect, SLM is better. Consequently, the spatial models (SLM and/or SEM) are better than the nonspatial model (OLS). Table 6 shows the test results of Wald and LR for SDM, SLM, and SEM. From the results of all four specific effects of pooled, individual fixed, time-period fixed, and two-way fixed, we find that both the Wald and LR test statistics are significant, which means that SDM is better than SLM and SEM. From the results of Tables 5 and 6, we can confirm that SLM and/or SEM are better than OLS, while SDM is more suitable than SLM and SEM, so SDM is the best model.
Given that SDM is the best model, we should compare the four specific effects (pooled, individual fixed, time fixed, and two-way fixed) of SDM to identify which is the best, and we take three steps to test it. First, we use two methods to compare the pooled effect and fixed effect. The first method is F statistic, and the F(29, 594) = 67.54, Prob > F = 0.00, so fixed effect is better than pooled effect. Then, another method is to evaluate the 29 individual effects shown as Appendix Tables 8, and 27 coefficients of individual are significant, only 2 are not. Therefore, the individual effects cannot be ignored, and the pooled effect is not the optimal model. Second, the fixed effect is compared with random effect with Hausman test, and statistic is 39.63, p value is 0.00, so we can confirm that the fixed effect is better than random effect. Finally, based on the two steps tests above, we can confirm that the fixed effect should be used, so we should choose the best effect from individual effect, time effect, and two-way effect. Using LR test method, we compare individual effect with two-way effect, and the LR statistic is 44.16, P value is 0.0014, so two-way effect is better than the individual effect. Similarly, the LR statistic between time effect and two-way effect is 954.10, and P value is 0.00, so two-way effect is also better than the individual effect. Therefore, based on the test above, we should take two-way fixed effect SDM model as our optimal model. After obtaining the optimal model specification, the following analytical logic is used. First, based on the partial differential approach, we calculate the direct, indirect, and total effects of all explanatory variables. Second, we perform a significance test for these three effects based on the variance-covariance matrices obtained from the ML results. Third, according to the direct, indirect, and total effects of Ln(GDP), (Ln(GDP)) 2 , and (Ln(GDP)) 3 , we can develop the direct, indirect, and total EKCs shown in Fig. 4.

Results of direct, indirect, and total EKCs
The traditional EKC theory considers the quadratic term of GDP and only considers whether there is an inverted "U" relationship between carbon emissions and income. On the basis of considering the quadratic term of GDP and the cubic term of GDP, this study aims to show the interaction process between income and carbon emissions in more detail.
First, we evaluated the direct EKC (a and c in Fig. 4). In 1997-2017, the short-term EKC also showed a flat upward stage in the long-term EKC. These characteristics showed that the CO 2 emissions in the same region also increased during the sample period with the economic growth of the same region in 1997-2017. China, as the fastest growing economy, continued to pursue rapid economic growth, which also led to a large release of carbon dioxide. With the regional economic growth, the carbon emissions caused by economic growth also increase. Similarly, Shan et al. (2016) also stated that provincial aggregated CO 2 emissions increased from 3160 million tons in 2000 to 8583 million tons in 2012. These findings indicate that the extensive economic growth mode has not been completely transformed and that the share of energy-intensive industries was still larger during 1997-2017.
We can deeply understand the relationship between CO 2 emissions and economic growth from the long-run curve (Fig. 4a), which looks like a lying-S shape. There are 3 obvious stages within the indirect EKC. The CO 2 emissions decrease at first, then have a steady transition stage at the second stage, and finally drop sharply. According to the long-run EKC (Fig. 4a), the CO 2 emissions in the local region decrease sharply with the economic growth of neighboring or related regions at an early stage and then enter a gentle rising stage; after the turning point, emissions will decrease rapidly. We analyze the reasons for this trend may be due to the actual situation in China. In the first stage, the level of economic development is low; economic growth mainly depends on agriculture and other primary industries. With the growth of economy and the improvement of economic level, more funds flow into the field of scientific and technological innovation, which makes the clean technology of pollution industry improve to a certain extent, and then the carbon emissions will decrease with the growth of economy. In the second stage, China continued to pursue rapid economic growth, which in turn brought China an unbalanced economic and social development (Pao et al. 2012).
With the acceleration of industrialization, China has become the fastest growing economy, and the carbon dioxide emissions released by industrialization have become the main source of China's carbon emissions. In the third stage, with the further development of the economy, China has begun to pay attention to the development of "quality" (Pao et al. 2012). The Chinese government has formulated a series of policies to reduce carbon emissions. The government and enterprises have paid more attention to the development and utilization of clean technology and green technology (Feng et al. 2022), which makes carbon emissions significantly reduced with further economic growth. Therefore, in the long run, the relationship between carbon dioxide emissions and economic growth is in the shape of lying-S. Second, we evaluated indirect EKC. The indirect EKC curve can show the relationship between economic growth and local environment in adjacent areas. Since the indirect effects of (Ln (GDP)), (Ln(GDP)) 2 , and (Ln(GDP)) 3 are not significant, there is no indirect EKC. The reason may be that economic development has a more obvious effect on carbon emissions in close areas than in far areas. Economic development will affect regional carbon emissions, and this promotion effect will weaken with the increase of distance. Therefore, the indirect effect of economic development on carbon emissions is not significant, indirect EKC does not exist.
Third, we evaluated the total EKC (b and d in Fig. 4). At present, the short-run total EKC, which is marked with a dashed line in the long-run total EKC, is in a gentle upward phase. The long-run total EKC is also similar to a lying-S, which is similar to the long-run indirect EKC. From the longrun view, after a sharp reduction phase, the CO 2 emissions gently increase until reaching the turning point and then drop rapidly again with further economic growth in China. From the total EKC perspective, the shape of the EKC is a lying-S. The EKC curve of the total effect is similar to the EKC curve of the direct effect. The reason may be that the direct effect of economic development on carbon emissions is closer and more effective, so the direct effect plays a leading role in the total effect. There are some deviations from the conventional inverted-U shape shown in Appendix Fig. 5, which suggests that CO 2 emissions do not have to increase all the time in the early economic development stage; namely, considering the implementation of high-quality economic development policies, emission reduction may occur in advance. Based on the results of direct, indirect, and total EKCs, we can obtain three interesting conclusions as follows.
First, the impact of economic growth on local carbon dioxide emissions is very different from that of adjacent regions. The effect of local economic growth on carbon emissions has always existed, while the effect of economic growth in surrounding areas on carbon emissions in the region is not significant. Therefore, the region's economic growth first limited the local carbon dioxide emissions, then gently stimulate local carbon dioxide emissions, and finally after the inflection point strongly curb local CO 2 emissions. This phenomenon is mainly due to the actual situation in China. Therefore, when formulating some carbon emission reduction policies, local economic growth and China's actual national conditions should be considered comprehensively.
Second, the shape of the total EKC curve is similar to that of the indirect EKC curve. The possible reason is that the direct effect reflects the impact of local economic development on local carbon emissions, while the indirect effect reflects the impact of neighboring economy on local economy, which is limited by geographical distance. Compared with the areas with long distance, the effect of economic development on the areas with close distance is more obvious. Thus, the impact of the local economy on local CO 2 emissions is far better than that of the neighboring economy on local CO 2 emissions, and the direct effect dominates. The shape of the direct EKC curve determines the shape of the total EKC curve. These important and interesting findings suggest that the direct impact of economic development and direct EKC are important and we should pay attention to them.
Finally, the relationship between economic growth and CO 2 emissions has a lying-S (similar to inverted-N) shape but not a regular inverted-U shape. This is slightly different from the results of the conventional EKC shown in Appendix  Fig. 5, which has a regular inverted-U shape. Similarly, Li et al. (2019a, b) confirmed that the relationship between the carbon intensity of human well-being (CIWB) and economic growth has an inverted-N shape in China. Another example came from Ana et al. (2014), who also reported that an inverted-N curve was estimated for the relationship between CO 2 emissions and GDP per capita with EU-27 panel data. It should be noted that the lying-S shape is not a denial of the EKC hypothesis but a useful supplement to it. The lying-S shaped curve, which indicates that the CO 2 emissions do not have to increase constantly during the early stage, contains more information than the conventional EKC hypothesis.
In summary, this study takes the spatial dependence of CO 2 emissions and economic growth into the analysis framework and investigates the relationship between CO 2 emissions and economic growth from direct, indirect, and total effects perspectives. Furthermore, the direct, indirect, and total EKCs can provide more information than the conventional EKC. Compared with the conventional EKC (Appendix Fig. 5), the turning point of the total EKC (c in Fig. 4) will occur much earlier due to the existence of an indirect EKC, and this finding was confirmed by Balado-Naves et al. (2018). However, note that we do not encourage the use of long-run accurate predictions, such as the turning point of the EKC, with these models. The value of these kinds of EKC analyses is the revelation of the developed trend of the relationship between CO 2 emissions and economic growth, as opposed to predicting the accurate value of dependent variables (CO 2 emissions) (Sarkodie and Strezov 2018).

The impacts of other control variables on CO 2 emissions
Based on the partial derivative approach, we can achieve the impacts of control variables, Ln(POP), Ln(URB), Ln(IND), Ln(ENER), and Ln(TRAD), from direct, indirect, and total perspectives (shown in Table 7). Considering the difference between the direction of direct and indirect effects, we take these variables into two types as follows.
The first type includes the variables Ln(POP). The direct and indirect effects of Ln(POP) are − 0.133 and − 0.667, respectively, and the effect was significant, which indicates that an average 1% rise in the local population will lead to decrease of 0.133% in CO 2 emissions in the local region, and a 1% growth in the population of neighboring regions will also result in a decrease of 0.667% in the local region. Thus, due to the significantly negative direct and indirect effects, the total effect is − 0.800, which is greater than that of the direct and indirect effects.
The second type contains Ln(URB), Ln(IND), Ln(ENER), and Ln(TRAD), which have different directions of direct and indirect effects. The direct effect coefficient of Ln (URB) is 0.121, and it is significant at the 1% level. The indirect effect coefficient is − 0.021, but the effect is not significant, which indicates that a 1% rise in the local urbanization rate will lead to a 0.121% increase in CO 2 emissions in the local region, but the impact of urbanization rate in neighboring areas on local CO 2 emissions is not obvious. For the variable Ln(IND), the direct effect coefficient is 0.107, and it is significant at the 1% level. The indirect effect coefficient is − 0.007, but the effect is not significant. It shows that every 1% increase in the proportion of local secondary industry will lead to 0.107% increase in CO 2 emissions in the region, while the proportion of secondary industry in neighboring areas has no obvious effect on local CO 2 emissions. Similarly, the direct effect of Ln(ENER) was 0.353, and the indirect effect coefficient was not significant through the 1% level significance test. It shows that every 1% increase in local energy intensity will increase carbon emissions by 0.353%. In addition, the direct effect coefficient of Ln(TRAD) was not significant, the indirect effect coefficient is − 0.057, and it was significant at the 1% level. The results show that for every 1% increase in the openness of adjacent regions, carbon emissions can be reduced by 0.057%, while the effect of regional openness on local carbon emissions is not obvious. The total effect is significantly negative, and its magnitude is greater than the indirect effect.
In summary, the analyses of direct and indirect effects can explore the impacts of explanatory variables on CO 2 emissions from direct (local) and indirect (neighboring) perspectives. The total effect will be reinforced with the same direction of direct and indirect effects, while it will be mitigated or even nonsignificant with different directions. Specifically, both local and neighboring population growth may reduce carbon emissions. Rising urbanization rates may spur more carbon emissions. The increase in the proportion of secondary industry in the GDP can result in the rise of CO 2 emissions in the local region, but the impact on the surrounding areas is not significant, and the total effect is not significant, which indicates that industrialization is a necessary process of economic development. Consequently, considering the development of the whole country in the long run, we should improve the industrial level in the appropriate region, although these may be harmful for the local region in the short run. In addition, increased local energy intensity can lead to increased carbon emissions. From the total effect of Ln(TRAD), trade openness may promote the reduction of emissions, so we should expand openness as much as possible to promote more technology transfer.

Conclusion
The EKC hypothesis for CO 2 emissions has always been the studied focus due to its importance in examining the relationship between income and the environment. However, few studies separate the direct and indirect effects from the total effect and introduce the third-order term of income for more detailed analysis. Therefore, this paper first studied the spatiotemporal evolutionary process of CO 2 emissions in China based SDE-GC model, and the spatial characteristics with GMI method, and then investigated the relationship between economic growth and CO 2 emissions from the three perspectives of direct, indirect, and total effects with SDM. According to the empirical results above, the following major conclusions can be drawn.
First, China's CO 2 emissions have obvious spatial and temporal evolution characteristics. China's carbon dioxide emissions are mainly concentrated in the economically developed eastern region. With the western regions in China development strategy, some industries, especially energy-intensive industries, gradually shift to the west of China, leading to carbon dioxide emissions also shift to the west. With the extension of the long axis and short axis of the ellipse, CO 2 emissions continue to increase, and the scope of increase is also expanding. High CO 2 emissions are mainly concentrated in China's most economically developed Beijing-Tianjin-Hebei region, the Yangtze River Delta region, and the Pearl River Delta region, and carbon emissions show significant spatial dependence.
Second, the greatest contribution to the existing literature is that we have obtained three types of relationship curves between CO 2 emissions and economic growth: direct EKC, indirect EKC, and total EKC. Compared with the traditional EKC, these three types of EKCs can provide more comprehensive information. In general, the indirect EKC does not exist, and the long-term direct EKC and the long-term total EKC are similar in shape with three obvious stages. From the long direct and total EKCs, CO 2 emissions decline first, then there is a stable transition phase in the second stage, and finally a sharp decline. Because the direct effect has strong pertinence and leading role, the EKC shape of the total effect and the direct effect is similar. Third, if the direct and indirect effects of the indicators are in the same direction, the total effect will strengthen along the same direction, while if not, the total effects will weaken or even not be significant. Local economic indicators may have the same directional impact as neighboring indicators (e.g., population, urbanization rate, energy intensity, trade openness), resulting in an enhanced overall impact on carbon emissions. However, the impact of local and neighboring economic indicators in different directions, such as the proportion of secondary industry, can reduce the impact on carbon emissions.
In addition, our findings provide new evidence that local carbon dioxide emissions are mainly closely related to local economic conditions and less correlated with economic growth in surrounding areas. Therefore, based on the empirical results of direct, indirect, and total effect EKCs, this paper provides more sufficient information for carbon market rule makers and policy-makers to formulate scientific CO 2 emission reduction strategies from the following four aspects.
First, local economic development will have an impact on regional carbon emissions, but the spillover effect of economic growth in neighboring regions on local carbon emissions is not obvious. Therefore, when formulating policies, policy-makers should combine China's basic national conditions and pay attention to guiding the development of local economy. When formulating policies for different regions, the actual situation of the region's geographical location, economic development level, and carbon emissions should be taken into account to adapt to local conditions. Second, since the total population plays a positive role in reducing carbon emissions in local and neighboring regions, appropriate carbon emission reduction can be achieved by changing the population of each region. Urbanization related to population structure is not conducive to the reduction of local carbon emissions. From the total effect, it can still significantly promote carbon emissions. Therefore, policymakers should appropriately reduce the level of urbanization on the basis of comprehensive consideration of the level of urbanization in various regions.
Third, the share of secondary industry represented by industrialization appears to be insignificant for overall emissions, similar to the assertion of Balado-Naves et al. (2018) that the level of service industries is almost insignificant for carbon emissions. Energy intensity is closely related to the level of industrialization, and raising the level of energy intensity will significantly promote carbon emissions in the region regions. Among these control variables, energy intensity has the greatest impact on carbon emissions. Therefore, policy-makers should focus on reducing the average level of energy intensity to reduce carbon emissions.
Fourth, from the perspective of indirect effects and total effects, expanding openness will reduce carbon emissions in the region. The expansion of openness will also bring new technologies while introducing a large number of polluting industries. Therefore, when expanding the degree of opening to the outside world, we should focus on encouraging technology transfer in, paying attention to the types of industries introduced, encouraging the introduction of clean, green and environmental protection industries, and increasing the introduction of green and clean technologies.
It should be noted that there is also some further work needed on this topic of this paper. For instance, under the direct, indirect, and total effect EKCs studied framework, it is meaningful to investigate the spatial heterogeneity of the three types of EKCs. In addition, it is another useful work to explore whether there are some differences in the three types of EKCs among countries with different income levels.