# Carbon Market Simulation with Adaptive Mechanism Design

Han Wang1, Wenhao Li1, Hongyuan Zha1, Baoxiang Wang1,2
1The Chinese University of Hong Kong, Shenzhen
2Vector Institute
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A carbon market is a market-based tool that incentivizes economic agents to
align individual profits with the global utility, i.e., reducing carbon
emissions to tackle climate change. Cap and trade stands as a critical
principle based on allocating and trading carbon allowances (carbon emission
credit), enabling economic agents to follow planned emissions and penalizing
excess emissions. A central authority is responsible for introducing and
allocating those allowances in cap and trade. However, the complexity of
carbon market dynamics makes accurate simulation intractable, which in turn
hinders the design of effective allocation strategies. To address this, we
propose an adaptive mechanism design framework, simulating the market using
hierarchical, model-free multi-agent reinforcement learning (MARL). Government
agents allocate carbon credits, while enterprises engage in economic
activities and carbon trading. This framework illustrates agents’ behavior
comprehensively. Numerical results show MARL enables government agents to
balance productivity, equality, and carbon emissions. Our project is available
at https://github.com/xwanghan/Carbon-Simulator.

_K_ eywords Carbon Market, Adaptive Mechanism Design, Multi-Agent
Reinforcement Learning, Cap and Trade

##  1 Introduction

Climate change has emerged as a pressing worldwide concern [1], significantly
imperiling global ecosystems, economic systems, and sociopolitical stability.
The United Nations reports that in developing regions, one in ten individuals
subsists on less than US$ $1.90$ daily [2], with $2.2$ billion people deprived
of access to safely managed potable water resources [3]. The burgeoning
climate crisis amplifies these challenges, as worldwide temperature
escalations provoke droughts and rising sea levels, exacerbating famines and
enhanced forced displacements [4].

In $2016$, $196$ nations endorsed the Paris Agreement to mitigate climate
change collaboratively. However, pursuing transnational environmental targets
often conflicts with short-term interests, requiring mechanisms to reconcile
local and international objectives [5]. Carbon markets exemplify such
mechanisms [6, 7, 8], incentivizing economic agents to curb emissions. The cap
and trade format, predominant in carbon markets, involves allocating and
trading allowances [9, 10, 11, 12]. Economic agents must possess sufficient
allowances to offset emissions or face penalties for surplus emissions. The
cap and trade system sets a predetermined limit on allowances within an
economy, with a central authority introducing and allocating allowances based
on specified objectives. While this policy helps balance efficiency and
fairness, determining the optimal allocation remains challenging in general
economic contexts. The high-dimensional dynamics of the carbon market,
influenced by rational, self-interested, and far-sighted economic agents, lead
to market simulation reliance on models like CGE (computable general
equilibrium) [13, 14, 15] or ABM (agent-based modeling) [16, 17, 18]
frameworks, employing simplifying assumptions that are arduous to validate,
such as production and trading behaviors.

Given the unique nature of the carbon market, we integrate the AI Economist
[19, 20] to simulate market dynamics. Our adaptive mechanism design framework,
employing hierarchical, model-free MARL, mimics the carbon market. Lower-level
enterprise agents engage in realistic economic activities, such as emitting
carbon dioxide, trading emission credits, and investing in emission reduction
projects. Higher-level government agents analyze diverse allocation strategies
to achieve balanced efficiency and fairness, leading to significant carbon
emission reductions. The framework demonstrates the conduct of rational, self-
interested, and far-sighted agents within the carbon market. We emphasize that
our approach is not a simple transfer of the AI-Economist from taxation to
carbon credit allocation. Simulating the carbon market is challenging due to
limited data, fluctuating regulations, and non-market factors.

To validate our simulator, we conducted comparisons with several widely
adopted indicator allocation approaches at the firm level [11]. The simulation
results indicate reasonable action responses by enterprise agents to these
allocation policies. Additionally, numerical findings demonstrate that
government agents, through MARL, effectively discover allocation policies
capable of balancing productivity, equity, and carbon emissions. Our primary
contributions encompass: 1) We propose a systematic carbon market simulator
featuring carbon credits allocation and trading, and achieve realistic carbon
economy simulation based on hierarchical, model-free MARL. 2) We implement
several widely adopted indicator allocation approaches at the firm level as
baselines. 3) We observe that learning-based allocation policies possess the
potential to effectively balance productivity, equity, and carbon emissions.

Figure 1: Simulator structure. Left: One episode is divided into several
periods, and in each periods, the government firstly acts to allocate carbon
credits; the remaining of time, enterprises do their economics activities.
Right: Enterprises’ economics activities are modeled in a Gather-Trade-Build
game; in this grid map, they can produce properties (build) to get coins, can
move which can gather carbon credits and increase community’s total power
(green project investment), do carbon reduction invest action to reduce the
carbon emission level (carbon reduction investment), also trade carbon credits
and coins with each other.

##  2 Preliminaries and Related Works

###  2.1 Carbon Market and Cap and Trade

A carbon market can manifest at various scales, from local to global. This
paper primarily focuses on regional carbon markets, which involve the
participation of the government and many enterprises. Carbon market allowances
can be classified into mandatory and voluntary types, with the former being
strictly regulated by the government. At the same time, the latter encourages
enterprises to decrease their emissions through investment in carbon credits
generated by government-certified projects, facilitating emissions avoidance,
reduction, or removal from the atmosphere. The most prevalent form of carbon
markets is the cap and trade system, encompassing a predefined cap on total
allowances within the economy at a specific time. A central authority or
government manages the introduction of these allowances, with the cap value
being established to align with their objectives. This mechanism ensures an
accurately targeted level of emissions.

Carbon credits can be introduced into the market through $3$ primary means.
The initial method, termed free allocation, entails the government determining
the number of free credits to allocate to enterprises based on factors such as
the enterprise’s size and historical emissions. Another approach involves the
government selling allowances in an auction market. The final method
incorporates certified projects, enabling enterprises to garner additional
allowances by undertaking and completing energy conservation and emissions
reduction projects, typically facilitated by technological investment and
advancement. This paper exclusively focuses on free allocation and government-
certified projects.

In the context of free allocation, a pivotal decision confronting the
government is the establishment of both the overall volume of allowances for a
given period and the distribution of these allowances among enterprises within
distinct regions [21, 22, 23]. Allocation policies can be primarily divided
into four approaches: indicator, optimization, game theoretic, and hybrid
approaches. Each approach possesses strengths and weaknesses in terms of
productivity and equity. The indicator approach is the most frequently
employed in practice; however, consensus on indicator selection remains
elusive, as differing indicator methodologies can yield substantially varied
allocation outcomes and are challenging to validate [11]. In our simulations,
we utilized the indicator policy as a baseline to compare against the MARL-
based allocation policy.

The paramount component within the cap and trade is carbon trading [6, 24,
12]. Enterprises that exceed allotted emissions limits can elude penalties by
purchasing allowances from entities possessing surplus allowances. This paper
simulates carbon trading using a bid-and-ask auction mechanism, facilitating
allowance-based transactions between enterprises, analogous to the AI
Economist [19].

###  2.2 Carbon Market Simulated Models

The prevailing quantitative simulator encompass the computable general
equilibrium model(CGE) [13, 14, 15] and the agent-based model(ABM) [16, 17,
18, 25]. CGE models are economic frameworks that obtain empirical data as
input to emulate economic structures and the behavioral response of economic
entities as accurately as feasible to examine potential influences of varied
policies and other disturbances. CGE and dynamic CGE models have been employed
extensively to broadly explore climate policies, specifically, the carbon
market [26, 14]. Moreover, multi-agent-based models have been utilized to
replicate distinct industries’ behavior at a national scale [16, 25, 18]. [17]
propose a multi-agent simulation framework to emulate regional emissions
trading systems and assess numerous regulatory policies and carbon auction
regulations. Additionally, other multi-agent-oriented models have been applied
for exploring specific aspects of emissions trading schemes, such as
prospective carbon auction rules [27] or international transportation patterns
[28]. Rafieisakhaei et al. [29] introduces models targeting the EU ETS and the
global oil market, subsequently examining the interrelation between these
components.

###  2.3 Machine Learning for Economic Design

In the domain of economic design, contemporary investigations have chiefly
centered on automated mechanism design [30] and auction design [31], with the
primary aim being the augmentation of market efficacy and resource allocation
via the employment of machine learning methodologies. Historically, agents
within such research contexts have functioned as static entities devoid of
learning capabilities. However, recent scholarly contributions have begun to
incorporate MARL for the examination of security games [32] on Stackelberg
equilibrium, as well as the investigation of resource allocation games [33].
Furthermore, certain studies have extended the application of RL by enabling
agents to modify their conduct in adherence to novel regulatory frameworks,
thereby intensifying the refinement of market design [34, 35] and aligning
with the "Positronic Economist" paradigm [36]. Moreover, AI economists [37]
have successfully optimized taxation policies within a two-level simulation
environment, surpassing the effectiveness of extant real-world tax frameworks.

