Cooperative Strategies in Transboundary Water Pollution Control: A Differential Game Approach

Abstract

This paper, based on differential game theory, examines governance models and cooperative strategies for managing cross-border water pollution in regions with uneven economic development. To address cross-regional water pollution, three differential game models are constructed under different scenarios: the Nash noncooperative mechanism, the pollution control cost compensation mechanism, and the collaborative cooperation mechanism. This study analyzes the dynamic changes in pollution emissions, governance investments, and economic returns within each model. The results indicate that the collaborative cooperation mechanism is the most effective, as it significantly reduces pollution emissions, maximizes overall regional benefits, and achieves Pareto optimality. In comparison, the pollution control cost compensation mechanism is suboptimal under certain conditions, while the Nash noncooperative mechanism is the least efficient, resulting in the highest pollution emissions. Furthermore, the research explores the influence of cooperation costs on the selection of governance models. It finds that high cooperation costs reduce local governments’ willingness to engage in collaborative cooperation. However, an appropriate compensation mechanism can effectively encourage less-developed regions to participate. Numerical analysis confirms the dynamic evolution of pollution stocks and economic returns under different models, and provides corresponding policy recommendations. This paper offers theoretical insights and practical guidance for cross-regional water pollution management, highlighting the importance of regional cooperation and cost-sharing in environmental governance.

1. Introduction

Watershed water pollution exhibits complex characteristics, including unidirectional, bidirectional, and cross-flow patterns, affecting multiple administrative regions. It presents transboundary, negative externalities, and spatial spillover effects. Rivers, as carriers of water resources, not only facilitate daily living for residents but also support industrial production. Consequently, both household and industrial wastewater often enter water bodies. With economic development and population growth, excessive water consumption and the cross-border transmission of pollutants in various regions have increasingly harmed the livelihoods and production of residents along riverbanks. The challenge of maintaining water quality standards at cross-border sections has severely hindered ecological restoration and balanced regional development [1]. The United Nations Environment Programme (UNEP) has warned that “many countries around the world are facing the threat of cross-regional water pollution, and the problem is worsening” [2]. In recent years, the accumulation of pollutants, the reduced self-purification capacity of rivers, and damage to vegetation—mainly due to industrial pollution and human activities—have become prominent issues. Given the challenges in resolving transboundary water pollution disputes, it is critical to implement robust water resource management systems and advance research on transboundary pollution compensation and cooperation within watersheds.

In this context, this paper addresses the issue of water pollution control between two regions with unequal economic development. To explore this, three regional pollution control scenarios are established: the Nash noncooperative mechanism, the pollution control cost compensation mechanism, and the collaborative cooperation mechanism. Differential game models are developed to conduct comparative analyses. Finally, numerical simulations are performed to validate the models, and management strategies are proposed based on the results, providing valuable insights for selecting regional water pollution control strategies.

In research on cooperative governance of water pollution, part of the literature focuses on the causes and nature of water pollution. Some scholars argue that the root cause of transboundary water pollution lies in conflicts of interest. Huisman et al. [3], through their analysis of water pollution control in the Rhine River, concluded that conflicting interests between countries are the primary reason for the failure of environmental governance agreements. List and Mason [4] contend that the key to solving transboundary water pollution issues lies in the “cost–benefit” dilemma. Governments often lack the awareness and incentives to manage transferable pollution. Only by finding viable solutions that consider cost–benefit factors can transboundary pollution problems be effectively resolved. Another stream of the literature focuses on pollution control, with research hotspots centered on pollution permit trading among enterprises and government–enterprise cooperation in managing water pollution. For example, Montgomery [5] and Eheart et al. [6] developed deterministic models of river pollution permit trading, analyzing the role of market access permits in controlling pollution. Building on this, Hung et al. [7] introduced a trading ratio system (TRS) for unidirectional river flow, demonstrating that TRS could account for emission location effects and achieve environmental quality standards at the lowest total abatement cost. Extending Hung’s work, Mesbah et al. [8] proposed an enhanced TRS as a cost-effective decision-making tool for managing river water quality in real time, while also accounting for risks. Although the aforementioned studies provide valuable insights into the causes of water pollution, research on cooperative pollution governance tends to focus on collaboration between enterprises or between governments and enterprises. Few studies address pollution compensation or collaborative control across administrative regions. Moreover, most of the literature on cooperative governance examines the effects or conditions of a single cooperation model, lacking in-depth analysis of different forms of cooperative governance.

Game theory, as a vital tool for analyzing conflicts and cooperation in human society, has been widely applied to the study of cross-regional water pollution governance. Common game models include evolutionary games, which are used to examine ecological compensation in water bodies [9] and address time consistency issues in downstream water pollution control [10]; cooperative games, such as those utilizing Monte Carlo numerical simulations for experimental research [11]; and dynamic games, such as river pollution models based on network externalities [12]. These studies encompass both traditional game methods and emerging theories like evolutionary games.

Differential games, a type of game in continuous-time systems, involve multiple participants who optimize their independent and conflicting objectives through ongoing interactions, eventually arriving at strategies that evolve over time and reach Nash equilibrium. This theory originated from U.S. Air Force research on military conflicts in the 1950s, combining optimal control and game theory, and has since been widely applied to areas such as joint emissions reduction [13], cooperative disaster relief [14], inventory management [15], product pricing [16], and cooperative governance [17]. In the context of dynamic pollution control, Huang et al. [18] proposed a differential game model for cross-border pollution between two regions; Xu et al. [19] used differential game models to examine the impact of the “River Chief System” on water pollution control; Wei et al. [20] analyzed how to balance sustainable development with water resource protection through ecological compensation; and Li et al. [21] investigated cross-border pollution using differential game models.

Early explorations into the application of game theory for pollution control laid a crucial theoretical groundwork in this area. Dockner and Long [22], along with Long [23], investigated international pollution control through the Cournot–Nash model, employing differential game techniques. Their findings underscored a key flaw in noncooperative approaches: when nations act independently, prioritizing economic benefits over environmental preservation, the result is often an increase in global emissions. This issue catalyzed further research aimed at coordinating international pollution control efforts to enhance emission reduction results.

Following this foundation, Wirl [24] examined the use of tax and permit systems within the Cournot–Nash framework, shedding additional light on the limitations of noncooperative methods. His study highlighted the challenges these strategies face in achieving effective pollution control.

In later work, Insley and Forsyth [25,26] utilized the Stackelberg model to analyze pollution control dynamics, introducing a leader–follower hierarchy as a means of achieving partial coordination in noncooperative settings. Compared with the Cournot–Nash approach, the Stackelberg model offers some reduction in the tension between individual and collective interests. Nevertheless, it still falls short of realizing globally optimal emission targets, especially in economically diverse regions.

Recent studies have further developed these theories by examining game structures under varying environmental policies. For instance, Nkuiya and Plantinga [27] investigated strategic pollution control within free trade conditions, concluding that noncooperative countries face challenges in maximizing global welfare through emissions reductions. Conversely, the Stackelberg strategy’s hierarchical leader–follower approach helps to mitigate some negative outcomes, though its effectiveness remains limited by the noncooperative nature of such frameworks.

In the context of cross-regional pollution management, balancing cooperative and noncooperative strategies presents significant difficulties. While the Stackelberg model provides a pathway for partial coordination, realizing optimal outcomes in cross-regional pollution management requires overcoming the inherent constraints of noncooperative models like Cournot–Nash. This is particularly challenging in economically disparate regions, where less-wealthy nations may “free-ride”, contributing minimally to pollution control efforts and exacerbating environmental burdens on neighboring areas.

To address these issues, cooperative approaches, such as cost-sharing mechanisms and collaborative governance structures, are critical. Such strategies enable economically diverse regions to jointly bear the costs and responsibilities of pollution control, thereby enhancing overall welfare.

To address the limitations of previous research, this paper analyzes the effects of production capacity, pollution control costs, and environmental pollution across regions, considering economic disparities. It explores three models: the Nash noncooperative mechanism, the pollution control cost-sharing mechanism, and collaborative cooperation. To conduct an in-depth analysis of these issues, this study employs “differential game” methods as a mathematical tool. Finally, by analyzing the dynamic trajectories of pollution levels and regional instantaneous returns, this paper examines the sensitivity of optimal strategies to factors such as short-sighted behavior and the natural degradation rate of the environment. Through numerical analysis, this paper discusses the dynamic changes in total pollutant levels and overall returns, along with the effects of discount rates and natural pollutant degradation rates on returns. This paper offers valuable theoretical references for regional cooperative strategies in water pollution governance.

Additionally, part of the inspiration for this paper stems from the work of Carbonaro and Menale (2024) [28], who investigated the application of Markov chains and kinetic theory in socioeconomic contexts. Their research on dynamic systems offers crucial insights for modeling decision-making processes in differential games. In particular, their focus on the evolution and optimization of strategies within changing environments closely aligns with the objectives of this study. By incorporating concepts from their research, this paper seeks to deepen the understanding of participant behavior in game psychology, especially in terms of how strategies adapt over time in response to various factors and uncertainties.

2. Problem Description and Model Construction

Assume that two adjacent regions, which are not administratively connected, are dealing with a cross-regional water pollution issue. These regions will be referred to as Region A and Region B. Assume that Region A has a higher economic level than Region B. The pollutants emitted by each region not only cause damage within their own region but also affect neighboring regions. These damages are represented by coefficients ${\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{\mathrm{g}}_{\mathrm{A}}(\mathrm{t}),{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}{\mathrm{g}}_{\mathrm{B}}(\mathrm{t}),{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}{\mathrm{g}}_{\mathrm{A}}(\mathrm{t}),$ and ${\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\mathrm{g}}_{\mathrm{B}}(\mathrm{t})$, where ${\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}>0$ and ${\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}>0$ indicate the coefficients of environmental damage caused by emissions within each region. Similarly, ${\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}>0$ and ${\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}>0$ represent the coefficients of damage caused to neighboring regions by emissions from each region.

Assume that each region produces a specific product, and each unit of this product generates a certain amount of atmospheric pollutants. According to Breton’s research [29], the revenue of each region can be described by the following concave function:

$${\mathrm{R}}_{\mathrm{i}}(\mathrm{t})={\mathrm{\alpha }}_{\mathrm{i}}{\mathrm{g}}_{\mathrm{i}}(\mathrm{t})-{\displaystyle \frac{1}{2}}{\mathrm{g}}_{\mathrm{i}}^2(\mathrm{t}), \mathrm{i}=\mathrm{A},\mathrm{B}$$

(1)

Here, the parameter $\mathrm{\alpha }$ is a positive constant that measures the disparity in production efficiency between the two regions. It is assumed that Region A is more developed than Region B, meaning ${\mathrm{\alpha }}_{\mathrm{A}}>{\mathrm{\alpha }}_{\mathrm{B}}$. Additionally, $0<\mathrm{\alpha }<1$.

The pollutants emitted by both regions can be treated using end-of-pipe technology, with the pollution control investment denoted as $\mathrm{h}$. The pollution control costs for Regions A and B are represented as:

$${\mathrm{C}}_{\mathrm{i}}({\mathrm{h}}_{\mathrm{i}}(\mathrm{t}))={\displaystyle \frac{1}{2}}{\mathrm{c}}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{i}}^2(\mathrm{t}), \mathrm{i}=\mathrm{A},\mathrm{B}$$

(2)

where ${\mathrm{c}}_{\mathrm{i}}$ is the pollution control cost coefficient, which measures the differences in pollution control technology between the two regions. Specifically, $1<{\mathrm{c}}_{\mathrm{A}}<{\mathrm{c}}_{\mathrm{B}}$, indicating that the pollution treatment cost in Region B is higher than in Region A.

To protect and improve the natural environment, it is assumed that the higher-level government (the environmental regulatory authority) imposes an environmental tax on pollutants for local governments. This environmental tax is denoted as $\mathrm{\tau }$.

Due to the mobility and accumulation characteristics of regional water pollution, each region not only bears the negative externalities of its own emissions but is also affected by the emissions of neighboring regions, resulting in a transregional pollution interaction problem. Additionally, pollutants continuously accumulate as emissions persist in each region. The pollutant stock $\mathrm{x}(\mathrm{t})$ is assumed to be represented by the following dynamic differential equation:

$$\stackrel{·}{\mathrm{x}}(\mathrm{t})={\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{g}}_{\mathrm{i}}(\mathrm{t})}-{\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{\omega }}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{i}}(\mathrm{t})}-\mathrm{\zeta }\mathrm{x}(\mathrm{t})$$

(3)

And $\mathrm{x}(0)={\mathrm{x}}_0,\mathrm{x}(\mathrm{t})\ge 0$. Here, $\mathrm{\zeta }>0$ represents the natural decomposition rate of pollutants, or the environment’s self-purification capacity, while $\mathrm{\omega }>0$ represents the amount of pollution reduced through local pollution control investments.

Assume that pollutant emissions cause pollution in both regions, leading to social losses. As a result, both regions must account for the negative effects of environmental damage in their decision making. Let the negative effect be represented by a linear function ${\mathrm{d}}_{\mathrm{i}}(\mathrm{x}),\mathrm{i}=\mathrm{A}、\mathrm{B}$, where ${\mathrm{d}}_{\mathrm{A}}=\mathrm{d},{\mathrm{d}}_{\mathrm{B}}=\mathrm{\beta }\mathrm{d},(0<\mathrm{\beta }<1)$, indicating that Region B experiences less environmental damage than Region A. While Region B has weaker economic development capabilities than Region A, its natural environment is in a better condition. For instance, residents in developed regions are more sensitive to environmental pollution, requiring the government to invest more resources to mitigate the negative impact of pollution.

In cross-regional water pollution governance, cooperation costs arise mainly from two sources: first, addressing the environmental damage caused by cross-border pollution; and second, the expenses related to pollution control investments and resource coordination.

First, consider that a portion of the cooperation cost is allocated to mitigating the environmental damage caused by cross-border pollution emissions. This cost is directly related to the pollution levels in each region. The higher the emissions, the more severe the environmental damage, leading to a significant increase in the measures required to reduce these emissions, thereby raising cooperation costs. Second, assume that both regions need to coordinate technologies and share resources during the governance process. Another part of the cooperation cost involves the total expenditure of both regions on governance investments, technical coordination, and facility construction.

By combining these two factors, the cooperation cost ${\mathrm{c}}_{\mathrm{C}}$ can be expressed as the total expenses incurred by the two regions for addressing cross-border pollution emissions and governance investments. The specific formula is

$${\mathrm{c}}_{\mathrm{C}}={\mathrm{k}}_1{\mathrm{g}}_{\mathrm{A}}+{\mathrm{k}}_2{\mathrm{g}}_{\mathrm{B}}+{\mathrm{k}}_3{\mathrm{h}}_{\mathrm{A}}+{\mathrm{k}}_4{\mathrm{h}}_{\mathrm{B}}$$

(4)

where the coefficients ${\mathrm{k}}_1$ and ${\mathrm{k}}_2$ represent the cost weights for each region in reducing pollution emissions, and ${\mathrm{k}}_3$, ${\mathrm{k}}_4$ reflect the cost weights related to governance investment coordination and technical sharing. By linking cooperation costs with pollution emission levels and governance investments, a more accurate assessment of the actual economic burden of cooperation and its impact on overall utility can be achieved.

3. Materials and Methods

3.1. Differential Game Method

This study utilizes differential game theory, a mathematical framework designed to model and analyze decision making in dynamic systems involving multiple agents with conflicting objectives. Specifically, the method was applied to formulate and solve problems where each agent’s strategy evolves over time. The Nash equilibrium was employed to identify optimal strategies, ensuring that no agent has an incentive to deviate from their chosen approach. This method facilitates the analysis of complex dynamic interactions and provides insights into optimal control within competitive scenarios.

3.2. AI Tools Disclosure

The author acknowledges the use of AI tools in preparing this manuscript. In particular, ChatGPT (GPT-4, https://www.openai.com, accessed on 15 October 2024) was used to enhance language clarity and accuracy by correcting grammar, refining sentence structures, and improving word choice. It is important to note that all academic content, research data, and analyses were conducted solely by the authors, with the AI tool serving only to improve the linguistic quality of the text.

4. Differential Game Model

4.1. Nash Noncooperative Game

In this scenario, Region A and Region B make decisions independently. Under the constraint of dynamic changes in pollution capacity, each region selects its optimal emission and pollution control paths to maximize its instantaneous net benefits. Therefore, in the noncooperative case, the utility functions and constraints for both regions are as follows:

Each local government seeks to maximize its own utility. The utility functions for the local governments in Region A and Region B are expressed as follows:

$${\mathrm{U}}_{\mathrm{A}}=\underset{{{\displaystyle \mathrm{g}}}_{\mathrm{A}}{\displaystyle (\mathrm{t})}}{\max }{\displaystyle {\int }_0^{\infty }{\mathrm{e}}^{-\mathrm{p}\mathrm{t}}}\left\{\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{g}}_{\mathrm{A}}(\mathrm{t})-{\displaystyle \frac{1}{2}}{\mathrm{g}}_{\mathrm{A}}^2(\mathrm{t})-{\displaystyle \frac{1}{2}}{\mathrm{c}}_{\mathrm{A}}{\mathrm{h}}_{\mathrm{A}}^2(\mathrm{t})-\mathrm{\tau }[{\mathrm{g}}_{\mathrm{A}}(\mathrm{t})-{\mathrm{\omega }}_{\mathrm{A}}{\mathrm{h}}_{\mathrm{A}}(\mathrm{t})]\\ -{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{\mathrm{g}}_{\mathrm{A}}(\mathrm{t})-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\mathrm{g}}_{\mathrm{B}}(\mathrm{t})-\mathrm{d}\mathrm{x}(\mathrm{t})\end{array}\right\}\mathrm{d}\mathrm{t}$$

(5)

$${\mathrm{U}}_{\mathrm{B}}=\underset{{{\displaystyle \mathrm{g}}}_{\mathrm{B}}{\displaystyle (\mathrm{t})}}{\max }{\displaystyle {\int }_0^{\infty }{\mathrm{e}}^{-\mathrm{p}\mathrm{t}}}\left\{\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{g}}_{\mathrm{B}}(\mathrm{t})-{\displaystyle \frac{1}{2}}{\mathrm{g}}_{\mathrm{B}}^2(\mathrm{t})-{\displaystyle \frac{1}{2}}{\mathrm{c}}_{\mathrm{B}}{\mathrm{h}}_{\mathrm{B}}^2(\mathrm{t})-\mathrm{\tau }[{\mathrm{g}}_{\mathrm{B}}(\mathrm{t})-{\mathrm{\omega }}_{\mathrm{B}}{\mathrm{h}}_{\mathrm{B}}(\mathrm{t})]\\ -{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}{\mathrm{g}}_{\mathrm{B}}(\mathrm{t})-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}{\mathrm{g}}_{\mathrm{A}}(\mathrm{t})-\mathrm{\beta }\mathrm{d}\mathrm{x}(\mathrm{t})\end{array}\right\}\mathrm{d}\mathrm{t}$$

(6)

$$\mathrm{s}.\mathrm{t}. \stackrel{·}{\mathrm{x}}(\mathrm{t})={\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{g}}_{\mathrm{i}}(\mathrm{t})}-{\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{\omega }}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{i}}(\mathrm{t})}-\mathrm{\zeta }\mathrm{x}(\mathrm{t})$$

(7)

Here, $\mathrm{p}>0$ represents the discount rate or discount factor.

Since the equilibrium of open-loop strategies cannot ensure time consistency and fails to reflect the dynamic changes in decision structures caused by pollution emissions, the optimal strategies derived from the feedback Nash equilibrium provide consistency in both state and time. Therefore, this paper addresses the Markov feedback Nash equilibrium of the differential game using the following Hamilton–Jacobi–Bellman (HJB) equation. For simplicity, the time variable $\mathrm{t}$ is omitted in the subsequent analysis.

$$\mathrm{\rho }{\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}(\mathrm{x})=\underset{{{\displaystyle \mathrm{g}}}_{\mathrm{A}}{{\displaystyle ,\mathrm{h}}}_{\mathrm{A}}}{\max }\left\{\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{g}}_{\mathrm{A}}-{\displaystyle \frac{1}{2}}{\mathrm{g}}_{\mathrm{A}}^2-{\displaystyle \frac{1}{2}}{\mathrm{c}}_{\mathrm{A}}{\mathrm{h}}_{\mathrm{A}}^2-\mathrm{\tau }\left({\mathrm{g}}_{\mathrm{A}}-{\mathrm{\omega }}_{\mathrm{A}}{\mathrm{h}}_{\mathrm{A}}\right)-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{\mathrm{g}}_{\mathrm{A}}\\ -{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\mathrm{g}}_{\mathrm{B}}-\mathrm{d}\mathrm{x}+{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}({\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{g}}_{\mathrm{i}}}-{\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{\omega }}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{i}}}-\mathrm{\zeta }\mathrm{x})\end{array}\right\}$$

(8)

$$\mathrm{\rho }{\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}(\mathrm{x})=\underset{{{\displaystyle \mathrm{g}}}_{\mathrm{B}}{{\displaystyle ,\mathrm{h}}}_{\mathrm{B}}}{\max }\left\{\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{g}}_{\mathrm{B}}-{\displaystyle \frac{1}{2}}{\mathrm{g}}_{\mathrm{B}}^2-{\displaystyle \frac{1}{2}}{\mathrm{c}}_{\mathrm{B}}{\mathrm{h}}_{\mathrm{B}}^2-\mathrm{\tau }\left({\mathrm{g}}_{\mathrm{B}}-{\mathrm{\omega }}_{\mathrm{B}}{\mathrm{h}}_{\mathrm{B}}\right)-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}{\mathrm{g}}_{\mathrm{B}}\\ -{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}{\mathrm{g}}_{\mathrm{A}}-\mathrm{\beta }\mathrm{d}\mathrm{x}+{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}}({\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{g}}_{\mathrm{i}}}-{\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{\omega }}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{i}}}-\mathrm{\zeta }\mathrm{x})\end{array}\right\}$$

(9)

The right side of the HJB equation above represents the maximization function. To determine its specific form, take the first-order partial derivatives with respect to the variables ${\mathrm{g}}_{\mathrm{A}}^{\mathrm{N}},{\mathrm{g}}_{\mathrm{B}}^{\mathrm{N}},{\mathrm{h}}_{\mathrm{A}}^{\mathrm{N}},{\mathrm{h}}_{\mathrm{B}}^{\mathrm{N}}$, and set these derivatives to zero. This yields the following maximization conditions:

$$\begin{array}{l}{\mathrm{g}}_{\mathrm{A}}^{\mathrm{N}}={\mathrm{\alpha }}_{\mathrm{A}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}\hspace{1em}\hspace{1em}{\mathrm{g}}_{\mathrm{B}}^{\mathrm{N}}={\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}}\\ {\mathrm{h}}_{\mathrm{A}}^{\mathrm{N}}={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{A}}\left(\mathrm{\tau }-\frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}\right)}{{\mathrm{c}}_{\mathrm{A}}}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathrm{h}}_{\mathrm{B}}^{\mathrm{N}}={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}\left(\mathrm{\tau }-\frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}\right)}{{\mathrm{c}}_{\mathrm{B}}}}\end{array}$$

(10)

Substitute Equation (10) into the HJB equation and simplify to obtain the following:

$$\mathrm{\rho }{\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}(\mathrm{x})=-\left(\mathrm{d}+\mathrm{\zeta }{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}\right)\mathrm{x}+\left\{\begin{array}{l}{\displaystyle \frac{\left({\mathrm{c}}_{\mathrm{A}}+{{\mathrm{\omega }}_{\mathrm{A}} }^2\right){\left(\frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}\right)}^2}{2{\mathrm{c}}_{\mathrm{A}}}}\\ +{\displaystyle \frac{1}{{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}}}{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}\left\{\begin{array}{l}{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{c}}_{\mathrm{A}}-{\mathrm{c}}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)-\mathrm{\tau }{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)\\ -{\mathrm{c}}_{\mathrm{A}}\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2+{\mathrm{c}}_{\mathrm{A}}\left({\mathrm{c}}_{\mathrm{B}}+{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\left({\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}}\right)}^2\end{array}\right\}\\ +{\displaystyle \frac{1}{{2\mathrm{c}}_{\mathrm{A}}}}\left\{\begin{array}{l}{\mathrm{c}}_{\mathrm{A}}\left(\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{A}}{}^2+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{}^2-2{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+2{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+{\mathrm{\tau }}^2\\ +2\mathrm{\tau }\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\right)-2{\mathrm{\alpha }}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+\mathrm{\tau }\right)\end{array}\right)\\ +{\mathrm{\tau }}^2{\mathrm{\omega }}_{\mathrm{A}}{}^2-2{\mathrm{c}}_{\mathrm{A}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}}\end{array}\right\}\end{array}\right\}$$

