# Power Availability of PV plus Thermal Batteries in real-world electric power
grids

Odin Foldvik Eikeland1,2, Colin C. Kelsall1, Kyle Buznitsky1, Shomik Verma1,
Filippo Maria Bianchi3,
Matteo Chiesa1,2,4, & Asegun Henry1

1Department of Mechanical Engineering, Massachusetts Institute of Technology,
Cambridge, MA, USA.
2Department of Physics and Technology, UiT – the Arctic University of Norway,
9037 Tromsø, Norway.
3Department of Mathematics and Statistics, UiT – the Arctic University of
Norway Technology, 9037 Tromsø, Norway.
4Laboratory for Energy and NanoScience (LENS), Khalifa University of Science
and Technology, Masdar.
Corresponding author<EMAIL_ADDRESS>author<EMAIL_ADDRESS>
###### Abstract

As variable renewable energy sources comprise a growing share of total
electricity generation, energy storage technologies are becoming increasingly
critical for balancing energy generation and demand.

In this study, we modeled an existing thermal energy storage unit with
estimated capital costs that are sufficiently low to enable large-scale
deployment in the electric power system. Our analysis emphasizes the value of
using such units to cost-effectively improve renewable energy dispatchability.
This study modeled an existing real-world grid rather than simulating
hypothetical future electric power systems. The storage unit coupled with a
photovoltaic (PV) system was modeled with different storage capacities,
whereas each storage unit size had various discharge capacities.

The modeling was performed under a baseline case with no emission constraints
and under hypothetical scenarios in which CO2 emissions were reduced. The
results show that power availability increases with increasing storage size
and vastly increases in the hypothetical CO2 reduction scenarios, as the
storage unit is utilized differently. When CO2 emissions are reduced, the
power system must be less dependent on fossil fuel technologies that currently
serve the grid, and thus rely more on the power that is served from the PV +
storage unit.

The proposed approach can provide increased knowledge to power system planners
regarding how adding PV + storage systems to existing grids can contribute to
the efficient stepwise decarbonization of electric power systems.

_K_ eywords Thermal Energy Grid Storage $\cdot$ Power System Modeling $\cdot$
Decarbonizing Electric Grids

##  1 Introduction

The use of variable renewable energy (VRE) resources, such as wind power and
solar photovoltaics (PV), is expanding rapidly as a share of total power
generation and is critical to the decarbonization of electrical power systems
[1, 2, 3]. The weather-dependent intermittency of VRE sources complicates the
planning and management of power systems as the electric power generation can
no longer be directly modulated to match the electricity demand. Energy
storage will therefore be an increasingly critical component of future energy
systems with high penetrations of VRE sources. Energy storage can charge
excess electricity in periods with high generation and low demand, and then
discharge the electricity in periods with low generation from the VRE sources
to match the load in periods with high electricity demand [4, 5].

The need for inexpensive storage over periods with different lengths, from
seconds to days and even seasonal storage, has accelerated in accordance with
the increasing share of VRE technologies in electric power systems [6, 7].
Pumped hydropower storage (PHS) is an established technique for large-scale
energy storage but can be used only in certain geographical areas. Lithium-ion
batteries have been the state-of-the-art technology for short-term storage.
However, capital costs between US$80 and US$100 kWh-1 make them unaffordable
for the multi-day storage objectives required to completely decarbonize the
grid [4, 5, 8]. Concentrated solar power with thermal energy storage (CSP-TES)
has been seen as a promising option, but major projects around the world have
been plagued by delays, cost overruns and mechanical issues, and interest has
waned in recent years [9, 10]. Studies suggest that achieving cost-efficient
multi-day storage requires a capital cost reduction to US$3-30 kWh-1 [5, 11].
Resolving this issue could enable more rapid decarbonization of the power
system, resulting in a 25% reduction in global GHG emissions [12, 13].
Therefore, one of the most significant issues that needs to be resolved to
achieve the GHG emission reduction targets is the energy storage challenge.

A storage concept based on Thermal Energy Storage (TES) has shown promising
potential to achieve sufficiently low capital cost in the multi-day storage
regime. TES stores the electricity as heat rather than electrochemically, and
then converts it back to electricity when needed [14]. The Thermal Energy Grid
Storage (TEGS) concept, detailed in [15], stores electricity as sensible heat
in graphite storage blocks and uses thermophotovoltaics (TPV) to convert heat
back to electricity on demand [15, 16]. While the conversion of heat to
electricity results in significant efficiency penalties, storing energy as
heat instead of electrochemically can be vastly cheaper, and thus the round-
trip efficiency (RTE) penalty compared to electrochemical batteries ($\sim$
90%) can potentially be a worthwhile tradeoff [15]. To maximize the conversion
efficiency from heat to electricity, the heat is stored at extremely high
temperatures ( 2400 °C). In a recent work by [17], the authors demonstrated a
world-record high conversion efficiency of 41% using TPV, and reported a
projected conversion efficiency of 50% in the future. As such this technology
can achieve a projected cost below US$ 20 kWh-1 at gigawatt scales.

In addition to the projected low cost, a unique property of TEGS compared to
Li-ion battery technology is the fact that, since energy is stored as heat in
graphite blocks and thereafter converted to electricity using TPV, it enables
the possibility of fully decoupling the charge and discharge capacities of the
storage unit. This allows the TEGS to charge (i.e., store heat) at a much
higher capacity than that required for discharging. The benefit of such a
mechanism is that a large amount of energy can be charged in a short amount of
time when generation surpluses exist and discharged over a longer period to
cover the electricity load in periods where demand exceeds supply.

The body of existing literature counts several studies that have employed
different approaches to evaluate the value of using storage to increase the
dispatchability of VRE sources, and the different studies have highlighted the
storage requirements (capital cost and storage duration) to enable the full
decarbonization of the power system. Tab. 1 shows an overview of some relevant
studies, where the key findings in each work are highlighted.

Despite the vast number of previous studies on modeling the value of energy
storage in emerging power systems, there is a lack of research that analyzes
specific real-world grids. Previous studies have largely focused on modeling
hypothetical future electric power systems starting from scratch (i.e.,
“greenfield” models) [7, 18, 19, 20, 21, 22, 23]. However, such studies can,
in many cases, lead to vague results as this requires a complete change of the
current electricity mix. This is in many cases challenging to implement due to
policy considerations. In addition, by modeling greenfield cases,
infrastructure that already exists is neglected.

Some studies have modeled existing grids (i.e., brownfield studies), but lack
in modeling the value of using specific storage technologies [5, 8, 24, 25].
These studies model the storage requirements in general, whereas all studies
show that the studies show that the capital cost must be below US $20 kWh-1.
No previous work has modeled the potential of using an emerging storage
concept based on TES that already exists on a lab scale [15, 16].

In this study, we propose a framework for addressing the value of using TEGS
that has sufficiently low capital cost to be economically used in an electric
grid. Using a Capacity Expansion Model (CEM) [26], we consider a hypothetical
PV + TEGS system that is interconnected to an existing real-world grid in the
Northeastern US. The TEGS unit charges excess electricity from PV during
periods of surplus generation. When the grid demands electricity and the PV
plant cannot deliver sufficient power due to a lack of solar availability, the
stored energy is discharged. Different storage sizes with varying discharge
capacities connected to the PV plant are modeled to optimize energy
availability. To investigate how emission constraints affect the energy
availability of PV + storage systems, a hypothetical future scenario is
modeled for the existing power system where CO2 emissions are reduced.

The main contributions of this study are: Rather than modeling the power
system as a greenfield case study, we analyze an abstract representation of an
existing grid, i.e., a “brownfield” model, and address how adding a PV +
storage system can contribute to decarbonizing the grid. Instead of modeling
general requirements of storage to enable the full decarbonization of the
power system, we here model a TES unit that currently exists at lab-scale and
has promising cost projections that are well-documented in the literature
[15]. This study can provide increased knowledge to power system planners on
how coupling emerging storage technologies and PV systems to existing grids
can contribute to stepwise decarbonizing of the grid in a more short-term
horizon.

The remainder of this paper is organized as follows. Section 2 presents the
share of technologies in the modeled existing electric power system, and the
methodology for modeling the grid in which the hypothetical PV + storage
system is added to interact with the existing power system. Section 3 presents
the metrics for computing the power availability of electric power systems. In
Section 4, the results are presented in terms of power availability under
different modeling scenarios. Finally, the conclusions are presented in
Section 5.

Table 1: Overview of relevant work addressing the value of energy storage Ref. | Year | Key findings
---|---|---
[19] | 2016 | Large-scale deployment of available battery technologies requires cost reductions.
[22] | 2016 | Pathways to fully renewable systems are feasible with high cost and overgeneration.
[4] | 2016 | Cost reduction for storage technologies is required to reach widespread profitability.
[20] | 2017 | The availability of how low-carbon technologies impact the optimal capacity mix and generation patterns were demonstrated.
[6] | 2018 | To reliably meet 100% of total annual electricity demand, weeks of energy storage are required to support with electricity.
[18] | 2018 |  | The role of energy storage units in power systems with high shares of VRE was analyzed.
---
The importance of storage increases with the increasing share of renewable-
based power technologies.
[21] | 2018 |  | Firm low-carbon resources consistently lower decarbonized electricity system costs,
---
and the availability of firm low-carbon resources reduces costs 10%–62% in
zero-CO2 cases.
[23] | 2019 |  | The benefits of hydropower and storage units were analyzed.
---
Three decarbonized power systems with distinct grid expansion strategies were
compared.
Cutting transmission volume does not increase the total costs.
[27] | 2019 | Curtailment of renewable energy generation can be avoided using energy storage.
[5] | 2019 | Energy storage cost below $20 kWh-1 can enable cost-competitive baseload power.
[28] | 2020 | Hydrogen storage with up to 2 weeks of discharge duration is expected to be cost-effective in future power systems.
[29] | 2020 | Hydrogen storage enable for sector coupling in real-world power systems.
[30] | 2020 | Electricity triangle assures a consistent framework for the energy transition.
[8] | 2020 | Current Li-ion capital cost exceeds storage value in many instances.
[24] | 2020 | Decarbonization is less expensive with Energy Storage Systems, given sufficient low-cost assumptions.
[7] | 2021 | Energy capacity costs must be $\leq$ US $ 20/kWh to reduce the electricity price by $\geq$ 10%.
[31] | 2021 | Power systems with 100% RE is possible using existing technologies.
[32] | 2021 | There is a need for analytic tool development to model how to achieve a power system that are 100% decarbonized.
[33] | 2021 | Clean firm resources are cost-effective in decarbonizing the grid.
[34] | 2022 | Green hydrogen cost between $0.79/kg and $1.94/kg in 2030 can be achieved.
[25] | 2022 | The demand for power capacity will drive future adoption of higher battery power capacity.