###  2.4 Multi-Agent Reinforcement Learning

In the simulated environment under investigation, numerous enterprise agents
adapt in response to distinct government policy interventions. These agents’
training process is grounded in the principles of MARL. Owing to the
complexity of environment state transitions, which encompass not only
stochastic environmental factors but also the actions of all agents involved,
a single agent’s perspective may inadvertently conflate environmental
randomness with the behavior of other agents, ultimately leading to improper
agent updates [38]. A growing body of recent research in the domain of MARL
has been devoted to improving multi-agent learning in non-stationary
environments [39], employing strategies such as reward sharing [40, 41] or
state information sharing [42, 43] amidst agents. This paper adopts an
alternative approach to bolster the learning process: sharing agent parameters
[40] while precluding the exchange of individual agent information.

##  3 Carbon Market Modeling

In this section, we will present an approach to conduct simulations of carbon
markets within the context of model-free MARL. The MARL framework exhibits a
hierarchical structure, consisting of the higher-level government RL agent and
lower-level enterprise RL agents. Consequently, it is referred to as a
hierarchical model-free MARL framework, which has also been known as a
manager-worker architecture [44, 45] in existing literature. Section 3.1 will
provide an in-depth mathematical representation of this framework, while
Sections 3.3 and 3.2 will elucidate the instantiated models for government and
enterprise RL agents, respectively. We use notations reference from RL and
economics, see Table 1.

###  3.1 Hierarchical Multi-Agent RL Framework

$t$ | timestep | $\pi_{g}$ | government allocation policy |
---|---|---|---|---
$i,j,k$ | enterprise agent indices | $p$ | punishment |
$g$ | government agent index | $z$ | income |
$\theta$, $\phi$ | model weights | $l$ | labor |
$h$ | hidden state | $c^{l}$ | labor-income coefficient |
$s$ | state | $u$ | utility |
$o$ | observation | $swf$ | social welfare |
$a$ | action | $prod$ | productivity |
$r$ | reward | $eq$ | equality |
$\pi$ | policy | $ee$ | excess carbon emissions |
$\gamma$ | discount factor | $c^{e}$ | economy-climate coefficient |
$\mathscr{T}$ | state transition |  |  |
x | enterprise attributes | $x^{q}$ | carbon emission credit |
$x^{c}$ | coin | $x^{l}$ | carbon emission level |

Table 1: Notations for our Hierarchical Multi-Agent RL Framework.

On the one hand, the issues encountered by higher-level government agents in
the carbon market (primarily referring to carbon allowances allocation within
this paper) can be modeled as a standard Markov decision process(MDP) [46]
$\mathcal{M}_{h}=(\mathcal{S},\mathcal{A},\mathcal{P},\mathcal{R},\gamma)$,
given the fixed lower-level enterprise agents. Given states
$s,s^{\prime}\in\mathcal{S}$ and an action $a\in\mathcal{A}$, the transition
probability function is expressed as $\mathcal{P}\left(s^{\prime}\mid
s,a\right):\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow[0,1]$ and
the reward function is defined by
$\mathcal{R}(s,a,s^{\prime}):\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\mathbb{R}$.
The discount factor is denoted as $\gamma\in(0,1]$ and the policy is
represented as $\pi:\mathcal{S}\times\mathcal{A}\rightarrow[0,1]$. For the
timestep $t\in[1,T]$, the discounted return is identified by
$G_{t}=\sum_{t^{\prime}=t}^{T}\gamma^{t^{\prime}-t}r_{t^{\prime}}$. The
objective of the government agent is to determine a policy $\pi$ that
maximizes $J=\mathbb{E}_{a_{t}\sim\pi(\cdot\mid
s_{t}),s_{t+1}\sim\mathcal{P}(\cdot\mid
s_{t},a_{t})}[\sum_{t=1}^{T}\gamma^{t-1}r_{t}(s_{t},a_{t},s_{t+1})]$.

On the other hand, the problems (production, trade, etc.) encountered by the
lower-level enterprise agents in carbon markets can formally defined by a
partially observable stochastic game(POSG) [47]
$\mathcal{G}_{l}:(\mathcal{S},\mathcal{A}_{i=1}^{|\mathcal{I}|},\mathcal{R}_{i=1}^{|\mathcal{I}|},\mathcal{P},\mathcal{O}_{i=1}^{|\mathcal{I}|},\mathcal{I})$,
given a fixed higher-level government agent. At each timestep $t\in T$, lower-
level enterprise agent $i\in\mathcal{I}$ obtains an observation,
$o_{t}^{i}=\mathcal{O}_{i}(s_{t})$, where $\mathcal{O}_{i}(s_{t})$ is emission
or observation function. Then each agent $i$ execute an action
$a_{t}^{i}\in\mathcal{A}_{i}$ simultaneously according to their historical-
dependent policy
$\pi_{i}:\mathcal{O}^{1}_{i}\times\mathcal{A}^{1}_{i}\times\cdots\times\mathcal{O}^{t}_{i}\times\mathcal{A}^{t}_{i}\rightarrow[0,1]$.
After that, environment will return rewards
$r_{t}^{i}\in\mathcal{R}_{i}(s_{t},\boldsymbol{a}_{t})$ to agent $i$, where
$\boldsymbol{a}_{t}$ denotes the joint action of all agents. Then state
$s_{t}$ will transitions to the next state $s_{t+1}$ according to transition
function $\mathcal{P}(s_{t+1}|s_{t},\boldsymbol{a}_{t})$. Each agent’s
objective is find an optimal policy $\pi$ to maximize the $\gamma$-discounted
expected return
$\max_{\pi_{i}}\mathbb{E}_{a^{i}\sim\pi_{i},\textbf{a}^{-i}\sim\pi_{-i},s^{\prime}\sim\mathcal{P}}\left[\sum_{t}\gamma^{t}r_{t}^{i}\right]$,
where $\pi_{-i}$ means the joint policy of other agents. Note that government
and enterprises are make decisions in different timescales.

Subsequently, we instantiate $\mathcal{M}_{h}$ and $\mathcal{G}_{l}$,
representing government and enterprises within the carbon market context, and
provide definitions of state spaces, action spaces, and reward functions.

###  3.2 Lower-Level: Enterprise Modeling

The Gather-Trade-Build game, as proposed in AI-Economist [20], is employed to
emulate enterprise economic behavior within the context of a carbon market.
This two-dimensional grid world enables agents to traverse the grid, acquire
resources, amass coins (representing profit) through resource usage in house
construction endeavors, and partake in trading activities with other agents by
exchanging resources for coins. More specifically, enterprises are
stochastically initialized on grid cells at the onset of each episode.
Subsequently, these entities initiate actions at every timestep throughout an
episode, excluding when the government intervenes (to be detailed later on).
Upon the completion of each episode, a penalty equal to $ee*p$ is imposed on
enterprises’ coin holdings in response to surpassing their allotted carbon
emission credits.

#### Reward Function.

The objective of each enterprise is to maximize its utility which is defined
as:

$\displaystyle
u(z,l)=\frac{z^{1-\eta}-1}{1-\eta}-c^{l}*l,\quad\quad\eta=0.23,$ (1)

where $z$ is income that the sum of enterprise‘s coins, and $l$ is also the
cumulative labor during all previous time steps. And for $\eta$ is isoelastic
coefficient, we use isoelastic utility function [48] to model the income
component of reward, it is a nonlinear curve with higher $\eta$ lower
sensitivity to the income change. During an episode, labor-income coefficient
$c^{l}$ is changing, it follows equation

$\displaystyle\alpha*(1-e^{t/\beta}),$ (2)

where $\alpha,\beta$ is constant and $t$ refers to time.

#### Skill.

Enterprise agents’ skills include $S$ and $Rc$, representing enterprise size
and research and development (R&D) capabilities.

#### Action Space.

An episode is divided into several periods. At the start of each period, each
enterprise is allocated a certain number of carbon emission credits by the
government agent. At each timestep within a period, enterprises must choose
from a set of actions, including Move, Produce, Invest, Trade or NO-OP (no
operation).

Trade action consumes one unit of labor. Enterprise agents can move to any of
the four neighboring grids. Still, this action will be blocked by the edge of
the grid world or other enterprises’ properties (or houses). Suppose an
enterprise agent moves to a grid already allocated a government-certified
project. The enterprise will consume coins and labor and receive carbon
emission credits. Once completed, the project will remain on the grid and
influence the carbon emission process of all agents.

Produce action consumes one unit of labor, $x^{l}$ carbon emission credits,
and provides the agent with $10S$ coins. This action is available if the
current position in the grid environment is not occupied by the enterprise’s
property or a government-certified project. Additionally, when a large amount
of carbon emission is created, it will randomly pollute the value of nearby
blank grids. This will introduce a discount on the coins received when the
enterprise produces in such grids.

