This study seeks to discuss the relationship between REP and CO₂, to determine the relationship, the subsequent empirical model is considered while controlling for other explanatory variables (Equation 1): ln ⁡ CO 2 it = β 0 + β 1 ⁡ ln ⁡ R E P i t + ∑ I = 2 k β I ⁡ ln ⁡ x I , i , t + ε i t (1) Where ε i t is the error term, ln ⁡ CO 2 it is the logarithm of the dependent variable, CO₂ emissions per capita, ln ⁡ R E P i t is the logarithm of the independent variable, REP, ln ⁡ x I , i , t represents the logarithm of other control variables ln ⁡ x i , t = ln ⁡ G D P , ⁡ ln ⁡ E C S , ⁡ ln ⁡ I S , i denotes the index of individual units, and t denotes the index of time.

The long- and short-term dynamics between REP, PGDP, ESC, IS, and per capita CO₂ emissions were investigated through the PMG-ARDL method developed by Pesaran and Yamagata (2008). This method is theoretically based on dynamic panel data theory, cointegration theory, and the ECM. The PMG-ARDL model integrates the ideas of the PMG estimator (Pesaran et al., 1999) in panel data analysis and the ARDL model (Pesaran et al., 2001) in time series analysis, which is a novel tool for panel analysis. The model comprehensively considers individual heterogeneity and the overall characteristics of panel data, and by introducing the lag term, it establishes a dynamic panel model that can effectively address the endogeneity problem. Several empirical studies (Wu et al., 2022; Wang et al., 2021) have employed PMG-ARDL to examine the energy-carbon nexus. A key point to consider is that the PMG-ARDL model presupposes that the order of integration of the variables cannot exceed the first order. However, considering the possibility of CSD in panel data, traditional unit root tests (URT) [LLC (Levin et al., 2002), IPS (Im et al., 2003)] may lead to biased results (Pesaran, 2006). Therefore, the CIPS test (Pesaran, 2007) is chosen for this study, which can effectively deal with the CSD problem and deliver more reliable test results. Meanwhile, to ensure the validity of the cointegration relationship between variables, we successively use the Pedroni cointegration test (CT) (Pedroni, 1999) (for individual heterogeneity) and the Westerlund ECM cointegration test (Westerlund, 2007) (more robust when CSD exists).

On the confirmation of cointegration between variables, the subsequent analysis can be carried out through the PMG-ARDL model. The model is as Equation 2: ∆ Y 1 it = α Ii + γ 1 i Y 1 , i , t − l + ∑ I = 2 k γ li X I , i , t − 1 + ∑ j = 1 p − 1 ∅ 1 ij ∆ Y 1 , i , t − j + ∑ j = 0 q − 1 ∑ I = 2 k ∅ I i j ∆ X I , i , t − j + ε I i t (2)

In Equation 2, Y 1 represents the dependent variable, X l is the explanatory variable, I = 1 , 2 , 3 , 4 , ∆ stands for the difference operator, p, q are the lag orders, and j is the index for the lag order.

The reparameterized model becomes a form of error correction that arranges the variables used in this study and can be constructed as follows Equation 3: ∆ ln ⁡ CO 2 it = φ i E C T i , t + ∑ j = 1 p − 1 δ i , j * ∆ ln ⁡ C O 2 i , t − j + ∑ j = 0 q − 1 ϑ i , j * ∆ ln ⁡ X i , t − j + α i + ε i t (3) where, E C T i , t = ln ⁡ CO 2 i , t − 1 − θ i X i , t , φ i = − 1 − ∑ j = 1 p δ i , j , θ i = − ∑ j = 0 q ϑ i , j /   φ i , δ i , j *   = ∑ s   = j + 1 p δ i , s , ϑ i , j *   = ∑ s = j + 1 q ϑ i , s

In Equation 3, E C T i , t represents the error correction term, φ i represents the speed of adjustment coefficient. θ i represents the long-run coefficients, δ i , j   * and ϑ i , j * represents the short-run dynamic coefficients.

To verify the reliability of the estimation results of the PMG-ARDL approach, FMOLS, and DOLS were applied for robustness testing following Dam et al. (2024).