##  2 Methods

###  2.1 Modeled electric power system

In this study, an idealized single node representing the electric grid region
in the New England power system in North America is modeled. This system
considers one grid zone that represents a simplified power system topology of
the states of Massachusetts, New Hampshire, and Rhode Island. The electricity
demand, capital cost and performance data for the different generation
technologies in these regions were collected from the NREL annual technology
baseline (ATB) and EIA. The Github library PowerGenome 111The PowerGenome
Github library collects source data from EIA, NREL, and EPA and formats the
input files for the CEM model. The GitHub library could with associated
documentation could be found here https://github.com/PowerGenome/PowerGenome
was used to collect the input data and shape them to the required format for
the CEM. The weather year for modeling VRE availability was 2020. Fig. 1 shows
the share of different technologies in the current electric power system.
Natural Gas (NG) is the dominant power supply technology, accounting for 59%
of the installed capacity. In the existing power grid, VRE sources such as
wind and solar PV represent a smaller share (14%) of the overall electricity
generation mix.

Figure 1: The modeled New England grid zone is interconnected with a
hypothetical PV + storage system. The existing power system is dominated by
electricity generation from the NG. Solar and wind power represent a small
share of electric power systems. To charge the storage unit, electricity (from
the grid or PV) is converted to heat using resistance heaters. Heat is stored
in graphite blocks at high temperatures (illustrated as gray blocks). To
discharge the system, heat is transferred from the graphite blocks using
liquid tin to the TPV cells (illustrated as yellow/red pipes to the left of
the graphite blocs), where the transported heat is converted to electricity.

The hypothetical PV + storage power system is connected to the existing power
system through a transmission grid network. Fig. 1 shows how the hypothetical
system is connected to the existing grid, where the combined system can fully
participate in the power system by exchanging electricity on demand. In this
study, we chose to model a storage unit in conjunction with a PV plant, which
is believed to be the most dominant source of electricity generation in the
future power market [1, 2, 35]. In addition, the normal profile of the daily
generation from PV plants is believed to be a good match with storage
technologies, as it can store electricity when the PV plant power generates a
large amount of electricity during mid-day (and the demand is often low during
mid-day) and discharge the stored electricity when the sun is set (during
early morning and afternoon/evening). In this study, a PV plants with
installed peak power capacity of 100 MW and 1 GW were analyzed (see the
Supplementary material for the GW scale modeling), while the storage unit was
modeled with different energy storage capacities. This system is small
compared to the rest of the electric power system and is considered a price
taker in this grid.

###  2.2 Capacity Expansion Model (CEM) configuration and storage modeling

####  2.2.1 Capacity Expansion Model (CEM)

The analysis utilizes GenX [26], an electric power system CEM that evaluates
the cost-optimal electricity mix of the generation, storage, and transmission
infrastructure. The cost-optimization is subject to several constraints, such
as operational (electricity demand and generation) and policy (CO2 emission)
constraints. The dominant constraint used in this study was the maximum limit
on the allowed CO2 emissions in the grid. Here, the power system was modeled
with a baseline scenario without any CO2 constraints (i.e., the model finds
the cost-optimized electricity mix regardless of CO2 emissions) and with a
scenario where the CO2 emissions are reduced by 50%. Constraining the maximum
allowed CO2 emissions in the grid will change the dynamics of how the power
system is operated as the grid cannot use fossil-fuel-based technologies such
as NG at the same scale anymore. To fully capture high-resolution temporal
dependencies in the grid, we modeled the grid operation for each hour of the
year. All scenarios were evaluated with the Gurobi optimization solver [36]
using 16 cores with 128 GB RAM.

####  2.2.2 Storage modeling

Three different TEGS sizes coupled with the 100 MW PV plant were modeled. The
modeled energy storage capacities were as follows: 1) 400 MWh, 2) 600 MWh, and
3) 800 MWh. For each storage size, the charging capacity (i.e., the amount of
power that can be charged within one hour) is 100 MW, and the discharging
capacity is in the range of [5, 100] MW. These storage configurations were
modeled under the baseline case (with no CO2 reduction constraint) and 50% CO2
reduction scenarios. In total, 66 scenarios were modeled to address PV + TEGS
energy availability with different storage sizes, discharge capacities, and
CO2 constraints. Because future grid scenarios are modeled in this study, the
TEGS system is assumed to have a 50% RTE.

##  3 Experimental evaluations

###  3.1 Computing the availability of electric power systems

In this study, we are interested in assessing the amount of time our
hypothetical PV + TEGS system is available to the grid on demand. The Power
availability factor (PAF) was computed to describe power availability. The PAF
is computed as the fraction of time that the PV + storage system can deliver a
given power to the grid. Moreover, PAF allows the examination of the
availability of the combined power plant, as it measures how often a power
plant supplies the rated power to the grid. A power plant with a 100% PAF
constantly supplies the grid with a given rated power. In the case of solar
PV, the electricity generation suddenly drops when the sun no longer shines
(owing to cloud cover or when the sun sets). This sudden drop in electricity
generation reduces the PAF because the PV plant no longer generates at the
rated power. Here, storage units can be used to charge whenever there is a low
net demand for power in the grid and to discharge when there is a higher
demand for power. Fig. 2 illustrates an example of how the storage unit can be
used to shape the output to provide constant baseload power. Once the PV plant
starts generating electricity over derated power (e.g., 20 MW), excess
electricity is used to charge the storage unit. When the solar plant generates
less electricity than the derated power, the storage unit starts discharging
to satisfy the demand for electricity.

Figure 2: Example illustration where storage is used to provide a constant
baseload power to the grid by charging when there is excess electricity and
discharging when the PV plant does not generate electricity. The figure also
illustrates the large variability in the daily generation profile of solar PV
systems.

Fig. 2 illustrates how the combined PV and storage system can provide constant
baseload power, disregarding the electricity demand in the grid. However, when
the system is connected to an electricity grid, it becomes significantly more
complex. The system should not provide constant baseload power to the grid but
should be able to supply the requested power to the grid system operator. In
this study, the periods where the hypothetical PV + storage system cannot
provide the requested electricity to the grid on demand are detected. These
unwanted periods arise for the following reasons: 1) The PV system does not
deliver the required power, and 2) The TEGS system cannot discharge the
requested power as the State of Charge (SOC) is already zero. Such critical
periods can be calculated as:

$\textrm{Percentage Not
Available}=\frac{\textrm{When}((PV_{gen}+TEGS_{discharge}<\textrm{Derated
power})+(SOC=0))}{\textrm{8760 hours}},$ (1)

The system Percentage Not Available (PNA) gives information about the
percentage of time during a full year the system cannot provide the requested
electricity to the grid. On contrary, the PAF will be calculated as:

$PAF=1-PNA$ (2)

##  4 Results and discussions

###  4.1 Power availability from PV and energy storage

In Fig. 3, the yearly PAF from the hypothetical PV + TEGS system as a function
of different discharge powers is provided. The more the system is derated, the
more often it can deliver the required power. This is reasonable because the
more the system is derated, the more the storage system can discharge at a
lower-rated output for a longer period. If the discharge capacity of the
storage system has a discharge capacity of 5 MW, the PV+TEGS system can
deliver 5 MW to the grid approximately 95% of the hours during the year. If
the discharge capacity is 100 MW, the PV + TEGS system can deliver 100 MW to
the grid approximately 55-60% of the hours during the year. The percentage of
time the system cannot deliver the required power also changes with different
storage sizes, and there are fewer times during the year that the system
cannot deliver the required power if the storage size is larger.
Interestingly, there is a large difference between the scenarios with and
without CO2 reduction constraints. The electric power system dynamics change
completely once CO2 emissions were reduced by 50%. In this case, the increased
retirement of the NG makes the grid more dependent on the hypothetical PV +
TEGS system, which results in the system supplying the necessary power 100% of
the time for a derated power between 5 MW and 20 MW for the TEGS unit with a
storage capacity of 600 MWh and 800 MWh. This is remarkably higher than the
baseline case, where the system cannot deliver the required power 5-15% of the
time for such derated powers. When modeling the PAF for the PV + TEGS system
at the GW scale (see the Supplementary material), similar results are
obtained. The PAF increases when lowering the derated power in both the
baseline and the CO2 reduction scenarios. In addition, the PAF is
significantly higher when modeling the CO2 reduction scenario compared to the
baseline scenario. Although there are large differences in the PAF for the
scenarios, the modelled PV + storage system has a significantly higher PAF
than the average Capacity Factor (CF) (i.e., how often a power plant generates
electricity at rated power) for single PV plants, which is approximately 24%
[37]. Thus, the results show that combining storage with PV plants increases
the dispatchability of the VRE technology and provides more reliable power
output to the grid.

Figure 3: Example illustration where storage is used to provide a constant
baseload power to the grid by charging when there is excess electricity and
discharging when the PV plant does not generate electricity. The figure also
illustrates the large variability in the daily generation profile of solar PV
systems.

The large difference between the CO2 reduction scenarios indicates that the
electric power system dynamic changes significantly if the grid must be less
dependent on NG. Now the hypothetical PV + storage system plays a more
important role in the grid because it does not emit any CO2, and as such, the
cost-minimization schedule of the CEM optimizes the grid to ensure that the PV
+ storage system can deliver the requested power to the grid more often during
the year.

