Regional Differences, Dynamic Evolution and Convergence of Carbon Emissions from Rural Residents’ Living Consumption: Evidence from China

Abstract

Actively exploring a reduction in carbon emissions from rural residents’ living consumption (RRLC) is necessary to address climate change and achieve high-quality development of the rural economy. Based on the measurement of the carbon emissions from RRLC in China between the years 2000 and 2021, and it uncovers regional differences, dynamic evolution and convergence. The main findings are as follows: (1) Using the Dagum Gini coefficient, it was found that the differences in carbon emissions from RRLC in the nationwide and low-income level group (LLLG), low-middle-income level group (LMLG), upper-middle-income level group (UMLG), and high-income level group (HHLG) are all significantly decreasing, and the intensity of transvariation is the primary source of the overall difference. (2) Using the kernel density estimation, it was found that the level of carbon emissions from RRLC in the nationwide and the four major regions have generally gone upward, as well as a polarisation phenomenon. (3) Using the Markov chain, it was shown that there is an instability in the carbon emissions from RRLC, which can be transferred downward to the ideal state, but there is also a risk of increasing the upward shift of carbon emissions. (4) The nationwide level and the four regions showed typical σ convergence characteristics and absolute β convergence. After considering the influence of socio-economic and natural climatic factors, conditions β convergence trend is shown. And there are significant regional differences in spatial β convergence. The limitation of this study is that the data on carbon emissions from RRLC are only obtained at the macro level, which cannot accurately reflect the micro and individual impact on RRLC. On this basis, the paper puts forward policy recommendations to reduce the spatial imbalance of carbon emissions from RRLC.

1. Introduction

Water, energy, and environmental protection are currently hot-button issues of global concern, and they are relevant to human life and sustainable development [1,2,3]. Since the mid-20th century, anthropogenic factors have been the major driver of the dramatic increase in global CO₂ emissions [4,5]. Global warming directly endangers the safety of human beings’ health [6,7], biodiversity, and ecosystem balance [8], and the damage to the environment is also irreversible on a time scale of hundreds to thousands of years [9]. Therefore, global coordinated control of human-induced CO₂ emissions is an inevitable choice to deal with climate change. Responding to climate change and advancing the process of sustainable global economic and social development actively, China has always kept its promises, demonstrating the role of a responsible and major country, and is an important contributor and guide in the establishment of a universal eco-civilization [10,11,12].

China has exceeded the United States as the world’s largest CO₂ emissions [13,14]. But China is also one of the important countries committed to taking action to reduce industrial CO₂ emissions [15,16]. Furthermore, China has made many worthwhile inroads into mitigating climate change. In September 2020, China stated its strategic target of “carbon peaking by 2030 and carbon neutrality by 2060”. Centered around the “dual-carbon” goal, the construction of ecological civilization during the “14th Five-Year Plan” period concentrates on achieving synergistic pollution reduction, carbon reduction, green expansion, and economic growth. This will facilitate a comprehensive green transformation of the economy and society, with the ecological environment improvement progressing from quantitative to qualitative changes [17].

According to research, the direct and indirect CO₂ emissions from residential life have exceeded those of the industrial sectors [18]. According to the division based on urban and rural attributes, CO₂ emissions are usually divided into urban residents’ living CO₂ and rural residents’ CO₂ emissions [19]. What is more, the proportion of CO₂ emissions from urban and rural consumption in total is on an upward trend in China [20,21,22]. As the economy grows, CO₂ emissions due to residents’ energy generation maintain a high trend [23], leading to a high vigilance between the environment and human activities [24]. The CO₂ emissions from rural residents’ living consumption (RRLC) exhibit a progressive growth accompanied by spatially clustered higher emissions, which is also an essential area for future CO₂ reduction [25]. Without exaggeration, reducing CO₂ emissions from RRLC is not only an effective method of accelerating the improvement of environmental quality, but also an essential engine to enhance China’s rural productivity growth; moreover, it is an inherent requirement for achieving sustainable development of the Chinese economy.

Consequently, it is imperative to consider multiple factors contributing to carbon emission reduction and include the carbon emissions of RRLC within the research scope and institutional design agenda. Compared with urban residents, CO₂ emissions from RRLC could be more easily ignored. So, what is the level of carbon emissions from RRLC in China? Are there spatial spillover effects in the adjacent areas? How do regional differences evolve, and what is the evolution trend? What factors may cause CO₂ emissions from RRLC? The responses to the above questions will help to scientifically identify the space pattern and regional differences of CO₂ emissions from RRLC and provide a factual reference and empirical support for promoting rural revitalization and achieving high-quality rural economic progression. The framing of the study is presented in Figure 1.

The remaining sections in the study are as follows: Section 2 is a literature review. The materials and methods are presented in Section 3. The results are presented in Section 4. Section 5 is a discussion. Section 6 is devoted to conclusions and policy recommendations.

2. Literature Review

Carbon emissions have always been a broad concern for scholars, and many relevant studies have been carried out at domestic locations and overseas. These studies are mainly reflected in the following three aspects.

2.1. Transitioning from the Production Side to the Consumption Side

In recent years, a significant amount of studies on CO₂ emissions at home and abroad have primarily focused on energy-related industries [26,27,28]. Specifically, these industries are cement production [29], the chemical industry [30], electricity [31], transportation [32,33], heat [34], construction [35], and several other sectors. However, due to shifting living standards and consumption patterns, the yearly CO₂ emissions from residents’ daily lives have steadily risen. As a considerable source of CO₂ emissions, this has become a growing concern for researchers. There have been fruitful research outcomes regarding CO₂ emissions from household consumption, mainly from the international community [36], countries [37], provinces and cities [38,39]. However, there appears to be relatively less emphasis on studying the CO₂ emissions associated with RRLC. With the increasingly prominent rural ecological and environmental problems, a few scholars have begun to pay attention to the CO₂ emissions from RRLC and analyze the contribution of CO₂ emissions from carbon sources and regions [25]. Scholars have utilized various data and methods, including age structure [40,41], gender structure [42], family size [43], and the number of households [44], among other perspectives, to investigate the CO₂ emissions resulting from residents’ living consumption. The demonstration results showed that consumption demand generated by household activities affects CO₂ emissions, which have an important impact on residents’ living consumption [45,46,47]. In addition, Dou et al. [48] presented an analysis of the effects of energy inequality on China’s CO₂ emissions at the household level. Another study by Chen et al. [49] revealed that the contributions of Beijing’s urban and rural residents’ energy usage and carbon footprint to the total household energy consumption were on the rise.

2.2. The Factors Affecting Carbon Emissions

The IPAT model was formulated by Ehrlich and Holdren in 1971 [50], and it incorporated population, affluence, and technology within a unified analytical framework. Rosa et al. [51], on the other hand, addressed the limitations of the IPAT model by introducing the elastic coefficient and extending it to the STIRPAT model, which has gained significant popularity [52,53]. In recent years, the LMDI decomposition method and the DEA model have been widely used in carbon emission measurement [54,55]. Scholars have also proposed many corresponding improved forms of LMDI according to specific research objects [56,57]. In order to conduct a comprehensive analysis of the impact of various factors on CO₂ emissions, extensive discussions have been launched on whether and how a reduction in CO₂ emissions can achieve the established target [58,59]. Considering the primary driving force of CO₂ emission reduction, the current consensus is centered around technological progress [60,61], efficiency improvement [62], economic output [63], and energy structure adjustment [64]. Moreover, scholars have also examined on various aspects such as trade opening [65], business investment [66], urbanization [67], digital industry [68], and other aspects of studying the impact of CO₂ emissions. With the deepening of research, there is an increasing attention given to effects of policy implementation for CO₂ emission reduction. Studies have shown that the carbon market is the most cost-effective means of achieving CO₂ emission reduction. China’s carbon market, in particular, exhibits significant potential for both cost savings and CO₂ emission reductions [69,70].

2.3. Regional and Dynamic Evolution Trend of Carbon Emissions

Based on an examination of urban–rural structure and regional disparities, scholars have observed that CO₂ emissions from households in urban areas are higher than those in rural areas [71]. Furthermore, substantial disparities in household CO₂ emissions exist between the northern and southern regions of China [72]. In the realm of carbon peaking and achieving carbon neutrality, conducting a comprehensive analysis of regional disparities in CO₂ emissions from residents’ consumption can facilitate the coordination of regional low-carbon development efforts and a reduction in such discrepancies [73].

Given the escalating influence of CO₂ on environmental pollution, as well as human physical and mental health, reducing CO₂ emissions has emerged as a worldwide consensus. Consequently, scholars have started focusing on the phenomenon of CO₂ convergence [74,75]. As the levels of CO₂ continued to increase, scholarly attention shifted towards the socio-economic factors influencing the convergence of CO₂ emissions. Xu [76] discovered that China’s per capita CO₂ emissions did not exhibit absolute $\beta$ convergence, but rather demonstrated conditional $\beta$ convergence. Moreover, other scholars have conducted studies on the convergence of CO₂ emissions considering industry heterogeneity and the specific characteristics of CO₂ emission convergence within various industries, including manufacturing [77], transportation [78], construction [79], and finance [80]. Additionally, adopting a dynamic evolutionary approach, Cui et al. [81] utilized the Theil index as an analytical tool to decompose the differences in agricultural regions and examine their spatial–temporal evolution trends. Mi et al. [82] used the Gini coefficient to study carbon footprints to measure the inequality in various provinces. Several studies have employed various methodologies such as kernel density estimation and spatial Markov chains to examine the CO₂ emission intensity of 283 cities in China within the time span of 1992 to 2013 [83]. Others mostly describes the spatial–temporal characteristics of CO₂ emissions in terms of distribution dynamics, convergence characteristics, and influencing factors [84,85]. However, there is a lack of comprehensive measurement of CO₂ emissions from RRLC.