Invest action can decrease an enterprise’s carbon emission level. This action
is available when an enterprise has positive coins. It consumes one unit of
labor and $5/Rc$ coins. The carbon emission level, denoted as $x^{l}$, is a
continuously changing enterprise attribute, with a floating range between $0$
and $1$. At the beginning of the episode, all enterprises have their $x^{l}$
set to $1$. Through investments, an enterprise can change its $x^{l}$. The
carbon emission level is defined as follows:

$\displaystyle x^{l}_{i}=Pc_{i}*(1-Gr),$ (3)

where $Pc_{i}$ is defined as $Pc_{i}=exp(-\delta*Rc_{i}*n_{i}^{r})$ and $Gr$
is:

$\displaystyle Gr=\frac{n^{p}}{\sum_{j}^{I}(Pc_{j}*S_{j})+n^{p}}.$ (4)

We use power consumption $Pc$ to represent the power consumed by an enterprise
when producing $10$ coins. This value depends on the number of investments
$n_{i}^{r}$, the enterprise’s R & D capabilities, and constant $\delta$. The
green rate $Gr$ is used to represent the proportion of green energy generated
from government-certified projects in total energy. It decreases with a higher
number of completed government-certified projects located in the grid world,
denoted as $n^{p}$, or a lower total power consumption
$\sum_{j}^{I}(Pc_{j}*S)$. Additionally, we have set lower bounds to ensure
that $x^{l}_{i}$ remains within a reasonable range and $Gr$ is non-negative.
We have introduced delay, forgetting, and random failures in the count of
investment actions $n_{i}^{r}$ to model actual carbon emission reduction
investments.

Trade encompasses bidding and asking in carbon trading, similar to AI
Economist. The action space for Trade is twice the number of specified price
levels. When an enterprise publishes a bid or ask request, that request
remains active for a predetermined number of timesteps until it is matched
with a higher-priced ask or a lower-priced bid. The range for the lifetime of
a request should be smaller than the length of an episode because, at the
start of each new episode, market requests are cleared. Publishing a request
in the market incurs a lower labor cost than other actions, encouraging trade
activity. However, this action is subject to limitations. Enterprises can only
publish a maximum allowed number of requests, and they must have a sufficient
amount of coins and carbon emission credit to participate in trading.

#### Observation Space.

The range of observation for enterprises is significantly smaller than that of
the government. Regarding position information, enterprises can gather
information about their properties’ positions and a limited square range of
the grid world. Within this limited square range, enterprises can observe
which grid areas are polluted by carbon emissions resulting from their
production actions or purified by completing government-certified projects.
Enterprises will only know their skills and attributes but can access
information about the average carbon emission level. In carbon trading,
enterprises can observe the embedded price history and all published requests,
similar to what the government can monitor. Additionally, enterprises’
observations will include self-published requests.

###  3.3 Higher-Level: Government Modeling

Reward Function. Government’s objective is to maximize social welfare, which
is defines as:

$\displaystyle swf=prod*eq*exp(-c^{e}*ee),$ (5)

where $prod$ is defined as:

$\displaystyle prod=\sum_{i}(x_{i}^{c}),$ (6)

and $eq$ defined from concept of $gini$ index:

$\displaystyle
eq=1-\frac{\sum_{i=1}^{N}\sum_{j=1}^{N}|x_{ci}-x_{cj}|}{2(N-1)\sum_{i=1}^{N}x_{ci}},\quad\quad
0<eq<1.$ (7)

#### Action Space.

The government agent makes decisions at the first timestep of each period. Its
action is represented by a multi-dimensional discrete vector. The first
dimension specifies the proportion of the total credit for this period that
will be allocated to the remaining unallocated carbon emission credits for all
periods. The next $|\mathcal{I}|$ dimensions, where $|\mathcal{I}|$ refers to
the number of enterprise agents, determine the allocation weight of carbon
emission credits for each enterprise agent in each period.

The total credit available for the period is calculated based on the specified
proportion of the period’s full credit. The government then allocates 10% of
this credit to a government-certified project randomly placed on an empty
grid. The remaining 90% of the total credit is distributed to enterprises
according to their respective allocation weights.

Each dimension of the action space has a size of 101, with the remaining
action 0 used for masking at the rest of the timesteps in each period.

To ensure the proper functioning of a cap and trade, it is essential for the
government agent not only to allocate carbon credits to each period and
enterprise agent but also to establish the severity of penalties for exceeding
emissions. Thus, the $|\mathcal{I}|+2$-th dimension determines the severity of
punishments.

Regarding different facets of government policy, the parameter that exerts the
most significant influence on enterprise behavior is the punishment setting,
as it directly impacts the utility of enterprises. When the punishment is
excessively high, enterprises will avoid exceeding the carbon emission credit
at all costs, negatively affecting productivity and trade. Conversely, suppose
the punishment is set too low. In that case, production actions will go
unchecked, resulting in low equality and exceedingly high carbon emissions.
Therefore, we have maintained a constant punishment for most of our
experiments. After each period, the carbon emission credits allocated to
enterprises are reset while the government project remains intact.

#### Observation Space.

The government possesses access to comprehensive information. Firstly, it can
observe the entire grid world, including the positions of all enterprise
agents, their properties’ locations, and the government-certified project’s
location. Next, concerning the private information of enterprises, the
government can collect data on their skills and attributes. Additionally,
real-time market information, including embedded price history and all
published requests, is readily available. Furthermore, the government has
access to information related to carbon emission credits, such as excess
records and the total remaining credit.

#### Remark:

Notably, in developing a carbon market model similar to AI-Economist, we opt
for a selection of hyperparameters without referencing empirical data, relying
solely on fundamental economic principles. This endows the simulation model
with several advantages. Firstly, government and enterprise agents possess no
prior knowledge of external markets or economic theory and lack any
understanding of each other’s behavioral patterns. Consequently, the simulator
can optimize arbitrary social outcomes through hyperparameter configuration.
Secondly, carbon market data is scarce and typically encompasses relatively
small market scales. Our proposed simulator, by using appropriate
hyperparameter settings, transcends the limitations of real-world data,
enabling the examination of allocation strategies and the impact of economic
activities under varying carbon market scales. Lastly, the simulator can
validate the rationality of different economic theories by adjusting
hyperparameters following various theoretical perspectives.

Figure 2: Quantitative results of different allocation policies.

##  4 Enterprise Behavior Simulation

Upon the completion of carbon market modeling, we can employ (MA)RL to train
the government (enterprise) agent(s), enabling them to exhibit behavior that
closely resembles real-world scenarios, thereby achieving a realistic
simulation of the carbon market. This section primarily focuses on the
simulation of the critical entities in the carbon market, namely the behavior
of enterprise agents. Based on reward function we defined in Equation (1), we
can obtain desired policy $\pi_{i}$ by maximizing total discounted reward:

$\displaystyle\max_{\pi_{i}}\mathbb{E}_{a^{i}\sim\pi_{i},\textbf{a}^{-i}\sim\pi_{-i},s^{\prime}\sim\mathscr{T}}\left[\sum_{t}^{T}\gamma^{t}(u_{i,t}-u_{i,t-1})+u_{i,0}\right].$
(8)

Like AI Economist, we use a strong baseline in MARL, the independent PPO [49],
to train adaptive response enterprise agents for government policy. We
incorporate $5$ enterprise agents into the carbon market in the simulation.

###  4.1 Training Details

Figure 3: Simulation details for MARL and other approaches with SI, which
include the social welfare, cumulative carbon emissions, the total number of
trade activities, properties, government projects, and the investments.

To enhance agents’ performance, the observation of enterprises includes
position information, which can be represented as an image. To process this
position information, a convolutional neural network (CNN) is employed. Other
information is added to the CNN’s output, and the combined data is passed
through a pipeline consisting of fully connected (FC) layers, a long short-
term memory (LSTM) network, and additional FC layers. This process yields
action logits. Agents’ policies are obtained by applying an action mask and
performing a softmax operation on these logits.

MARL training aims to discover an allocation policy that maximizes the
Government’s reward $r$, while also finding a balance between productivity,
equality, and carbon emissions under the designated economy-climate
coefficient $c^{e}$.

For the joint optimization of enterprise and government policies, we first
initialize the parameters of enterprise agents to those trained under
government policies based on Flat and Enterprise size as the indicator (SI)
scenario. Subsequently, the parameters of government policies are randomly
initialized. During training, we utilize the PPO algorithm [49] under the
RLlib framework [50]. Additionally, we experiment with various hyperparameters
for both enterprise and government agents, including learning rate and entropy
regularization. Following training with $400$ million samples, we find both
enterprise and government agents to converge to stable policies, which
effectively balance productivity, equity, and carbon emissions. (Figure 2)

###  4.2 Baseline Allocation Policies