Fig. 4 illustrates an example of how the PAF changes with different CO2
reduction scenarios for a storage unit of 600 MWh that discharges 20 MW to the
grid. Clearly, reducing CO2 emissions results in a higher PAF. The maximum PAF
was achieved at 50% CO2 reduction. Reducing CO2 emissions requires the power
system to be less dependent on fossil fuel technologies, such as NG, and thus
must rely more on the power served by the PV + storage system. Modeling
scenarios with more than 50% CO2 reduction results in an infeasible solution
with the CEM optimization. Therefore, to further reduce the CO2 emissions in
the modeled grid, more renewable energy sources need to be installed in the
current grid to replace NG sources.

Figure 4: Percentage of time the hypothetical system can deliver the requested
20 MW power to the grid as a function of CO2 reduction scenarios. The PAF
increases as a function of CO2 reduction as the power system must rely more on
the hypothetical PV + storage system that does not emit any CO2.

Fig. 5 shows how the hypothetical system operates in the different emissions
scenarios (baseline and 50% CO2 reduction). Here, the hourly operation with a
600 MWh storage unit and 20 MW discharge capacity is shown during a typical
winter week in January. The uppermost and lowermost figure show how the system
operates in scenarios with and without constraining CO2 emissions. From both
graphs, TEGS is used frequently to discharge power to the grid whenever there
is low electricity generation from the solar PV, which reduces the
intermittency problem of PVs by increasing the number of hours the
hypothetical system can deliver the required power to the grid. In addition,
TEGS also charges power from the grid whenever there is a drop in the demand
to increase the SOC, which illustrates the benefit of using storage units that
are coupled with VRE resources. In the uppermost graph, a critical period is
highlighted by the yellow squared region. Clearly, the grid has an increasing
demand for electricity (blue line); however, the system cannot deliver the
requested power because there is no generation from the PV system, and the
storage system cannot discharge the required power to the grid because the SOC
is already zero. However, considering the 50% CO2 reduction scenario, the
system interacts differently with the grid, and it is clear that the CEM
optimizes the hypothetical system to have more energy available at more times
because the grid now is more dependent on the hypothetical system (because the
grid can use less NG). For the particular example week, the PV + storage
system can deliver the requested power to the grid at all times for the 50%
CO2 reduction scenario.

Figure 5: Winter week illustrations for a 100 MW PV plus a 600 MWh TEGS system
that can discharge 20 MW. Baseline case without constraining CO2 emissions
(uppermost graph) and future scenario with reduced emissions lowermost. The
baseline case shows a critical period during which the hypothetical system
cannot deliver the requested power to the grid.

Like the winter week example shown in Fig. 5, Fig. 6 shows how the
hypothetical system operates with the grid during a typical summer week. Here,
thanks to the higher solar availability, the system can deliver most of the
required power to the grid using only the PV plant, and the storage system is
used less frequently. However, for the baseline case, there are several
periods in which the grid requests energy and the storage unit cannot provide
sufficient power to the grid because the SOC is already zero. This is because
the grid is mainly dependent on NG, which can supply power whenever there is
low solar availability, and the TEGS system is only used to provide additional
peak power when the demand suddenly increases. Under the 50% CO2 reduction
scenario, the grid is much more dependent on the power from the storage unit
whenever there is no PV generation. It is clear that instead of providing peak
power to the grid, the PV system is used to charge the TEGS unit to a higher
degree to ensure that the SOC is never zero and thus can discharge the derated
power to the grid at all times when the PV system does not generate
electricity.

Figure 6: Illustrations of a typical summer week for a 100 MV PV plus 600 MWh
TEGS system that can discharge 20 MW. Baseline case without constraining the
CO2 emissions (uppermost graph), and future scenario with reduced emissions
lowermost. The PAF is significantly higher in the future scenarios.

##  5 Conclusion

In this study, we analyzed how storage can be used to tackle the intermittency
problem of VRE sources by increasing the dispatchability of a hypothetical PV
plant. We modeled an existing electricity grid region in North America using a
CEM and investigated how different storage configurations can reduce the
number of periods in which a hypothetical PV + storage system cannot provide
the required power to the grid.

Because of the high capital cost of electrochemical batteries, a TES
technology with a projected capital cost that fulfills the requirements (< US$
20 kWh-1) to enable full decarbonization of the grid was considered. The power
availability of the hypothetical system was modeled under different storage
sizes and discharge capacities. Additionally, the optimization schedule was
repeated under a hypothetical future scenario in which CO2 emissions were
constrained to be reduced by 50%. In total, 66 different scenarios were
modeled. To capture the high-resolution dependencies in the electricity
generation balance, a full year with hourly resolution was optimized using the
CEM.

The results support the added value of using storage to increase the
dispatchability of PV, as it significantly increases the PAF compared to PV
systems alone. The percentage of time during the year the system could not
deliver the required power to the grid decreased proportionally to the
discharge capacity of the storage unit. In addition, increasing the storage
size increases the energy availability, as more energy can be stored and
thereafter discharged over a longer period when there is a demand for
electricity in the grid. The findings were consistent when the PV plus TEGS
system were evaluated at both the megawatt and gigawatt scale.

Interestingly, there was a significant change in the electricity grid
generation dynamics when the CO2 emissions were reduced by 50%. Here, because
the grid can no longer rely on the same share of NG technology, the most cost-
efficient grid is achieved when the PV + TEGS system is utilized to a higher
degree, as these technologies do not emit any CO2. This shows that decreasing
the maximum allowed GHG emissions in the grid significantly increases the
value of using storage to increase the dispatchability of PV systems. We
believe that our findings will provide increased knowledge to power system
planners regarding how adding PV + storage systems to existing grids can
contribute to the efficient stepwise decarbonization of power systems.

###  5.1 Limitations and suggested future research

This study presented an idealized representation of an existing grid (i.e.,
“brownfield” CEM approach) in the New England grid region. However, we are
fully aware that the grid representation might not fully capture all details
of the existing grid, and there can be differences between our abstract grid
representation and the current real-world grid that is operated by ISO New
England. In addition to the transmission line between the existing grid and
hypothetical PV + storage system, we modeled the current grid as a single-zone
grid region without considering transmission losses or congestion between
generators and demand. This study’s major objectives are to evaluate a
hypothetical PV + storage system’s power availability and discuss the
significance of integrating such technologies to assure a successful step-by-
step decarbonization of the electric grid. Therefore, modeling the
transmission lines between the existing generators is not considered, as it is
outside the scope of this study and will significantly increase the
computational intensity of the CEM. The CEM is fully deterministic, assumes
perfect foresight in planning and operational decisions, and does not account
for uncertainty in VRE generation [38, 39, 40, 41]. Therefore, this study does
not aim to be used as a power planning tool for ISOs to assess the PAF of PV +
storage systems in the day-ahead electricity market, but shows how storage, in
general, will be a valuable technology to address the intermittency issue of
VRE. A suggested future study is to frame the CEM to account for the
uncertainty regarding the expected electricity generation from VRE sources.
This will allow the use of the CEM as a decision-making tool for optimizing
the management of the electricity grid in the day-ahead market.

##  6 Acknowledgement

O.F.E. and M.C. acknowledge the support from the research project
“Transformation to a Renewable & Smart Rural Power System Community (RENEW)”,
connected to the Arctic Centre for Sustainable Energy (ARC) at UiT-the Arctic
University of Norway through Grant No. 310026. We thank Maritsa Kissamitaki
for designing Figure 1 and Figure 2.