In summary, the previous literature primarily focuses on CO₂ emissions from industrial sources, ignoring the changing patterns in CO₂ emissions from rural residents’ consumption. China possesses a vast land area with variations in resource endowments, economic development levels, and climatic conditions across different regions. Additionally, previous research has failed to conduct a comprehensive analysis of inter-regional and intra-regional dynamics. In terms of spatial disparities, limited scholarly attention has been given to the comprehensive utilization of the Dagum Gini coefficient, kernel density estimation, Markov chain, and convergence for investigating the CO₂ emissions from RRLC in China. Therefore, a deeper analysis of the regionwide variations and sources of RRLC in China will accurately identify the current situation and the underlying reasons for unbalanced development and provide a factual basis for reducing CO₂ emissions from RRLC.

Given this, this paper’s contribution margin is mainly reflected in three main aspects. First, this paper makes up for the deficiencies of existing literature research perspectives and focuses on CO₂ emissions from RRLC in China. It aims to propose emission reduction programs aligned with the carbon peaking and carbon neutrality goals. Thus, the primary focus of this study is to analyze the changes and driving factors of CO₂ emissions from RRLC while accurately quantifying their levels. Second, in view of research, current studies primarily focus on describing the spatial distribution of CO₂ emissions from RRLC but lack analysis of the spatio-temporal probability matrix and the convergence mechanism of the spatio-temporal distribution. This paper combines the Dagum Gini coefficient method, kernel density estimation, and Markov chain to examine the regional differences and dynamic evolution of CO₂ emissions from RRLC. By doing so, it addresses the limitations of previous studies that lacked dynamic analysis. Third, this paper provides strong empirical evidence to support the promotion of rural revitalization and the achievement of CO₂ emission reduction goals. To overcome regional limitations, the paper divides the full sample into four income level groups: low-income level group (LLLG), lower-middle-income level group (LMLG), upper-middle-income level group (UMLG), and high-income level group (HHLG) based on rural GDP. Additionally, the paper proposes measures to reduce the CO₂ emissions from RRLC, offering factual reference and empirical support for promoting rural revitalization, achieving high-quality rural economic development, and helping China achieve its goals of “dual-carbon” development.

3. Materials and Methods

3.1. Carbon Emissions Measurement

The two primary sources of CO₂ emissions from RRLC are direct and indirect living consumption CO₂ emissions. Direct living consumption CO₂ primarily results from the energy directly consumed by rural residents, while indirect living consumption CO₂ mainly arises from the CO₂ emitted during the production of RRLC products. However, China lacks corresponding statistics for indirect living consumption CO₂ emissions due to the involvement of input–output tables in their calculation. Therefore, this research primarily focuses on the CO₂ emissions resulting from rural residents’ direct living consumption. Among them, three primary sources of carbon emissions from rural residents’ living consumption include fossil energy, heat, and electricity. Although rural residents do not directly consume fossil energy in the use of heat and electricity, the production of these two involves a significant amount of fossil energy consumption. They are also included in CO₂ emissions from direct RRLC. Equation (1) for measuring CO₂ emissions is as follows:

$$C=C_a+C_e+C_h$$

(1)

In the Equation (1), $C$ represents the total amount of CO₂ emissions from direct consumption by rural residents; $C_a$ represents the CO₂ emissions from rural fossil energy consumption; $C_e$ represents the CO₂ emissions from rural electricity energy consumption; and $C_h$ represents the CO₂ emissions from rural heat consumption.

Based on references such as the 2006 IPCC National Greenhouse Gas Inventory and the work of Wang et al. [86], the direct CO₂ emissions from the energy consumption of rural residents are estimated by using the consumption of various fossil fuels. The calculation method is demonstrated in Equation (2):

$$C_a=\sum _i^n{E}_i\times NC{V}_i\times CE{F}_i\times CO{F}_i\times \left(44/12\right)$$

(2)

In Equation (2), $E,NCV,CEF$, and $COF$ represent various energy consumption, average low calorific value, carbon content per unit calorific value, carbon oxidation rate, respectively; 44/12 is the coefficient of carbon conversion to carbon dioxide, and $i$ represents the type of fossil energy.

This article includes 14 kinds of fossil energy, namely raw coal, cleaned coal, other washed coal, briquettes coal, coke, coke oven gas, crude oil, gasoline, kerosene, diesel oil, fuel oil, liquefied petroleum gas, refinery gas, natural gas.

Referring to related studies [87,88,89], the CO₂ emissions resulting from energy consumption in the electricity and heating sectors can be calculated by multiplying the electricity consumption and heat supply with their respective emissions factors, as Equations (3) and (4), and then added together:

$$C_e=E_{ele}\times F_{ele}$$

(3)

$$C_h=E_{heat}\times F_{heat}$$

(4)

Here, $C_e$ and $C_h$ represent the CO₂ emissions generated from electricity and heat consumption for each province, respectively; $E_{ele}$ represents the electricity consumption for each province; $E_{heat}$ represents the heat consumption for each province; $F_{ele}$ represents the carbon emission coefficient for electricity consumption; $F_{heat}$ represents the carbon emission coefficient of heat consumption.

3.2. Dagum Gini Coefficient and Decomposition Method

Superior to the traditional Gini coefficient, the Dagum Gini coefficient has continuity, which can better describe the distribution differences. It is more sensitive to changes in both higher-level and lower-level groups and can better reflect overall distribution inequality. However, there are also some issues, such as subjectivity in parameter selection. The Dagum Gini coefficient can be decomposed into three aspects, namely intra-group difference $G_w$, inter-group difference $G_{nb}$, and intensity of transvariation $G_t$, and satisfy $G=G_w+G_{nb}+G_t$ [90]. In this paper, we adopt the Dagum Gini coefficient and its decomposition approach to study the regional differences in CO₂ emissions from RRLC in China [91]. The definition of the Dagum Gini coefficient is shown in Equation (5), and it satisfies condition (6).

$$G=\frac{\sum _{j=1}^k\sum _{h=1}^k\sum _{i=1}^{n_j}\sum _{r=1}^{n_h}\left|{RC}_{ji}-{RC}_{hr}\right|}{2\overline{y}{n}^2}$$

(5)

$$\overline{R{C}_1}\le \overline{R{C}_h}\le \cdot \cdot \cdot \overline{R{C}_j}\le \cdot \cdot \cdot \overline{R{C}_k}$$

(6)

In Equation (5), ${RC}_{ji}$($R{C}_{hr}$) represents CO₂ emissions from RRLC in any province i(r) within region j(h), $\overline{y}$ denotes the overall average level of CO₂ from RRLC in the provincial area, $k$ is the number of regions, $n$ is the number of provinces, $j$ and $h$ are the regional subscripts, $n_j$ and $n_h$ represent the numbers of provinces contained within the regions $j$ and $h$, and $i$ and $r$ are the provincial subscripts. $\overline{R{C}_j}$ denotes the average value of carbon emissions from RRLC in region j. The higher the Gini coefficient is, the more unbalanced the distribution of CO₂ emissions from RRLC is.

3.3. Kernel Density Estimation

The kernel density estimation is a non-parametric estimation method that can adapt to various complex data distributions. It can reflect the detailed patterns of the data and visualize the data distribution with a smooth density curve, with the characteristics of strong robustness and weak model dependency. However, compared to parametric density estimation methods, it lacks interpretability and does not provide explicit parameters to describe the data generation process. This method has gained popularity among many scholars [92,93]. The kernel density estimation can be obtained according to the empirical distribution function, which is depicted in (7):

$$f\left(x\right)=\frac{1}{Nh}\sum _{i=1}^{N}K\left(\frac{X_i-x}{h}\right)$$

(7)

$$K\left(x\right)=\frac{1}{\sqrt{2\pi }}exp\left(-\frac{x^2}{2}\right)$$

(8)

In Equation (8), $K\left(x\right)$ denotes the kernel density function; N stands for the number of observations; $X_i$ denote the independently and identically distributed observations, and $x$ denotes the average of observations. $h$ represents the bandwidth. Generally, we choose among the triangular kernel, quartic kernel, Epanechnikov kernel, and Gaussian kernel. In this paper, we propose adopting the Gaussian kernel to evaluate the dynamic distribution characteristics of the carbon emissions from RRLC.

3.4. Markov Chain Model

By constructing transition matrices, the Markov chain classifies continuous, discrete values into K types and describes the state distribution and evolutionary trend of a random variable sequence within a given state space [94,95,96]. However, the Markov chain method has limitations in the assumptions, stationarity requirements, and data demands. This paper investigates the internal dynamic evolution trend of the CO₂ emission distribution from RRLC in China over time by using the Markov chain transition matrix. The transition probability matrix formula is shown in Equation (9):