(11a)

$$\mathrm{\rho }{\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}(\mathrm{x})=-\left(\mathrm{\beta }\mathrm{d}+\mathrm{\zeta }{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}}\right)\mathrm{x}+\left\{\begin{array}{l}{\displaystyle \frac{\left({\mathrm{c}}_{\mathrm{B}}+{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\left(\frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}\right)}^2}{2{\mathrm{c}}_{\mathrm{B}}}}\\ +{\displaystyle \frac{1}{2}}\left\{\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{B}}{}^2-2\left({\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right){\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+2\mathrm{\tau }{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\tau }}^2{{\mathrm{\omega }}_{\mathrm{B}} }^2\\ -2{\mathrm{\alpha }}_{\mathrm{B}}\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)+{\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)}^2-2{\mathrm{c}}_{\mathrm{B}}{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}\end{array}\right\}\\ +{\displaystyle \frac{\frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}}{{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}}}\left\{\begin{array}{l}{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{c}}_{\mathrm{A}}-{\mathrm{c}}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)-\mathrm{\tau }{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)\\ -{\mathrm{c}}_{\mathrm{A}}\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2+{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{c}}_{\mathrm{A}}+{{\mathrm{\omega }}_{\mathrm{A}} }^2\right){\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}\end{array}\right\}\end{array}\right\}$$

(11b)

Based on the structure of Equation (11a,b), it can be concluded that the linear optimal value function of $\mathrm{X}$ is the solution to the HJB equation. Let

$${\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}(\mathrm{x})={\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}\mathrm{x}+{\mathrm{b}}_{\mathrm{A}}^{\mathrm{N}}, {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}(\mathrm{x})={\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}\mathrm{x}+{\mathrm{b}}_{\mathrm{B}}^{\mathrm{N}}$$

(12)

Here, ${\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}$, ${\mathrm{b}}_{\mathrm{A}}^{\mathrm{N}}$, ${\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}$, and ${\mathrm{b}}_{\mathrm{B}}^{\mathrm{N}}$ are constants to be determined. By substituting Equation (12) and its derivatives into Equation (11a,b), we obtain the following:

$$\mathrm{\rho }({\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}\mathrm{x}+{\mathrm{b}}_{\mathrm{A}}^{\mathrm{N}})=-\left(\mathrm{d}+\mathrm{\zeta }{\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}\right)\mathrm{x}+\left\{\begin{array}{l}{\displaystyle \frac{\left({\mathrm{c}}_{\mathrm{A}}+{{\mathrm{\omega }}_{\mathrm{A}} }^2\right){\left({\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}\right)}^2}{2{\mathrm{c}}_{\mathrm{A}}}}\\ +{\displaystyle \frac{{\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}}{{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}}}\left\{\begin{array}{l}{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{c}}_{\mathrm{A}}-{\mathrm{c}}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)-\mathrm{\tau }{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)\\ -{\mathrm{c}}_{\mathrm{A}}\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2+{\mathrm{c}}_{\mathrm{A}}\left({\mathrm{c}}_{\mathrm{B}}+{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\left({\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}\right)}^2\end{array}\right\}\\ +{\displaystyle \frac{1}{2}}\left\{\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{A}}{}^2+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{}^2-2{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+2{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\\ +2\mathrm{\tau }\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\right)+{\mathrm{\tau }}^2-2{\mathrm{\alpha }}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+\mathrm{\tau }\right)\end{array}\right\}\\ +{\displaystyle \frac{{\mathrm{\tau }}^2{\mathrm{\omega }}_{\mathrm{A}}{}^2-2{\mathrm{c}}_{\mathrm{A}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}}{2{\mathrm{c}}_{\mathrm{A}}}}\end{array}\right\}$$

(13a)

$$\mathrm{\rho }({\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}\mathrm{x}+{\mathrm{b}}_{\mathrm{B}}^{\mathrm{N}})=-\left(\mathrm{\beta }\mathrm{d}+\mathrm{\zeta }{\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}\right)\mathrm{x}+\left\{\begin{array}{l}{\displaystyle \frac{\left({\mathrm{c}}_{\mathrm{B}}+{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\left({\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}\right)}^2}{2{\mathrm{c}}_{\mathrm{B}}}}\\ +{\displaystyle \frac{1}{2}}\left\{\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{B}}{}^2-2\left({\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right){\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+2\mathrm{\tau }{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\tau }}^2{{\mathrm{\omega }}_{\mathrm{B}} }^2\\ -2{\mathrm{\alpha }}_{\mathrm{B}}\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)+{\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)}^2-2{\mathrm{c}}_{\mathrm{B}}{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}{\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}\end{array}\right\}\\ +{\displaystyle \frac{{\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}}{{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}}}\left\{\begin{array}{l}{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{c}}_{\mathrm{A}}-{\mathrm{c}}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)-\mathrm{\tau }{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)\\ -{\mathrm{c}}_{\mathrm{A}}\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2+{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{c}}_{\mathrm{A}}+{{\mathrm{\omega }}_{\mathrm{A}} }^2\right){\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}\end{array}\right\}\end{array}\right\}$$

(13b)

Based on the previous assumptions, Equations ${\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}(\mathrm{x})$ and ${\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}(\mathrm{x})$ must hold for all $\mathrm{x}\ge 0$. Consequently, the values of ${\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}$, ${\mathrm{b}}_{\mathrm{A}}^{\mathrm{N}}$ and ${\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}$, ${\mathrm{b}}_{\mathrm{B}}^{\mathrm{N}}$ can be determined as follows:

$${\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}=-{\displaystyle \frac{\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}$$

(14a)

$${\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}=-{\displaystyle \frac{\mathrm{\beta }\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}$$

(14b)

$${\mathrm{b}}_{\mathrm{A}}^{\mathrm{N}}=\left\{\begin{array}{l}{\displaystyle \frac{{\mathrm{c}}_{\mathrm{B}}{\left(\mathrm{d}+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)}^2{\mathrm{\omega }}_{\mathrm{A}}{}^2+2{\mathrm{c}}_{\mathrm{A}}\mathrm{d}\left(\mathrm{d}\mathrm{\beta }+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right){{\mathrm{\omega }}_{\mathrm{B}} }^2}{2\mathrm{\rho }{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{d}}^2\left(1+2\mathrm{\beta }\right)-2\mathrm{d}\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-\mathrm{\beta }{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-2\mathrm{\tau }\right)}{2\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\left({\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)}^2+2{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-{\mathrm{\alpha }}_{\mathrm{B}}\right)+2\mathrm{\tau }\left(-{\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\right)+{\mathrm{\tau }}^2}{2\mathrm{\rho }}}\end{array}\right\}$$

(15a)

$${\mathrm{b}}_{\mathrm{B}}^{\mathrm{N}}=\left\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}\mathrm{\beta }\left(\mathrm{d}+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right){{\mathrm{\omega }}_{\mathrm{A}} }^2}{\mathrm{\rho }{\mathrm{c}}_{\mathrm{A}}{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}+{\displaystyle \frac{+{\left(\mathrm{d}\mathrm{\beta }+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)}^2{{\mathrm{\omega }}_{\mathrm{B}} }^2}{2\mathrm{\rho }{\mathrm{c}}_{\mathrm{B}}{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{d}}^2\mathrm{\beta }\left(2+\mathrm{\beta }\right)+2\mathrm{d}\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\beta }\left(-{\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\alpha }}_{\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)\right)}{2\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{\alpha }}_{\mathrm{B}}{}^2+2\left(-{\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right){\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+2\mathrm{\tau }{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}-2{\mathrm{\alpha }}_{\mathrm{B}}{\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)}^2}{2\mathrm{\rho }}}\end{array}\right\}$$

(15b)

By substituting ${\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}$, ${\mathrm{b}}_{\mathrm{A}}^{\mathrm{N}}$, ${\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}$, and ${\mathrm{b}}_{\mathrm{B}}^{\mathrm{N}}$ into Equations (10) and (12), we derive Proposition 1:

Proposition 1.

Under the Nash equilibrium between the two local governments, the pollution emissions ${\mathrm{g}}_{\mathrm{A}}^{{\mathrm{N}}^{\ast }}$and ${\mathrm{g}}_{\mathrm{B}}^{{\mathrm{N}}^{\ast }}$, the pollution control investments ${\mathrm{h}}_{\mathrm{A}}^{{\mathrm{N}}^{\ast }}$and ${\mathrm{h}}_{\mathrm{B}}^{{\mathrm{N}}^{\ast }}$, and the total benefits ${\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})$ and ${\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})$ are as follows:

$${\mathrm{g}}_{\mathrm{A}}^{\mathrm{N}\ast }={\displaystyle \frac{\left({\mathrm{\alpha }}_{\mathrm{A}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}$$

(16a)

$${\mathrm{g}}_{\mathrm{B}}^{\mathrm{N}\ast }={\displaystyle \frac{\left({\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\mathrm{\beta }\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}$$

(16b)

$${\mathrm{h}}_{\mathrm{A}}^{\mathrm{N}\ast }={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{A}}\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\mathrm{d}\right]}{{\mathrm{c}}_{\mathrm{A}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}$$

(17a)

$${\mathrm{h}}_{\mathrm{B}}^{\mathrm{N}\ast }={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\mathrm{\beta }\mathrm{d}\right]}{{\mathrm{c}}_{\mathrm{B}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}$$

(17b)

$${\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})=-{\displaystyle \frac{\mathrm{d}\mathrm{x}}{\mathrm{\rho }+\mathrm{\zeta }}}+\left\{\begin{array}{l}{\displaystyle \frac{{\mathrm{c}}_{\mathrm{B}}{\left(\mathrm{d}+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)}^2{\mathrm{\omega }}_{\mathrm{A}}{}^2+2{\mathrm{c}}_{\mathrm{A}}\mathrm{d}\left(\mathrm{d}\mathrm{\beta }+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right){{\mathrm{\omega }}_{\mathrm{B}} }^2}{2\mathrm{\rho }{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{d}}^2\left(1+2\mathrm{\beta }\right)-2\mathrm{d}\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-\mathrm{\beta }{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-2\mathrm{\tau }\right)}{2\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\left({\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)}^2+2{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-{\mathrm{\alpha }}_{\mathrm{B}}\right)+2\mathrm{\tau }\left(-{\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\right)+{\mathrm{\tau }}^2}{2\mathrm{\rho }}}\end{array}\right\}$$

(18a)

$${\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})=-{\displaystyle \frac{\mathrm{\beta }\mathrm{dx}}{\mathrm{\rho }+\mathrm{\zeta }}}+\left\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}\mathrm{\beta }\left(\mathrm{d}+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right){{\mathrm{\omega }}_{\mathrm{A}} }^2}{\mathrm{\rho }{\mathrm{c}}_{\mathrm{A}}{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}+{\displaystyle \frac{+{\left(\mathrm{d}\mathrm{\beta }+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)}^2{{\mathrm{\omega }}_{\mathrm{B}} }^2}{2\mathrm{\rho }{\mathrm{c}}_{\mathrm{B}}{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{d}}^2\mathrm{\beta }\left(2+\mathrm{\beta }\right)+2\mathrm{d}\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\beta }\left(-{\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\alpha }}_{\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)\right)}{2\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{\alpha }}_{\mathrm{B}}{}^2+2\left(-{\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right){\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+2\mathrm{\tau }{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}-2{\mathrm{\alpha }}_{\mathrm{B}}{\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)}^2}{2\mathrm{\rho }}}\end{array}\right\}$$

(18b)