## References

  * [1] D. E. Remi Eriksen, Ulrike Haugen, Trond Hodne, Liv Hovem, Jin James Huang, Bent Erik Bakken, Kaveh Dianati, Thomas Horschig, Anne Louise Koefoed, Erica McConnell, Mats, S. S. Rinaldo, and Adrien Zambon.  Energy transition outlook 2022, a global and regional forecast to 2050, 2022.
  * [2] Dolf Gielen, Francisco Boshell, Deger Saygin, Morgan D Bazilian, Nicholas Wagner, and Ricardo Gorini.  The role of renewable energy in the global energy transformation.  Energy Strategy Reviews, 24:38–50, 2019.
  * [3] Harry Apostoleris, Sgouris Sgouridis, Marco Stefancich, and Matteo Chiesa.  Utility solar prices will continue to drop all over the world even without subsidies.  Nature Energy, 4(10):833–834, 2019.
  * [4] William A Braff, Joshua M Mueller, and Jessika E Trancik.  Value of storage technologies for wind and solar energy.  Nature Climate Change, 6(10):964–969, 2016.
  * [5] Micah S Ziegler, Joshua M Mueller, Gonçalo D Pereira, Juhyun Song, Marco Ferrara, Yet-Ming Chiang, and Jessika E Trancik.  Storage requirements and costs of shaping renewable energy toward grid decarbonization.  Joule, 3(9):2134–2153, 2019.
  * [6] Matthew R Shaner, Steven J Davis, Nathan S Lewis, and Ken Caldeira.  Geophysical constraints on the reliability of solar and wind power in the united states.  Energy & Environmental Science, 11(4):914–925, 2018.
  * [7] Nestor A Sepulveda, Jesse D Jenkins, Aurora Edington, Dharik S Mallapragada, and Richard K Lester.  The design space for long-duration energy storage in decarbonized power systems.  Nature Energy, 6(5):506–516, 2021.
  * [8] Dharik S Mallapragada, Nestor A Sepulveda, and Jesse D Jenkins.  Long-run system value of battery energy storage in future grids with increasing wind and solar generation.  Applied Energy, 275:115390, 2020.
  * [9] Xinhai Xu, K Vignarooban, Ben Xu, Keng Hsu, and Arunachala Mada Kannan.  Prospects and problems of concentrating solar power technologies for power generation in the desert regions.  Renewable and Sustainable Energy Reviews, 53:1106–1131, 2016.
  * [10] Kathleen M Kennedy, Tyler H Ruggles, Katherine Rinaldi, Jacqueline A Dowling, Lei Duan, Ken Caldeira, and Nathan S Lewis.  The role of concentrated solar power with thermal energy storage in least-cost highly reliable electricity systems fully powered by variable renewable energy.  Advances in Applied Energy, 6:100091, 2022.
  * [11] Paul Albertus, Joseph S Manser, and Scott Litzelman.  Long-duration electricity storage applications, economics, and technologies.  Joule, 4(1):21–32, 2020.
  * [12] R. K. & Meyer L. A. P. eds Core Writing Team.  Ipcc climate change 2014: Synthesis report, 2014.
  * [13] U. E. P. A. (EPA).  Inventory of us greenhouse gas emissions and sinks: 1990–2009,, 2011\.
  * [14] Asegun Henry, Ravi Prasher, and Arun Majumdar.  Five thermal energy grand challenges for decarbonization.  Nature Energy, 5(9):635–637, 2020.
  * [15] Caleb Amy, Hamid Reza Seyf, Myles A Steiner, Daniel J Friedman, and Asegun Henry.  Thermal energy grid storage using multi-junction photovoltaics.  Energy & Environmental Science, 12(1):334–343, 2019.
  * [16] Colin C Kelsall, Kyle Buznitsky, and Asegun Henry.  Technoeconomic analysis of thermal energy grid storage using graphite and tin.  arXiv preprint arXiv:2106.07624, 2021.
  * [17] Alina LaPotin, Kevin L Schulte, Myles A Steiner, Kyle Buznitsky, Colin C Kelsall, Daniel J Friedman, Eric J Tervo, Ryan M France, Michelle R Young, Andrew Rohskopf, et al.  Thermophotovoltaic efficiency of 40%.  Nature, 604(7905):287–291, 2022.
  * [18] Wolf-Peter Schill and Alexander Zerrahn.  Long-run power storage requirements for high shares of renewables: Results and sensitivities.  Renewable and Sustainable Energy Reviews, 83:156–171, 2018.
  * [19] Fernando J De Sisternes, Jesse D Jenkins, and Audun Botterud.  The value of energy storage in decarbonizing the electricity sector.  Applied Energy, 175:368–379, 2016.
  * [20] Clara F Heuberger, Iain Staffell, Nilay Shah, and Niall Mac Dowell.  A systems approach to quantifying the value of power generation and energy storage technologies in future electricity networks.  Computers & Chemical Engineering, 107:247–256, 2017.
  * [21] Nestor A Sepulveda, Jesse D Jenkins, Fernando J de Sisternes, and Richard K Lester.  The role of firm low-carbon electricity resources in deep decarbonization of power generation.  Joule, 2(11):2403–2420, 2018.
  * [22] Bethany A Frew, Sarah Becker, Michael J Dvorak, Gorm B Andresen, and Mark Z Jacobson.  Flexibility mechanisms and pathways to a highly renewable us electricity future.  Energy, 101:65–78, 2016.
  * [23] Hailiang Liu, Tom Brown, Gorm Bruun Andresen, David P Schlachtberger, and Martin Greiner.  The role of hydro power, storage and transmission in the decarbonization of the chinese power system.  Applied Energy, 239:1308–1321, 2019.
  * [24] Mehdi Jafari, Magnus Korpås, and Audun Botterud.  Power system decarbonization: Impacts of energy storage duration and interannual renewables variability.  Renewable Energy, 156:1171–1185, 2020.
  * [25] Anna H Schleifer, Caitlin A Murphy, Wesley J Cole, and Paul Denholm.  Exploring the design space of pv-plus-battery system configurations under evolving grid conditions.  Applied Energy, 308:118339, 2022.
  * [26] Jesse D Jenkins and Nestor A Sepulveda.  Enhanced decision support for a changing electricity landscape: the genx configurable electricity resource capacity expansion model.  2017\.
  * [27] Paul Denholm and Trieu Mai.  Timescales of energy storage needed for reducing renewable energy curtailment.  Renewable energy, 130:388–399, 2019.
  * [28] Omar J Guerra, Jiazi Zhang, Joshua Eichman, Paul Denholm, Jennifer Kurtz, and Bri-Mathias Hodge.  The value of seasonal energy storage technologies for the integration of wind and solar power.  Energy & Environmental Science, 13(7):1909–1922, 2020.
  * [29] Michael Schulthoff, Ivan Rudnick, Abhishek Bose, and Emre Gençer.  Role of hydrogen in a low-carbon electric power system: a case study.  Frontiers in Energy Research, 8:585461, 2021.
  * [30] Ettore Bompard, Audun Botterud, Stefano Corgnati, Tao Huang, Mehdi Jafari, Pierluigi Leone, Stefano Mauro, Giuseppe Montesano, Carlo Papa, and Francesco Profumo.  An electricity triangle for energy transition: Application to italy.  Applied Energy, 277:115525, 2020.
  * [31] Wesley J Cole, Danny Greer, Paul Denholm, A Will Frazier, Scott Machen, Trieu Mai, Nina Vincent, and Samuel F Baldwin.  Quantifying the challenge of reaching a 100% renewable energy power system for the united states.  Joule, 5(7):1732–1748, 2021.
  * [32] Paul Denholm, Douglas J Arent, Samuel F Baldwin, Daniel E Bilello, Gregory L Brinkman, Jaquelin M Cochran, Wesley J Cole, Bethany Frew, Vahan Gevorgian, Jenny Heeter, et al.  The challenges of achieving a 100% renewable electricity system in the united states.  Joule, 5(6):1331–1352, 2021.
  * [33] Ejeong Baik, Kiran P Chawla, Jesse D Jenkins, Clea Kolster, Neha S Patankar, Arne Olson, Sally M Benson, and Jane CS Long.  What is different about different net-zero carbon electricity systems?  Energy and Climate Change, 2:100046, 2021.
  * [34] Krishiv Gandhi, Harry Apostoleris, and Sgouris Sgouridis.  Catching the hydrogen train: economics-driven green hydrogen adoption potential in the united arab emirates.  International Journal of Hydrogen Energy, 2022.
  * [35] Felix Creutzig, Peter Agoston, Jan Christoph Goldschmidt, Gunnar Luderer, Gregory Nemet, and Robert C Pietzcker.  The underestimated potential of solar energy to mitigate climate change.  Nature Energy, 2(9):1–9, 2017.
  * [36] Miles Lubin and Iain Dunning.  Computing in operations research using julia.  INFORMS Journal on Computing, 27(2):238–248, 2015.
  * [37] NREL.  Electricity annual technology baseline (atb) data download., 2022.
  * [38] Yao Zhang, Jianxue Wang, and Xifan Wang.  Review on probabilistic forecasting of wind power generation.  Renewable and Sustainable Energy Reviews, 32:255–270, 2014.
  * [39] David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski.  Deepar: Probabilistic forecasting with autoregressive recurrent networks.  International Journal of Forecasting, 36(3):1181–1191, 2020.
  * [40] Aleksei Mashlakov, Toni Kuronen, Lasse Lensu, Arto Kaarna, and Samuli Honkapuro.  Assessing the performance of deep learning models for multivariate probabilistic energy forecasting.  Applied Energy, 285:116405, 2021.
  * [41] Odin Foldvik Eikeland, Finn Dag Hovem, Tom Eirik Olsen, Matteo Chiesa, and Filippo Maria Bianchi.  Probabilistic forecasts of wind power generation in regions with complex topography using deep learning methods: An arctic case.  Energy Conversion and Management: X, page 100239, 2022.


# Stability and strength of monolayer polymeric C60

Bo Peng<EMAIL_ADDRESS>Theory of Condensed Matter Group, Cavendish
Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE,
United Kingdom

###### Abstract

Two-dimensional fullerene networks have been synthesized in several forms [Hou
et al., Nature 2022, 606, 507], and it is unknown which monolayer form is
stable at ambient condition. Using first principles calculations, I show that
the believed stability of the quasi-tetragonal phases is challenged by
mechanical, dynamic or thermodynamic stability. For all temperatures, the
quasi-hexagonal phase is thermodynamically least stable. However, the
relatively high dynamic and mechanical stabilities suggest that the quasi-
hexagonal phase is intrinsically stronger than the other phases under various
strains. The origin of the high stability and strength of the quasi-hexagonal
phase can be attributed to the strong covalent C$-$C bonds that strongly hold
the linked C60 clusters together, enabling the closely packed hexagonal
network. These results rationalize the experimental observations that so far
only the quasi-hexagonal phase has been exfoliated experimentally as
monolayers.

Recent attempts to synthesize layers of connected buckyballs, i.e. C60
molecules linked by carbon$-$carbon bonds, have obtained different
arrangements of cluster cages through the formation of bonds between
neighboring C60 molecules Hou et al. (2022). The obtained allotropes include a
few-layer rectangular structure in which each C60 molecule has four
neighboring buckyballs and a monolayer hexagonal structure in which each C60
cage binds to six neighbors, namely, a few-layer quasi-tetragonal phase (qTP)
and a monolayer quasi-hexagonal phase (qHP) respectively. Great efforts have
been devoted to stabilizing the linking bonds between neighboring cluster
cages by introducing magnesium atoms to form a quasi-2D fullerene network with
strong intralayer covalent bonds Hou et al. (2022) because Mg atoms tend to
promote covalent bonds Pontiroli et al. (2013); Tanaka and Yamanaka (2018). To
aid exfoliation, the Mg ions that hold the C60 cages together can be then
replaced by large organic ions, which can be removed afterwards by hydrogen
peroxide, leading to pure, charge neutral fullerene networks in 2D Hou et al.
(2022); Gottfried (2022). Unfortunately, only qHP C60 has been obtained as
monolayers, while all the qTP C60 flakes are few-layer Hou et al. (2022).
These results raise doubts regarding the stability of monolayer fullerene
networks.