$$P_{i,j}^{t,t+d}=P\left\{X_{t+d}=j|X_t=i\right\}=\frac{n_{i,j}^{t,t+d}}{n_i^t}$$

(9)

In Equation (9), $P_{i,j}^{t,t+d}$ represents the probability that CO₂ emissions from RRLC in a particular area will shift from level i in year t to level j in year t + d. $n_{i,j}^{t,t+d}$ is the sum of the number of provinces with the level of CO₂ emission from RRLC transferred from level i in year t to level j in year t + d. $n_i^t$ represents the sum of the number of provinces with the i level of CO₂ emissions from RRLC in year t. In order to fully consider the spatial effects of carbon emissions from RRLC, and to compensate for the traditional Markov chain ignoring the interaction of geospatial factors. Spatial Markov chains are also introduced in this study, which indicates the shifting possibilities of the CO₂ emissions from neighboring areas with different levels to the RRLC in this area.

3.5. Convergence Model

3.5.1. σ Convergence Model

In distribution estimation methods, σ convergence usually applies statistical indicators such as standard deviation and the coefficient of variation to measure the gap between different regions, thus reflecting convergence in a better visualization way. However, σ convergence only focuses on the overall gap and cannot comprehensively evaluate internal variations within a specific region. Supposing the gap in the CO₂ emissions from RRLC at different levels shows a trend of gradually narrowing over time, there is σ convergence. Otherwise, σ convergence does not exist. In this paper, we use the coefficient of variation measure of σ convergence. The calculation equation is as follows (10):

$$\sigma =\frac{\sqrt{\raisebox{1ex}{${\sum }_i^{nj}({RC}_{ij}-\overline{{RC}_{ij}}{)}^2$}\!\left/ \!\raisebox{-1ex}{$n_j$}\right.}}{\overline{R{C}_{ij}}}$$

(10)

Here, $j$ ($j$ = 1, 2, 3, …) represents nationwide and the four regions LLLG, LMLG, UMLG, and HHLG; $i$ ($i$ = 1, 2, 3, …) indicates provinces in the region; $n_j$ is the number of provinces in region $j$; $R{C}_{ij}$ denotes the CO₂ emissions from RRLC in the $j$ region’s $i$ province; and $\overline{R{C}_{ij}}$ denotes the average value of CO₂ emissions from RRLC in region $j$.

3.5.2. β Convergence Model

β convergence is an improved convergence method, which has advantages over the traditional σ convergence in terms of increasing convergence speed and accuracy. However, in some cases, β convergence may result in suboptimal solutions and fail to find the global optimum. In this study, β convergence test is conducted on the CO₂ emissions from RRLC [97,98]. β convergence refers to the state in which the carbon emission from RRLC finally reaches a convergent growth rate with time, including absolute β convergence, conditional β convergence, spatial absolute β convergence, and spatial conditional β convergence. The absolute β convergence model is presented in Equation (11):

$$\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)=\alpha +\beta \mathit{ln}\left(R{C}_{i,t}\right)+u_i+v_t+{\epsilon }_{it}$$

(11)

where $i$ represents the province, $t$ represents the time, $R{C}_{i,t}$ denotes the CO₂ emissions from RRLC in province i during period t, $R{C}_{i,t+1}$ represents the CO₂ emissions from RRLC in province i during period t + 1, and $\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)$ represents the growth rate of province i during period t + 1. β stands for the convergence coefficient, β < 0, and passes the testing of significance, which indicates the existence of a convergence trend in the CO₂ emissions from RRLC in China. $u_i$ and $v_t$ indicate the individual and time effects, respectively, and ${\epsilon }_{it}$ denotes the stochastic disturbance term.

Considering the obvious spatial effect of CO₂ emissions from RRLC, there is a growing level of inter-regional spatial dependence. Due to the presence of various spatial effects, there will be endogenous, exogenous, and error-like spatial interaction effects. Consequently, this paper employs three spatial econometric models: the spatial autoregressive model (SAR), the spatial error model (SEM), and the spatial Durbin model (SDM). The spatial absolute β convergence models are given in Equations (12)–(14):

$$\mathrm{S}\mathrm{A}\mathrm{R}:\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)=\alpha +\beta \mathit{ln}\left(R{C}_{i,t}\right)+\rho W_{\mathrm{ij}}\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)+u_i+v_t+{\epsilon }_{it}$$

(12)

$$\mathrm{S}\mathrm{E}\mathrm{M}:\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)=\alpha +\beta \mathit{ln}\left(R{C}_{i,t}\right)+u_i+v_t+{\phi }_{\mathrm{it}},{\phi }_{\mathrm{it}}=\lambda W_{\mathrm{ij}}{\phi }_{\mathrm{it}}{+\epsilon }_{it}$$

(13)

$$\mathrm{S}\mathrm{D}\mathrm{M}:\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)=\alpha +\beta \mathit{ln}\left(R{C}_{i,t}\right)+\rho W_{\mathrm{ij}}\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)+\theta W_{\mathrm{ij}}\mathit{ln}\left(R{C}_{i,t}\right)+u_i+v_t+{\epsilon }_{it}$$

(14)

where $\rho$ is the spatial lag coefficient, indicating the impact of the growth rate of CO₂ emissions from RRLC in neighboring provinces on the province; $\lambda$ is the spatial error coefficient, indicating the presence of spatial effects in the randomized perturbation term; $\theta$ is the spatial lag coefficient of the independent variable, indicating its influence on the CO₂ emissions from RRLC in adjacent provinces. $W_{\mathrm{ij}}$ is the spatial weight matrix. This paper utilizes the reciprocal of the squared geographic distance as the spatial weight matrix; that is to say, with the shortening of geographical distance, the CO₂ emissions from RRLC are more closely related. The weights are set as follows in Equation (15):

$$W_{\mathrm{ij}}=\left\{\begin{array}{ll}\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$d_{ij}^2$}\right.& (i\ne j)\\ 0& (i=j)\end{array}\right.$$

(15)

However, there are many factors affecting the convergence of CO₂ emissions from RRLC. Due to significant differences in the social development, economic conditions, geographical locations, and natural climate of different regions, the conclusions of convergence will vary after the incorporation of control variables that affect convergence [99]. With reference to relevant studies, this paper includes the rural residents’ disposable income (Cjcsr), urban–rural income gap (Cxcjp), per capita education years (Educa), average household size (Hjgmb), total dependency ratio (Zfyba), and urbanization level (Urban); average temperature (Tempa); average annual relative humidity (Humid); annual precipitation (Rainf); annual sunshine hours (Sunsh); and total sown area of crops (Agric) as the control variables. Cjcsr is represented by the income level of the rural population. Cxcjp is calculated by dividing urban disposable income by rural disposable income. Educa is calculated as years of schooling per rural resident. Hjgmb refers to the average number of individuals per household. Zfyba is expressed as the ratio of the number of persons of non-working age to the number of persons of working age. Urban is characterized by the share of the urban population in the total population. The rest of the data are raw data from statistical sources. The conditional β convergence model is shown in Equation (16):

$$ln\left(\frac{{RC}_{i,t+1}}{{RC}_{i,t}}\right)=\alpha +\beta ln\left({RC}_{i,t}\right)+\gamma X_{i,t+1}+u_i+v_t+{\epsilon }_{it}$$

(16)

The spatial conditional $\beta$ convergence models are described in Equations (17)–(19):

$$\mathrm{S}\mathrm{A}\mathrm{R}:\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)=\alpha +\beta \mathit{ln}\left(R{C}_{i,t}\right)+\rho W_{\mathrm{ij}}\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)+\gamma X_{i,t+1}+u_i+v_t+{\epsilon }_{it}$$

(17)

$$\mathrm{S}\mathrm{E}\mathrm{M}:\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)=\alpha +\beta \mathit{ln}\left(R{C}_{i,t}\right)+\gamma X_{i,t+1}+u_i+v_t+{\phi }_{\mathrm{it}},{\phi }_{\mathrm{it}}=\lambda W_{\mathrm{ij}}{\phi }_{\mathrm{it}}{+\epsilon }_{it}$$

(18)

$$\mathrm{S}\mathrm{D}\mathrm{M}:\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)=\alpha +\beta \mathit{ln}\left(R{C}_{i,t}\right)+\rho W_{\mathrm{ij}}\mathit{ln}\left(\frac{R{C}_{i,t+1}}{R{C}_{i,t}}\right)+\theta W_{\mathrm{ij}}\mathit{ln}\left(R{C}_{i,t}\right)+\gamma X_{i,t+1}+\delta W_{\mathrm{ij}}{X}_{i,t+1}+u_i+v_t+{\epsilon }_{it}$$

(19)

$X_{i,t+1}$ is the control variable that reflects the characteristics of CO₂ emissions from RRLC, and $\gamma$ is the corresponding coefficient matrix.

3.6. Data Sources