The WCA method is based on wavelet transform theory. This method is derived from mathematics and has been studied in the field of energy economics, such as by Kuşkaya and Bilgili (2020), Abban et al. (2022), and AlNemer et al. (2023), who have applied the WCA to the United States, African oil-producing countries, and Saudi Arabia, respectively, to study the time-frequency relationship between RE and CO₂. These studies provide a strong reference for this study. This paper uses the WCA developed by Goupillaud et al. (1984) in 1984 to capture the time-frequency dependence among the variables. The breakthrough of WCA is the multiscale decomposition of the time series data, thus mining the interaction between the two at different frequencies, as described in Equations 4–9.

The wavelets used in this study are from the Morlet wavelet series with the following equations: χ t = π − 1 4 e − i ω 0 t e − 1 2 t 2 (4) Where χ is applied to finite time series analysis.

By varying the scale parameter f (determining frequency) and translation parameter g (determining position on the time axis), different wavelet functions with varying frequencies and positions can be obtained. The results after the wavelet transform are as follows: χ g , f t = 1 f χ t − g f , g , f ∈ R , f ≠ 0 (5)

The continuous wavelet transform is the inner product operation of the original time series n(t) with wavelet functions of different scales and positions as follows: W n g , f = ∫ − ∞ ∞ n t 1 f χ t − g ¯ f d t (6)

A cross-wavelet transform is applied to the two-time series: W n m g , f = W n g , f W m g , f ¯ (7) Where W n g , f and W m g , f denote the wavelet transform of two time series n(t) and m(t). In turn, the square formula for wavelet coherence is derived as follows: R 2 g , f = C f − 1 W n m g , f 2 C f − 1 W n g , f 2 C f − 1 W m g , f 2 (8) Where C is the smoothing factor, indicating the time smoothing process. R 2 reflects the degree of correlation, with a value between 0 and 1. The larger R 2 value, the stronger the correlation between the two-time series, conversely, the weaker the correlation.

Since R 2 is a non-negative number, it fails to provide information about the direction of interaction. Therefore, Torrence and Compo (1998) and Pal and Mitra (2017) proposed a phase-based method to detect differences in wavelet coherence and reveal the lead-lag relationship between two-time series. The formula is as Equation 9: Π n m = tan − 1 L C f − 1 W n m g , f O C f − 1 W n m g , f (9) Where L and O represent the Hilbert transform’s imaginary and real part operators, respectively.

Exploratory spatial data analysis (ESDA) is a set of technical methods for exploring the characteristics of spatial data (Hou et al., 2022), the core of which is to delve into the distributional characteristics of spatial objects (Overmars et al., 2003), including global SAC and local SAC analysis (Sun et al., 2013). ESTDA (Rey et al., 2011) systematically analyzes temporal behavior based on ESDA, which compensates for the shortcomings of ESDA ephemeral measures. ESTDA is mainly based on spatial econometrics and the theory of spatial dynamic change, which can reveal the spatial dependence patterns of variables and their evolution paths. Already studies (Kao et al., 2023) have analyzed the characteristics of the spatial relationship between EC and CO₂ using the ESDA method, and then Liu et al. (2023b) explored the spatio-temporal dynamics of CO₂ and the role of influencing factors on CO₂ based on ESTDA. Therefore, the application of ESTDA to the study of the spatial and temporal dynamics of REP and CO₂ in this study is a justifiable extension.

SAC is usually quantified by calculating Moran’s I value. The range of Moran’s index is when the calculation result falls in the interval of (0,1], it indicates a positive autocorrelation and the positive correlation is enhanced with the increase of the value; when the calculation result falls in the interval of, the formula is as follows Equation 10: Moran ′ s   I = ∑ i = 1 n ∑ j = 1 n w i j x i − x ¯ x j − x ¯ ∑ i = 1 n ∑ j = 1 n w i j , S 2 = 1 n ∑ i = 1 n x i − x ¯ 2 (10)