Ever since the discovery of C60 Kroto et al. (1985), the formation mechanism
and stability of the fullerene molecules are far from completely understood
Kroto (1987, 1990); Goroff (1996); Bernal et al. (2019). When forming
structural units of C60 clusters in a 2D plane, it is unclear whether ordered
structures of monolayer polymeric C60 are stable under ambient conditions such
as strain and temperature. Recent first principles calculations have
investigated various structural phases of monolayer C60 with different bonding
characters Tromer et al. (2022); Yu et al. (2022); Yuan et al. (2022); Peng
(2022); Ying et al. (2023); Dong et al. (2022). The mechanical stability of
several phases has been confirmed Yu et al. (2022); Tromer et al. (2022); Ying
et al. (2023). More recently, the thermal stability of monolayer C60 has been
addressed using molecular dynamics simulations, showing that both qTP and qHP
C60 monolayers can remain stable at temperatures near 800 K Ribeiro et al.
(2022), which is partially consistent with the experimental result that
monolayer qHP C60 does not decompose at 600 K Hou et al. (2022). However,
previous analysis based on mechanical and thermal stability cannot explain why
the qTP monolayers have not yet been exfoliated experimentally. Furthermore,
the dynamic stability of monolayer fullerene networks with respect to lattice
vibrations, which indicates whether the crystal structure is in a local
minimum of the potential energy surface Peng et al. (2017); Malyi et al.
(2019); Luo et al. (2022); Pallikara et al. (2022), is still unexplored.
Additionally, the thermodynamic stability of different phases, which
energetically classifies the stability (especially at finite temperatures)
Pavone et al. (1998); Pavone (2001); Masago et al. (2006); van Setten et al.
(2007); Stoffel et al. (2010); Zhang et al. (2011); Deringer et al. (2014,
2014); White et al. (2015); Nyman et al. (2016); Skelton et al. (2017);
Pallikara and Skelton (2021); Bartel (2022), remains unknown. Because of such
knowledge gaps, several questions need to be answered to understand the phase
stability of monolayer fullerene networks: ($i$) Are qTP and qHP C60, as pure
carbon monolayers without extra Mg or organic ions to bind the C60 cages
together, dynamically stable? ($ii$) What is their relative stability in a
thermodynamic aspect? ($iii$) Can the calculated mechanical strength support
their phase stability?

Figure 1: Bond structures for (a) qTP1, (b) qTP2 and (c) qHP C60. The default
isosurface level in vesta Momma and Izumi (2011) is used. Maps of the ELF on
the (010) plane are also present.

In this work, I investigate the mechanical, dynamic and thermodynamic
stability of monolayer qTP and qHP C60 by using first principles calculations.
Structural relaxation obtains two crystal structures of the quasi-tetragonal
phase (denoted as qTP1 and qTP2 respectively). I show that the qTP1 monolayer,
albeit being thermodynamically stable at all temperatures above 380 K,
possesses low dynamic and mechanical stability due to its weak bonding
perpendicular to the straight chains of C60 buckyballs. On the other hand,
although qTP2 fullerene might be the ground-state structure with the lowest
Gibbs free energy at 0 K and exhibits good dynamic and mechanical stability,
it is only thermodynamically stable with respect to qTP1 C60 at low
temperatures. Instead, monolayer qHP C60 should be experimentally accessible
owing to their dynamic and mechanical stability, in spite of its lowest
thermodynamic stability among all three phases. In addition, qHP C60 has the
highest strength under various strains (hydrostatic, uniaxial and shear)
because of the closely packed crystal structures.

First principles calculations are performed using the Vienna ab-initio
simulation package (vasp) Kresse and Furthmüller (1996, 1996). The projector
augmented wave (PAW) potential is used with C $2s^{2}2p^{2}$ valence states
Blöchl (1994); Kresse and Joubert (1999), under the generalized gradient
approximation (GGA) with the Perdew-Burke-Ernzerhof parameterization revised
for solids (PBEsol) as the exchange-correlation functional Perdew et al.
(2008). The crystal structures are optimized by fully relaxing the lattice
constants and internal atomic coordination (for computational details, see the
Supporting Information). Geometry optimization starting from the quasi-
tetragonal phase consisting of only carbon atoms leads to a quasi-1D qTP
monolayer (qTP1), as shown in Fig. 1(a). On the other hand, a two-step
structural relaxation, starting from monolayer qTP Mg2C60 and then removing
the Mg ions before the second relaxation, obtains a tightly bound qTP
monolayer (qTP2), as shown in Fig. 1(b). The two-step structure relaxation
mimics the experimental procedure to remove the charged ions introduced during
synthesis Hou et al. (2022); Gottfried (2022). The computed lattice constants
for all three phases are listed in Table 1, which are in good agreement with
previous results Yu et al. (2022); Tromer et al. (2022); Dong et al. (2022),
therefore confirming the reliability of the present calculations.

Table 1: Calculated static lattice constants (in Å) and cohesive energy $E_{\mathrm{c}}$ (in eV/atom) of qTP1, qTP2 and qHP C60 monolayers, 1D qTP C60 chain and 0D C60 molecule. The cohesive energy is defined as $E_{\mathrm{c}}=E_{\mathrm{tot}}/N-E_{\mathrm{isolated}}$, where $E_{\mathrm{tot}}$ is the total energy of the crystal, $N$ is the number of atoms in the unit cell and $E_{\mathrm{isolated}}$ is the total energy of an isolated carbon atom. The room-temperature lattice constants calculated under the quasi-harmonic approximation are also listed in parentheses for comparison.  phase | $a$ | $b$ | $E_{\mathrm{c}}$
---|---|---|---
2D qTP1 | 10.491 | 9.063 | –9.2582
| (10.522) | (9.090) |
2D qTP2 | 9.097 | 9.001 | –9.2587
| (9.132) | (9.031) |
2D qHP | 15.848 | 9.131 | –9.2465
| (15.896) | (9.162) |
1D | 9.062 | - | –9.2579
| (9.098) | - |
0D | - | - | –9.2564
Table 2: Elastic properties for qTP1, qTP2 and qHP C60, with the elastic constants $C_{ij}$, shear modulus $G^{\mathrm{2D}}$, layer modulus $\gamma$, Young’s modulus $Y^{\mathrm{2D}}$ in N/m, and Poisson’s ratio $\nu$ dimensionless. The elastic constants $C_{11}$, $C_{22}$ and $C_{66}$ calculated from phonon speed of sound are also listed in parentheses for comparison. phase | $C_{11}$ | $C_{22}$ | $C_{12}$ | $C_{66}=G^{\mathrm{2D}}$ | $\gamma$ | $Y^{\mathrm{2D}}_{a}$ | $Y^{\mathrm{2D}}_{b}$ | $\nu_{a}$ | $\nu_{b}$
---|---|---|---|---|---|---|---|---|---
qTP1 | 5.4 | 123.7 | –1.2 | –0.2 | 31.7 | 5.4 | 123.5 | –0.010 | –0.225
| (2.5) | (121.3) | - |  |  |  |  |  |
qTP2 | 149.9 | 148.7 | 22.9 | 53.4 | 86.1 | 146.4 | 145.2 | 0.154 | 0.153
| (150.5) | (141.2) | - | (54.5) |  |  |  |  |
qHP | 150.8 | 186.8 | 22.5 | 60.6 | 95.6 | 148.1 | 183.4 | 0.120 | 0.149
| (142.4) | (172.7) | - | (61.7) |  |  |  |  |

The bond structures at equilibrium are examined in Fig. 1. The relaxed
structure for qTP1 fullerene can be regarded as one-dimensional chains of C60
cages along the $b$ direction that are linked by the nearly in-plane [2+2]
cycloaddition bonds. In comparison, qTP2 fullerene is a two-dimensional
network of C60 cages connected by the out-of-plane vertical [2+2]
cycloaddition bonds along the $a$ direction and the in-plane [2+2]
cycloaddition bonds along the $b$ direction. The major difference between qTP1
and qTP2 is the absence of the vertical [2+2] cycloaddition bonds along $a$ in
the former. Regarding the qHP monolayer, the C60 cages form a hexagonal
network through the similar planar [2+2] cycloaddition bonds along the $b$
direction and C$-$C single bonds along the other two directions diagonal to
the rectangular unit cell.

Figure 1 also shows the maps of the electron localization function (ELF) on
the (010) plane. A high value of ELF indicates strong electron localization
Becke and Edgecombe (1990); Savin et al. (1992); Gatti (2005); Chen and Kamran
(2013). As shown in Fig. 1(a), the covalent [2+2] bonds along $b$ in qTP1
fullerene lead to high electron localization there (plane 1), whereas no bonds
are formed between neighboring C60 cages along $a$ (plane 2). In contrast, the
vertical [2+2] bonds along $a$ in qTP2 fullerene result in high electron
localization between neighboring C60 cages, as demonstrated in plane 2 of Fig.
1(b). For qHP C60, the hexagonal network has higher electron localization in
both directions, as one can see from Fig. 1(c). As a result, one can expect
that the hexagonal networks should stabilize and strengthen the structure of
qHP C60, making it slightly more stable than qTP2 C60 while much more stable
than qTP1 C60. However, as shown below, although the mechanical and dynamic
stabilities are consistent with the ELF picture, high electron localization in
qHP C60 does not guarantee their thermodynamic stability.

To confirm the mechanical stability, the elastic constants are calculated by
finite differences through finite distortions of the lattice Le Page and Saxe
(2002); Wu et al. (2005). There are different ways to define the 2D elastic
constants from the computed 3D coefficients Blonsky et al. (2015); Peng et al.
(2017); Pashartis et al. (2022). Here the 2D coefficients
$C_{ij}^{\mathrm{2D}}$ are renormalized by the $c$ lattice constant (the
spacing between 2D layers) Blonsky et al. (2015); Peng et al. (2017), i.e.
$C_{ij}^{\mathrm{2D}}=c\times C_{ij}^{\mathrm{3D}}$. The obtained 2D elastic
constants (including ionic relaxations) are listed in Table 2 using the Voigt
notation: $1-xx$, $2-yy$, $6-xy$, and the present results agree well with
previous calculations Yu et al. (2022); Tromer et al. (2022). According to
Born-Huang’s lattice dynamical theory Born and Huang (1954); Wu et al. (2007),
in monoclinic crystals (qTP1 and qHP with space group $P2/m$ and $Pc$
respectively), the mechanical stability criteria is given by