Limited by the accessibility of data, this paper focuses on the CO₂ emissions of RRLC in 30 provinces (cities and regions) in China (excluding Tibet, Hong Kong, Macao, and Taiwan) from 2000 to 2021 as a sample of the study. This paper breaks through the regional restrictions by quartile division of the full sample into a low-income level group (LLLG), lower-middle-income level group (LMLG), upper-middle-income level group (UMLG), and high-income level group (HHLG) based on rural GDP. The basic data of rural residents’ energy consumption and CO₂ emissions are mainly from the annual “China Energy Statistical Yearbook” for multiple years; other data mainly primarily derive from “China Statistical Yearbook”, “China Rural Statistical Yearbook”, “China Population and Employment Statistical Yearbook”, provincial statistical yearbooks, the China Economic Net statistical database, and the CSMAR database. The disposable income of rural residents was deflated by the consumer price index, and missing values were filled in by interpolation.

4. Results

4.1. Spatial–Temporal Evolution Patterns of Carbon Emissions from RRLC in China

In order to specifically depict the present status of carbon emissions from RRLC in China, the per capita carbon emissions from RRLC from 2000 to 2021 were measured. Using the partition method in ArcGIS software (v. 10.8.2), the carbon emissions from RRLC in the years 2000 and 2021 were divided into five levels to explore the characteristics of their spatial patterns, and the results are presented in Figure 2.

In 2000, the provinces with low carbon emissions from RRLC areas were Hainan, Guangxi, Jiangxi, Anhui, and Shandong, basically in an “adjacent or scattered distribution”. Provinces with high carbon emissions are Beijing, Shanghai, Guizhou, Tianjin, Shanxi, and Xinjiang, with a decentralized distribution. Intermediate-level CO₂ emission areas, such as Ningxia, Qinghai, Gansu, Hebei, Heilongjiang, and Liaoning, show the spatial evolution of “scattered distribution–clustered distribution”. In 2021, the provinces with low carbon emissions from RRLC were Guangxi, Hainan, Sichuan, and Chongqing, which are characterized by decentralized distribution. The provinces with high carbon emissions were Tianjin, Inner Mongolia, Zhejiang, Beijing, Hebei, and Fujian, evolving with a pattern of “decentralized distribution–clustered distribution”.

Comparative analysis between 2000 and 2021 shows that there is an overall upward trend of carbon emissions from RRLC, accompanied by significant changes in spatial distribution patterns. The main reasons for this can be attributed to two factors. On the one hand, with the improvement of rural residents’ living standards, their lifestyles and consumption habits have been changed, with an example being the widespread use of household appliances. On the other hand, rural residents rely mainly on traditional energy sources, such as coal and oil, leading to an increase in carbon emissions from RRLC. Further analysis reveals that Shanghai, which is located in the LLLG region, transformed from a higher carbon-emission area in 2000 to a low carbon emission area in 2021, indicating a decrease in CO₂ emissions from RRLC. On the other hand, Inner Mongolia, which is in the HHLG region, transformed from being close to a low-carbon-emission area in 2000 to a high-carbon-emission area. This may be attributed to the government’s policies in rural areas of Shanghai that encourage the use of clean energy, at the same time, gradually popularize clean energy sources such as solar power, wind power, and geothermal energy. As we all know, the energy structure of rural areas in Inner Mongolia, which are rich in mineral resources, relies mainly on coal resources, which is compounded by the impact of energy consumption under unique climate conditions.

4.2. Regional Disparities and Sources of Carbon Emissions from RRLC in China

This article utilizes means of the Gini coefficient to assess the overall differences, intra-regional differences, inter-regional differences, and contribution rates of carbon emissions from RRLC in China. The aim is to uncover the spatial variations and sources of carbon emissions from RRLC. The estimation results are presented in

Table 1. Nationwide and intra-regional Gini coefficients of carbon emissions from RRLC in China.

Year	G_Nationwide	G_sub(1) (LLLG)	G_sub(2) (LMLG)	G_sub(3) (UMLG)	G_sub(4) (HHLG)
2000	0.4212	0.2354	0.5360	0.3303	0.2554
2001	0.4052	0.2230	0.5384	0.2762	0.1826
2002	0.4043	0.2377	0.5138	0.2913	0.1915
2003	0.4074	0.2282	0.5355	0.2794	0.2178
2004	0.4127	0.2308	0.5472	0.2572	0.3498
2005	0.3824	0.2984	0.5057	0.2152	0.3013
2006	0.3777	0.2636	0.5329	0.2112	0.2615
2007	0.3692	0.2661	0.5167	0.2295	0.2518
2008	0.3297	0.2598	0.4609	0.2280	0.2256
2009	0.3408	0.2214	0.4843	0.2337	0.2653
2010	0.3428	0.1707	0.5075	0.2415	0.2548
2011	0.3275	0.1681	0.4746	0.2456	0.2480
2012	0.2875	0.1482	0.3990	0.2382	0.2441
2013	0.2708	0.1622	0.3398	0.2427	0.1788
2014	0.2565	0.1404	0.3195	0.2256	0.1789
2015	0.2446	0.1212	0.3114	0.2195	0.1573
2016	0.2524	0.1038	0.3026	0.2639	0.1646
2017	0.2261	0.1057	0.2934	0.1918	0.1562
2018	0.2261	0.1016	0.2611	0.1776	0.1865
2019	0.2166	0.0922	0.2453	0.1672	0.1832
2020	0.2222	0.1115	0.2271	0.1902	0.1876
2021	0.2372	0.1203	0.2196	0.1988	0.2237

,

Table 2. Inter-regional Gini coefficients of carbon emissions from RRLC in China.

Year	G_jh (LLLG–LMLG)	G_jh (LLLG–UMLG)	G_jh (LLLG–HHLG)	G_jh (LMLG–UMLG)	G_jh (LMLG–HHLG)	G_jh (UMLG–HHLG)
2000	0.5254	0.3525	0.4688	0.4891	0.4701	0.3581
2001	0.5159	0.3204	0.4835	0.4707	0.4606	0.3254
2002	0.5082	0.3647	0.4823	0.4512	0.4384	0.3091
2003	0.5135	0.3361	0.4561	0.4694	0.4706	0.3094
2004	0.5187	0.3514	0.4126	0.4703	0.4967	0.3230
2005	0.5057	0.2927	0.3704	0.4499	0.4726	0.2759
2006	0.5120	0.2784	0.3191	0.4627	0.4854	0.2478
2007	0.4946	0.2720	0.3055	0.4558	0.4689	0.2523
2008	0.4288	0.2523	0.2636	0.4096	0.4124	0.2382
2009	0.4332	0.2372	0.2727	0.4290	0.4321	0.2649
2010	0.4553	0.2185	0.2405	0.4556	0.4537	0.2580
2011	0.4187	0.2232	0.2271	0.4313	0.4234	0.2561
2012	0.3580	0.2056	0.2109	0.3736	0.3683	0.2470
2013	0.3280	0.2253	0.1764	0.3712	0.3355	0.2422
2014	0.3157	0.2096	0.1636	0.3603	0.3208	0.2275
2015	0.2831	0.2182	0.1468	0.3486	0.3010	0.2195
2016	0.2763	0.2337	0.1528	0.3622	0.3071	0.2350
2017	0.2770	0.1851	0.1393	0.3362	0.2842	0.1995
2018	0.2157	0.2234	0.1945	0.3110	0.2829	0.1976
2019	0.2011	0.2149	0.1952	0.2949	0.2751	0.1900
2020	0.1958	0.2252	0.2052	0.2944	0.2688	0.2070
2021	0.1855	0.2531	0.2508	0.2945	0.2865	0.2214

and

Table 3. Contribution values and rates of regional differences in carbon emissions from RRLC.

Year	Contribution Value			Contribution Rate (%)
	Intra-Regional	Inter-Regional	Intensity of Transvariation	Intra-Regional	Inter-Regional	Intensity of Transvariation
2000	0.0863	0.1690	0.1659	20.50	40.13	39.37
2001	0.0814	0.1679	0.1559	20.08	41.44	38.48
2002	0.0807	0.1778	0.1458	19.97	43.96	36.07
2003	0.0841	0.1651	0.1581	20.66	40.53	38.81
2004	0.0872	0.1405	0.1850	21.14	34.03	44.83
2005	0.0842	0.0990	0.1992	22.01	25.90	52.09
2006	0.0839	0.0889	0.2050	22.21	23.52	54.27
2007	0.0832	0.0818	0.2042	22.54	22.15	55.31
2008	0.0762	0.0441	0.2094	23.12	13.37	63.51
2009	0.0787	0.0610	0.2011	23.08	17.90	59.02
2010	0.0776	0.0545	0.2107	22.64	15.90	61.46
2011	0.0751	0.0484	0.2040	22.94	14.76	62.30
2012	0.0656	0.0303	0.1916	22.83	10.53	66.64
2013	0.0594	0.0570	0.1544	21.93	21.05	57.02
2014	0.0554	0.0604	0.1407	21.59	23.54	54.87
2015	0.0526	0.0699	0.1221	21.50	28.58	49.92
2016	0.0545	0.0646	0.1333	21.58	25.61	52.81
2017	0.0475	0.0780	0.1007	20.99	34.48	44.53
2018	0.0466	0.1001	0.0794	20.60	44.27	35.13
2019	0.0441	0.0993	0.0732	20.35	45.84	33.81
2020	0.0462	0.1030	0.0730	20.79	46.34	32.87
2021	0.0493	0.1135	0.0743	20.80	47.86	31.34

.

4.2.1. Nationwide Differences in Carbon Emissions from RRLC