We utilize the indicator approach [11] to allocate the proportion of total
carbon credits for each enterprise annually, with Emission (also called
grandfathering or GF [51]), Emission intensity (also called benchmarking or BM
[52]), and Enterprise size selected as indicators (shorten as SI). Due to the
large temporal scale of our simulation spanning $10$ years, the government
agent needs to allocate the total carbon emission credits for each year over
the $10$-year period. We refer to the global carbon emission historical data
and future forecasts provided by the IPCC [53] to establish emission scenarios
on a large temporal scale: we call this scenario Convex. Additionally, we
provide scenarios where the annual emission targets decrease gradually over
time (Decreasing), and scenarios where the annual emission targets remain
constant over time (Flat).

###  4.3 Simulation Details

In this section, we delve into the specifics of our simulation setup and the
observed behaviors of enterprise agents under various government policies. We
analyze key activities such as trade, production, and the construction of
government-certified projects, assessing their impact on overall productivity,
equality, and carbon emissions. Through these simulations, we aim to
illustrate how different allocation policies influence enterprise behavior and
the resulting economic and environmental outcomes.

#### Trade Activities.

In each trade activity, a carbon emission credit circulates among the
enterprises, along with coins corresponding to the current price of the carbon
emission credit. The price of emission credit fluctuates significantly with
changes in enterprise carbon emission level, making it challenging to analyze.
Therefore, we analyze enterprise trading activities by averaging trading
volumes over multiple episodes. As shown in the graph in the middle-upper part
of Figure 3, under various government policies, the circulation of coins among
enterprises due to trading aligns closely with the equality of government
policies.

#### Production Activities.

Produce action is the only action that generates coins, and the coins
generated in a single production are solely related to the enterprise’s size.
Therefore, the total income of all enterprises in an episode is closely
related to the total number of produce actions, which is essentially the
number of properties.

#### Government Certified project.

Government-certified projects in our setup represent government-led public
green energy initiatives, aligning with real-world scenarios. In our
simulator, these projects can store carbon credits across periods. When
constructed, they initially reduce the carbon emissions of all enterprises
(defined in Equation 3). They also pay $1$ unit of carbon emission credit to
the constructor (i.e., enterprise) but charge high coins and labor costs. As
observed in the bottom-left part of Figure 3, all baseline policies exhibit
lower total government-certified project construction. In contrast, the MARL-
based allocation policy stands out with significantly higher government-
certified project construction. This difference may stem from the MARL-based
policy’s more diversified allocation strategies, resulting in higher social
welfare.

#### Carbon emission.

According to Equation (3), invest actions can reduce the carbon emissions of
the investing enterprise while also having a minor impact on the carbon
emissions of other enterprises. The overall carbon emissions are related to
the overall production activities and the carbon emissions of individual
enterprises. From the lower-middle and lower-right parts of Figure 3, in
conjunction with the previous analyses of production activities and
government-certified projects, we can observe different behavioral patterns
associated with the $4$ different policies. Under the flat policy, enterprises
engage in more total production and investment activities, resulting in higher
carbon emissions. However, due to the lack of differentiation in this policy,
some enterprises exceed their allocated carbon credits significantly. In the
decreasing policy, economic activities of enterprises are subdued, leading to
lower productivity but higher equality and fewer excess carbon emissions. For
the convex policy, there are fewer total production activities but more
investment activities, resulting in lower and fewer excess carbon emissions.
Under the MARL-based policy, enterprises construct more government-certified
projects that can enhance overall income. Due to differentiated allocation
strategies, even with lower investment activity, there are more total
production activities and higher carbon emissions but fewer excess carbon
emissions.

##  5 Visualization

In this section, we present the visual tools developed to illustrate the
simulation results and analyze the behaviors of the agents in our carbon
market model. Visualization plays a crucial role in understanding the complex
interactions and outcomes within the multi-agent reinforcement learning (MARL)
environment. Our visual tools provide insights into the dynamics of enterprise
and government agents, their economic activities, and the impact of different
policies.

Figure 4: Simulator dashboard, it presents detailed information encompassing
enterprises’ attributes, assets, and actions within a single time step across
various example episodes under different policies. Additionally, it provides
visual representations of the average carbon prices over different periods and
presents rewards for both enterprises and government.

###  5.1 Dashboard Overview

The simulator dashboard (Figure 4) offers a comprehensive view of the
simulation at any given time step within an episode. It includes detailed
information about the attributes, assets, and actions of enterprises, as well
as visual representations of the average carbon prices, trade volumes, and
rewards for both enterprises and the government. This allows for an intuitive
comparison between different baseline and MARL strategies.

The dashboard is divided into several key sections:

  * •

Skills: Displays the research ability and manufacturing volume of each
enterprise.

  * •

Grid World: Visualizes the positions and movements of enterprises within the
grid world.

  * •

Carbon Price: Tracks the fluctuations in carbon emission credit prices.

  * •

Rewards: Shows the cumulative rewards for both enterprises and the government.

  * •

Carbon Emission Credits: Illustrates the remaining carbon credits for each
enterprise.

  * •

Coins: Displays the accumulation of coins by each enterprise.

  * •

Carbon Emission Level: Monitors the changes in carbon emission levels over
time.

  * •

Labor: Tracks the labor usage of each enterprise.

By leveraging these visualizations, we gain valuable insights into the trade-
offs between productivity, equality, and carbon emissions, ultimately guiding
the development of more effective carbon market policies.

##  6 Closing Remarks

In this paper, enterprise and government agents participate in carbon market
simulations via MARL-based adaptive mechanism design. By fine-tuning the
government’s reward function, we can exploit the adaptability to strike a
balance between various economic and climate objectives. Unlike the commonly
used indicator approach, MARL-based agents can incorporate more comprehensive
information, enabling them to formulate more personalized and diversified
allocation strategies. We also illustrates the practicality of employing
hierarchical model-free MARL for carbon market simulation. It envisions the
potential of machine learning to contribute to global emission reduction
endeavors. However, the proposed simulator still needs to be improved, notably
the absence of empirical modeling for emissions reduction investments made by
enterprises. Consequently, future simulations can enhance their realism by
integrating more real-world data.

## Ethical Statement

In our development of the carbon market simulator, we adhere to principles of
transparency, integrity, and fairness, ensuring compliance with the highest
ethical standards while advancing understanding in environmental economics. We
prioritize privacy, equity, and social responsibility throughout our research
and development process. However, it’s important to acknowledge that our
simulator may not encompass all aspects of the real world. As such, we do not
endorse the use of learned policies for actual policy making.

## Acknowledgements

Baoxiang Wang is partially supported by National Natural Science Foundation of
China (62106213, 72394361).