$C_{11}>0,\ C_{22}>0,\ C_{66}>0,\ C_{11}+C_{22}+2C_{12}>0.$ (1)

In orthorhombic crystals (qTP2 C60 with space group $Pmmm$), the Born
stability criteria has an extra requirement in addition to Eq. (1)

$C_{11}+C_{22}-2C_{12}>0.$ (2)

The elastic constants of qTP2 and qHP C60 satisfy their corresponding
criteria, indicating that they are mechanically stable. Interestingly, the
shear strength $G^{\mathrm{2D}}$ of qTP1 C60 is negative, demonstrating its
low shear resistance. The 1D chains in qTP1 C60 are prone to bending under
shear deformation, which may lead to a sliding of C60 chains and even lattice
instability. In addition, $C_{11}$ in qTP1 fullerene is more than one order of
magnitude lower than $C_{22}$. Such weak stiffness is correlated to weak
interchain bonding effect along $a$ as discussed above and the weak dynamic
stability as will be demonstrated below. In contrast, $C_{11}$ and $C_{22}$
are nearly the same in qTP2 fullerene because of the similar [2+2]
cycloaddition bonds along both $a$ and $b$. For qHP fullerene, the elastic
constants $C_{11}$, $C_{22}$ and $C_{66}$ are the highest among the three
phases, in consistent with the high electron localization in the hexagonal
networks that strengthens the crystal structure.

The strength of monolayer fullerene networks is obtained from the computed
elastic constants, as summarized in Table 2. The layer modulus $\gamma$ is the
2D equivalent to the bulk modulus, which measures the resistance to
hydrostatic stretching in 2D materials Andrew et al. (2012). The layer moduli
show an increasing trend from qTP1 to qTP2 and to qHP C60. The $\gamma$ for
qTP2 C60 is more than twice that of qTP1 C60, while slightly lower than that
of qHP C60, which concurs with the bonding structures of three fullerene
networks. In general, qTP2 and qHP C60 have comparable moduli and therefore
similar hardness properties, whereas qTP1 C60 has less resilience to both
shear and hydrostatic strains.

Figure 2: Low-frequency phonons of (a) qTP1, (b) qTP2 and (c) qHP C60 using
the static and room-temperature lattice constants. Entire phonon spectra for
(d) qTP1, (e) qTP2 and (f) qHP C60 using the room-temperature lattice
constants, with the phonon occupation number $N_{\mathrm{ph}}$ determined from
the Bose-Einstein distribution function at 300 K.

The anisotropy of strength is also investigated by calculating the Young’s
modulus and Poisson’s ratio. In qTP1 fullerene, the Young’s modulus
$Y^{\mathrm{2D}}$ along $a$ is more than 22 times lower than that along $b$,
indicating that qTP1 C60 is much less structurally rigid to elongations along
$a$. In qTP2 C60, the Young’s moduli along $a$ and $b$ are nearly the same due
to similar [2+2] cycloaddition bonds along both directions, showing that they
have the same resilience to linear strain. Regarding the qHP monolayers, the
$Y^{\mathrm{2D}}_{a}$ has a value 80% the $Y^{\mathrm{2D}}_{b}$, indicating
slightly weaker stiffness of the C$-$C single bonds in the presence of strain
along $a$. The Poisson’s ratio $\nu$ for qTP1 C60 is negative, i.e. the qTP1
fullerene monolayers expand laterally when stretched, and the $|\nu_{a}|$ is
significantly lower than $|\nu_{b}|$ because of much less bond stretching
under uniaxial strain. Monolayer qTP2 C60 has a nearly isotropic $\nu$ of
$0.153-0.154$, while the $\nu_{a}$ in qHP fullerene is slightly lower than
$\nu_{b}$. These results indicate that qTP1 C60 is unable to withstand greater
strains along $a$ than those along $b$, which is the origin of its overall low
strength.

To evaluate the dynamic stability of monolayer fullerene networks, lattice
dynamical properties are calculated within the harmonic approximation based on
density functional perturbation theory Baroni et al. (2001); Gonze (1995,
1995). The phonon spectra of all three phases are gathered in Fig. 2. As shown
in Fig. 2(a), the phonon dispersion of qTP1 C60 using the static lattice
constants exhibits small imaginary frequency (¡ 0.6i THz) along the entire
$\Gamma$-X high-symmetry line. An imaginary frequency indicates a decrease in
potential energy when the atoms are displaced away from their equilibrium
positions, corresponding to a non-restorative force Pallikara et al. (2022).
Therefore, the imaginary frequency along $\Gamma$-X implies that monolayer
qTP1 C60 can be split into individual 1D chains in the presence of interchain
(out-of-plane) vibrations, demonstrating its weak dynamic stability along the
$a$ direction. There is a fourth mode at $\Gamma$ with nearly zero frequency
in qTP1 C60, which, sometimes known as the torsional acoustic mode, is a
strong indication of the (quasi-)1D nature Peng et al. (2018); Lange et al.
(2022). The thermal expansion is included by computing the Gibbs free energy
under the quasi-harmonic approximation Huang et al. (2015, 2016); Peng et al.
(2019), and the room-temperature lattice constants are listed in Table 1. At
300 K, the imaginary mode in qTP1 C60 remains along $\Gamma$-X, though the
imaginary frequency becomes smaller (¡ 0.2i THz). In contrast, qTP2 and qHP
fullerene are dynamically more stable as there is no imaginary mode in Fig.
2(b) and (c) using both the static and room-temperature lattice constants,
indicating these structures are a local minimum on the potential energy
surface and the atoms vibrate harmonically around their equilibrium positions.

From the phonon speed of sound, the elastic constants $C_{11}$, $C_{22}$ and
$C_{66}$ can be calculated Feng et al. (2012) (for details on the phonon group
velocity, see the Supporting Information). As shown in Table 2, the calculated
elastic constants are in reasonably good agreement with those computed from
the finite difference method Le Page and Saxe (2002); Wu et al. (2005).
Moreover, qHP fullerene has the highest speed of sound along $b$ and the
highest phonon frequency with four phonon branches higher than 48.4 THz
throughout the entire Brillouin zone, whereas the highest phonon frequencies
in qTP1 and qTP2 fullerene are lower than 47.7 THz, which is in line with the
high mechanical strength in qHP C60.

To clarify the thermodynamic stability of C60 monolayers, the cohesive energy
is calculated, as listed in Table 1. The resulting $E_{\mathrm{c}}$ of qTP2
fullerene is 0.5 meV/atom (30 meV per formula unit) lower than qTP1 and 12.2
meV/atom (732 meV per formula unit) lower than qHP, suggesting its
thermodynamic stability. However, because the energy difference between the
qTP1 and qTP2 monolayers is quite small, phonons can play an important role in
determining the thermodynamic stability at both 0 K and finite temperatures
Pavone et al. (1998); Pavone (2001); Masago et al. (2006); van Setten et al.
(2007). The contribution of phonons can be examined by calculating the Gibbs
free energy $F$ Dove (1993); Togo et al. (2008); Togo and Tanaka (2015)

$F=\mathop{\min}_{a,b}\bigg{[}E_{\mathrm{tot}}+\frac{1}{2}\sum_{\lambda}\hbar\omega_{\lambda}+k_{B}T\sum_{\lambda}\ln{\big{(}1-\mathrm{e}^{-\hbar\omega_{\lambda}/k_{B}T}\big{)}}\bigg{]},$
(3)

where $\mathop{\min}_{a,b}[\ ]$ finds unique minimum value in the brackets by
changing the lattice constants $a$ and $b$ to include thermal expansion,
$E_{\mathrm{tot}}$ is the total energy of the crystal, $\hbar$ is the reduced
Planck constant, $\omega_{\lambda}$ is the phonon frequency at mode $\lambda$,
$k_{B}$ is the Boltzmann constant, and $T$ is the temperature. The second term
in Eq. (3) is temperature-independent, corresponding to the zero point energy
(ZPE) of phonons; and the last term refers to the thermally excited population
of phonons, as demonstrated by the Bose-Einstein distribution
$N_{\mathrm{ph}}$ at 300 K in Fig. 2(d)-(f).

Figure 3: Relative thermodynamic stability of monolayer fullerene networks,
one-dimensional fullerene chain and zero-dimensional fullerene molecule, with
the Gibbs free energy $F$ of monolayer qTP2 C60 set to zero to compare the
relative stability.

To quantify the relative thermodynamic stability at finite temperatures, the
difference in Gibbs free energy $\Delta F$ with respect to the free energy of
monolayer qTP2 C60 is plotted as a function of temperature $T$ for all three
phases, as illustrated in Fig. 3. With increasing temperature, the free energy
of qTP1 C60 drops faster than that of qTP2 C60 owing to their smaller
vibrational frequencies. According to Eq. (3), smaller vibrational frequencies
give rise to lower free energy and higher entropy (for details on the phonon
density of states and entropy, see the Supporting Information), which is
similar to the case in $\alpha$\- and $\beta$-tin Pavone et al. (1998); Pavone
(2001). At 150 K, monolayer qTP1 C60 becomes thermodynamically more stable
than the other two phases. At 300 K, the free energy of qTP1 C60 lies 47 meV
per formula unit (meV/f.u.) below that of qTP2 C60. At higher temperatures,
the difference becomes even larger, further stabilizing the qTP1 structure in
a thermodynamic perspective.