Figure 3 presents the overall evolutionary trend of CO₂ emissions of China’s RRLC from 2000 to 2021. During the observed period, nationwide differences in CO₂ emissions of RRLC display a significant downward trend, and the Gini coefficient dropped from 0.4212 in 2000 to 0.2372 in 2021, with an average annual decline rate of approximately 2.70% (Table 1). Although there was a slight upward trend in 2004, 2010, 2016, and 2021, it did not affect the overall downward trend. This indicates a gradual reduction in the spatial disparity of carbon emissions from RRLC. The primary factors contributing to the phenomenon may be as follows: On the one hand, China places significant emphasis on addressing the rural energy consumption structure. Efforts are made to promote energy-saving technologies and low-consumption products and gradually improve the backward rural areas’ energy structure and RRLC structure. Since 2015, on the other hand, with the guidance of accurate poverty alleviation strategies, rural living standards and life patterns in poor areas have been gradually improved, with the difference in rural life quality gradually being narrowed. This has contributed to the adoption of sustainable development characteristics in agriculture, the upgrading of agricultural practices, and a reduction in the overall difference in carbon emissions from RRLC.

4.2.2. Intra-Regional Differences in Carbon Emissions from RRLC

From the evolution trend of the source of difference within the region described in Figure 4, the LLLG can be roughly divided into the change process of “slow decline–rapid rise–obvious decline–steady rise” during the sample investigation period. From 2000 to 2004, it showed a steady downward trend, with the range of change being small. From 0.2354 in 2000 to 0.2308 in 2004, it shows a sharp increase after 2005 and reaches a sample maximum of 0.2984 in 2005. The years 2006–2019 showed a clear downward trend, from 0.2636 in 2006 to 0.0922 in 2019, and the minimum value for the sample examination period was reached in 2019 (Table 1). During the period, the Gini coefficient displayed a trend of significant decline, from 0.2354 in 2000 to 0.1203 in 2021, and the average annual decrease rate reached 3.15%. The LMLG overall fluctuating downward trend. Although the Gini coefficient rebounded in 2004 and reached a maximum value of 0.5472 in 2004, the overall downward trend of the Gini coefficient has not changed. And, there was a sharp downward trend from 0.4746 in 2011 to 0.3398 in 2013; followed by a steady downward momentum from 2014-2021 and reaches a minimum value of 0.2196 for the sample period in 2021. In 2021, compared with 2000, there was a general decline, down by 0.3164. The UMLG generally showed the evolution trend of “obvious decline–weak rise–sharp decline–rise”. Specific manifestations include the following: 2000–2006 showed a significant decline, from 0.3303 in 2000 to 0.2112 in 2006, reaching the maximum for the observation period in 2000; a slight upward trend in 2006–2016; a sharp decline in 2016–2019, followed by an upward trend of fluctuation and lasting until 2021, with the minimum value of 0.1672 for sample investigation period being reached in 2019. During the observation period, the Gini coefficient of carbon emissions from RRLC in UMLG regions showed a significant downward trend, from 0.3303 in 2000 to 0.1988 in 2021, with an average yearly decrease rate of approximately 2.39%. HHLG regions exhibited irregular fluctuations. The main performance had the change characteristics of “obvious decline–sharp increase–rapid decline–fluctuation rise”. In 2000–2001 showed a significant decline, then there was a sharp increase in 2004 to 0.3498, reaching the maximum value during the observation period; in 2004–2017 was a fluctuating decline phase; 2017–2021 showed a steady rise, and in 2017 the minimum value of 0.1562 was reached; the Gini coefficient value in 2021 was 0.2237, compared to 0.2554 in 2000, and its coefficient changes showed less variation. However, judging from its current development trend, there is an overall decline.

On the whole, the intra-regional difference values for the LLLG, LMLG, UMLG and LMLG, all showed a downward trend of fluctuation during the period under examination. The Gini coefficient values of the four were ranked, in descending order, LMLG, UMLG, HHLG, and LLLG, and the mean values were 0.4124, 0.2343, 0.2212, and 0.1823, respectively. The mean of the intra-regional difference in the LMLG is the largest, which indicates that the imbalance in the LMLG is the most prominent, followed by the UMLG, and then the HHLG, and the difference in the LLLG is the smallest. Overall, the differences within the UMLG are relatively close to those of the HHLG. The main reason for this may be that there are significant differences in the level of rural economic development among the provinces within the LMLG regions, which in turn lead to large differences in CO₂ emissions from RRLC. For example, in Beijing, because of the regional basic conditions, resource endowment, and the first mover advantage, farmers’ income not only includes the primary industry but also other property income, and the overall income level is higher; meanwhile, the demand for household appliances, private cars, housing, and other energy-intensive products is larger, so energy consumption correspondingly increases rapidly, leading to increased carbon emissions. However, the rural development in Guangxi and Yunnan is relatively backward, the economic development level is relatively concentrated, the farmers’ income source is single, the consumption level is low, and the consumption structure is converging, so the carbon emissions from RRLC mainly come from daily necessities.

4.2.3. Inter-Regional Differences in Carbon Emissions from RRLC

Figure 5 reflects the evolution of inter-regional differences in carbon emissions from the living consumption of rural residents. It can be found that LLLG–LMLG roughly shows the evolution trend of “steady decline–weak rise–obvious decline”. The period from 2000 to 2008 was a steady decline, from 0.5254 in 2000 to 0.4288 in 2008 (Table 2), followed by a weak upward phase in 2009 to 2010, then a significant downward in 2010–2021. During the sampling examination period, the Gini coefficient of the LLLG–LMLG regions reached the maximum of 0.5254 in 2000, and reached the minimum in 2021 was 0.1855, with a decrease of nearly 0.3399. The overall the LLLG–UMLG regions showed a decreasing trend. Although there were small floating fluctuations in 2000–2002 and 2003–2005, there was no impact on the overall evolution trend. The Gini coefficient value generally declined, reaching 0.2531 in 2021, compared with 0.3525 in 2000, and there was an average annual decline rate of about 1.57%. The overall LLLG–HHLG inter-regional difference also showed a decreasing trend. During the sample survey period, there was a slightly rising trend in 2000–2002 and 2017–2021, but the overall evolution trend showed a steady downward trend. It declined from 0.4688 in 2000 to 0.2508 in 2021, reached the minimum of 0.1393 in 2017, and decreased by 0.218 in 2021 compared to 2000. The LMLG–UMLG inter-regional difference is generally declining. From 2000 to 2011, there was a fluctuating downward trend from 0.4891 to 0.4313; then, there was a steep decline in 2012, followed by a steady decline after 2013 and continuing until 2021. The Gini coefficient in 2021 was 0.2945, a decrease of 0.1946 compared with 2000. The LMLG–HHLG inter-regional difference showed an irregular trend of fluctuating decline. The main performance was characterized by the change of “weak decline–sharp steep increase–sharp decline–rebound–steady decline”. The period from 2000 to 2002 showed a weak decline, with a sharp increase to 0.4967 in 2004, reaching the maximum; 2004–2008 was a rapid decline phase, and 2008–2010 was a rebound and recovery phase; then, it showed a steady decline that lasted until 2021. The Gini coefficient value decreased from 0.4701 in 2000 to 0.2865 in 2021, with an average annual decline rate of about 2.33%. From the perspective of its evolution trend, it was developing in the direction of steady decline. The UMLG–HHLG inter-regional differences are generally on a downward trend. Among them, 2008-2009, 2015-2016, 2017–2021 showed a rebound and recovery trend, and in 2019, the minimum of 0.1900 was reached; the Gini coefficient in 2021 declined by 0.1367 compared with 2000.

The inter-regional difference in CO₂ emissions from RRLC showed a downward trend of decline or volatility, which means that the inter-regional gap in CO₂ emissions from RRLC is narrowing and developing towards a balanced trend. The average value of the Gini coefficient between regions was ranked in the following order: LMLG–UMLG, LMLG–HHLG, LLLG–LMLG, LLLG–HHLG, LLLG–UMLG, and UMLG–HHLG, and the average values were 0.3996, 0.3871, 0.3848, 0.2790, 0.2588, and 0.2548, respectively. This indicates that the difference between LMLG–UMLG is the largest, ahead of the difference in other regions. The difference between UMLG–HHLG is relatively small, which reflects that the carbon emissions from RRLC in HHLG and UMLG are comparable. In addition, the LMLG–HHLG difference ranked second, indicating that under the trend of the overall narrowing of the differences among the regions, the development difference between LMLG–HHLG is still quite different.

4.2.4. The Contribution Rate of Regional Differences in Carbon Emissions from RRLC