To ensure that emissions remain positive, conditions $\left({\mathrm{\alpha }}_{\mathrm{A}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\mathrm{d}>0$ and $\left({\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\mathrm{\beta }\mathrm{d} > 0$ must be satisfied.

From Equation (18a,b), the optimal benefit functions for local governments’ collaborative environmental governance can be derived.

$${\mathrm{V}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})={\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})+{\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})=(18\mathrm{a})+(18\mathrm{b})$$

In this model, the dynamic changes in pollutant stock over time $\mathrm{t}$ for the two regions are represented by ${\mathrm{x}}^{\mathrm{N}}(\mathrm{t})$:

$$\left\{\begin{array}{l}\stackrel{·}{\mathrm{x}}(\mathrm{t})={\mathrm{A}}_{\mathrm{N}}-\mathrm{\zeta }\mathrm{x}(\mathrm{t})\\ {\mathrm{x}}^{\mathrm{N}}(0)={\mathrm{x}}_0\\ {\mathrm{A}}_{\mathrm{N}}=\left\{\begin{array}{l}{\displaystyle \frac{\left({\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\zeta }-2\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\mathrm{d}\left(1+\mathrm{\beta }\right)}{\mathrm{\rho }+\mathrm{\zeta }}}\\ -{\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{A}}{}^2\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\mathrm{d}\right]}{{\mathrm{c}}_{\mathrm{A}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}-{\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}{}^2\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\mathrm{\beta }\mathrm{d}\right]}{{\mathrm{c}}_{\mathrm{B}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\end{array}\right\}\\ {\mathrm{x}}^{\mathrm{N}}(\mathrm{t})={\displaystyle \frac{{\mathrm{A}}_{\mathrm{N}}}{\mathrm{\zeta }}}+\left({\mathrm{x}}_0-{\displaystyle \frac{{\mathrm{A}}_{\mathrm{N}}}{\mathrm{\zeta }}}\right){\mathrm{e}}^{-\mathrm{\zeta }\mathrm{t}}\end{array}\right.$$

(19)

4.2. Stackelberg Game

In this scenario, since Region A possesses more advanced pollution control technology and experiences more severe pollution, a pollution control cost compensation mechanism is implemented to reduce total pollution capacity. Under this mechanism, the more developed Region A will take on part of the pollution control costs. As a result, the two independent governments make decisions following a Stackelberg leader–follower game structure: Region A first determines the pollution emission level, pollution control investment, and compensation rate for pollution control costs that maximize its own value function. Based on this, Region B selects the pollution emission level and control investment that maximize its total benefits. The Stackelberg feedback Nash equilibrium is solved through backward induction, starting with the decision-making process of the less-developed Region B. Assuming that Region A’s compensation rate for pollution control costs is $\mathrm{\theta }$.

The optimal value function ${\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}$ for Region B satisfies the following HJB equation:

$$\mathrm{\rho }{\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}(\mathrm{x})=\underset{{{\displaystyle \mathrm{g}}}_{\mathrm{B}}{{\displaystyle ,\mathrm{h}}}_{\mathrm{B}}}{\max }\left\{\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{g}}_{\mathrm{B}}-{\displaystyle \frac{1}{2}}{\mathrm{g}}_{\mathrm{B}}^2-{\displaystyle \frac{1}{2}}\left(1-\mathrm{\theta }\right){\mathrm{c}}_{\mathrm{B}}{\mathrm{h}}_{\mathrm{B}}^2-\mathrm{\tau }\left({\mathrm{g}}_{\mathrm{B}}-{\mathrm{\omega }}_{\mathrm{B}}{\mathrm{h}}_{\mathrm{B}}\right)-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}{\mathrm{g}}_{\mathrm{B}}\\ -{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}{\mathrm{g}}_{\mathrm{A}}-\mathrm{\beta }\mathrm{d}\mathrm{x}+{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}}{\partial \mathrm{x}}}({\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{g}}_{\mathrm{i}}}-{\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{\omega }}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{i}}}-\mathrm{\zeta }\mathrm{x})\end{array}\right\}$$

(20)

By taking the derivatives of Equation (20) with respect to pollution emissions and pollution control investment, and maximizing the first-order optimal conditions, we obtain the following:

$${\mathrm{g}}_{\mathrm{B}}={\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+{{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}}{\partial \mathrm{x}}}}_{,} {\mathrm{h}}_{\mathrm{B}}={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}\left(\mathrm{\tau }-\frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}}{\partial \mathrm{x}}\right)}{\left(1-\mathrm{\theta }\right){\mathrm{c}}_{\mathrm{B}}}}$$

(21)

Since Region A is more developed, it will rationally predict Region B’s choices and select its optimal pollution emissions and pollution control investment based on Region B’s reaction function. Substituting this reaction function into Region A’s HJB equation and simplifying, we obtain the following:

$$\mathrm{\rho }{\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}(\mathrm{x})=\underset{{{\displaystyle \mathrm{g}}}_{\mathrm{A}}{{\displaystyle ,\mathrm{h}}}_{\mathrm{A}}{\displaystyle ,\mathrm{\theta }}}{\max }\left\{\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{g}}_{\mathrm{A}}-{\displaystyle \frac{1}{2}}{\mathrm{g}}_{\mathrm{A}}^2-{\displaystyle \frac{1}{2}}{\mathrm{c}}_{\mathrm{A}}{\mathrm{h}}_{\mathrm{A}}^2-{\displaystyle \frac{1}{2}}\mathrm{\theta }{\mathrm{c}}_{\mathrm{B}}{\left[{\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}\left(\mathrm{\tau }-\frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}}{\partial \mathrm{x}}\right)}{\left(1-\mathrm{\theta }\right){\mathrm{c}}_{\mathrm{B}}}}\right]}^2\\ -\mathrm{\tau }\left({\mathrm{g}}_{\mathrm{A}}-{\mathrm{\omega }}_{\mathrm{A}}{\mathrm{h}}_{\mathrm{A}}\right)-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{\mathrm{g}}_{\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\left({\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}}{\partial \mathrm{x}}}\right)-\mathrm{d}\mathrm{x}\\ +{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}}{\partial \mathrm{x}}}\left[\begin{array}{l}{\mathrm{g}}_{\mathrm{A}}+\left({\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}}{\partial \mathrm{x}}}\right)-{\mathrm{\omega }}_{\mathrm{A}}{\mathrm{h}}_{\mathrm{A}}\\ -{\mathrm{\omega }}_{\mathrm{B}}{\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}\left(\mathrm{\tau }-\frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}}{\partial \mathrm{x}}\right)}{\left(1-\mathrm{\theta }\right){\mathrm{c}}_{\mathrm{B}}}}-\mathrm{\zeta }\mathrm{x}\end{array}\right]\end{array}\right\}$$

(22)

From this, the first-order optimal conditions for maximizing Region A’s pollution emissions, pollution control investments ${\mathrm{g}}_{\mathrm{A}}$ and ${\mathrm{h}}_{\mathrm{A}}$, and the pollution control cost compensation rate $\mathrm{\theta }$ are derived as follows:

$${\mathrm{g}}_{\mathrm{A}}={\mathrm{\alpha }}_{\mathrm{A}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}}{\partial \mathrm{x}}}}_{,} {\mathrm{h}}_{\mathrm{A}}={{\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{A}}\left(\mathrm{\tau }-\frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}}{\partial \mathrm{x}}\right)}{{\mathrm{c}}_{\mathrm{A}}}}}_{,} \mathrm{\theta }={\displaystyle \frac{2\frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}}{\partial \mathrm{x}}-\frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}}{\partial \mathrm{x}}+\mathrm{\tau }}{2\frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}}{\partial \mathrm{x}}-\mathrm{\tau }}}$$

(23)

Substituting the optimal decision conditions for both Region A and Region B into their respective HJB equations and simplifying, we obtain

$$\mathrm{\rho }{\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}(\mathrm{x})=-\left(\mathrm{d}+\mathrm{\zeta }{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}}{\partial \mathrm{x}}}\right)\mathrm{x}+\left\{\begin{array}{l}{\displaystyle \frac{{\mathrm{c}}_{\mathrm{B}}{\mathrm{\tau }}^2{{\mathrm{\omega }}_{\mathrm{A}} }^2}{2{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}}}+{\displaystyle \frac{\left(-2{\mathrm{c}}_{\mathrm{B}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+\left(2{\mathrm{c}}_{\mathrm{B}}+{{\mathrm{\omega }}_{\mathrm{B}} }^2\right)\frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}\right)\frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}}{2{\mathrm{c}}_{\mathrm{B}}}}\\ +{\displaystyle \frac{\frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}{2{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}}}\left\{\begin{array}{l}2{\mathrm{c}}_{\mathrm{B}}\left(\begin{array}{l}{\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{c}}_{\mathrm{A}}-\mathrm{\tau }{{\mathrm{\omega }}_{\mathrm{A}} }^2\\ -{\mathrm{c}}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)\end{array}\right)\\ -{\mathrm{c}}_{\mathrm{A}}\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2+\left({\mathrm{c}}_{\mathrm{B}}\left({\mathrm{c}}_{\mathrm{A}}+{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)+2{\mathrm{c}}_{\mathrm{A}}{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}\end{array}\right\}\\ +{\displaystyle \frac{1}{2}}\left({\mathrm{\alpha }}_{\mathrm{A}}{}^2+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{}^2+{\mathrm{\tau }}^2\right)-{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\\ +\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\right)\mathrm{\tau }-{\mathrm{\alpha }}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+\mathrm{\tau }\right)\end{array}\right\}$$

(24a)

$$\mathrm{\rho }{\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}(\mathrm{x})=-\left(\mathrm{\beta }\mathrm{d}+\mathrm{\zeta }{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}}{\partial \mathrm{x}}}\right)\mathrm{x}+\left\{\begin{array}{l}{\displaystyle \frac{1}{2}}\left\{{\mathrm{\alpha }}_{\mathrm{B}}{}^2+{\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)}^2\right\}-\left(\left({\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right){\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}\right)+\mathrm{\tau }{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}-{\mathrm{\alpha }}_{\mathrm{B}}\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)\\ +{\displaystyle \frac{1}{8{\mathrm{c}}_{\mathrm{B}}}}\left\{\begin{array}{l}3{\mathrm{\tau }}^2{\mathrm{\omega }}_{\mathrm{B}}{}^2+\left(4{\mathrm{c}}_{\mathrm{B}}+3{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\left({\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}}\right)}^2\\ -4{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}\left(2{\mathrm{c}}_{\mathrm{B}}{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2+{{\mathrm{\omega }}_{\mathrm{B}} }^2{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}\right)\end{array}\right\}\\ +{\displaystyle \frac{\frac{\partial {\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}}{\partial \mathrm{x}}}{{4\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}}}\left\{\left(\begin{array}{l}4{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{c}}_{\mathrm{A}}-{\mathrm{c}}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)-\mathrm{\tau }{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)\\ -3{\mathrm{c}}_{\mathrm{A}}\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2+2\left(2{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{c}}_{\mathrm{A}}+{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)+{\mathrm{c}}_{\mathrm{A}}{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}}{\partial \mathrm{x}}}\end{array}\right)\right\}\end{array}\right\}$$

(24b)

Similarly, based on the structure of Equation (24a,b), it can be deduced that the linear structure of $\mathrm{X}$ is the solution to the value function, i.e., let

$${\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}(\mathrm{x})={\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}\mathrm{x}+{\mathrm{b}}_{\mathrm{A}}^{\mathrm{S}}, {\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}(\mathrm{x})={\mathrm{k}}_{\mathrm{B}}^{\mathrm{S}}\mathrm{x}+{\mathrm{b}}_{\mathrm{B}}^{\mathrm{S}}$$

(25)

Here, ${\mathrm{k}}_{\mathrm{A}}^{\mathrm{N}}$, ${\mathrm{b}}_{\mathrm{A}}^{\mathrm{N}}$, ${\mathrm{k}}_{\mathrm{B}}^{\mathrm{N}}$, and ${\mathrm{b}}_{\mathrm{B}}^{\mathrm{N}}$ are constants to be determined. By substituting Equation (25) and its first-order derivative with respect to $\mathrm{X}$ into Equation (24a,b), we obtain the following:

$$\mathrm{\rho }({\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}\mathrm{x}+{\mathrm{b}}_{\mathrm{A}}^{\mathrm{S}})=-\left(\mathrm{d}+\mathrm{\zeta }{\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}\right)\mathrm{x}+\left\{\begin{array}{l}{\displaystyle \frac{{\mathrm{c}}_{\mathrm{B}}{\mathrm{\tau }}^2{{\mathrm{\omega }}_{\mathrm{A}} }^2}{2{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}}}+{\displaystyle \frac{\left(-2{\mathrm{c}}_{\mathrm{B}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+\left(2{\mathrm{c}}_{\mathrm{B}}+{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}\right){\mathrm{k}}_{\mathrm{B}}^{\mathrm{S}}}{2{\mathrm{c}}_{\mathrm{B}}}}-{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\\ +{\displaystyle \frac{{\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}}{2{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}}}\left\{\begin{array}{l}2{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{c}}_{\mathrm{A}}-{\mathrm{c}}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)-\mathrm{\tau }{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)\\ -{\mathrm{c}}_{\mathrm{A}}\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2+\left({\mathrm{c}}_{\mathrm{B}}\left({\mathrm{c}}_{\mathrm{A}}+{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)+2{\mathrm{c}}_{\mathrm{A}}{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}\end{array}\right\}\\ +{\displaystyle \frac{1}{2}}\left({\mathrm{\alpha }}_{\mathrm{A}}{}^2+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{}^2+{\mathrm{\tau }}^2\right)+\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\right)\mathrm{\tau }-{\mathrm{\alpha }}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+\mathrm{\tau }\right)\end{array}\right\}$$

(26a)

$$\mathrm{\rho }({\mathrm{k}}_{\mathrm{B}}^{\mathrm{S}}\mathrm{x}+{\mathrm{b}}_{\mathrm{B}}^{\mathrm{S}})=-\left(\mathrm{\beta }\mathrm{d}+\mathrm{\zeta }{\mathrm{k}}_{\mathrm{B}}^{\mathrm{S}}\right)\mathrm{x}+\left\{\begin{array}{l}{\displaystyle \frac{1}{2}}\left({\mathrm{\alpha }}_{\mathrm{B}}{}^2+{\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)}^2\right)-\left({\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right){\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\tau }{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}-{\mathrm{\alpha }}_{\mathrm{B}}\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)\\ {\displaystyle \frac{1}{8{\mathrm{c}}_{\mathrm{B}}}}\left\{3{\mathrm{\tau }}^2{\mathrm{\omega }}_{\mathrm{B}}{}^2-4{\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}\left(2{\mathrm{c}}_{\mathrm{B}}{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2+{\mathrm{\omega }}_{\mathrm{B}}{}^2\mathrm{k}_{\mathrm{A}}^{\mathrm{S}}\right)+\left(4{\mathrm{c}}_{\mathrm{B}}+3{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\left({\mathrm{k}}_{\mathrm{B}}^{\mathrm{S}}\right)}^2\right\}\\ +{\displaystyle \frac{{\mathrm{k}}_{\mathrm{B}}^{\mathrm{S}}}{{4\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}}}\left\{\left(\begin{array}{l}4{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{c}}_{\mathrm{A}}-{\mathrm{c}}_{\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)-\mathrm{\tau }{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)\\ -3{\mathrm{c}}_{\mathrm{A}}\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2+2\left(2{\mathrm{c}}_{\mathrm{B}}\left({\mathrm{c}}_{\mathrm{A}}+{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)+{\mathrm{c}}_{\mathrm{A}}{{\mathrm{\omega }}_{\mathrm{B}} }^2\right){\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}\end{array}\right)\right\}\end{array}\right\}$$

(26b)

Based on the previous assumptions, Equations ${\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}(\mathrm{x})$ and ${\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}(\mathrm{x})$ must hold for all $\mathrm{x}\ge 0$; thus, the values of ${\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}$, ${\mathrm{b}}_{\mathrm{A}}^{\mathrm{S}}$, ${\mathrm{k}}_{\mathrm{B}}^{\mathrm{S}}$, and ${\mathrm{b}}_{\mathrm{B}}^{\mathrm{S}}$ can be determined as follows:

$${\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}=-{\displaystyle \frac{\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}$$

(27a)

$${\mathrm{k}}_{\mathrm{B}}^{\mathrm{S}}=-{\displaystyle \frac{\mathrm{\beta }\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}$$

(27b)

$${\mathrm{b}}_{\mathrm{A}}^{\mathrm{S}}=\left\{\begin{array}{l}{\displaystyle \frac{2\left({\mathrm{d}}^2\left(1+2\mathrm{\beta }\right)-2\mathrm{d}\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-\mathrm{\beta }{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-2\mathrm{\tau }\right)\right)}{\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{2\left({\left({\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)}^2+2{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-{\mathrm{\alpha }}_{\mathrm{B}}\right)+2\mathrm{\tau }\left(-{\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\right)+{\mathrm{\tau }}^2\right)}{\mathrm{\rho }}}\\ +{\displaystyle \frac{4{\mathrm{c}}_{\mathrm{B}}{\left(\mathrm{d}+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)}^2{\mathrm{\omega }}_{\mathrm{A}}{}^2+{\mathrm{c}}_{\mathrm{A}}{\left(\mathrm{d}\left(2+\mathrm{\beta }\right)+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)}^2{{\mathrm{\omega }}_{\mathrm{B}} }^2}{2{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\end{array}\right\}$$

(28a)

$${\mathrm{b}}_{\mathrm{B}}^{\mathrm{S}}=\left\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}\mathrm{\beta }\left(\mathrm{d}+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right){{\mathrm{\omega }}_{\mathrm{A}} }^2}{{\mathrm{c}}_{\mathrm{A}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}+{\displaystyle \frac{\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2\left(\mathrm{d}\mathrm{\beta }+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)\left(\mathrm{d}\left(2+\mathrm{\beta }\right)+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\right)}{4{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}{}^2\left({\mathrm{d}}^2\mathrm{\beta }\left(2+\mathrm{\beta }\right)+2\mathrm{d}\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\beta }\left(-{\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\alpha }}_{\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)\right)\right)}{2\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{{\mathrm{\omega }}_{\mathrm{B}} }^2\left({\mathrm{\alpha }}_{\mathrm{B}}{}^2+2\left(-{\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right){\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+2{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}\mathrm{\tau }-2{\mathrm{\alpha }}_{\mathrm{B}}\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)+{\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)}^2\right)}{2\mathrm{\rho }}}\end{array}\right\}$$

(28b)

By substituting ${\mathrm{k}}_{\mathrm{A}}^{\mathrm{S}}$, ${\mathrm{b}}_{\mathrm{A}}^{\mathrm{S}}$, ${\mathrm{k}}_{\mathrm{B}}^{\mathrm{S}}$, and ${\mathrm{b}}_{\mathrm{B}}^{\mathrm{S}}$ into Equations (21) and (23), we derive Proposition 2:

Proposition 2.

Under the Nash equilibrium between the two local governments, the pollution emissions ${\mathrm{g}}_{\mathrm{A}}^{{\mathrm{S}}^{\ast }}$ and ${\mathrm{g}}_{\mathrm{B}}^{{\mathrm{S}}^{\ast }}$, pollution control investments ${\mathrm{h}}_{\mathrm{A}}^{{\mathrm{S}}^{\ast }}$and ${\mathrm{h}}_{\mathrm{B}}^{{\mathrm{S}}^{\ast }}$, pollution control cost compensation rate ${\mathrm{\theta }}^{\ast }$, and total benefits ${\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})$ and ${\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})$ are as follows:

$$\left\{\begin{array}{l}{\mathrm{g}}_{\mathrm{A}}^{\mathrm{S}\ast }={\displaystyle \frac{\left({\mathrm{\alpha }}_{\mathrm{A}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}\\ {\mathrm{g}}_{\mathrm{B}}^{\mathrm{S}\ast }={\displaystyle \frac{\left({\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\mathrm{\beta }\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}\end{array}\right.$$

(29)

$$\left\{\begin{array}{l}{\mathrm{h}}_{\mathrm{A}}^{\mathrm{S}\ast }={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{A}}\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\mathrm{d}\right]}{{\mathrm{c}}_{\mathrm{A}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\\ {\mathrm{h}}_{\mathrm{B}}^{\mathrm{S}\ast }={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}\left[\mathrm{d}\left(2+\mathrm{\beta }\right)+\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)\right]}{{2\mathrm{c}}_{\mathrm{B}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\end{array}\right.$$

(30)

$${\mathrm{\theta }}^{\ast }={\displaystyle \frac{\mathrm{d}\left(2-\mathrm{\beta }\right)-\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}{\mathrm{d}\left(2+\mathrm{\beta }\right)+\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}},\mathrm{d}>{\displaystyle \frac{\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}{2{\mathrm{c}}_{\mathrm{B}}-\mathrm{\beta }}}$$

(31)

$${\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})=-{\displaystyle \frac{\mathrm{d}\mathrm{x}}{\mathrm{\rho }+\mathrm{\zeta }}}+\left\{\begin{array}{l}{\displaystyle \frac{2\left({\mathrm{d}}^2\left(1+2\mathrm{\beta }\right)-2\mathrm{d}\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-\mathrm{\beta }{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-2\mathrm{\tau }\right)\right)}{\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{2\left({\left({\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)}^2+2{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-{\mathrm{\alpha }}_{\mathrm{B}}\right)+2\mathrm{\tau }\left(-{\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\right)+{\mathrm{\tau }}^2\right)}{\mathrm{\rho }}}\\ +{\displaystyle \frac{4{\mathrm{c}}_{\mathrm{B}}{\left(\mathrm{d}+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)}^2{\mathrm{\omega }}_{\mathrm{A}}{}^2+{\mathrm{c}}_{\mathrm{A}}{\left(\mathrm{d}\left(2+\mathrm{\beta }\right)+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)}^2{{\mathrm{\omega }}_{\mathrm{B}} }^2}{2{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\end{array}\right\}$$

(32a)

$${\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})=-{\displaystyle \frac{\mathrm{\beta }\mathrm{dx}}{\mathrm{\rho }+\mathrm{\zeta }}}+\left\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}\mathrm{\beta }\left(\mathrm{d}+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right){{\mathrm{\omega }}_{\mathrm{A}} }^2}{{\mathrm{c}}_{\mathrm{A}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}+{\displaystyle \frac{\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}{}^2\left(\mathrm{d}\mathrm{\beta }+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)\left(\mathrm{d}\left(2+\mathrm{\beta }\right)+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\right)}{4{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}{}^2\left({\mathrm{d}}^2\mathrm{\beta }\left(2+\mathrm{\beta }\right)+2\mathrm{d}\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\beta }\left(-{\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{\alpha }}_{\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+2\mathrm{\tau }\right)\right)\right)}{2\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{{\mathrm{\omega }}_{\mathrm{B}} }^2\left({\mathrm{\alpha }}_{\mathrm{B}}{}^2+2\left(-{\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right){\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+2{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}\mathrm{\tau }-2{\mathrm{\alpha }}_{\mathrm{B}}\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)+{\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)}^2\right)}{2\mathrm{\rho }}}\end{array}\right\}$$

(32b)

To ensure that emissions remain positive, conditions $\left({\mathrm{\alpha }}_{\mathrm{A}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\mathrm{d}>0$ and $\left({\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\mathrm{\beta }\mathrm{d} > 0$ must be satisfied.

From Equation (32a,b), the optimal benefit functions for local governments’ collaborative environmental governance can be derived.

$${\mathrm{V}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})={\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})+{\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})=(32\mathrm{a})+(32\mathrm{b})$$

In this model, the dynamic changes in pollutant stock over time $\mathrm{t}$ for the two regions are represented by ${\mathrm{x}}^{\mathrm{S}}(\mathrm{t})$:

$$\left\{\begin{array}{l}\stackrel{·}{\mathrm{x}}(\mathrm{t})={\mathrm{A}}_{\mathrm{S}}-\mathrm{\zeta }\mathrm{x}(\mathrm{t})\\ {\mathrm{x}}^{\mathrm{S}}(0)={\mathrm{x}}_0\\ {\mathrm{A}}_{\mathrm{S}}=\left\{\begin{array}{l}{\displaystyle \frac{\left({\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\zeta }-2\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\mathrm{d}\left(1+\mathrm{\beta }\right)}{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\\ -{\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{A}}{}^2\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\mathrm{d}\right]}{{\mathrm{c}}_{\mathrm{A}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}-{\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}{}^2\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\left(2+\mathrm{\beta }\right)\mathrm{d}\right]}{{2\mathrm{c}}_{\mathrm{B}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\end{array}\right\}\\ {\mathrm{x}}^{\mathrm{S}}(\mathrm{t})={\displaystyle \frac{{\mathrm{A}}_{\mathrm{S}}}{\mathrm{\zeta }}}+\left({\mathrm{x}}_0-{\displaystyle \frac{{\mathrm{A}}_{\mathrm{S}}}{\mathrm{\zeta }}}\right){\mathrm{e}}^{-\mathrm{\zeta }\mathrm{t}}\end{array}\right.$$

(33)

4.3. Cooperative Game

In this scenario, the two regions collaborate to form a unified ecological environment system, aiming to maximize the overall system benefit. By determining the optimal pollution emissions and pollution control investments, the entire system reaches its optimal state. Therefore, the total value function for both regions satisfies the following HJB equation:

$$\mathrm{\rho }{\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}(\mathrm{x})=\underset{{{\displaystyle \mathrm{g}}}_{\mathrm{A}}{{\displaystyle ,\mathrm{h}}}_{\mathrm{A}}{{\displaystyle ,\mathrm{g}}}_{\mathrm{B}}{{\displaystyle ,\mathrm{h}}}_{\mathrm{B}}}{\max }\left\{\begin{array}{l}-\left(1+\mathrm{\beta }\right)\mathrm{d}\mathrm{x}+{\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{g}}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{g}}_{\mathrm{B}}-{\displaystyle \frac{1}{2}}{\mathrm{g}}_{\mathrm{A}}^2-{\displaystyle \frac{1}{2}}{\mathrm{g}}_{\mathrm{B}}^2-{\displaystyle \frac{1}{2}}{\mathrm{c}}_{\mathrm{A}}{\mathrm{h}}_{\mathrm{A}}^2-{\displaystyle \frac{1}{2}}{\mathrm{c}}_{\mathrm{B}}{\mathrm{h}}_{\mathrm{B}}^2\\ -\mathrm{\tau }\left({\mathrm{g}}_{\mathrm{A}}+{\mathrm{g}}_{\mathrm{B}}-{\mathrm{\omega }}_{\mathrm{A}}{\mathrm{h}}_{\mathrm{A}}-{\mathrm{\omega }}_{\mathrm{B}}{\mathrm{h}}_{\mathrm{B}}\right)-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{\mathrm{g}}_{\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\mathrm{g}}_{\mathrm{B}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}{\mathrm{g}}_{\mathrm{B}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}{\mathrm{g}}_{\mathrm{A}}\\ +{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}({\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{g}}_{\mathrm{i}}}-{\displaystyle \sum _{\mathrm{i}=\mathrm{A},\mathrm{B}}{\mathrm{\omega }}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{i}}}-\mathrm{\zeta }\mathrm{x})-\left({\mathrm{k}}_1{\mathrm{g}}_{\mathrm{A}}+{\mathrm{k}}_2{\mathrm{g}}_{\mathrm{B}}+{\mathrm{k}}_3{\mathrm{h}}_{\mathrm{A}}+{\mathrm{k}}_4{\mathrm{h}}_{\mathrm{B}}\right)\end{array}\right\}$$

(34)

The right-hand side of this HJB equation represents the maximization function. To determine its specific form, take the first-order partial derivatives with respect to the variables ${\mathrm{g}}_{\mathrm{A}}^{\mathrm{C}}$, ${\mathrm{g}}_{\mathrm{B}}^{\mathrm{C}}$, ${\mathrm{h}}_{\mathrm{A}}^{\mathrm{C}}$, and ${\mathrm{h}}_{\mathrm{B}}^{\mathrm{C}}$, and set them to zero. This results in the following maximization conditions:

$$\begin{array}{l}{\mathrm{g}}_{\mathrm{A}}^{\mathrm{C}}={\mathrm{\alpha }}_{\mathrm{A}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}-{\mathrm{k}}_1\hspace{1em}\hspace{1em}{\mathrm{g}}_{\mathrm{B}}^{\mathrm{C}}={\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}-{\mathrm{k}}_2\\ {\mathrm{h}}_{\mathrm{A}}^{\mathrm{C}}={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{A}}\left(\mathrm{\tau }-\frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}\right)-{\mathrm{k}}_3}{{\mathrm{c}}_{\mathrm{A}}}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathrm{h}}_{\mathrm{B}}^{\mathrm{C}}={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}\left(\mathrm{\tau }-\frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}\right)-{\mathrm{k}}_4}{{\mathrm{c}}_{\mathrm{B}}}}\end{array}$$

(35)

Substituting Equation (35) into the HJB Equation (34) and simplifying, we obtain the following:

$$\mathrm{\rho }{\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}(\mathrm{x})=-\left(1+\mathrm{\beta }\right)\mathrm{d}\mathrm{x}+\left\{\begin{array}{l}{\displaystyle \frac{{\left({\mathrm{\alpha }}_{\mathrm{A}}-{\mathrm{k}}_1\right)}^2+{\left({\mathrm{\alpha }}_{\mathrm{B}}-{\mathrm{k}}_2\right)}^2+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{}^2+{{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}} }^2}{2}}+{\displaystyle \frac{{{\mathrm{k}}_3 }^2}{2{\mathrm{c}}_{\mathrm{A}}}}+{\displaystyle \frac{{{\mathrm{k}}_4 }^2}{2{\mathrm{c}}_{\mathrm{B}}}}\\ +{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+{\mathrm{\tau }}^2\\ +\mathrm{\tau }\left({\mathrm{k}}_1+{\mathrm{k}}_2-2{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)\\ -{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)+{\left({\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}\right)}^2\\ +{\mathrm{k}}_1\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}-{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}\right)+{\mathrm{k}}_2\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}-{\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}\right)\\ +{\mathrm{\alpha }}_{\mathrm{A}}\left({\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\tau }\right)+{\mathrm{\alpha }}_{\mathrm{B}}\left({\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-\mathrm{\tau }\right)\\ +{\displaystyle \frac{1}{2{\mathrm{c}}_{\mathrm{A}}}}\left\{\left({\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}-\mathrm{\tau }\right){\mathrm{\omega }}_{\mathrm{A}}\left(2{\mathrm{k}}_3+\left({\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}-\mathrm{\tau }\right){\mathrm{\omega }}_{\mathrm{A}}\right)\right\}\\ +{\displaystyle \frac{1}{2{\mathrm{c}}_{\mathrm{B}}}}\left\{\left({\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}-\mathrm{\tau }\right){\mathrm{\omega }}_{\mathrm{B}}\left(2{\mathrm{k}}_4+\left({\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}-\mathrm{\tau }\right){\mathrm{\omega }}_{\mathrm{B}}\right)\right\}\end{array}\right\}$$

(36)

Based on the structure of Equation (36), it can be concluded that the linear optimal value function of $\mathrm{X}$ is the solution to the HJB equation. Let

$${\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}(\mathrm{x})={\mathrm{k}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}\mathrm{x}+{\mathrm{b}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}$$

(37)

Substitute Equation (37) and its first-order derivative with respect to ${\displaystyle \frac{\partial {\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}}{\partial \mathrm{x}}}$ into Equation (36), and simplify to obtain the following:

$${\mathrm{k}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}=-{\displaystyle \frac{\left(1+\mathrm{\beta }\right)\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}$$

(38)

$${\mathrm{b}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}=\left\{\begin{array}{l}{\displaystyle \frac{{\mathrm{\alpha }}_{\mathrm{A}}{}^2+{\mathrm{\alpha }}_{\mathrm{B}}{}^2+\left({\mathrm{k}}_1+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)\left({\mathrm{k}}_1+2{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)+\left({\mathrm{k}}_2+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)\left({\mathrm{k}}_2+2{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)}{2\mathrm{\rho }}}\\ +{\displaystyle \frac{{\left(\mathrm{d}\left(1+\mathrm{\beta }\right){\mathrm{\omega }}_{\mathrm{A}}-\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{k}}_3-\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{A}}\right)\right)}^2}{2{\mathrm{c}}_{\mathrm{A}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}+{\displaystyle \frac{{\left(\mathrm{d}\left(1+\mathrm{\beta }\right){\mathrm{\omega }}_{\mathrm{B}}-\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{k}}_4-\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}\right)\right)}^2}{2{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{1}{\mathrm{\rho }}}\left\{\begin{array}{l}\mathrm{\tau }\left({\mathrm{k}}_1+{\mathrm{k}}_2+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)\\ -{\mathrm{\alpha }}_{\mathrm{A}}{}^2\left({\mathrm{k}}_1+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\tau }\right)-{\mathrm{\alpha }}_{\mathrm{B}}{}^2\left({\mathrm{k}}_2+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)\end{array}\right\}\\ +{\displaystyle \frac{{\mathrm{d}}^2{\left(1+\mathrm{\beta }\right)}^2-\mathrm{d}\left(1+\mathrm{\beta }\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}-{\mathrm{k}}_1-{\mathrm{k}}_2-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-2\mathrm{\tau }\right)}{\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\end{array}\right\}$$

(39)

By substituting ${\mathrm{k}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}$ and ${\mathrm{b}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}$ into Equations (35) and (37), we derive Proposition 3:

Proposition 3.

Under the Nash equilibrium between the two local governments, the pollution emissions ${\mathrm{g}}_{\mathrm{A}}^{{\mathrm{C}}^{\ast }}$and ${\mathrm{g}}_{\mathrm{B}}^{{\mathrm{C}}^{\ast }}$, pollution control investments ${\mathrm{h}}_{\mathrm{A}}^{{\mathrm{C}}^{\ast }}$and ${\mathrm{h}}_{\mathrm{B}}^{{\mathrm{C}}^{\ast }}$, and total benefits ${\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}{}^{\ast }(\mathrm{x})$ for the two regions are as follows:

$$\left\{\begin{array}{l}{\mathrm{g}}_{\mathrm{A}}^{{\mathrm{C}}^{\ast }}={\displaystyle \frac{\left({\mathrm{\alpha }}_{\mathrm{A}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-{\mathrm{k}}_1\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\left(1+\mathrm{\beta }\right)\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}\\ {\mathrm{g}}_{\mathrm{B}}^{{\mathrm{C}}^{\ast }}={\displaystyle \frac{\left({\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-{\mathrm{k}}_2\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-\left(1+\mathrm{\beta }\right)\mathrm{d}}{\mathrm{\rho }+\mathrm{\zeta }}}\end{array}\right.$$

(40)

$$\left\{\begin{array}{l}{\mathrm{h}}_{\mathrm{A}}^{\mathrm{C}\ast }={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{A}}\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\left(1+\mathrm{\beta }\right)\mathrm{d}\right]-{\mathrm{k}}_3\left(\mathrm{\rho }+\mathrm{\zeta }\right)}{{\mathrm{c}}_{\mathrm{A}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\\ {\mathrm{h}}_{\mathrm{B}}^{\mathrm{C}\ast }={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\left(1+\mathrm{\beta }\right)\mathrm{d}\right]-{\mathrm{k}}_4\left(\mathrm{\rho }+\mathrm{\zeta }\right)}{{\mathrm{c}}_{\mathrm{B}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\end{array}\right.$$

(41)

$${\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}{}^{\ast }(\mathrm{x})=-{\displaystyle \frac{\left(1+\mathrm{\beta }\right)\mathrm{d}\mathrm{x}}{\mathrm{\rho }+\mathrm{\zeta }}}+\left\{\begin{array}{l}{\displaystyle \frac{{\mathrm{d}}^2{\left(1+\mathrm{\beta }\right)}^2}{\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}+{\displaystyle \frac{{\left(\mathrm{d}\left(1+\mathrm{\beta }\right){\mathrm{\omega }}_{\mathrm{A}}-\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{k}}_3-\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{A}}\right)\right)}^2}{2{\mathrm{c}}_{\mathrm{A}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{\alpha }}_{\mathrm{A}}{}^2+{\mathrm{\alpha }}_{\mathrm{B}}{}^2+\left({\mathrm{k}}_1+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)\left({\mathrm{k}}_1+2{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}\right)+\left({\mathrm{k}}_2+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)\left({\mathrm{k}}_2+2{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}\right)}{2\mathrm{\rho }}}\\ +{\displaystyle \frac{1}{\mathrm{\rho }}}\left\{\begin{array}{l}\mathrm{\tau }\left({\mathrm{k}}_1+{\mathrm{k}}_2+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)\\ -{\mathrm{\alpha }}_{\mathrm{A}}{}^2\left({\mathrm{k}}_1+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\tau }\right)-{\mathrm{\alpha }}_{\mathrm{B}}{}^2\left({\mathrm{k}}_2+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)\end{array}\right\}\\ +{\displaystyle \frac{{\left(\mathrm{d}\left(1+\mathrm{\beta }\right){\mathrm{\omega }}_{\mathrm{B}}-\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{k}}_4-\mathrm{\tau }{\mathrm{\omega }}_{\mathrm{B}}\right)\right)}^2}{2{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{d}}^2{\left(1+\mathrm{\beta }\right)}^2-\mathrm{d}\left(1+\mathrm{\beta }\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)\left({\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}-{\mathrm{k}}_1-{\mathrm{k}}_2-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-2\mathrm{\tau }\right)}{\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\end{array}\right\}$$

(42)

In this model, the dynamic changes in pollutant stock over time $\mathrm{t}$ for the two regions are represented by ${\mathrm{x}}^{\mathrm{C}}(\mathrm{t})$:

$$\left\{\begin{array}{l}\stackrel{·}{\mathrm{x}}(\mathrm{t})={\mathrm{A}}_{\mathrm{C}}-\mathrm{\zeta }\mathrm{x}(\mathrm{t})\\ {\mathrm{x}}^{\mathrm{C}}(0)={\mathrm{x}}_0\\ {\mathrm{A}}_{\mathrm{C}}=\left\{\begin{array}{l}{\displaystyle \frac{\left({\mathrm{\alpha }}_{\mathrm{A}}+{\mathrm{\alpha }}_{\mathrm{B}}-\mathrm{\zeta }-2\mathrm{\tau }-{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}-{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}-{\mathrm{k}}_1--{\mathrm{k}}_2\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)-2\mathrm{d}\left(1+\mathrm{\beta }\right)}{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\\ -{\displaystyle \frac{\left({\mathrm{c}}_{\mathrm{A}}{\mathrm{\omega }}_{\mathrm{B}}{}^2+{\mathrm{c}}_{\mathrm{B}}{{\mathrm{\omega }}_{\mathrm{A}} }^2\right)\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\left(1+\mathrm{\beta }\right)\mathrm{d}\right]-{\mathrm{c}}_{\mathrm{A}}{\mathrm{k}}_4\left(\mathrm{\rho }+\mathrm{\zeta }\right)-{\mathrm{c}}_{\mathrm{B}}{\mathrm{k}}_3\left(\mathrm{\rho }+\mathrm{\zeta }\right)}{{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\end{array}\right\}\\ {\mathrm{x}}^{\mathrm{C}}(\mathrm{t})={\displaystyle \frac{{\mathrm{A}}_{\mathrm{C}}}{\mathrm{\zeta }}}+\left({\mathrm{x}}_0-{\displaystyle \frac{{\mathrm{A}}_{\mathrm{C}}}{\mathrm{\zeta }}}\right){\mathrm{e}}^{-\mathrm{\zeta }\mathrm{t}}\end{array}\right.$$

(43)

5. Comparison of Equilibrium Results

The following is a comparison of the pollution emissions between the two regions under three different scenarios. Based on the analysis of Propositions 1 through 3:

$${\mathrm{g}}_{\mathrm{A}}^{\mathrm{N}\ast }={\mathrm{g}}_{\mathrm{A}}^{\mathrm{S}\ast }>{\mathrm{g}}_{\mathrm{A}}^{\mathrm{C}\ast }, {\mathrm{g}}_{\mathrm{B}}^{\mathrm{N}\ast }={\mathrm{g}}_{\mathrm{B}}^{\mathrm{S}\ast }>{\mathrm{g}}_{\mathrm{B}}^{\mathrm{C}\ast }, {\mathrm{h}}_{\mathrm{A}}^{\mathrm{N}\ast }={\mathrm{h}}_{\mathrm{A}}^{\mathrm{S}\ast }<{\mathrm{h}}_{\mathrm{A}}^{\mathrm{C}\ast }, {\mathrm{h}}_{\mathrm{B}}^{\mathrm{C}\ast }>{\mathrm{h}}_{\mathrm{B}}^{\mathrm{S}\ast }>{\mathrm{h}}_{\mathrm{B}}^{\mathrm{N}\ast }$$

The results indicate that when the developed region (Region A) compensates the less-developed region (Region B) for pollution control costs, the total pollution emissions of both regions remain unchanged. The pollution control investment of the developed region (Region A) also remains the same. However, the pollution control investment of the less-developed region (Region B) increases.

Based on Equations (19), (33), and (43), the differences in pollution stock between the two regions under noncooperative and cooperative games are

$${\mathrm{x}}^{\mathrm{N}}(\mathrm{t})-{\mathrm{x}}^{\mathrm{C}}(\mathrm{t})={\displaystyle \frac{\mathrm{d}\left[{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}\left(\mathrm{\beta }+1\right)+{\mathrm{c}}_{\mathrm{B}}{\mathrm{\omega }}_{\mathrm{A}}{}^2\mathrm{\beta }+{\mathrm{c}}_{\mathrm{A}}{{\mathrm{\omega }}_{\mathrm{B}} }^2\right]}{{\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}>0$$

(44)

The difference in pollution stock between the two regions under the pollution control cost compensation game and the cooperative game is

$${\mathrm{x}}^{\mathrm{S}}(\mathrm{t})-{\mathrm{x}}^{\mathrm{C}}(\mathrm{t})={\displaystyle \frac{{2\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}\mathrm{d}\left(\mathrm{\beta }+1\right)+{2\mathrm{c}}_{\mathrm{B}}{\mathrm{\omega }}_{\mathrm{A}}{}^2\mathrm{\beta }\mathrm{d}+{\mathrm{c}}_{\mathrm{A}}{\mathrm{\omega }}_{\mathrm{B}}{}^2\left[\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)+\mathrm{d}\mathrm{\beta }\right]}{{2\mathrm{c}}_{\mathrm{A}}{\mathrm{c}}_{\mathrm{B}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}>0$$

(45)

The difference in pollution stock between the two regions under the noncooperative game and the pollution control cost compensation game is

$${\mathrm{x}}^{\mathrm{N}}(\mathrm{t})-{\mathrm{x}}^{\mathrm{S}}(\mathrm{t})={\displaystyle \frac{{\mathrm{\omega }}_{\mathrm{B}}{}^2\left[2\mathrm{d}-\mathrm{d}\mathrm{\beta }-\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)\right]}{{2\mathrm{c}}_{\mathrm{B}}\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}$$

(46)

From the above analysis, it is clear that when the two regions engage in collaborative cooperation, the pollution capacity reaches its lowest level. When $\mathrm{d}>{\displaystyle \frac{\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}{2-\mathrm{\beta }}}$, indicating that environmental damage exceeds a certain threshold, Region A compensates Region B for pollution control costs. As a result, the pollution capacity is lower than when the two regions act independently. Therefore, from an environmental perspective, joint pollution control represents the Pareto optimal state, while pollution control cost compensation is the second-best Pareto state. In contrast, the Nash noncooperative game produces the least efficient outcome.

Based on Equations (18a,b), (32a,b), and (42), we obtain the following:

(1)

The difference in Region A’s benefits between the pollution control cost compensation game and the noncooperative game is

$${\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})={\displaystyle \frac{{\left(\mathrm{d}\left(\mathrm{\beta }-2\right)+\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)\right)}^2{{\mathrm{\omega }}_{\mathrm{B}} }^2}{8{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}>0$$

(47)

(2)

The difference in Region B’s benefits between the pollution control cost compensation game and the noncooperative game is

$${\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})={\displaystyle \frac{\left(2\mathrm{d}-\mathrm{d}\mathrm{\beta }-\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)\right)\left(\mathrm{d}\mathrm{\beta }+\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)\right){{\mathrm{\omega }}_{\mathrm{B}} }^2}{4{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}>0$$

(48)

(3)

The difference in the benefits of both regions between the cooperative game and the noncooperative game is

$${\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})=\left\{\begin{array}{l}{\displaystyle \frac{{\mathrm{d}}^2\left(1+{\mathrm{\beta }}^2\right)+2\mathrm{d}\left({\mathrm{k}}_1+{\mathrm{k}}_2+{\mathrm{k}}_1\mathrm{\beta }+{\mathrm{k}}_2\mathrm{\beta }+\mathrm{\beta }{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\right)\left(\mathrm{\rho }+\mathrm{\zeta }\right)}{2\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{{\mathrm{k}}_3 }^2}{2{\mathrm{c}}_{\mathrm{A}}\mathrm{\rho }}}+{\displaystyle \frac{{{\mathrm{k}}_4 }^2}{2{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }}}+{\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\beta }}^2{{\mathrm{\omega }}_{\mathrm{A}} }^2}{2{\mathrm{c}}_{\mathrm{A}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}-{\displaystyle \frac{{\mathrm{d}}^2{{\mathrm{\omega }}_{\mathrm{B}} }^2}{2{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}-{\displaystyle \frac{{\mathrm{\alpha }}_{\mathrm{A}}{{\mathrm{k}}_1 }^2}{\mathrm{\rho }}}-{\displaystyle \frac{{\mathrm{\alpha }}_{\mathrm{B}}{\mathrm{k}}_2}{\mathrm{\rho }}}\\ +{\displaystyle \frac{{\mathrm{k}}_1{}^2+{\mathrm{k}}_2{}^2+2{\mathrm{k}}_1\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\tau }\right)+2{\mathrm{k}}_2\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)}{2\mathrm{\rho }}}\\ -{\displaystyle \frac{{\mathrm{k}}_3\left(\mathrm{d}\left(1+\mathrm{\beta }\right)+\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)\right){\mathrm{\omega }}_{\mathrm{A}}}{{\mathrm{c}}_{\mathrm{A}}\mathrm{\rho }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}-{\displaystyle \frac{{\mathrm{k}}_4\left(\mathrm{d}\left(1+\mathrm{\beta }\right)+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right){\mathrm{\omega }}_{\mathrm{B}}}{{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\end{array}\right\}$$

(49)

(4)

The difference in the benefits of both regions between the cooperative game and the pollution control cost compensation game is

$${\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})=\left\{\begin{array}{l}{\displaystyle \frac{{{\mathrm{k}}_3 }^2}{2{\mathrm{c}}_{\mathrm{A}}\mathrm{\rho }}}+{\displaystyle \frac{{{\mathrm{k}}_4 }^2}{{2\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }}}+{\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\beta }}^2{{\mathrm{\omega }}_{\mathrm{A}} }^2}{2{\mathrm{c}}_{\mathrm{A}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}+{\displaystyle \frac{{\left(\mathrm{d}\mathrm{\beta }+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right)}^2{{\mathrm{\omega }}_{\mathrm{B}} }^2}{8{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}\\ +{\displaystyle \frac{{\mathrm{d}}^2\left(1+{\mathrm{\beta }}^2\right)}{2\mathrm{\rho }{\left(\mathrm{\rho }+\mathrm{\zeta }\right)}^2}}-{\displaystyle \frac{2{\mathrm{\alpha }}_{\mathrm{A}}{\mathrm{k}}_1}{2\mathrm{\rho }}}-{\displaystyle \frac{{\mathrm{k}}_3\left(\mathrm{d}\left(1+\mathrm{\beta }\right)+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right){\mathrm{\omega }}_{\mathrm{A}}}{{\mathrm{c}}_{\mathrm{A}}\mathrm{\rho }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\\ +{\displaystyle \frac{{\mathrm{k}}_1{}^2+2{\mathrm{k}}_1\left({\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+\mathrm{\tau }\right)+{\mathrm{k}}_2\left({\mathrm{k}}_2+2\left({\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}+\mathrm{\tau }\right)-2{\mathrm{\alpha }}_{\mathrm{B}}\right)}{2\mathrm{\rho }}}\\ +{\displaystyle \frac{\mathrm{d}\left({\mathrm{k}}_1+{\mathrm{k}}_2+{\mathrm{k}}_1\mathrm{\beta }+{\mathrm{k}}_2\mathrm{\beta }+\mathrm{\beta }{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}+{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}\right)}{\mathrm{\rho }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\\ -{\displaystyle \frac{{\mathrm{k}}_4\left(\mathrm{d}\left(1+\mathrm{\beta }\right)+\left(\mathrm{\rho }+\mathrm{\zeta }\right)\mathrm{\tau }\right){\mathrm{\omega }}_{\mathrm{B}}}{{\mathrm{c}}_{\mathrm{B}}\mathrm{\rho }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}}\end{array}\right\}$$

(50)

The above analysis leads to Proposition 4:

Proposition 4:

 When the negative effects of pollution $\mathrm{d}$ exceed a certain threshold, i.e., $\mathrm{d}>{\displaystyle \frac{\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}{2-\mathrm{\beta }}}$, the environmental pollution capacity under the pollution control cost compensation mechanism is lower than that under the noncooperative mechanism, i.e., ${\mathrm{x}}^{\mathrm{N}}(\mathrm{t})-{\mathrm{x}}^{\mathrm{S}}(\mathrm{t})>0$. At the same time, the instantaneous returns of both regions are greater than under the noncooperative mechanism, i.e., ${\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})>0$  and ${\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})>0$. In the collaborative cooperation scenario, the environmental pollution capacity is better than in both the noncooperative mechanism and the pollution control cost compensation mechanism, i.e., ${\mathrm{x}}^{\mathrm{N}}(\mathrm{t})-{\mathrm{x}}^{\mathrm{C}}(\mathrm{t})>0$ and ${\mathrm{x}}^{\mathrm{S}}(\mathrm{t})-{\mathrm{x}}^{\mathrm{C}}(\mathrm{t})>0$. However, due to the existence of cooperation costs ${\mathrm{c}}_{\mathrm{C}}$, the governments of both regions will only choose collaborative cooperation if ${\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})>0$. Otherwise, due to cost constraints, local governments will still opt for the noncooperative approach based on their self-interests. Similarly, the governments will only adopt cooperative play if ${\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})>0$ is satisfied. Otherwise, they will revert to the cost compensation game mode. If both conditions ${\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{A}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{B}}^{\mathrm{N}}{}^{\ast }(\mathrm{x})<0$ and ${\mathrm{V}}_{\mathrm{A}\mathrm{B}}^{\mathrm{C}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{A}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})-{\mathrm{V}}_{\mathrm{B}}^{\mathrm{S}}{}^{\ast }(\mathrm{x})<0$ are met, the choice of game mode depends on the value of $\mathrm{d}$. When $\mathrm{d}>{\displaystyle \frac{\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}{2-\mathrm{\beta }}}$, the pollution control cost compensation mechanism is superior to the noncooperative mechanism.

From both environmental protection and economic benefit perspectives, regional cooperative governance of pollution represents the ideal Pareto optimal state. On one hand, regional cooperation in pollution control allows both parties to share the environmental burden, resulting in individual costs that are lower than social costs, thereby reducing pollution levels. On the other hand, collaborative cooperation prevents economic losses caused by “free-rider” behavior, enhances overall returns, and achieves Pareto improvements. However, cooperation costs have a significant influence on the feasibility of cooperative governance. These costs mainly include expenses related to addressing cross-border pollution, technical coordination, and resource integration. High cooperation costs may reduce the willingness of governments to engage in collaborative governance, especially when the costs are excessively high, leading them to prefer noncooperative governance models. Therefore, reducing cooperation costs—particularly by optimizing pollution control cost structures through technological innovation and resource sharing—is essential to achieving long-term sustainable cooperation. Compared to full collaboration, when pollution exceeds a certain threshold due to its mobility and accumulation, the developed region can reduce the pollution control costs of the less-developed region by providing ecological subsidies for pollution control technologies. This can result in a Pareto improvement in both environmental pollution levels and the instantaneous returns for both regions.

6. Numerical Analysis

In this section, the optimal decisions outlined above are analyzed numerically using Mathematica. Referring to the studies by Li [30] and Shoude [31], the benchmark parameters are set as ${\mathrm{\alpha }}_{\mathrm{A}}=140,{\mathrm{\alpha }}_{\mathrm{B}}=120,\mathrm{d}=5,\mathrm{\beta }=0.85,{\mathrm{c}}_{\mathrm{A}}=1,{\mathrm{c}}_{\mathrm{B}}=1.2,\mathrm{\rho }=0.1,\mathrm{\zeta }=0.2,\mathrm{\tau }=4,$ ${\mathrm{x}}_0=350,{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{A}}=0.1,{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{B}}=0.1,{\mathrm{\epsilon }}_{\mathrm{A}\mathrm{B}}=0.1,{\mathrm{\epsilon }}_{\mathrm{B}\mathrm{A}}=0.1,{\mathrm{\omega }}_{\mathrm{A}}=1,{\mathrm{\omega }}_{\mathrm{B}}=2,{\mathrm{k}}_1=1.8,{\mathrm{k}}_2=1.8,{\mathrm{k}}_3=2.2,$${\mathrm{k}}_4=2.2$. Based on this, a numerical analysis is conducted.

6.1. Optimal Pollution Capacity and Instantaneous Benefit Trajectory

Figure 1 illustrates the time evolution trajectories of pollutant capacity and benefits for the two regions under the feedback Nash equilibrium. Figure 1a,c show that the pollution capacity in both regions during collaborative cooperation is significantly lower than in the other two scenarios and gradually converges to a stable state over time. Additionally, the total benefits for both regions in collaborative cooperation are clearly superior to those in the other two mechanisms. For the noncooperative mechanism and the pollution control cost compensation mechanism, when the damage coefficient borne by Region A exceeds a certain threshold, $\mathrm{d}>{\displaystyle \frac{\mathrm{\tau }\left(\mathrm{\rho }+\mathrm{\zeta }\right)}{2-\mathrm{\beta }}}$, Region A chooses to share pollution control costs with Region B. Figure 1b,c show that when Region A provides pollution control cost compensation, the returns for both regions are higher than those under the noncooperative mechanism. However, when Region A cooperates with Region B in pollution control, the initial returns are lower compared to the cost compensation mechanism. Over time, though, the returns from cooperative pollution control gradually exceed those of the cost compensation mechanism.

6.2. Numerical Analysis in Steady State

In this section, we analyze the impact of factors such as the discount rate and the natural decomposition rate of pollutants on the steady-state equilibrium results.

Figure 2a–c illustrate the effect of the discount factor on the pollution capacity and benefits of the two regions. It can be observed that the discount rate is positively correlated with pollution capacity and negatively correlated with benefits. As the discount rate $\mathrm{\rho }$ increases, it indicates that both regions exhibit “short-sighted” behavior, with the slopes of the pollution capacity and benefits curves decreasing. This suggests that both regions focus more on maximizing short-term benefits and reduce their investment in pollution control technologies. Additionally, due to this “short-sighted” behavior, Region A reduces its pollution control cost compensation to Region B, causing the pollution capacity under the “cost compensation mechanism” to gradually approach that under the “noncooperative mechanism”.

Figure 3a–c illustrate the relationship between the natural decomposition rate of pollutants $\mathrm{\zeta }$, environmental pollution capacity, and the benefits of the two regions. When adopting the “noncooperative mechanism” and the “pollution control cost compensation mechanism”, an increase in the natural decomposition rate leads to a significant reduction in pollution capacity, thereby reducing environmental damage in both regions and increasing their benefits.

Furthermore, as shown in Figure 3b, with the rise in the natural decomposition rate, the slope of Region A’s benefit curve is steeper compared to Region B, and the benefit curves for the noncooperative mechanism and the cost compensation mechanism tend to converge. This suggests that Region A is more sensitive to environmental damage than Region B. As the natural decomposition rate $\mathrm{\zeta }$ increases, pollution capacity initially grows at a decreasing rate (concave growth) and then declines; correspondingly, the benefits of both regions first decrease and then increase, eventually stabilizing. Additionally, this study found that if the natural degradation rate of pollutants continues to rise, the total returns from cooperative pollution control will gradually decline and eventually fall below the returns of the other two models. This occurs because, as the ecological environment improves and pollutants can self-degrade at a higher rate, the need for cooperation diminishes.

When the natural decomposition rate $\mathrm{\zeta }$ is low, it implies that the ecological environment of the region is fragile and has poor recovery capability, making it more sensitive to environmental damage, and vice versa. Figure 3b also shows that the slope of Region A’s benefit curve is greater than that of Region B, indicating that Region A is more sensitive to environmental damage and, thus, requires higher pollution control costs, while Region B exhibits “free-rider” behavior. As the natural decomposition rate increases, the impact of environmental damage on Region A gradually diminishes, and the cooperation between the two regions becomes more stable.

7. Conclusions and Recommendations

7.1. Conclusions

This paper, based on a differential game model, investigated three governance models for controlling cross-border water pollution between two regions with asymmetric economic development: the Nash noncooperative mechanism, the pollution control cost compensation mechanism, and the collaborative cooperation mechanism. Comparative analysis of these models reveals significant differences in terms of pollution emissions, pollution control investments, economic returns, and cooperation costs.

The key research findings are as follows:

Nash noncooperative mechanism: Each region makes independent decisions, resulting in higher pollution emissions, poor governance outcomes, and lower total returns. Due to the absence of coordinated action, each region bears the full cost of governance, leading to reduced efficiency.

Pollution control cost compensation mechanism: The developed region takes on a larger share of the pollution control costs, which reduces overall pollution levels, while the less-developed region increases its pollution control investment. However, the overall impact remains limited. Despite the developed region’s greater investment, the less-developed region’s pollution control capacity does not reach optimal levels.

Collaborative cooperation mechanism: This is the optimal model. Cooperation not only reduces pollution emissions but also maximizes the total returns for both regions. Through resource sharing and technical coordination, cooperation costs are distributed equitably, leading to Pareto improvement. Collaborative cooperation effectively prevents the “free-rider” problem and significantly enhances long-term returns for both regions.

This study also highlights that cooperation costs significantly affect the feasibility of the cooperative mechanism. Local governments will only opt for collaborative cooperation if cooperation costs remain below a certain threshold. If these costs are too high, governments are more likely to choose the pollution control cost compensation mechanism or revert to the noncooperative model. Therefore, controlling cooperation costs is crucial for sustaining long-term cooperation.

7.2. Recommendations

Strengthen regional cooperation mechanisms and equitably distribute costs: Governments should establish stable cooperation frameworks to ensure that costs are fairly shared. Developed regions can provide technical and financial support to compensate for the pollution control investments of less-developed regions, thereby preventing a decline in cooperation due to high costs.

Optimize cooperation cost structures: Reducing cooperation costs through technological innovation and resource integration can minimize expenses related to technical coordination and infrastructure. Market mechanisms such as emissions trading systems can further optimize the allocation of pollution control costs.

Develop flexible compensation mechanisms: Developed regions should bear a larger share of cooperation costs, particularly in heavily polluted areas. Flexible compensation mechanisms can encourage less-developed regions to participate actively in governance and ensure effective control of cross-border water pollution. The compensation mechanism should be dynamically adjusted to prevent cooperation breakdowns caused by excessive costs.

Establish a long-term evaluation system: A system should be implemented to regularly assess cooperation costs and governance outcomes. Periodic evaluation of costs and environmental benefits will help adjust strategies to ensure the long-term success of the cooperation model.

Improve pollution control technologies and reduce cooperation costs: Promoting technological innovation in pollution control can enhance governance efficiency and reduce the costs of technical coordination. By fostering government–enterprise partnerships to develop new technologies, pollution control investments can be lowered, creating a win–win situation for both the environment and the economy.

By controlling and optimizing cooperation costs, regional cooperation mechanisms for cross-border water pollution governance can be stabilized, leading to effective pollution control and balanced economic development between regions.