Traditional univariate Moran’s I index primarily measures the SAC of a single variable. Anselin (1995) extended the traditional Moran’s I to a bivariate Moran’s I, enabling the exploration of spatial relationships between two variables and effectively capturing their spatial covariation. This study employs bivariate SA analysis to investigate the spatial response pattern between REP and per capita CO₂ emissions. The formula is as Equations 11, 12: I = ∑ i = 1 n ∑ j = 1 n w i j x i − x ¯ y j − y ¯ S 2 ∑ i = 1 n ∑ j = 1 n w i j (11) I i = z i ∑ j w i j z j (12)

In Equations 11, 12, I and I i are respectively bivariate global Moran’s I index and bivariate local Moran’s I index, n the number of the study area, x i , and y j denote the observed values of two variables in regions i , j , respectively, x ¯ and y ¯ are the corresponding means, S 2 is the sample variance, w i j represent the spatial weight matrix, z i and z j are the variance normalization of different attribute values of region i , j , respectively. Based on I i four types of clusters can be identified: H-H (high-high), L-L (low-low), L-H (low-high), and H-L (high-low).

LISA time paths (Zhang et al., 2022) can describe the dynamic migration changes realized by LISA coordinates in Moran scatter plots in the time dimension, i.e., it can reveal the degree of stability of LISA coordinates over time, and the geometrical characteristics of the LISA time paths are expressed in terms of the relative lengths, and tortuosity, the formula is shown in Equations 13–14. Γ i = n ∑ t = 1 T − 1   d L i , t , L i , t + 1 ∑ i = 1 n   ∑ t = 1 T − 1   d L i , t , L i , t + 1 (13) D i = ∑ t = 1 T − 1   d L i , t , L i , t + 1 d L i , 1 , L i , T (14) Where Γ i is the relative length; D i is the degree of curvature; d L i , t , L i , t + 1 denotes the moving distance of province i between years t and t + 1; the larger Γ i is, the more dynamic the local spatial structure is; Γ i >1 indicates that the moving distance of province i is more than the national average distance, or else it is the other way around. The larger value of D i indicates that the dynamic path of province i is more curved, showing a more dynamic local spatial dependence process, while the smaller value of D i indicates a more stable local spatial dependence direction.

LISA space-time transition (Liu et al., 2023b) is used to reveal the changes of local spatial correlation types in Moran’s I scatterplot, combining local Markov transfer and spatio-temporal leaps to classify the leaps into four types, as shown in Table 3. Among them, “Type III” can be subdivided into two types according to whether the direction of its own and neighboring leaps are consistent or not ( T y p e I I I A : HH_t→LL_t+1, LL_t→HH_t+1, T y p e I I I B : HL_t→LH_t+1, LH_t→HL_t+1), and on the basis of which, further compute the characteristics of the spatiotemporal flow(SF) and spatiotemporal cohesion(SC) of the object of study, the formula is shown in Equations 15–16. m is the total number of leaping units. m is the total number of transition units. S F = Type   I + Type   I I m (15) S C = T y p e   I I I A + Type   I V m (16)

To mitigate the potential estimation bias caused by neglecting CSD in existing studies, this paper examines the CSD in the panel data by using four methods, including the Breusch-Pagan LM test (Breusch and Pagan, 1980) and the Pesaran CD test. The results in Table 4 provide strong evidence of CSD in the sample data, as all tests reject the null hypothesis at the 1% significance level.

Considering the potential impact of CSD on the results of the first-generation panel URT, this paper further employs Pesaran’s CIPS test to conduct more robust tests of stationarity while controlling for CSD, thereby providing a more accurate assessment of the integration properties of the variables. Table 5 provides strong evidence against the occurrence of a unit root at first differences for all variables. Regardless of whether the model contains a constant term or both a constant and trend term. Given that variables are smooth either at I (0) or I (1), there is no evidence of I (2) processes. Therefore, the PMG-ARDL model is appropriate for this paper.

In this study, Pedroni and Westerlund CT were applied and the results are presented in Tables 5, 6. Four out of seven statistics in Table 6 and three out of four statistics in Table 7 reject the null hypothesis at the 1% significance level. Thus, the outcome provides evidence for a long-term association between the variables.