## References

  * [1] Richard Schmalensee, Thomas M. Stoker, and Ruth A. Judson.  World carbon dioxide emissions: 1950-2050.  The Review of Economics and Statistics, 80(1):15–27, 1998.
  * [2] United Nations.  Ending Poverty.  https://www.un.org/en/global-issues/ending-poverty/, 2023.  Accessed on .
  * [3] United Nations.  Water.  https://www.un.org/en/global-issues/water/, 2023.  Accessed on .
  * [4] United Nations.  Climate Change.  https://www.un.org/en/global-issues/climate-change/, 2023.  Accessed on .
  * [5] UNFCCC.  The paris agreement, 2015.
  * [6] UNFCCC.  Kyoto protocol.  1997\.  Accessed on .
  * [7] Michael Wara.  Is the global carbon market working?  Nature, 445(7128):595–596, 2007.
  * [8] Kaile Zhou and Yiwen Li.  Carbon finance and carbon market in china: Progress and challenges.  Journal of Cleaner Production, 214:536–549, 2019.
  * [9] Lawrence H Goulder and Andrew R Schein.  Carbon taxes versus cap and trade: a critical review.  Climate Change Economics, 4(03):1350010, 2013.
  * [10] Richard Schmalensee and Robert N Stavins.  The design of environmental markets: What have we learned from experience with cap and trade?  Oxford Review of Economic Policy, 33(4):572–588, 2017.
  * [11] P. Zhou and M. Wang.  Carbon dioxide emissions allocation: A review.  Ecological Economics, 125:47–59, 2016.
  * [12] Cameron Hepburn.  Carbon trading: A review of the kyoto mechanisms.  Annual Review of Environment and Resources, 32(1):375–393, 2007.
  * [13] Michael Hübler, Sebastian Voigt, and Andreas Löschel.  Designing an emissions trading scheme for china—an up-to-date climate policy assessment.  Energy policy, 75:57–72, 2014.
  * [14] Ling Tang, Jiarui Shi, and Qin Bao.  Designing an emissions trading scheme for china with a dynamic computable general equilibrium model.  Energy Policy, 97:507–520, 2016.
  * [15] Huimin Bi, Hao Xiao, and Kejuan Sun.  The impact of carbon market and carbon tax on green growth pathway in china: a dynamic cge model approach.  Emerging Markets Finance and Trade, 55(6):1312–1325, 2019.
  * [16] Ling Tang, Jiaqian Wu, Lean Yu, and Qin Bao.  Carbon allowance auction design of china’s emissions trading scheme: A multi-agent-based approach.  Energy Policy, 102:30–40, 2017.
  * [17] Meirong Zhou, Ming Zhou, Yanchun Pan, Zhimin Chen, and Jun Zeng.  Multi-agent-based simulation for policy evaluation of carbon emissions.  In Theory, Methodology, Tools and Applications for Modeling and Simulation of Complex Systems: 16th Asia Simulation Conference and SCS Autumn Simulation Multi-Conference, AsiaSim/SCS AutumnSim 2016, Beijing, China, October 8-11, 2016, Proceedings, Part III 16, pages 265–272. Springer, 2016.
  * [18] João Bernardo Narciso de Sousa.  Carbon market multi-agent simulation metamodel: a testbed for coordination mechanisms in sociotechnical systems.  2021\.
  * [19] Stephan Zheng, Alexander Trott, Sunil Srinivasa, David C Parkes, and Richard Socher.  The ai economist: Taxation policy design via two-level deep multiagent reinforcement learning.  Science advances, 8(18):eabk2607, 2022.
  * [20] Alexander Trott, Sunil Srinivasa, Douwe van der Wal, Sebastien Haneuse, and Stephan Zheng.  Building a foundation for data-driven, interpretable, and robust policy design using the ai economist.  arXiv preprint arXiv:2108.02904, 2021.
  * [21] Paul Baer, Glenn Fieldman, Tom Athanasiou, and Sivan Kartha.  Greenhouse development rights: towards an equitable framework for global climate policy.  In The Politics of Climate Change, pages 192–212. Routledge, 2013.
  * [22] Christoph Böhringer and Andreas Lange.  On the design of optimal grandfathering schemes for emission allowances.  European Economic Review, 49(8):2041–2055, 2005.
  * [23] Wen-Jing Yi, Le-Le Zou, Jie Guo, Kai Wang, and Yi-Ming Wei.  How can china reach its co2 intensity reduction targets by 2020? a regional allocation based on equity and development.  Energy Policy, 39(5):2407–2415, 2011.
  * [24] Mikel González-Eguino.  The importance of the design of market-based instruments for co2 mitigation: An age analysis for spain.  Ecological Economics, 70(12):2292–2302, 2011.
  * [25] Song-min Yu, Ying Fan, Lei Zhu, and Wolfgang Eichhammer.  Modeling the emission trading scheme from an agent-based perspective: System dynamics emerging from firms’ coordination among abatement options.  European Journal of Operational Research, 286(3):1113–1128, 2020.
  * [26] T Huw Edwards and John P Hutton.  Allocation of carbon permits within a country: a general equilibrium analysis of the united kingdom.  Energy Economics, 23(4):371–386, 2001.
  * [27] Akram Esmaeili Avval, Farzad Dehghanian, and Mohammadali Pirayesh.  Auction design for the allocation of carbon emission allowances to supply chains via multi-agent-based model and q-learning.  Computational and Applied Mathematics, 41(4):170, 2022.
  * [28] Hideyuki Mizuta, Sei Kato, and Hideki Tai.  Agent-based simulation of co2 emissions trading and large-scale traffic flow.  In SICE, 2008.
  * [29] Mohammadhussein Rafieisakhaei, Babak Barazandeh, Amirbahador Moosavi, Masoud Fekri, and Kaveh Bastani.  Modeling dynamics of the carbon market: A system dynamics approach on the co2 emissions and its connections to the oil market.  In ICSDS, 2016.
  * [30] Vincent Conitzer and Tuomas Sandholm.  Complexity of mechanism design.  arXiv preprint cs/0205075, 2002.
  * [31] Paul Dütting, Felix Fischer, Pichayut Jirapinyo, John K Lai, Benjamin Lubin, and David C Parkes.  Payment rules through discriminant-based classifiers, 2015.
  * [32] Yufei Wang, Zheyuan Ryan Shi, Lantao Yu, Yi Wu, Rohit Singh, Lucas Joppa, and Fei Fang.  Deep reinforcement learning for green security games with real-time information.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 1401–1408, 2019.
  * [33] David Balduzzi, Marta Garnelo, Yoram Bachrach, Wojciech Czarnecki, Julien Perolat, Max Jaderberg, and Thore Graepel.  Open-ended learning in symmetric zero-sum games.  In International Conference on Machine Learning, pages 434–443. PMLR, 2019.
  * [34] Pingzhong Tang.  Reinforcement mechanism design.  In IJCAI, pages 5146–5150, 2017.
  * [35] Weiran Shen, Binghui Peng, Hanpeng Liu, Michael Zhang, Ruohan Qian, Yan Hong, Zhi Guo, Zongyao Ding, Pengjun Lu, and Pingzhong Tang.  Reinforcement mechanism design: With applications to dynamic pricing in sponsored search auctions.  In Proceedings of the AAAI conference on artificial intelligence, volume 34, pages 2236–2243, 2020.
  * [36] David Thompson, Neil Newman, and Kevin Leyton-Brown.  The positronic economist: A computational system for analyzing economic mechanisms.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 31, 2017.
  * [37] Stephan Zheng, Alexander Trott, Sunil Srinivasa, Nikhil Naik, Melvin Gruesbeck, David C Parkes, and Richard Socher.  The ai economist: Improving equality and productivity with ai-driven tax policies.  arXiv preprint arXiv:2004.13332, 2020.
  * [38] Caroline Claus and Craig Boutilier.  The dynamics of reinforcement learning in cooperative multiagent systems.  AAAI/IAAI, 1998(746-752):2, 1998.
  * [39] Guillaume J Laurent, Laëtitia Matignon, Le Fort-Piat, et al.  The world of independent learners is not markovian.  International Journal of Knowledge-based and Intelligent Engineering Systems, 15(1):55–64, 2011.
  * [40] Tabish Rashid, Mikayel Samvelyan, Christian Schroeder De Witt, Gregory Farquhar, Jakob Foerster, and Shimon Whiteson.  Monotonic value function factorisation for deep multi-agent reinforcement learning.  The Journal of Machine Learning Research, 21(1):7234–7284, 2020.
  * [41] Peter Sunehag, Guy Lever, Audrunas Gruslys, Wojciech Marian Czarnecki, Vinicius Zambaldi, Max Jaderberg, Marc Lanctot, Nicolas Sonnerat, Joel Z Leibo, Karl Tuyls, et al.  Value-decomposition networks for cooperative multi-agent learning.  arXiv preprint arXiv:1706.05296, 2017.
  * [42] Ryan Lowe, Yi I Wu, Aviv Tamar, Jean Harb, OpenAI Pieter Abbeel, and Igor Mordatch.  Multi-agent actor-critic for mixed cooperative-competitive environments.  Advances in neural information processing systems, 30, 2017.
  * [43] Chao Yu, Akash Velu, Eugene Vinitsky, Jiaxuan Gao, Yu Wang, Alexandre Bayen, and Yi Wu.  The surprising effectiveness of ppo in cooperative multi-agent games.  Advances in Neural Information Processing Systems, 35:24611–24624, 2022.
  * [44] Tianmin Shu and Yuandong Tian.  M3RL: Mind-aware multi-agent management reinforcement learning.  In ICLR, 2019.
  * [45] Jinming Ma and Feng Wu.  Feudal multi-agent deep reinforcement learning for traffic signal control.  In AAMAS, 2020.
  * [46] Richard Bellman.  A markovian decision process.  Journal of mathematics and mechanics, pages 679–684, 1957.
  * [47] Eric A Hansen, Daniel S Bernstein, and Shlomo Zilberstein.  Dynamic programming for partially observable stochastic games.  In AAAI, 2004.
  * [48] Gerard Debreu et al.  Representation of a preference ordering by a numerical function.  Decision processes, 3:159–165, 1954.
  * [49] Christian Schroeder de Witt, Tarun Gupta, Denys Makoviichuk, Viktor Makoviychuk, Philip HS Torr, Mingfei Sun, and Shimon Whiteson.  Is independent learning all you need in the starcraft multi-agent challenge?  arXiv preprint arXiv:2011.09533, 2020.
  * [50] Eric Liang, Richard Liaw, Robert Nishihara, Philipp Moritz, Roy Fox, Ken Goldberg, Joseph Gonzalez, Michael Jordan, and Ion Stoica.  Rllib: Abstractions for distributed reinforcement learning.  In International conference on machine learning, pages 3053–3062. PMLR, 2018.
  * [51] Lars Zetterberg, Markus Wråke, Thomas Sterner, Carolyn Fischer, and Dallas Burtraw.  Short-run allocation of emissions allowances and long-term goals for climate policy.  Ambio, 41:23–32, 2012.
  * [52] Heleen Groenenberg and Kornelis Blok.  Benchmark-based emission allocation in a cap-and-trade system.  Climate Policy, 2(1):105–109, 2002.
  * [53] Climate.  Climate Change 2023: Synthesis Report. Contribution of Working Groups I, II and III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change.  IPCC, Geneva, 2023.