Figure 6 reveals the trends in the sources of the regional differences in carbon emissions from RRLC. In terms of the changing process, the intra-regional difference during the survey period showed a relatively flat straight line; that is, an upward trend from 2000 to 2008, with a maximum intra-regional contribution of 23.12% in 2008. And, followed by a steady downward trend until 2021, when the contribution rate was 20.80% (Table 3). The contribution rate between inter-regional generally showed a fluctuating upward trend, experiencing a process of “fluctuation rise–obvious decline–steady rise”. Specifically, the main trend was characterized by a fluctuating rising in 2000–2002, followed by a clear downward trend in 2003–2012 and a steady upward trend in 2013–2021, with a contribution rate 47.86% in 2021, compared with 40.13% in 2000, which is 0.84% per year on average. The intensity of transvariation shows the evolution trend of “fluctuation rise–steady decline”. The period of 2000–2012 showed a trend of fluctuating rising, and the maximum of 66.64% was reached in 2012; then, the momentum of a steady decline was shown in 2021. In 2021, there was a decrease by 8.03% compared to 2000. In terms of the magnitude of contribution rate, the average contribution rates of the intra-regional variation, inter-regional variation, and intensity of transvariation were 21.54%, 30.08%, and 48.38%, respectively. This demonstrates that the intensity of transvariation is the main source of the difference in carbon emissions from RRLC, followed by the inter-regional variation. This shows that the intensity of transvariation rate should be used as a breakthrough to reduce carbon emissions from RRLC.

4.3. The Kernel Density Estimation of Carbon Emissions from RRLC

In order to explore the dynamic evolution of carbon emissions from RRLC in China, this study uses kernel density estimation to analyze the distribution, shape, spread, and polarization phenomena of carbon emissions from RRLC at the nationwide level and in the LLLG, LMLG, UMLG, and HHLG regions. The results are shown below.

4.3.1. Kernel Density Estimates at the Nationwide Level

Figure 7 depicts the shape of the distribution and the dynamic evolution trend of carbon emissions from RRLC in the nationwide during the observation period. In terms of distributional location, the distribution curve’s center position gradually shifts to the right; based on this, it is inferred that the carbon emissions from RRLC continue to increase. From the distribution pattern, the height of the main peak shows a decreasing trend of evolution. On the whole, the distribution curve shows an increasingly gentle trend, and the width of the main peak has broadened, which indicates that the carbon emission from RRLC has a decentralized trend. In terms of ductility, the distribution curve does not display an apparent left-trailing symptom; still, there is an apparent right-trailing symptom, and it tends to shorten over time, which means that the carbon emissions from RRLC in most of the provinces nationwide are low, but some provinces are still in a phase of high-carbon development. From the view of polarization tendency, there are main and side peaks in the distribution curve, with an obvious height difference between these main peaks and the side peaks, which means that there is a polarization phenomenon of CO₂ emissions from RRLC nationwide. In summary, it can be seen that the overall trend of CO₂ emissions from RRLC is gradually rising, with a decentralized distribution pattern, and the overall development gap is narrowing.

4.3.2. Kernel Density Estimation for Four Level Regions

Figure 8 describe the shape and evolution trend of CO₂ emissions from RRLC in the LLLG, LMLG, UMLG, and HHLG during the investigation period, respectively.

First, from the viewpoint of distribution, the peaks of the location curves for the LLLG, LMLG, UMLG, and HHLG regions generally show the trend of a rightward shift, but the processes of change differ from each other, indicating that the level of carbon emissions from RRLC in the four major regions has generally gone upward. Second, in terms of the evolution of the distribution pattern, the heights of the main peaks of the distribution curves of the LLLG and LMLG regions are roughly rising, while the heights of the main peaks of the UMLG and HHLG regions have a general the evolution trend of declining. This manifests that the difference in the development of carbon emissions from RRLC in the LLLG and LMLG regions is gradually narrowing, while the development of the UMLG and HHLG regions is expanding. Again, in terms of the perspective of distribution extension, there is a clear right-trailing phenomenon in the distribution curve of LMLG, which indicates the existence of provinces with higher levels of development in carbon emissions from RRLC within the region. Finally, regarding the phenomenon of distribution polarization, the LLLG, LMLG, UMLG, and HHLG regions all show clear bimodal or multimodal peaks, indicating those regions the existence of the multi-polarization phenomenon.

Taken together, whether at the nationwide level or in the four level regions, the development of carbon emissions from RRLC has its own characteristics. It should be pointed out that the absolute amount of carbon emissions from RRLC in China are large. This shows that the problem of carbon emissions from RRLC should not be ignored and that rural areas have great potential for energy saving and carbon reduction. In the face of the opportunity period of the “dual-carbon” development strategy, we are facing unprecedented challenges. Considering the different levels of economic development in various regions, we should break through the development bottleneck by changing the energy consumption structure in rural areas and aim to build a beautiful and livable countryside.

4.4. Markov Chain Analysis of Carbon Emissions from RRLC in China

The Dagum Gini coefficient and kernel density estimation method have been used to study the carbon emissions from RRLC. In the following subsections, Markov chain analysis is applied to explore the internal dynamic evolution of the distribution of carbon emissions from RRLC in China over time.

4.4.1. Markov Chain Analysis

Kernel density estimation generally grasped the distribution dynamics of carbon emissions from RRLC. We use the Markov chain analysis method to further investigate the direction and transfer probability of carbon emissions from RRLC in China. Based on the carbon emissions from RRLC, the total sample quartiles are divided into the low carbon emission level group (LLCE), low-middle carbon emission level group (LMCE), upper-middle carbon emission level group (UMCE), and high carbon emission level group (HHCE); then, the transition probability matrix of the carbon emissions of China’s RRLC from 2000 to 2021 was obtained. The results are shown in

Table 4. Traditional Markov chain transition probability matrix of carbon emissions from RRLC.

t/t + 1	I	II	III	IV
I	0.8363	0.1515	0.0061	0.0061
II	0.0305	0.7683	0.2012	0.0000
III	0.0064	0.0577	0.7821	0.1538
IV	0.0069	0.0000	0.0345	0.9586

Note: I indicates low carbon emission level group (LLCE), II indicates low-middle carbon emission level group (LMCE), III indicates upper-middle carbon emission level group (UMCE), and IV indicates high carbon emission level group (HHCE).

.

The probability values on the diagonal are always greater than those on the off-diagonal. The probability of provinces in LLCE, LMCE, UMCE, and HHCE maintaining the original level after one year is 83.63%, 76.83%, 78.21%, and 95.86%. It can be inferred that the development of carbon emissions from RRLC is influenced by path dependence and lock-in. In addition, the probability values at both diagonal ends are higher than those in the middle of the diagonal. In other words, the convergence probability of LLCE and HHCE is higher than that of LMCE and UMCE provinces, and the club effect of LLCE and HHCE convergence is more obvious. This shows that there is a certain “Matthew effect” in the carbon emissions from RRLC. In further discussion, there is an ease of transfer between adjacent types of classes, which suggests that the carbon emission from RRLC is a gradual process and that it is difficult to achieve a “leapfrog” transfer. At the same time, the probabilities of LLCE, LMCE, and UMCE being transferred to an upward level after one year are 15.15%, 20.12%, and 15.38%, respectively. It is evident that the carbon emission from RRLC is a dynamic process, tortuous and fluctuating. The probability of a downward shift of the LMCE, UMCE, and HHCE is 3.05%, 5.77%, and 3.45%, respectively, which indicates that there is some space for decreasing carbon emissions from RRLC. Generally, there is a greater probability of an upward shift. Therefore, attention should be paid to preventing an increase in carbon emissions from RRLC, maintaining the current level of development and striving for a downward shift.

4.4.2. Spatial Markov Chain Analysis

In this paper, spatial factors are taken into consideration, and a spatial Markov transition probability matrix is established. The obtained results are displayed in

Table 5. Spatial Markov chain transition probability matrix of RRLC carbon emissions.

Type of Lag	t/t + 1	I	Ⅱ	III	IV
I	I	0.9216	0.0784	0.0000	0.0000
	II	0.4000	0.6000	0.0000	0.0000
	III	0.0000	0.0769	0.8462	0.0769
	IV	0.0000	0.0000	0.0000	0.0000
II	I	0.8632	0.1368	0.0000	0.0000
	II	0.0364	0.8182	0.1454	0.0000
	III	0.0400	0.0800	0.8800	0.0000
	IV	0.0000	0.0000	0.0000	1.0000
III	I	0.4737	0.4211	0.0526	0.0526
	II	0.0118	0.7647	0.2235	0.0000
	III	0.0000	0.1017	0.6949	0.2034
	IV	0.0256	0.0000	0.0513	0.9231
IV	I	0.0000	0.0000	0.0000	0.0000
	II	0.0000	0.6842	0.3158	0.0000
	III	0.0000	0.0000	0.8136	0.1864
	IV	0.0000	0.0000	0.0312	0.9688

Note: I indicates low carbon emission level group (LLCE), II indicates low-middle carbon emission level group (LMCE), III indicates upper-middle carbon emission level group (UMCE), and IV indicates high carbon emission level group (HHCE).

.

First, there are differences in the transition probability matrices that remained stable or shifted upward or downward under the four different levels of spatial lag types. Taking the spatial lag in LLCE neighborhoods as an example, the probabilities of remaining stable in LLCE, LMCE, UMCE, and HHCE areas are 92.16%, 60.00%, 84.62%, and 0%, respectively; the LLCE probability value in Table 4 without considering the spatial factors is 83.63%, with an obvious difference. Second, the probabilities on the LLCE, LMCE, UMCE and HHCE types of lag diagonal lines are greater than those on the nondiagonal lines. This shows that under the spatial spillover effect, the carbon emissions from RRLC are still stable, and there is a phenomenon of “club convergence”. Third, there are non-zero probabilities on both sides of the diagonal line, indicating that the carbon emissions from RRLC are unstable. Although they can shift downward toward the desired state, there is also a risk of upward shifts in carbon emissions. Fourth, the same type of lag has different effects on different levels. For example, when spatial lag is considered for UMCE neighborhoods, the probability of downward transfer is 1.18%, 10.17%, and 5.13% for LMCE, UMCE, and HHCE areas, respectively. This indicates that the transfer is affected by both the lag type and the original level of carbon emissions from RRLC.