Table 8 displays the PMG-ARDL model results, which demonstrate a robust negative correlation between REP and CO₂ in both the long and short term. Specifically, every 1% increase in REP leads to a reduction in long-term and short-term per capita CO₂ emissions by 0.05% and 0.26%, respectively. Consistent with the previous findings of Zheng et al. (2021), and Wang (2022), which highlights the suppressive effect of RE on CO₂ emissions and emphasizes the importance of increasing the share of RE. A noteworthy finding is that the coefficient of influence of REP on per capita CO₂ emissions in the results is low compared to other variables. Although RE is regarded as an important means of reducing CO₂, its effectiveness is influenced by numerous factors, insufficient capacity of the power system and inefficient transmission are the main bottlenecks limiting the role of RE (Chen and Lei, 2018), and the development of RE has not yet led to a fundamental shift in the structure of power production, the resulting in its relatively limited impact on CO₂ (Yu et al., 2020). To maximize the carbon reduction potential of RE, the policy level should focus on the construction of power grid infrastructure, accelerate the implementation of the strategy of building a smart grid, and synchronize the enhancement of the power system’s acceptance capacity and transmission efficiency.

According to the long-term estimation results, among other explanatory variables, PGDP and ECS exert a powerful positive influence on CO₂, with an increase of 0.75% and 0.53% in CO₂ for every 1% increase in both. This suggests that economic growth and the proportion of coal consumption impede progress toward carbon emission reduction goals in China in the long run. Numerous studies have corroborated this finding (Chen et al., 2019), economic growth is often accompanied by an increase in production activities, which will directly or indirectly increase CO₂ emissions; coal combustion is the main source of CO₂, and the higher the proportion of coal in the ECS, the more CO₂ will be emitted. In contrast, the improvement of IS has a significant contribution to the reduction of per capita CO₂ emissions, and every 1% increase in IS index will reduce CO₂ by 0.32%. Similar findings are found in the literature by Zhao et al. (2022) and Yu et al. (2018). Industrial restructuring is a substantial initiative to promote low-carbon development by improving energy efficiency (Wang and Yang, 2024), which in turn mitigates CO₂ emissions. Based on the above findings it is necessary to promote a society-wide low-carbon transition through market mechanisms and policy regulation. Green investment can be encouraged through tax incentives and green financial instruments to reduce the dependence of economic growth on high-carbon energy. Gradually reduce the share of coal in the ECS, promote coal consumption to reach the peak, and accelerate the formulation and implementation of a stricter plan for the orderly rollout of coal. Improve the decarbonization of the industrial structure, develop strategic emerging industries, encourage the development of green industries such as smart manufacturing and the digital economy, and implement stricter energy-efficiency standards for electricity, a high-energy-consuming and high-emission industry.

The PMG-ARDL model also provides results for short-term dynamics, and as expected the error correction term is significantly negative, confirming the presence of a long-run equilibrium relationship; in other words, shocks in the short run do not permanently change the path of the system, which will always ultimately return to the long-run equilibrium. The ECS has a similar short-term relationship to its long-run relationship, increasing per capita CO₂ emissions, and the short-term relationships for the other control variables are less clear.

This paper employs FMOLS and DOLS to examine the robustness of the estimation outcomes, the outcomes are tabulated in Table 9. The observations of FMOLS and DOLS are consistent with the results of the PMG-ARDL model.

As Figure 3 illustrates, the global Moran’s I value of REP and per capita CO₂ emission in 30 provinces cluster below the zero value and pass the 1% significance test, demonstrating that there is a negative spatial correlation between the two, i.e., the regions with high REP tend to be neighboring to the regions with low values of CO₂, which corroborates with the results of the PMG-ARDL model.