11institutetext: Nuclear and Chemical Sciences Division, Lawrence Livermore
National Laboratory, Livermore, California 94551, USA

# Numerical Convergence of Electromagnetic Responses with the Finite-Amplitude
Method

Tong Li 11<EMAIL_ADDRESS>Nicolas Schunck 11<EMAIL_ADDRESS>
###### Abstract

The response of a nucleus to an electromagnetic probe is a key quantity to
simulate photabsorption or photodeexcitation processes. For large calculations
at the scale of the entire mass table, this response can be estimated by
linear response theory. Thanks to the introduction of the finite-amplitude
method (FAM), calculations are computationally efficient. In this paper, we
investigate in more details the convergence of FAM calculations of the
response function as a function of the parameters controlling the numerical
implementation of the theory. We show that the response is much less sensitive
to the details of the single-particle basis than, e.g., Hartree-Fock-
Bogoliubov calculations.

##  1 Introduction

Photonuclear processes play central roles as probes to nuclear structure and
in applications of nuclear physics to nuclear technology, isotope production
or astrophysics. Much work is thus devoted to building reliable and
comprehensive data libraries of properties such as $\gamma$-strength functions
goriely2019reference . Notwithstanding fully phenomenological approaches based
on explicit fit to experimental data, large-scale theoretical predictions of
such quantities typically rely on the quasiparticle random-phase approximation
(QRPA) and linear response theory goriely2002largescale ;
goriely2004microscopic ; daoutidis2012largescale ; martini2016largescale ;
goriely2018gognyhfb . Because of their unfavorable computational scaling, most
of these calculations have often been restricted to spherical symmetry and/or
even-even nuclei.

With the introduction of the finite-amplitude method (FAM), the computational
cost of the linear response in deformed nuclei became considerably lower and
is now somewhat comparable to the one of a deformed Hartree-Fock-Bogoliubov
(HFB) calculation nakatsukasa2007finite ; avogadro2011finite . This has
enabled studies of nuclear multipole responses in heavy, deformed nuclei
inakura2009selfconsistent ; stoitsov2011monopole ; kortelainen2015multipole ;
oishi2016finite  as well as more accurate estimates of nuclear moments of
inertia and the collective inertia tensor petrik2018thoulessvalatin ;
washiyama2021finiteamplitude . The FAM was originally introduced and tested
with non-relativistic Skyrme functionals and later also implemented with
relativistic functionals niksic2013implementation .

In most of these early applications, the uncertainties of the results induced
by numerical inputs such as the basis parameters, quasiparticle cutoff and
smearing width were not fully discussed, with the exception of Refs.
stoitsov2011monopole  and oishi2016finite  which presented results for two
different basis sizes. Such uncertainties should be quantified before
embarking on any large-scale calculations of the response function. The goal
of this paper is thus to perform a more systematical study of the convergence
of FAM calculations with Skyrme functionals as a function of various numerical
features of the implementation. We briefly recall the theoretical formalism in
Section 2 before presenting a selection of the results in Section 3.

##  2 Theoretical Framework

We use the general framework of the energy density functional (EDF) to compute
the electromagnetic response of nuclei. The theory is presented in great
details in textbooks and review articles, most notably valatin1961generalized
; mang1975selfconsistent ; bender2003selfconsistent ; ring2004nuclear ;
schunck2019energy . In the following, we only recall the most essential
formulas.

###  2.1 Hartree-Fock-Bogoliubov Theory with Skyrme Functionals

The ground-state many-body wavefunction $\ket{\Phi}$ of the nucleus is
computed at the Hartree-Fock-Bogoliubov (HFB) approximation. We recall that
$\ket{\Phi}$ is then a product state of quasiparticle annihilation operators
that are related to any single-particle basis by a (unitary) Bogoliubov, or
canonical, transformation $\mathcal{W}$ characterized by the matrices $U$ and
$V$,

$\displaystyle\beta_{\mu}$ $\displaystyle=\sum_{m}\quantity[U^{\dagger}_{\mu
m}\,c_{m}+V^{\dagger}_{\mu m}\,c_{m}^{\dagger}]\,,\vskip 6.0pt plus 2.0pt
minus 2.0pt$ (1) $\displaystyle\beta_{\mu}^{\dagger}$
$\displaystyle=\sum_{m}\quantity[V^{T}_{\mu m}\,c_{m}+U^{T}_{\mu
m}\,c_{m}^{\dagger}]\,,$ (2)

where $\quantity{\beta,\beta^{\dagger}}$ are the quasiparticle operators and
$\quantity{c,c^{\dagger}}$ the single-particle operators. The one-body density
matrix $\rho=V^{*}V^{T}$ and pairing tensor $\kappa=V^{*}U^{T}$ can both be
expressed as functions of the Bogoliubov transformation, and the total energy
is a functional of $\rho$, $\kappa$ and $\kappa^{*}$. In practice, the
matrices $U$ and $V$ are determined by requiring that the total energy
$E[\rho,\kappa,\kappa^{*}]$ be a minimum with respect to all independent
matrix elements $\rho_{ij}$, $\kappa_{ij}$, and $\kappa^{*}_{ij}$. The
solution to the resulting HFB equation fully determines the Bogoliubov
matrices $U$ and $V$.

In this paper, we choose the functional form $E[\rho,\kappa,\kappa^{*}]$ to be
a Skyrme functional in the particle-hole channel (terms proportional to
$\rho$) and a simple zero-range, density-dependent functional in the particle-
particle channel (terms proportional to $\kappa$ and $\kappa^{*}$). As is well
known, the zero range of the pairing functional requires applying a cutoff on
the number of quasiparticles used to compute the density matrix and pairing
tensors.

###  2.2 Linear Response Theory with the Finite-Amplitude Method

Linear response theory quantifies the effect of applying a small one-body
perturbation operator $\hat{F}(t)$ to the ground state described by the many-
body wavefunction $\ket{\Phi}$ and is most naturally derived as the small-
amplitude limit of time-dependent Hartree-Fock-Bogoliubov theory (TDHFB)

$i\hbar\partialderivative{\hat{\mathcal{R}}}{t}=\quantity[\hat{\mathcal{H}}(t)+\hat{F}(t),\hat{\mathcal{R}}(t)]\,,$
(3)

where $\hat{\mathcal{R}}(t)$ is the time-dependent generalized density and
$\hat{\mathcal{H}}(t)$ the time-dependent HFB matrix. Here, one will assume
that

$\hat{F}(t)=\hat{F}(\omega)e^{-i\omega t}+\hat{F}^{\dagger}(\omega)e^{+i\omega
t}\,,$ (4)

where the frequency $\omega$ plays the role of the excitation energy in linear
response theory. Under the effect of the perturbation operator, the
generalized density $\hat{\mathcal{R}}_{0}$ associated with the HFB ground
state $\ket{\Phi}$ becomes, to the first order,
$\hat{\mathcal{R}}(t)=\hat{\mathcal{R}}_{0}+\delta\hat{\mathcal{R}}(\omega)e^{-i\omega
t}+\delta\hat{\mathcal{R}}(\omega)^{\dagger}e^{i\omega t}$, and
correspondingly
$\hat{\mathcal{H}}(t)=\hat{\mathcal{H}}_{0}+\delta\hat{\mathcal{H}}(\omega)e^{-i\omega
t}+\delta\hat{\mathcal{H}}(\omega)^{\dagger}e^{i\omega t}$. In the Bogoliubov
(quasiparticle) basis, the perturbed density
$\delta\hat{\mathcal{R}}(\omega)$, perturbed mean field
$\delta\hat{\mathcal{H}}$, and perturbation operator $\hat{F}(\omega)$ take
the form

$\delta\mathcal{R}(\omega)=\matrixquantity(0&X\\\
Y^{T}&0)\,,\qquad\delta\mathcal{H}(\omega)=\matrixquantity(0&\delta H^{20}\\\
-\delta H^{02}&0)\,,\qquad\mathcal{F}(\omega)=\matrixquantity(0&F^{20}\\\
-F^{02}&0)\,.$ (5)

The standard, matrix equation of linear response reads

$\left[\matrixquantity(A&B\\\ B^{*}&A^{*})-\omega\matrixquantity(1&0\\\
0&-1)\right]\matrixquantity(X(\omega)\\\
Y(\omega))=\matrixquantity(F^{20}(\omega)\\\ F^{02}(\omega))\,,$ (6)

where the matrices $A$ and $B$ together form the QRPA matrix and are obtained
from the second derivatives of the energy with respect to the generalized
density $\mathcal{R}$. Solving (6) directly can become computationally
prohibitive in the case of heavy deformed nuclei where the number of relevant
quasiparticles can be of the order of $10^{3}$, leading $A$ and $B$ to be
complex-valued dense matrices of size of about $10^{6}$.

####  2.2.1 Finite-Amplitude Method