4.5. Convergence Analysis of Carbon Emissions from RRLC in China

4.5.1. σ Convergence Analysis

In order to more accurately examine the evolution trajectory and trend characteristics of the differences in carbon emissions from RRLC, σ convergence model is adopted to test the convergence of carbon emissions from RRLC in China. Specifically speaking, σ convergence, absolute β convergence, conditional β convergence, spatial absolute β convergence, and spatial conditional β convergence are applied, so as to further comprehensively reveal the convergence feature. The outcomes are listed in

Table 6. The σ coefficients of carbon emissions from RRLC at the nationwide level and four regions.

Year	Nationwide	LLLG	LMLG	UMLG	HHLG
2000	0.8818	0.4534	1.4465	0.6375	0.5055
2001	0.8866	0.4390	1.4392	0.5288	0.3949
2002	0.8577	0.4668	1.3593	0.5764	0.4052
2003	0.9409	0.4567	1.4797	0.5415	0.4484
2004	0.9208	0.4519	1.4971	0.5191	0.7685
2005	0.8588	0.6031	1.4303	0.4076	0.5737
2006	0.9334	0.5271	1.5377	0.4132	0.5010
2007	0.9057	0.5470	1.4936	0.4527	0.4802
2008	0.7628	0.5375	1.2856	0.4448	0.4368
2009	0.8291	0.4679	1.3769	0.4619	0.5162
2010	0.8614	0.3687	1.4686	0.4672	0.5019
2011	0.8152	0.3508	1.3677	0.4657	0.4950
2012	0.6109	0.2962	1.0750	0.4523	0.4791
2013	0.5437	0.3132	0.9613	0.4619	0.3430
2014	0.5094	0.2671	0.9105	0.4271	0.3400
2015	0.4832	0.2335	0.8370	0.4158	0.3000
2016	0.4783	0.2000	0.7893	0.4957	0.3159
2017	0.4261	0.2144	0.7348	0.3625	0.3058
2018	0.4210	0.1965	0.6231	0.3356	0.3606
2019	0.4019	0.1841	0.5728	0.3214	0.3522
2020	0.4146	0.2143	0.5211	0.3747	0.3614
2021	0.4474	0.2376	0.4842	0.3923	0.4261

.

To measure the convergence trend of the nationwide level and four regions of carbon emissions from RRLC, this study adopts coefficient of variation conduct the σ convergence test, as presented in Figure 9. From the tendency of change, the coefficient of variation at the nationwide level shows a fluctuating downward evolution process, fluctuating down from 0.8818 in 2000 to 0.4474 in 2021, and in 2019 reached the minimum of 0.4019 within the sample observation period. The coefficients of variation for LLLG areas show a “slow rise–continuous decline” process, with the coefficients of variation slowly rising from 0.4534 in 2000 to 0.6031 in 2005, reaching a maximum value in 2005, and then decreasing continuously to 0.2376 in 2021, reaching a minimum of 0.1841 in 2019. The coefficients of variation of LMLG regions follow a trend of “fluctuating downward–transient upward–continuous downward”, with the coefficient fluctuating downward trend from 1.4465 in 2000 to 1.2856 in 2008, showing an upward trend until 2010, and then showing a continuous decline until 2021. The UMLG regions show a fluctuating downward trend, from 0.6375 in 2000 to 0.3923 in 2021, and the minimum of 0.3214 for the sample period is reached in 2019. The HHLG regions show a fluctuating upward and then fluctuating downward process of change, increasing steeply to the maximum of 0.7685 in 2004, and then exhibit a fluctuating downward process.

It can be clearly observed that the σ coefficients of carbon emissions from RRLC at the nationwide level and the four regions generally show a decreasing trend, indicating that the development trends of carbon emissions from RRLC at the nationwide level and the four regions are characterized by a typical σ convergence. The decrease in the LMLG regions is the largest and closest to the nationwide trend of convergence; the decrease in the σ convergence coefficient of the LLLG and UMLG regions is the next largest and slightly weaker convergence; the fluctuation in the σ convergence coefficient of the HHLG region is larger, and the rate of convergence is slower.

4.5.2. β Convergence Analysis

The absolute β convergence test is a method with a relatively high degree of quantization.

Table 7. Absolute $\beta$ convergence test results for carbon emissions from RRLC.

Variable	Nationwide	LLLG	LMLG	UMLG	HHLG
$\beta$	−0.1634 ***	−0.0680 *	−0.0904 **	−0.3268 ***	−0.3038 ***
	(0.0210)	(0.0374)	(0.0363)	(0.0541)	(0.0566)
_cons	−0.2382 ***	−0.1310 **	−0.0994	−0.4127 ***	−0.6143 ***
	(0.0433)	(0.0600)	(0.0783)	(0.1021)	(0.1258)
control variables	No	No	No	No	No
time fixed effect	Yes	Yes	Yes	Yes	Yes
regional fixed effect	Yes	Yes	Yes	Yes	Yes
N	630	147	168	168	147
R2	0.1945	0.2391	0.2886	0.4084	0.3088

Note: *, **, and *** indicate significance at the 10%, 5%, and 1%, respectively; standard errors are in parentheses.

quantifies the test results of absolute β convergence in nationwide and the four major regions.

As can be noticed, the estimated index of β is markedly negative, passing the 1% significance test at the nationwide level and the UMLG and HHLG regions. The coefficients are negative, even though the LLLG and LMLG regions only pass the 10% and 5% significance levels, respectively. This reveals that there is an absolute β convergence tendency in the CO₂ emissions from RRLC across the nationwide and the LLLG, LMLG, UMLG, and HHLG regions. The carbon emission level from RRLC will finally converge to reach a stable level as time goes by. Differences exist in the rate of convergence within each region without controlling for extraneous influences. From the absolute β value of the estimation coefficient, the speed of convergence is the highest in the HHLG and second highest in the HHLG, nationwide moderate, and the LLLG and the nationwide level have relatively lower speeds. The possible reason is that widespread emphasis in UMLG and HHLG regions on the adoption of clean energy and energy efficiency improvements to reduce carbon emissions from RRLC, resulting in a faster convergence of emission levels.

The results of conditional β convergence for RRLC are shown in

Table 8. Conditional β convergence test results for carbon emissions from RRLC.

Variable	Nationwide	LLLG	LMLG	UMLG	HHLG
β	−0.1732 ***	−0.1718 ***	−0.1970 ***	−0.5456 ***	−0.5154 ***
	(0.0217)	(0.0465)	(0.0453)	(0.0759)	(0.0770)
_cons	−4.6334 *	−3.2425	−11.6904 **	4.5491	−22.0068 ***
	(2.4118)	(3.8253)	(5.3930)	(7.4132)	(6.5120)
control variables	Yes	Yes	Yes	Yes	Yes
time fixed effect	Yes	Yes	Yes	Yes	Yes
regional fixed effect	Yes	Yes	Yes	Yes	Yes
N	630	147	168	168	147
R2	0.2190	0.4181	0.3845	0.5080	0.4247

Note: (1) *, **, and *** indicate significance at the 10%, 5%, and 1%, respectively, with standard errors in parentheses. (2) The control variables include the effects of the disposable income of rural residents, urban–rural income gap, per capita education years, average household size, total dependency ratio, urbanization level, average temperature, average annual relative humidity, annual precipitation, annual sunshine hours, and total sown area of crops on the carbon emissions from RRLC.

. It can be observed that the estimation values are always negative, and all of them pass the 1% significance test; in other words, the CO₂ emissions from RRLC at the nationwide level and for the four major regions have the characteristics of conditional β convergence. Also, the results indicate that after considering the influence of socio-economic and natural climatic factors such as the disposable income of rural residents, per capita education years, urban–rural income gap, average household size, total dependency ratio, urbanization level, and total sown area of crops, average temperature, average annual relative humidity, annual precipitation, and annual sunshine hours, the convergence trend of the CO₂ emission index from RRLC in the research areas of the paper remains stable, and the results are robust. The absolute β values of the estimation coefficient for the nationwide level, LLLG, LMLG, UMLG, and HHLG have been raised, demonstrating that the relevant socio-economic and natural climatic characteristics have accelerated the convergence of the CO₂ emissions from RRLC.

Considering the apparent spatial effect of carbon emissions, the inter-regional spatial interdependence of carbon emissions from RRLC has been strengthening. To show the stability and validity of the β convergence spatial model, for this paper, a spatial dynamic panel data regression model was constructed to further examine the spatial absolute β convergence and spatial conditional β convergence characteristics of carbon emissions from RRLC. The results are shown in

Table 9. Spatial absolute β convergence test and spatial conditional β convergence results for carbon emissions from RRLC.

Variable/Type	Spatial Absolute β Convergence	Spatial Conditional β Convergence
	(1) LMLG	(2) LLLG	(3) LMLG	(4) UMLG
β	−0.0447 ***	−0.1542 ***	−0.0935 ***	−0.2396 ***
	(0.0120)	(0.0373)	(0.0265)	(0.0567)