To further explore the thermodynamic stability of C60 in different dimensions,
the monolayer qTP1 network is further isolated into 1D qTP C60 chain and 0D
C60 molecule (for details about 1D and 3D C60, see the Supporting
Information). The cohesive energy $E_{\mathrm{c}}$ of 1D qTP C60 and 0D C60 is
listed in Table 1. Interestingly, the $E_{\mathrm{c}}$ of 1D qTP C60 chain is
higher than monolayer qTP1 and qTP2 C60 network by merely 0.3 and 0.8 meV/atom
respectively (18 and 48 meV/f.u.), whereas 11.4 meV/atom (684 meV/f.u.) lower
than the $E_{\mathrm{c}}$ of 2D qHP C60. The difference in $E_{\mathrm{c}}$
between 1D qTP C60 chain and monolayer qTP1 C60 network is even lower than the
thermal fluctuation energy $k_{\mathrm{B}}T$ at room temperature (26 meV),
implying that 2D qTP1 C60 can be transformed into 1D chains in the presence of
thermal fluctuations. Taking the finite temperature effects into account, the
Gibbs free energy of 1D qTP C60 chain and 0D C60 molecule is shown in Fig. 3.
The Gibbs free energy of 1D qTP C60 chain is higher than that of 2D qTP1 C60
in the entire temperature range ($0-900$ K), and their free energy difference
is 22 meV/f.u. at 300 K. On the other hand, the free energy of 1D qTP C60
chain becomes lower than that of 2D qTP2 C60 at temperatures above 220 K. Most
interestingly, the Gibbs free energy of 0D C60 molecule drops faster than all
the other phases below room temperature, and becomes lower than 2D qTP1 and
qTP2 C60 at 120 and 150 K respectively. However, the free energy of 2D qTP1
C60 decreases the fastest above room temperature, and consequently 2D qTP1 C60
is energetically more favored than all the other phases at temperatures above
380 K. As a result, monolayer qTP2 C60 is thermodynamically most stable at
temperatures below 150 K, 0D C60 molecule has the lowest energy for
temperatures between 150 and 380 K, and 2D qTP1 C60 is thermodynamically
favored above 380 K.

Looking back at the calculated mechanical properties and stabilities, they
seem in line with the experimental findings. It has been reported that
fullerene monolayers can only be isolated experimentally for the honeycomb
structure qHP, whereas the obtained rectangular structure qTP is few-layered
Hou et al. (2022). Although qTP2 C60 is thermodynamically favored over qTP1
C60 at low temperatures, clearly qTP1 C60 is thermodynamically more stable
than the other two phases at all temperatures above 150 K. However,
thermodynamic stability of qTP1 C60 does not guarantee high dynamic stability
in the presence of interchain (out-of-plane) vibrations perpendicular to the
quasi-1D chains. In addition, the low moduli and strength of qTP1 C60
originated from the chain crystal structures, in addition to its low shear
resistance, indicate that qTP1 C60 cannot be intrinsically resilient under
deformation. Moreover, monolayer qTP1 C60 is thermodynamically less stable
than 0D C60 molecule for temperatures between 120 and 380 K. These results
indicate the plausibility that monolayer qTP1 fullerene network can be further
split into individual 1D chains or 0D molecules by thermal fluctuations,
interchain acoustic vibrations, or external strains. In contrast, qHP C60 is
both dynamically and mechanically more stable with respect to qTP1 C60.
Therefore, monolayer polymeric C60 has so far only been exfoliated from the
quasi-hexagonal bulk single crystals. These results indicate that a systematic
analysis of mechanical, dynamic and thermodynamic stabilities is necessary to
rationalize the experimental data.

In conclusion, I carry out first principles calculations to evaluate the
mechanical, dynamic and thermodynamic stabilities of qTP1, qTP2 and qHP C60
monolayers, which have been so far believed to be stable. Electron
localization analysis reveals that the low mechanical and dynamic stabilities
in qTP1 fullerene are associated with the lack of C$-$C bonds connecting the
adjacent C60 chains, which also limits its achievable strength. Monolayer qTP2
C60 is thermodynamically more stable at temperatures below 150 K, while
thermally populated phonons hinder its thermodynamic stability with increasing
temperature. The relatively high moduli of qHP fullerene indicate that it has
a high strength because of the closely packed hexagonal fullerene network
linked through both [2+2] cycloaddition bonds and C$-$C single bonds. This, in
combination with the phonon stability, endows monolayer qHP C60 with high
stability and strength.

## Acknowledgements

I thank Prof. Bartomeu Monserrat and Dr. Ivona Bravić at the University of
Cambridge for helpful discussions. I acknowledge support from the Winton
Programme for the Physics of Sustainability, and from Magdalene College
Cambridge for a Nevile Research Fellowship. The calculations were performed
using resources provided by the Cambridge Tier-2 system, operated by the
University of Cambridge Research Computing Service (www.hpc.cam.ac.uk) and
funded by EPSRC Tier-2 capital grant EP/P020259/1, as well as with
computational support from the U.K. Materials and Molecular Modelling Hub,
which is partially funded by EPSRC (EP/P020194), for which access is obtained
via the UKCP consortium and funded by EPSRC grant ref. EP/P022561/1.