The finite-amplitude method (FAM) provides an alternative way to determine the
amplitudes $X$ and $Y$ characterizing the perturbed density and scales much
more favorably than the direct matrix approach nakatsukasa2007finite ;
avogadro2011finite . The FAM yields a system of coupled equations

$\displaystyle\quantity(E_{\mu}+E_{\nu}-\omega)X_{\mu\nu}(\omega)+\delta
H^{20}_{\mu\nu}(\omega)$ $\displaystyle=-F^{20}_{\mu\nu}(\omega)\,,$ (7)
$\displaystyle\quantity(E_{\mu}+E_{\nu}+\omega)Y_{\mu\nu}(\omega)+\delta
H^{02}_{\mu\nu}(\omega)$ $\displaystyle=-F^{02}_{\mu\nu}(\omega)\,,$ (8)

where $\delta H^{20}$ and $\delta H^{02}$ are functionals of the perturbed
densities

$\displaystyle\delta\rho$ $\displaystyle=UXV^{T}+V^{*}Y^{T}U^{\dagger}\,,$ (9)
$\displaystyle\delta\kappa_{+}$
$\displaystyle=UXU^{T}+V^{*}Y^{T}V^{\dagger}\,,$ (10)
$\displaystyle\delta\kappa_{-}$
$\displaystyle=V^{*}X^{\dagger}V^{\dagger}+UY^{*}U^{T}\,,$ (11)

given by the Bogoliubov transformation of $\delta\mathcal{R}(\omega)$ as

$\matrixquantity(\delta\rho&\delta\kappa_{+}\\\
\delta\kappa_{-}^{\dagger}&-\delta\rho^{*})=\matrixquantity(U&V^{*}\\\
V&U^{*})\matrixquantity(0&X\\\
Y^{T}&0)\matrixquantity(U^{\dagger}&V^{\dagger}\\\ V^{T}&U^{T})\,.$ (12)

Given the initial HFB solution characterized by the matrices $U$ and $V$ and
the matrix elements $F^{20}$ and $F^{02}$ of the perturbation operator
$\hat{F}$, equations (7)-(8) can be solved iteratively.

After the $X$ and $Y$ amplitudes have been determined, the actual excitation
strength can be readily obtained as

$\frac{dB(\omega)}{d\omega}=-\frac{1}{\pi}\text{Im}\,S(\omega)\,,$ (13)

with the response function $S(\omega)$ given by

$S(\omega)=\sum_{\mu<\nu}\quantity[F^{20*}_{\mu\nu}X_{\mu\nu}(\omega)+F^{02*}_{\mu\nu}Y_{\mu\nu}(\omega)]\,.$
(14)

In practice, calculations are performed at complex values of the frequency
$\omega$ as $\omega\rightarrow\omega+i\tfrac{\Gamma}{2}$, where $\Gamma$ is
the smearing width. This is not only mathematically motivated by the need to
avoid the poles of the response function on the real axis, but also
corresponds, physically, to folding the solutions to (6) with a Lorentzian
distribution. In the following we will choose $\Gamma=1.0$ MeV if it is not
explicitly specified.

####  2.2.2 Quasiparticle Cutoff

Additional complications emerge when performing practical calculations with
zero-range pairing forces because of the ultra-violet divergence of the
pairing tensor bulgac2002renormalization ; borycki2006pairing ;
schunck2019energy . At the HFB level, all quasiparticles $\mu$ such that
$E_{\mu}>E_{\rm cut}$ are discarded. Depending on the value of the actual
cutoff energy and of the parameters of the HO basis, this implies that the $U$
and $V$ matrices may not be square but of dimension $(N,N_{\rm qp})$ with $N$
the size of the basis and $N_{\rm qp}\leq N$ the actual number of active
quasiparticles. In such a case, the one-body density matrix reads

$\rho_{ij}=\sum_{\mu\leq N_{\rm qp}}V_{i\mu}V^{*}_{\mu j}\,.$ (15)

How to handle quasiparticle truncation at the FAM level is less clear. Writing
(9) explicitly, we have for example

$\delta\rho_{ij}=\sum_{\mu,\nu}U_{i\mu}X_{\mu\nu}V_{\nu
j}+\sum_{\mu,\nu}V^{*}_{i\mu}Y_{\mu\nu}U^{*}_{\nu j}\,.$ (16)

The question here is whether summations over both $\mu$ and $\nu$ should be
truncated up to $N_{\rm qp}$, only one of them, or none. Some insight can be
gained by examining the time-dependent evolution of the quasiparticle operator
$\beta_{\mu}(t)$ of TDHFB,

$\beta_{\mu}(t)=\quantity[\beta_{\mu}+\delta\beta_{\mu}(t)]e^{iE_{\mu}t}\,,\qquad\delta\beta_{\mu}(t)=\sum_{\nu}\quantity[X_{\nu\mu}(\omega)e^{-i\omega
t}+Y^{*}_{\nu\mu}(\omega)e^{i\omega t}]\beta^{\dagger}_{\nu}\,.$ (17)

Quasiparticle index $\mu$ should obey $\mu\leq N_{\rm qp}$ to be consistent
with the definition of the one-body density matrix (15). However, the
expression for $\delta\beta_{\mu}(t)$ is not informative about any additional
truncation on $\nu$. Therefore, one may be allowed to let the row indices of
the matrices $X$ and $Y$ span the whole quasiparticle space while their column
indices must be truncated. We realize that this is an a minima argument that
should call for further clarifications, perhaps by extending the analysis the
behavior of the HFB densities $\rho$ and $\kappa$ in homogeneous nuclear
matter at the limit of large momenta; see discussion in Chapter 4 of
schunck2019energy . In the rest of this proceeding, we will adopt the “1/2”
truncation (only column indices are truncated).

##  3 Benchmark Results

In this section, we present a series of calculations aimed at benchmarking the
convergence of FAM calculations of the response function as a function of
different numerical input parameters. All calculations were performed in the
240Pu nucleus with the SLy4 parameterization of the Skyrme force
chabanat1997skyrme . In the pairing channel, we use $V_{0}^{(n)}=-219.973$ MeV
and $V_{0}^{(p)}=-285.524$ MeV as in navarroperez2022controllinga . HFB
calculations were performed with the code HFBTHO marevic2022axiallydeformed
(axial symmetry) using a quadrature grid of $N_{\rm GH}=N_{\rm GL}=40$ Gauss-
Hermite and Gauss-Laguerre integration points and $N_{\rm Leg}=80$ Gauss-
Legendre integration points stoitsov2013axially .

In this work, we considered both electric ($EL$) and magnetic ($ML$)
operators. We looked at the convergence patterns of both isoscalar and
isovector operators and, for each operator, at the $K=0$ and $K>0$ components.
Results were very similar in all cases, so we will only show results for the
isoscalar, electric quadrupole ($E2$) transition with $K=0$ defined by the
operator

$\hat{F}=e\sum_{i=1}^{A}r_{i}^{2}Y_{20}(\theta_{i},\varphi_{i})\,,$ (18)

where the $Y_{lm}(\theta,\varphi)$ are spherical harmonics
varshalovich1988quantum  and we took $e$ as the regular electric charge.

###  3.1 Convergence of the Response Function with Respect to the HO Basis

Both the HFB and FAM equations are solved in the harmonic oscillator (HO)
basis. Basis functions are characterized by an oscillator length $b_{0}$,
which is equivalent to an oscillator frequency
$\omega_{0}=\hbar/(mb_{0}^{2})$. In realistic calculations of deformed nuclei,
the basis itself is deformed, or stretched. For axially-deformed basis, this
means that the two frequencies $\omega_{z}$ and $\omega_{\perp}$ are
different, leading to different oscillator lengths $b_{z}\neq b_{\perp}$. In
HFBTHO, the ratios of the two frequencies is related to the basis deformation
parameter $\beta_{2}$; see stoitsov2013axially  for details. The oscillator is
thus characterized by two parameters, $b_{0}$ and $\beta_{2}$. The total
number of basis functions is determined by the number of shells
$N_{\mathrm{sh}}$.

Figure 1: Left: Isoscalar electric quadrupole excitation strength
$dB(\omega)/d\omega$ in 240Pu as a function of the excitation energy $\omega$,
for different values of the total number of shells in the basis. Right:
Variation of the HFB energy across the same range of oscillator shells.

In the left panel of Fig.1, we show how the isoscalar, electric quadrupole
excitation strength changes with the size of the basis as determined by the
total number of shells $N_{\mathrm{sh}}$. For each calculation, the oscillator
length is $b_{0}=2.3$ fm and the basis is spherical, $\beta_{2}=0$. It is
somewhat surprising to see that the first peak at $\omega\approx 5$ MeV is
hardly changed, the height of the second peak at $\omega\approx 10$ MeV varies
by less than 15%, and most of the variability of the results is in the small
third peak at around $\omega\approx 12$ MeV. This rather weak dependency of
the results with the number of shells should be contrasted with the rather
strong dependency of the HFB energy; see the right panel of Fig.1.