$\rho$/$\lambda$	0.2252 **	−0.4690 ***	0.2324 **	−0.2511 **
	(0.0913)	(0.1361)	(0.0955)	(0.1023)
control variables	No	Yes	Yes	Yes
time fixed effect	No	Yes	Yes	Yes
regional fixed effect	No	Yes	Yes	Yes
model selection	SEM	SAR	SEM	SAR
N	168	147	168	168
R2	0.0725	0.0099	0.1145	0.0195

Note: (1) **, and *** indicate significance at the 5%, and 1% levels, respectively, with standard errors in parentheses. (2) In the absolute β convergence model, the nationwide level, LLLG, UMLG, and HHLG did not pass the test of the spatial econometric model, so the traditional OLS regression model was directly used. Because the results have been shown in the table, they are not listed repeatedly. (3) In the conditional β convergence model, the nationwide and HHLG did not pass the test of the spatial econometric model, and the traditional OLS regression model was directly used. The results have been shown in the table, so they are not listed repeatedly. (4) The spatial econometric models were selected after correlation tests such as LM, Wald, LR, and Hausman tests.

.

Column (1) is the spatial absolute β convergence test results of the LMLG region, where the β convergence coefficient is negative and passes the 1% significance test. This demonstrates that there is a spatial absolute β convergence characteristic of the CO₂ emissions from RRLC in this region. Similarly, columns (2) to (4) are spatial conditional β convergence tests for the LLLG, LMLG and UMLG regions, and it is not difficult to find that the conditional β convergence coefficients are always significantly negative at 1% level in these regions. This also indicates the existence of spatial conditional β convergence characteristics of CO₂ emissions from RRLC in these regions. Through the comparison of spatial absolute β convergence and spatial conditional β convergence tests, as it can be observed, upon the inclusion of pertinent variables, the magnitude of the estimated coefficient for spatial conditional β convergence is greater than that of spatial absolute β convergence in the LMLG regions. This illustrates the influence of related socio-economic and natural climatic factors expediting the rate of convergence in carbon emissions from RRLC. That is to say, the heterogeneity of socio-economic factors and natural climatic characteristics significantly influences spatial conditional β convergence.

5. Discussion

In the existing literature, both local and overseas studies on carbon emissions have achieved abundant achievements, but most of them concentrate on industrial sectors [28,29,30]. Insufficient attention has been paid to CO₂ emissions from RRLC. However, the contribution of CO₂ emissions from RRLC to the overall emissions should not be underestimated. Today’s global push to mitigate carbon emissions has penetrated all spheres.

Comparative analysis between 2000 and 2021 shows that the overall trend of carbon emissions from RRLC has been upward, accompanied by significant changes in spatial distribution patterns. Overall, CO₂ emissions from RRLC exhibited a trend of progressive growth, and high carbon emissions exhibited a spatial pattern characterized by “decentralized distribution–clustered distribution”. Consequently, reducing CO₂ emissions from RRLC has emerged as a key focus for future emission reduction initiatives. To this end, our findings present rich evidence on the regional differences and sources, dynamic evolutionary trends, and factors affecting the CO₂ emissions from RRLC in China in depth. More importantly, the methodology and focus of our study are different compared to previous studies. Our findings align with the research conclusion of Zhang and Li [25], namely that carbon CO₂ from RRLC exhibits regional heterogeneity characteristics and a spatial clustering pattern. This may be associated with the prevailing energy mix in rural areas in China, which is based on traditional forms of energy such as coal and diesel. Those generate large carbon emissions into the environment. Because of the limited supply of energy in rural areas, rural residents are limited in their choice of energy sources, making it difficult for them to use cleaner, low-carbon energy sources. The governments, rural residents, and relevant stakeholders must work together, create a virtuous cycle, and contribute to mitigating climate change and protecting the environment. Our findings are conducive to a more precise investigation into the current state and underlying causes of the imbalanced progression of carbon emissions from RRLC. This will accelerate the achievement of the “double carbon” goal and present novel solutions for sustainable human development.

However, it is also noticeable that several limitations remain in this article. First, the data on carbon emissions from RRLC are only obtained at the macro level, which cannot accurately reflect the micro and individual impact on RRLC. In coming studies, the evolution characteristics of carbon emissions from RRLC should be portrayed from microscopic perspectives, such as township enterprises, individuals, or households. Second, the pathways which influence the carbon emissions from RRLC are diverse, and it would be meaningful to estimate how these routes impact the environmental behavior of rural residents. Moreover, for the perspective of convergence, this study only explores the socio-economic and natural climatic factors affecting the carbon emissions from RRLC, such as the disposable income of rural residents, average years of education, urban–rural income gap, average annual relative humidity, annual precipitation, and annual sunshine, by using a convergence model. In future research, more empirical studies are needed to identify the elements influencing the CO₂ emissions from RRLC, such as green finance, digital economy, and high-quality development. Third, this paper considers the three main sources of carbon emissions from rural residents' direct living consumption include fossil energy, heat, and electricity. Therefore, some indirect consumption data, such as processing energy consumption, traffic, and manufacturing in rural areas, should be included in future studies.

6. Conclusions and Policy Implications

6.1. Conclusions

In this article, the carbon emissions from RRLC in the years 2000–2021 were measured. The regional disparities and sources were analyzed by using the Dagum Gini coefficient; then, kernel density estimation and the Markov chain method were used to characterize dynamic evolution; finally, conditional σ convergence and β convergence were applied to examine the convergence characteristics and driving factors of carbon emissions from RRLC. The findings were revealed the following:

First, the Dagum Gini coefficient was adopted to examine the spatial disparities and sources of carbon emissions from RRLC in China. In general, intra-regional and inter-regional differences overall exhibit a declining or fluctuating downward trend. The intensity of transvariation is the primary source of differences in carbon emissions from RRLC.

Second, the dynamic evolution tendency of CO₂ emissions from RRLC was validly evaluated based on kernel density estimation. The level of carbon emissions from RRLC in the nationwide and the four regions have generally gone upward. The distribution curves of the LLLG and LMLG regions show an increasingly steep trend. The nationwide and the four regions all show clear bimodal or multimodal peaks, indicating those regions the existence of the multi-polarization phenomenon.

Third, spatial Markov chains were further utilized to explore the internal evolution of carbon emissions from RRLC. The carbon emissions from RRLC are gradual, and it is hard to achieve a “leapfrog” transition. Generally, the probability of upward transfer is higher than the probability of a descending transition. After considering the space factors, it was found that the conclusion still holds.

Finally, convergence was used to inspect the evolutionary trajectory and drivers of carbon emissions from RRLC. There is $\sigma$ convergence and absolute β convergence for the national level and four regional levels of carbon emissions from RRLC. After considering the influence of socio-economic and natural climatic factors, a stable convergence trend was shown. Upon the inclusion of pertinent variables, the magnitude of the estimated coefficient for spatial conditional β convergence is greater than that of spatial absolute β convergence.

6.2. Policy Implications

On the basis of the above research results, this paper presents the following policy implications:

Firstly, it is important to acknowledge the high carbon emissions from RRLC in China. During the sample study period, there has been an overall increase in carbon emissions from RRLC. The improved living standards of rural residents have led to changes in their lifestyles and consumption habits, including the widespread adoption of household appliances. Moreover, rural residents predominantly depend on traditional energy sources, such as coal and oil, leading to an increase in carbon emissions from RRLC. However, there is enormous potential for energy conservation and emission reduction. In considering the “dual-carbon” strategy, China should prioritize various factors including rural residents’ disposable income, average education level, urban–rural income gap, average household size, total dependency ratio, and urbanization level. It is crucial to intensify efforts in reducing carbon emissions from RRLC.

Secondly, we should pay more attention to the uneven spatial distribution characteristics of carbon emissions from RRLC in China. The intensity of transvariation is the primary contributor to the overall differences in RRLC. In other words, the issue of cross-overlap between different regions has a significant effect on the carbon emissions from RRLC. Based on this, it is necessary to regulate the carbon emissions from RRLC and adopt a systemic perspective in planning socio-economic developments, regional development policies, and resource allocation strategies, to achieve a more equitable and sustainable development. Cooperation and synergistic governance between major regions should be enhanced, clearly delineating the key jurisdiction and responsibility for a reduction in carbon emissions from RRLC at a high speed, high intensity, and high quality. Furthermore, there is a need to identify the core factors for reducing carbon emissions from RRLC in distinct regions, explore the models of carbon emission reduction in different regions, and accelerate the transition to a low-carbon transformation and upgrading system.

Thirdly, it is important to formulate progressive regional development policies. Efforts should be made to gradually reduce carbon emissions from RRLC towards an optimal state, with prioritization given to rural areas that have the conditions to achieve significant progress such as Shanghai. Given the significant variations in social development, economic conditions, geographical locations, and resource endowments among different regions in China, when formulating regional development plans, it is crucial to consider the current levels of carbon emissions from RRLC, as well as the comparative advantages. It is essential to develop tailored regional strategies for carbon emission reduction based on local conditions, enhance market-oriented policies. This approach can play a fundamental role in unlocking the inherent potential for emission reduction in rural areas and effectively narrowing the regional disparities in carbon emissions from RRLC.