The bivariate LISA agglomeration map (Figure 4) reflects the clustering of the six-time points, and during the study period, the regions with significant correlations (5%) are mainly clustered in northern and southern, and the regions with negative correlations of L-H clustering and H-L clustering are widely distributed. The former spatial distribution is concentrated in the northern part of China, especially in Inner Mongolia, Liaoning Province, and Hebei Province. The latter are primarily situated in the southward regions such as Yunnan, Guizhou, and Hunan. This may be due to the strong industrial base in northern China and its high dependence on coal. Although wind power has been vigorously developed in recent years, the heavy historical burden has led to a relatively slow adjustment of the energy structure, and the pressure to reduce carbon emissions remains significant. In contrast, southern regions are rich in hydropower and solar resources, providing ample conditions for hydropower and photovoltaic power generation. Moreover, the construction of the southern power grid is relatively complete, which is conducive to the grid connection of renewable electricity, the combined effect has led to significant carbon reduction in these southern regions. This spatial difference reflects the unbalanced regional development in China and also provides important insights for regional synergistic emission reduction. In response to the above phenomenon of differentiation between the North and the South and its analysis, the government should implement differentiated low-carbon transformation policies, accelerate coal substitution, industrial upgrading and power grid construction and upgrading in the North, and optimize green power consumption and market-based trading in the south, to enhance the synergistic effect of the national carbon emission reduction.

Interestingly, even within the same region with similar renewable resource conditions, the carbon reduction effects of REP may exhibit different characteristics across provinces, suggesting that intra-regional spatial heterogeneity also deserves attention. Although both Yunnan and Sichuan are major hydropower provinces, Yunnan has performed more prominently in terms of RE emission reductions. The main reason may be that, the industry in Yunnan Province is dominated by tourism and services, which have relatively low fossil energy demands, thus facilitating the expansion of RE. In contrast, Sichuan, with a stronger industrial base and greater reliance on fossil fuels, has higher CO₂ emissions. Additionally, Sichuan as an important power exporter in China, has long undertaken the task of transmitting power to the eastern region, and the priority of power transmission is higher, which to some extent squeezes the space for local consumption of clean energy, whereas Yunnan has less pressure on power transmission, the highest power self-sufficiency rate, with less pressure on carbon emission reduction. This phenomenon suggests that the energy consumption structure and industrial structure are important resistance to the pressure of emission reduction in high RE production regions. Differentiated policies should be formulated according to the regional characteristics. Sichuan Province should optimize the local energy consumption structure and increase the proportion of local clean energy consumption while ensuring the transmission of electricity, and Yunnan Province should further take advantage of its service and tourism industries to promote the development of low-carbon industries and consolidate its renewable energy emission reduction results.

The spatial distribution of LISA time path length and curvature is shown in Figure 5. Provinces with LISA time paths less than 1 accounted for 70% of the total in the study period, and the whole shows a relatively stable local spatial pattern. The spatial structure of Northeast and Eastern China is more stable, but the causes are different, as Northeast China found it difficult to get rid of the fossil-energy-dominated energy structure quickly resulting in slower changes; Eastern China is due to the relatively mature economy and energy structure, and there is less room for advancement, and therefore shows a stable structure. Therefore, Northeast China should accelerate the adjustment of energy structure and reduce the dependence on fossil energy, and East China should further improve the efficiency of renewable energy utilization through technological innovation. The LISA time path tortuosity is generally greater than 1, while 66.7% of the provinces are lower than the national average, suggesting that the direction of dependence between REP and CO₂ local spatial patterns is stable. Regionally, Central and South China and southwestern China have more dynamic spatial variability processes, indicating that these regions have more complex and dynamic interactions between REP and CO₂. For these regions with large spatial variations, regional energy coordination and policy guidance should be strengthened to promote the synergistic development of renewable energy and carbon emission reduction, to adapt to the dynamically evolving energy structure.

The local Moran’s I transfer probability matrix is shown in Table 10. Overall, the structure of local spatial association between REP and CO₂ in China is relatively stable, and most of the provinces and regions have not been separated from their original clustering categories, and the probability of spatial and temporal mobility is significantly lower than the probability of spatial and temporal cohesion, presenting a certain degree of inertia in the transitions. The most prevalent type of transition is TypeⅣ (94.3%), indicating a strong path-locking characteristic, and the proportion of TypeⅠ, TypeⅡ, and TypeⅢ transition is only 3.5%, 2.2%, and 0.0%, respectively, which is a relatively small number of transition phenomena, suggesting that the reforms and changes in this area are more inclined to be adjusted gradually. Given this, policymaking should balance stability and breakthroughs, and while maintaining the continuity and stability of energy transition policies, it should increase reform efforts in key regions, break path dependency, and promote synergistic development of green energy among regions.