Figure 2: Left: Isoscalar electric quadrupole excitation strength
$dB(\omega)/d\omega$ in 240Pu as a function of the excitation energy $\omega$,
for different values of the oscillator length (in fm). Right: Variation of the
HFB energy across the same range of oscillator lengths.

In the left panel of Fig.2, we show a similar study for the variations of the
response with the oscillator length $b_{0}$. Here, in each calculation the
number of shells is $N_{\mathrm{sh}}=20$ and the basis is also spherical.
Although there are some small differences in the height of each of the peaks,
the response function is pretty insensitive to the value of $b_{0}$. The right
panel of the figure shows that the variations of the HFB energy with $b_{0}$
do not exceed 700 keV at $N_{\mathrm{sh}}=20$. Similar conclusions apply when
varying the deformation $\beta_{2}$ of the basis.

###  3.2 Impact of Quasiparticle Cutoff

As mentioned earlier, HFB and QRPA calculations with zero-range pairing forces
require introducing a cutoff on the number of quasiparticles used to define
the densities (15) and (16). In Fig.3, we illustrate the dependence of the
response function on the choice of the cutoff.

Given a cutoff value $E_{\rm cut}$, a consistent calculation requires first
calibrating the strengths $V_{0}^{(n)}$ and $V_{0}^{(p)}$ of the pairing force
on some experimental data. Here we perform this calibration at the HFB level
by employing the same prescription as in navarroperez2022controllinga . In
practice, for each $E_{\rm cut}$ we refitted the pairing strengths so that the
HFB pairing gaps match the three-point odd-even mass differences
$\Delta^{(3)}$ in 232Th. After this calibration, the pairing channel of the
FAM calculations is well defined and we can compute the perturbed densities
(16) according to the prescription discussed in Section 2.2.2.

Figure 3: Left: Isoscalar electric quadrupole excitation strength
$dB(\omega)/d\omega$ in 240Pu as a function of the excitation energy $\omega$,
for different quasiparticle cutoff energies. For each $E_{\rm cut}$, the
strengths $V_{0}^{(n)}$ and $V_{0}^{(p)}$ of the pairing force have been
refitted to keep the pairing gaps of 232Th at their reference values. Right:
Same as in the left panel, only a single HFB solution with $E_{\rm cut}=60$
MeV is used for all four FAM calculations; see text for additional details.

In the left panel of Fig. 3, we show results for the same isoscalar electric
quadrupole excitation strength as before for four values of the cutoff. Not
surprisingly, calculations are undistinguishable from one another since
recalibrating the pairing strengths imply that the underlying HFB solutions
are nearly identical by construction. For completeness, we also show in the
right panel the non-consistent case where the cutoff of the HFB solution is
fixed, here at $E_{\rm cut}=60$ MeV, while the cutoff of the FAM calculation
(16) itself varies. In other words a single identical HFB calculation is used
as the basis for all FAM calculations. Even in this not very realistic
scenario, the response depends only weakly on the cutoff.

###  3.3 Dependence of the Response Function on the Smearing Width

Finally, we discuss the dependency of the results on the smearing width
$\Gamma$ of the FAM response. As mentioned earlier, this parameter is
introduced because the response function (13) has poles corresponding to QRPA
eigenmodes: computing the response at complex $\omega$ values makes
calculations more efficient nakatsukasa2007finite . At the same time, it is
equivalent to computing the response (14) by summing explicitly over the
solutions of (6) and folding the result with a Lorentzian distribution of
width $\Gamma$. The net result is that $\Gamma$ becomes a parameter of the
calculation and its choice will have an impact on the results.

Figure 4: Isoscalar electric quadrupole excitation strength
$dB(\omega)/d\omega$ in 240Pu as a function of the excitation energy $\omega$,
for different smearing widths $\Gamma$.

From a computational perspective, large values of $\Gamma$ make calculations
faster and more stable. From a physics standpoint, semi-classical models of
the giant dipole response are typically based on a fit of the photoabsorption
cross section with a Lorentzian shape bohr1998nucleara . Surveys of
experimental results suggest a rather large variation of $\Gamma$ with the
mass number $A$ of the nucleus goriely2019reference . In Fig. 4 we show how
the isoscalar electric quadrupole excitation strength changes with $\Gamma$.
As expected, larger values of $\Gamma$ smooth out all the features of the
response while smaller values make its fine structure more visible. Out of all
the input parameters we considered in this work, $\Gamma$ is by far the one
with the largest impact on the results.

##  4 Conclusions

We have investigated the convergence patterns of the electromagnetic response
function as a function of the numerical parameters controlling the
calculation. Overall, the dependency of the results HO basis parameters is
weaker than that of the underlying HFB calculations themselves. The response
function is also insensitive to the quasiparticle cutoff – when the latter is
implemented consistently at both HFB and FAM levels. Only the smearing width
$\Gamma$ used to define the complex frequency $\omega$ of the response has a
visible impact on the results.

## Acknowledgments

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.
Computing support came from the Lawrence Livermore National Laboratory (LLNL)
Institutional Computing Grand Challenge program.

## References

  * (1) S. Goriely, P. Dimitriou, M. Wiedeking, T. Belgya, R. Firestone, J. Kopecky, M. Krtička, V. Plujko, R. Schwengner, S. Siem et al., Eur. Phys. J. A 55, 172 (2019)
  * (2) S. Goriely, E. Khan, Nucl. Phys. A 706, 217 (2002)
  * (3) S. Goriely, E. Khan, M. Samyn, Nucl. Phys. A 739, 331 (2004)
  * (4) I. Daoutidis, S. Goriely, Phys. Rev. C 86, 034328 (2012)
  * (5) M. Martini, S. Péru, S. Hilaire, S. Goriely, F. Lechaftois, Phys. Rev. C 94, 014304 (2016)
  * (6) S. Goriely, S. Hilaire, S. Péru, K. Sieja, Phys. Rev. C 98, 014327 (2018)
  * (7) T. Nakatsukasa, T. Inakura, K. Yabana, Phys. Rev. C 76, 024318 (2007)
  * (8) P. Avogadro, T. Nakatsukasa, Phys. Rev. C 84, 014314 (2011)
  * (9) T. Inakura, T. Nakatsukasa, K. Yabana, Phys. Rev. C 80, 044301 (2009)
  * (10) M. Stoitsov, M. Kortelainen, T. Nakatsukasa, C. Losa, W. Nazarewicz, Phys. Rev. C 84, 041305 (2011)
  * (11) M. Kortelainen, N. Hinohara, W. Nazarewicz, Phys. Rev. C 92, 051302(R) (2015)
  * (12) T. Oishi, M. Kortelainen, N. Hinohara, Phys. Rev. C 93, 034329 (2016)
  * (13) K. Petrík, M. Kortelainen, Phys. Rev. C 97, 034321 (2018)
  * (14) K. Washiyama, N. Hinohara, T. Nakatsukasa, Phys. Rev. C 103, 014306 (2021)
  * (15) T. Nikšić, N. Kralj, T. Tutiš, D. Vretenar, P. Ring, Phys. Rev. C 88, 044327 (2013)
  * (16) J.G. Valatin, Phys. Rev. 122, 1012 (1961)
  * (17) H.J. Mang, Phys. Rep. 18, 325 (1975)
  * (18) M. Bender, P.H. Heenen, P.G. Reinhard, Rev. Mod. Phys. 75, 121 (2003)
  * (19) P. Ring, P. Schuck, _The Nuclear Many-Body Problem_ , Texts and Monographs in Physics (Springer, 2004)
  * (20) N. Schunck, _Energy Density Functional Methods for Atomic Nuclei._ , IOP Expanding Physics (IOP Publishing, Bristol, UK, 2019), ISBN 978-0-7503-1423-7, https://doi.org/10.1088/2053-2563/aae0ed
  * (21) A. Bulgac, Y. Yu, Phys. Rev. Lett. 88, 042504 (2002)
  * (22) P. Borycki, J. Dobaczewski, W. Nazarewicz, M. Stoitsov, Phys. Rev. C 73, 044319 (2006)
  * (23) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, R. Schaeffer, Nucl. Phys. A 627, 710 (1997)
  * (24) R. Navarro Pérez, N. Schunck, Phys. Lett. B 833, 137336 (2022)
  * (25) P. Marević, N. Schunck, E.M. Ney, R. Navarro Pérez, M. Verriere, J. O’Neal, Comput. Phys. Commun. 276, 108367 (2022)
  * (26) M. Stoitsov, N. Schunck, M. Kortelainen, N. Michel, H. Nam, E. Olsen, J. Sarich, S. Wild, Comput. Phys. Commun. 184, 1592 (2013)
  * (27) D. Varshalovich, A. Moskalev, V. Khersonskii, _Quantum Theory of Angular Momentum_ (World Scientific, Singapore, 1988)
  * (28) A. Bohr, B. Mottelson, _Nuclear Structure_ , Vol. II, Nuclear Deformations (World Scientific, 1998)