## References

  * Hou et al. (2022) Hou, L.; Cui, X.; Guan, B.; Wang, S.; Li, R.; Liu, Y.; Zhu, D.; Zheng, J. Synthesis of a monolayer fullerene network. _Nature_ 2022, _606_ , 507–510.
  * Pontiroli et al. (2013) Pontiroli, D.; Aramini, M.; Gaboardi, M.; Mazzani, M.; Gorreri, A.; Riccò, M.; Margiolaki, I.; Sheptyakov, D. Ionic conductivity in the Mg intercalated fullerene polymer Mg2C60. _Carbon_ 2013, _51_ , 143–147.
  * Tanaka and Yamanaka (2018) Tanaka, M.; Yamanaka, S. Vapor-Phase Growth and Structural Characterization of Single Crystals of Magnesium Doped Two-Dimensional Fullerene Polymer Mg2C60. _Crystal Growth & Design_ 2018, _18_ , 3877–3882.
  * Gottfried (2022) Gottfried, J. M. Molecular soccer balls connected to make a 2D material. _Nature_ 2022, _606_ , 470–471.
  * Kroto et al. (1985) Kroto, H. W.; Heath, J. R.; O’Brien, S. C.; Curl, R. F.; Smalley, R. E. C60: Buckminsterfullerene. _Nature_ 1985, _318_ , 162–163.
  * Kroto (1987) Kroto, H. W. The stability of the fullerenes Cn, with $n$ = 24, 28, 32, 36, 50, 60 and 70. _Nature_ 1987, _329_ , 529–531.
  * Kroto (1990) Kroto, H. C60, fullerenes, giant fullerenes and soot. _Pure and Applied Chemistry_ 1990, _62_ , 407–415.
  * Goroff (1996) Goroff, N. S. Mechanism of Fullerene Formation. _Acc. Chem. Res._ 1996, _29_ , 77–83.
  * Bernal et al. (2019) Bernal, J. J.; Haenecour, P.; Howe, J.; Zega, T. J.; Amari, S.; Ziurys, L. M. Formation of Interstellar C60 from Silicon Carbide Circumstellar Grains. _The Astrophysical Journal Letters_ 2019, _883_ , L43.
  * Tromer et al. (2022) Tromer, R. M.; Ribeiro, L. A.; Galvão, D. S. A DFT study of the electronic, optical, and mechanical properties of a recently synthesized monolayer fullerene network. _Chemical Physics Letters_ 2022, _804_ , 139925\.
  * Yu et al. (2022) Yu, L.; Xu, J.; Peng, B.; Qin, G.; Su, G. Anisotropic Optical, Mechanical, and Thermoelectric Properties of Two-Dimensional Fullerene Networks. _J. Phys. Chem. Lett._ 2022, _13_ , 11622–11629.
  * Yuan et al. (2022) Yuan, D.; Pi, H.; Jiang, Y.; Zhou, Y. H. L.; Jia, Y.; Su, G.; Fang, Z.; Weng, H.; Ren, X.; Zhang, W. Highly in-plane anisotropic optical properties of fullerene monolayers. _arXiv:2207.11366_ 2022,
  * Peng (2022) Peng, B. Monolayer Fullerene Networks as Photocatalysts for Overall Water Splitting. _J. Am. Chem. Soc._ 2022, _144_ , 19921–19931.
  * Ying et al. (2023) Ying, P.; Dong, H.; Liang, T.; Fan, Z.; Zhong, Z.; Zhang, J. Atomistic insights into the mechanical anisotropy and fragility of monolayer fullerene networks using quantum mechanical calculations and machine-learning molecular dynamics simulations. _Extreme Mechanics Letters_ 2023, _58_ , 101929\.
  * Dong et al. (2022) Dong, H.; Cao, C.; Ying, P.; Fan, Z.; Qian, P.; Su, Y. Anisotropic and high thermal conductivity in monolayer quasi-hexagonal fullerene: A comparative study against bulk phase fullerene. _arXiv:2208.03982_ 2022,
  * Ribeiro et al. (2022) Ribeiro, L.; Pereira, M.; Giozza, W.; Tromer, R.; Galvão, D. S. Thermal stability and fracture patterns of a recently synthesized monolayer fullerene network: A reactive molecular dynamics study. _Chemical Physics Letters_ 2022, _807_ , 140075.
  * Peng et al. (2017) Peng, B.; Zhang, H.; Shao, H.; Ning, Z.; Xu, Y.; Ni, G.; Lu, H.; Zhang, D. W.; Zhu, H. Stability and strength of atomically thin borophene from first principles calculations. _Materials Research Letters_ 2017, _5_ , 399–407.
  * Malyi et al. (2019) Malyi, O. I.; Sopiha, K. V.; Persson, C. Energy, Phonon, and Dynamic Stability Criteria of Two-Dimensional Materials. _ACS Appl. Mater. Interfaces_ 2019, _11_ , 24876–24884.
  * Luo et al. (2022) Luo, D.; Yin, K.; Dronskowski, R. Existence of BeCN2 and Its First-Principles Phase Diagram: Be and C Introducing Structural Diversity. _J. Am. Chem. Soc._ 2022, _144_ , 5155–5162.
  * Pallikara et al. (2022) Pallikara, I.; Kayastha, P.; Skelton, J. M.; Whalley, L. D. The physical significance of imaginary phonon modes in crystals. _Electron. Struct._ 2022, _4_ , 033002.
  * Pavone et al. (1998) Pavone, P.; Baroni, S.; de Gironcoli, S. $\alpha$$\leftrightarrow$$\beta$ phase transition in tin: A theoretical study based on density-functional perturbation theory. _Phys. Rev. B_ 1998, _57_ , 10421–10423.
  * Pavone (2001) Pavone, P. Old and new aspects in lattice-dynamical theory. _Journal of Physics: Condensed Matter_ 2001, _13_ , 7593–7610.
  * Masago et al. (2006) Masago, A.; Shirai, K.; Katayama-Yoshida, H. Crystal stability of $\alpha$\- and $\beta$-boron. _Phys. Rev. B_ 2006, _73_ , 104102.
  * van Setten et al. (2007) van Setten, M. J.; Uijttewaal, M. A.; de Wijs, G. A.; de Groot, R. A. Thermodynamic Stability of Boron: The Role of Defects and Zero Point Motion. _J. Am. Chem. Soc._ 2007, _129_ , 2458–2465.
  * Stoffel et al. (2010) Stoffel, R. P.; Wessel, C.; Lumey, M.-W.; Dronskowski, R. Ab Initio Thermochemistry of Solid-State Materials. _Angew. Chem. Int. Ed._ 2010, _49_ , 5242–5266.
  * Zhang et al. (2011) Zhang, R. F.; Lin, Z. J.; Mao, H.-K.; Zhao, Y. Thermodynamic stability and unusual strength of ultra-incompressible rhenium nitrides. _Phys. Rev. B_ 2011, _83_ , 060101.
  * Deringer et al. (2014) Deringer, V. L.; Stoffel, R. P.; Dronskowski, R. Thermochemical Ranking and Dynamic Stability of TeO2 Polymorphs from Ab Initio Theory. _Crystal Growth & Design_ 2014, _14_ , 871–878.
  * Deringer et al. (2014) Deringer, V. L.; Stoffel, R. P.; Dronskowski, R. Vibrational and thermodynamic properties of GeSe in the quasiharmonic approximation. _Phys. Rev. B_ 2014, _89_ , 094303.
  * White et al. (2015) White, M. A.; Cerqueira, A. B.; Whitman, C. A.; Johnson, M. B.; Ogitsu, T. Determination of Phase Stability of Elemental Boron. _Angew. Chem. Int. Ed._ 2015, _54_ , 3626–3629.
  * Nyman et al. (2016) Nyman, J.; Pundyke, O. S.; Day, G. M. Accurate force fields and methods for modelling organic molecular crystals at finite temperatures. _Phys. Chem. Chem. Phys._ 2016, _18_ , 15828–15837.
  * Skelton et al. (2017) Skelton, J. M.; Burton, L. A.; Oba, F.; Walsh, A. Chemical and Lattice Stability of the Tin Sulfides. _J. Phys. Chem. C_ 2017, _121_ , 6446–6454.
  * Pallikara and Skelton (2021) Pallikara, I.; Skelton, J. M. Phase stability of the tin monochalcogenides SnS and SnSe: a quasi-harmonic lattice-dynamics study. _Phys. Chem. Chem. Phys._ 2021, _23_ , 19219–19236.
  * Bartel (2022) Bartel, C. J. Review of computational approaches to predict the thermodynamic stability of inorganic solids. _Journal of Materials Science_ 2022, _57_ , 10475–10498.
  * Kresse and Furthmüller (1996) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Phys. Rev. B_ 1996, _54_ , 11169–11186.
  * Kresse and Furthmüller (1996) Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. _Computational Materials Science_ 1996, _6_ , 15 – 50\.
  * Blöchl (1994) Blöchl, P. E. Projector augmented-wave method. _Phys. Rev. B_ 1994, _50_ , 17953–17979.
  * Kresse and Joubert (1999) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. _Phys. Rev. B_ 1999, _59_ , 1758–1775.
  * Perdew et al. (2008) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. _Phys. Rev. Lett._ 2008, _100_ , 136406.
  * Momma and Izumi (2011) Momma, K.; Izumi, F. VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data. _Journal of Applied Crystallography_ 2011, _44_ , 1272–1276.
  * Becke and Edgecombe (1990) Becke, A. D.; Edgecombe, K. E. A simple measure of electron localization in atomic and molecular systems. _The Journal of Chemical Physics_ 1990, _92_ , 5397–5403.
  * Savin et al. (1992) Savin, A.; Jepsen, O.; Flad, J.; Andersen, O. K.; Preuss, H.; von Schnering, H. G. Electron Localization in Solid-State Structures of the Elements: the Diamond Structure. _Angew. Chem. Int. Ed. Engl._ 1992, _31_ , 187–188.
  * Gatti (2005) Gatti, C. Chemical bonding in crystals: new directions. _Zeitschrift für Kristallographie_ 2005, _220_ , 399–457.
  * Chen and Kamran (2013) Chen, K.; Kamran, S. Bonding Characteristics of TiC and TiN. _Modeling and Numerical Simulation of Material Science_ 2013, _3_ , 7–11.
  * Le Page and Saxe (2002) Le Page, Y.; Saxe, P. Symmetry-general least-squares extraction of elastic data for strained materials from ab initio calculations of stress. _Phys. Rev. B_ 2002, _65_ , 104104.
  * Wu et al. (2005) Wu, X.; Vanderbilt, D.; Hamann, D. R. Systematic treatment of displacements, strains, and electric fields in density-functional perturbation theory. _Phys. Rev. B_ 2005, _72_ , 035105.
  * Blonsky et al. (2015) Blonsky, M. N.; Zhuang, H. L.; Singh, A. K.; Hennig, R. G. Ab Initio Prediction of Piezoelectricity in Two-Dimensional Materials. _ACS Nano_ 2015, _9_ , 9885–9891.
  * Pashartis et al. (2022) Pashartis, C.; van Setten, M.; Houssa, M.; Pourtois, G. On the elastic tensors of ultra-thin films: A study of ruthenium. _Applied Surface Science_ 2022, _592_ , 153194.
  * Born and Huang (1954) Born, M.; Huang, K. _Dynamical theory of crystal lattices_ ; Clarendon Press: Oxford, 1954.
  * Wu et al. (2007) Wu, Z.-j.; Zhao, E.-j.; Xiang, H.-p.; Hao, X.-f.; Liu, X.-j.; Meng, J. Crystal structures and elastic properties of superhard $\mathrm{Ir}{\mathrm{N}}_{2}$ and $\mathrm{Ir}{\mathrm{N}}_{3}$ from first principles. _Phys. Rev. B_ 2007, _76_ , 054115.
  * Andrew et al. (2012) Andrew, R. C.; Mapasha, R. E.; Ukpong, A. M.; Chetty, N. Mechanical properties of graphene and boronitrene. _Phys. Rev. B_ 2012, _85_ , 125428\.
  * Baroni et al. (2001) Baroni, S.; de Gironcoli, S.; Dal Corso, A.; Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. _Rev. Mod. Phys._ 2001, _73_ , 515–562.
  * Gonze (1995) Gonze, X. Perturbation expansion of variational principles at arbitrary order. _Phys. Rev. A_ 1995, _52_ , 1086–1095.
  * Gonze (1995) Gonze, X. Adiabatic density-functional perturbation theory. _Phys. Rev. A_ 1995, _52_ , 1096–1114.
  * Peng et al. (2018) Peng, B.; Xu, K.; Zhang, H.; Ning, Z.; Shao, H.; Ni, G.; Li, J.; Zhu, Y.; Zhu, H.; Soukoulis, C. M. 1D SbSeI, SbSI, and SbSBr With High Stability and Novel Properties for Microelectronic, Optoelectronic, and Thermoelectric Applications. _Adv. Theory Simul._ 2018, _1_ , 1700005\.
  * Lange et al. (2022) Lange, G. F.; Bouhon, A.; Monserrat, B.; Slager, R.-J. Topological continuum charges of acoustic phonons in two dimensions and the Nambu-Goldstone theorem. _Phys. Rev. B_ 2022, _105_ , 064301.
  * Huang et al. (2015) Huang, L.-F.; Gong, P.-L.; Zeng, Z. Phonon properties, thermal expansion, and thermomechanics of silicene and germanene. _Phys. Rev. B_ 2015, _91_ , 205433.
  * Huang et al. (2016) Huang, L.-F.; Lu, X.-Z.; Tennessen, E.; M.Rondinelli, J. An efficient ab-initio quasiharmonic approach for the thermodynamics of solids. _Computational Materials Science_ 2016, _120_ , 84–93.
  * Peng et al. (2019) Peng, B.; Bravić, I.; MacManus-Driscoll, J. L.; Monserrat, B. Topological semimetallic phase in ${\mathrm{PbO}}_{2}$ promoted by temperature. _Phys. Rev. B_ 2019, _100_ , 161101.
  * Feng et al. (2012) Feng, J.; Xiao, B.; Zhou, R.; Pan, W.; Clarke, D. R. Anisotropic elastic and thermal properties of the double perovskite slab-rock salt layer Ln2SrAl2O7 (Ln=La, Nd, Sm, Eu, Gd or Dy) natural superlattice structure. _Acta Materialia_ 2012, _60_ , 3380–3392.
  * Dove (1993) Dove, M. T. _Introduction to Lattice Dynamics_ ; Cambridge University Press, 1993.
  * Togo et al. (2008) Togo, A.; Oba, F.; Tanaka, I. First-principles calculations of the ferroelastic transition between rutile-type and CaCl2-type SiO2 at high pressures. _Phys. Rev. B_ 2008, _78_ , 134106.
  * Togo and Tanaka (2015) Togo, A.; Tanaka, I. First principles phonon calculations in materials science. _Scripta Materialia_ 2015, _108_ , 1–5.