# Circularly symmetric thin-shell wormholes in F(R) gravity
with (2+1)-dimensions

Cecilia Bejarano, Ernesto F. Eiroa, Griselda Figueroa-Aguirre
Instituto de Astronomía y Física del Espacio (IAFE, CONICET-UBA),
Casilla de Correo 67, Sucursal 28, 1428, Buenos Aires, Argentina e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Within the framework of $F(R)$ theories of gravity with (2+1)-dimensions and
constant scalar curvature $R$, we construct a family of thin-shell wormholes
with circular symmetry and we analyze the stability of the static
configurations under radial perturbations. We show an example of
asymptotically anti-de Sitter thin-shell wormholes with charge, finding that
stable configurations with normal matter are possible for a suitable range of
the parameters.

Keywords: Wormholes, F(R) gravity, Thin shells.

##  1 Introduction

One of the main goals of alternative theories of gravity is focused on some
fundamental issues in the current cosmological picture, from the early-time
inflation scenario to the late-time cosmic accelerated expansion. Basically,
the intention is to modify the geometrical side of the gravitational field
equations without adding any unknown matter-energy content. Certainly, any
alternative approach also intends to provide a geometrical explanation for the
spacetime singularities as well as a cornerstone for quantum version of
gravity theory. There is a plethora of different modified gravity models. In
particular, $F(R)$ gravity [1, 2, 3] is one of the most straightforward
modifications to General Relativity (GR) since the Einstein-Hilbert action is
replaced by

$S=\frac{1}{2\kappa}\int d^{4}x\sqrt{|g|}\left(F(R)+\mathcal{L}_{m}\right),$
(1)

with $\kappa=8\pi$, and $R$ is the scalar curvature. We adopt units such that
c = G = 1, with c the speed of light and G the gravitational constant. In the
metric formalism, in which the metric tensor is the only dynamical variable to
perform the variation of the action, the field equations read

$F^{\prime}(R)R_{\mu\nu}-\dfrac{g_{\mu\nu}}{2}F(R)-\nabla_{\mu}\nabla_{\nu}F^{\prime}(R)+g_{\mu\nu}\nabla_{\rho}\nabla^{\rho}F^{\prime}(R)=\kappa
T_{\mu\nu},$ (2)

where $\nabla$ is the covariant derivative and $T_{\mu\nu}$ is the energy-
momentum tensor.

Here, we present charged thin-shell wormholes in $F(R)$ gravity in (2+1)
dimensions, finding stable static solutions for certain values of the
parameters. The construction follows the junction formalism [4, 5, 6, 7] which
is also applied, for instance, to study thin layers of matter around (i)
vacuum (namely, bubbles) or (ii) black holes. Bubbles made of pure double
layer in quadratic $F(R)$ have been studied a few years ago [8]; and layers of
matter surrounding black holes have also been explored [9, 10, 11]. The
stability of (3+1)-dimensional thin-shell wormholes in $F(R)$ gravity has been
recently investigated [12, 13].

##  2 Junction conditions

The cut-and-paste formalism used to construct some wormhole geometries is
based on the junction conditions of the theory by allowing to match different
solutions across a hypersurface $\Sigma$, which represents a boundary if there
is not matter-energy content or a thin shell of matter otherwise. We adopt the
notation that the jump of any quantity $\Upsilon$ across $\Sigma$ is defined
by $[\Upsilon]\equiv(\Upsilon^{2}-\Upsilon^{1})|_{\Sigma}$.

###  2.1 General Relativity

In GR, the junction conditions are expressed by the Darmois-Israel formalism
[4, 5] which can be sketched as follows:

  * •

Two manifolds $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are glued at a
hypersurface $\Sigma$ to give ${\cal M}=\mathcal{M}_{1}\cup\mathcal{M}_{2}$.

  * •

The first fundamental form (the induced metric) at $\Sigma$ is

$h^{1,2}_{ij}=\left.g^{1,2}_{\mu\nu}\frac{\partial
X^{\mu}_{1,2}}{\partial\xi^{i}}\frac{\partial
X^{\nu}_{1,2}}{\partial\xi^{j}}\right|_{\Sigma},$ (3)

where ${X^{\mu}_{1,2}}$ are the coordinates that describe $\mathcal{M}_{1,2}$,
respectively; and ${\xi^{a}}$ are the intrinsic local coordinates at the
hypersurface, hence $X^{\mu}_{1,2}=X^{\mu}_{1,2}(\xi^{a})$.

  * •

The second fundamental form (the extrinsic curvature) at $\Sigma$ is given by

$K_{ij}^{1,2}=-n_{\gamma}^{1,2}\left.\left(\frac{\partial^{2}X^{\gamma}_{1,2}}{\partial\xi^{i}\partial\xi^{j}}+\Gamma_{\alpha\beta}^{\gamma}\frac{\partial
X^{\alpha}_{1,2}}{\partial\xi^{i}}\frac{\partial
X^{\beta}_{1,2}}{\partial\xi^{j}}\right)\right|_{\Sigma},$ (4)

with $n_{\gamma}^{1,2}$ the unit normal ($n^{\gamma}n_{\gamma}=1$) pointing
from $\mathcal{M}_{1}$ to $\mathcal{M}_{2}$.

  * •

The matching between $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ at $\Sigma$,
should fulfill the so-called junction conditions:

– The line element should be continuous across $\Sigma$:

$[h_{\mu\nu}]=0.$ (5)

– The field equations at $\Sigma$ (Lanczos equations) reduce to

$\kappa S_{\mu\nu}=-[K_{\mu\nu}]-h_{\mu\nu}K^{\gamma}_{\;\;\gamma},\quad
n^{\mu}S_{\mu\nu}=0$ (6)

where $S_{\mu\nu}$ is the energy-momentum tensor at the gluing hypersurface.

The hypersurface $\Sigma$ is called a boundary surface if $S_{\mu\nu}=0$, or
thin shell if $S_{\mu\nu}\neq 0$.

###  2.2 $F(R)$ gravity

The formalism for the junction conditions in $F(R)$ gravity was developed in
recent years [6, 7]. In this theory, there are two separate cases to take into
account.

#### General (non-quadratic) case:

$F^{\prime\prime\prime}(R)\neq 0$

At the matching hypersurface $\Sigma$, the continuity of the following
quantities should be accomplished

$[h_{\mu\nu}]=0,\quad[K^{\mu}{}_{\mu}]=0,\quad[R]=0,$ (7)

and the field equations result

$\kappa
S_{\mu\nu}=-F^{\prime}(R)[K_{\mu\nu}]+F^{\prime\prime}(R)[n^{\rho}\nabla_{\rho}R]h_{\mu\nu},\quad
n^{\mu}S_{\mu\nu}=0.$ (8)

#### Quadratic case:

$F^{\prime\prime\prime}(R)=0$

In this case, there is one less restriction and the discontinuity of $R$ is
allowed

$[h_{\mu\nu}]=0,\quad[K^{\mu}{}_{\mu}]=0.$ (9)

By taking $F(R)=\alpha R^{2}+R-2\Lambda$, the field equations read

$\kappa
S_{\mu\nu}=-[K_{\mu\nu}]+2\alpha([n^{\gamma}\nabla_{\gamma}R]h_{\mu\nu}-[RK_{\mu\nu}]),\quad
n^{\mu}S_{\mu\nu}=0.$ (10)

In order to guarantee the divergence-free energy-momentum tensor (locally
conservation) at $\Sigma$, it is necessary to contemplate the following extra
contributions [14, 15]

$-$ External energy flux vector:

$\kappa\mathcal{T}_{\mu}=-2\alpha\nabla_{\mu}[R],\quad
n^{\mu}\mathcal{T}_{\mu}=0,$

$-$ External scalar pressure/tension:

$\kappa\mathcal{T}=2\alpha[R]K^{\gamma}{}_{\gamma},$

$-$ Double layer contribution:

$\kappa\mathcal{T}_{\mu\nu}=\nabla_{\gamma}\left(2\alpha[R]h_{\mu\nu}n^{\gamma}\delta^{\Sigma}\right),$

with $\delta^{\Sigma}$ the Dirac delta on $\Sigma$ (which is analogous to the
dipole distributions in classical electrodynamics). It is easy to verify that
these contributions vanish if the scalar curvature is continuous at $\Sigma$.

##  3 Charged solutions with constant curvature

$F(R)$ gravity coupled to a nonlinear electromagnetic field is described by
the action[16]

$S=\frac{1}{2\kappa}\int
d^{3}x\sqrt{|g|}(F(R)+(-\mathcal{F}_{\mu\nu}\mathcal{F}^{\mu\nu})^{s}),$ (11)

where $s$ is a positive parameter ($s\neq 1/2$) and
$\mathcal{F}_{\mu\nu}=\partial_{\mu}\mathcal{A}_{\nu}-\partial_{\nu}\mathcal{A}_{\mu}$
is the electromagnetic tensor. The corresponding energy-momentum tensor
becomes traceless for $s=3/4$, so this value is adopted in what follows.

The dynamical equations given by Eq. (2) for the gravitational sector are
completed with the electromagnetic field part

$\partial_{\mu}\left(\sqrt{-g}F^{\mu\nu}(-F_{\alpha\beta}F^{\alpha\beta})^{-1/4}\right)=0.$
(12)

For constant $R=R_{0}$, the spherically symmetric geometry [16]

$ds^{2}=-A(r)dt^{2}+A(r)^{-1}dr^{2}+r^{2}d\theta^{2},$ (13)

with

$A(r)=-M-\frac{\left(2\mathcal{Q}^{2}\right)^{3/4}}{2F^{\prime}(R_{0})r}-\frac{r^{2}R_{0}}{6},$
(14)

is the solution of the theory represented by the action (11), where
$R_{0}=6\Lambda_{eff}$ can be thought in terms of an effective cosmological
constant, $M$ is the mass and $Q$ is the charge. If $\mathcal{Q}=0$, the
vacuum static BTZ geometry of GR is obtained. Therefore, this is also a
solution in $F(R)$ gravity. Notice that this spacetime is asymptotically anti-
de Sitter for $R_{0}<0$.

For future purpose, we define an effective charge with its sign determined by
the sign of $F^{\prime}(R_{0})$

$Z\equiv\frac{\left(2\mathcal{Q}^{2}\right)^{3/4}}{2F^{\prime}(R_{0})}.$ (15)

Different values of $Z$ establish the horizon structure of the geometry given
by (14). From this spacetime, in the next Section, we construct thin-shell
wormholes by applying the junction formalism in $F(R)$ gravity and we study
the stability by following the standard potential approach [17].

##  4 Wormhole geometry: construction and stability

The thin-shell wormhole spacetime is built by glueing the two side of the
spherically symmetric solution (13) on $\Sigma$, which is defined by
$G(r_{1,2})=r-a=0$, to create a geodesically complete manifold
$\mathcal{M}=\mathcal{M}_{1}\cup\mathcal{M}_{2}$ where each side is given by
$\mathcal{M}_{1,2}=\\{X^{\alpha}=(t,r,\theta)/r\geqslant a\\}$:

  * •

The area $4\pi r^{2}$ is minimal when $r=a$, so $\mathcal{M}$ represents a
wormhole connecting two regions with a throat of radius $a$ (flare-out
condition).

  * •

$\xi^{i}=(\tau,\theta)$ are the coordinates on $\Sigma$, being $\tau$ the
proper time.

  * •

The thin-shell radius is taken as a function of time $a=a(\tau)$, in order to
analyze the stability of the static shells.

Then, the metrics at the sides of the throat are given by

$ds_{1,2}^{2}=-A_{1,2}(r)dt^{2}+A_{1,2}(r)^{-1}dr^{2}+r^{2}d\theta^{2},$ (16)

with

$A_{i}(r)=1-\frac{2M_{i}}{r}+\frac{Q_{i}^{2}}{F^{\prime}(R_{i})r^{2}}-\frac{R_{i}r^{2}}{12},\quad
i={1,2}$ (17)

In the orthonormal basis
$\\{e_{\hat{\tau}}=e_{\tau},e_{\hat{\theta}}=a^{-1}e_{\theta}\\}$, at the
hypersurface $\Sigma$, the first fundamental form results

$h^{1,2}_{\hat{\imath}\hat{\jmath}}=\mathrm{diag}(-1,1),$ (18)

and the non-null components of the second fundamental form for static
configurations are

$K_{\hat{\tau}\hat{\tau}}^{1,2}=\pm\frac{A^{\prime}_{1,2}(a_{0})}{2\sqrt{A_{1,2}(a_{0})}},\quad
K_{\hat{\theta}\hat{\theta}}^{1,2}=\mp\frac{1}{a_{0}}\sqrt{A_{1,2}(a_{0})},$
(19)

where the prime represents derivative with respect to $r$. The condition
$[K^{\hat{\imath}}_{\;\;\hat{\imath}}]=0$, gives the following relation

$\frac{A_{2}^{\prime}(a_{0})}{2\sqrt{A_{2}(a_{0})}}+\frac{A_{1}^{\prime}(a_{0})}{2\sqrt{A_{1}(a_{0})}}+\frac{1}{a_{0}}\left(\sqrt{A_{1}(a_{0})}+\sqrt{A_{2}(a_{0})}\right)=0,$
(20)

that should be always fulfilled by the throat radius $a_{0}$.

We can obtain the dynamics of the thin shell from the time dependent version
of $[K^{\hat{\imath}}_{\;\;\hat{\imath}}]=0$. Then, we analyze the stability
of the configuration by rewriting the dynamical equation of the throat in
terms of the potential

$\dot{a}^{2}=-V(a),$ (21)

where

$V(a)=-\frac{a^{2}\left(A_{1}(a)-A_{2}(a)\right)^{2}}{4a_{0}^{2}\left(\sqrt{A_{1}\left(a_{0}\right)}+\sqrt{A_{2}\left(a_{0}\right)}\right)^{2}}-\frac{a_{0}^{2}\left(\sqrt{A_{1}\left(a_{0}\right)}+\sqrt{A_{2}\left(a_{0}\right)}\right)^{2}}{4a^{2}}+\frac{A_{1}(a)+A_{2}(a)}{2}.$
(22)

It is straightforward to see that $V(a_{0})=0$ and $V^{\prime}(a_{0})=0$.
Then, the static solutions are stable under radial perturbations when
$V^{\prime\prime}(a_{0})>0$, otherwise they are unstable.

Figure 1: Wormholes in general $F(R)$ gravity symmetric across the throat.
Solid lines represent stable solutions while dashed ones unstable solutions.
Meshed zones represent matter satisfying WEC. Grey zones have no physical
meaning. Left plot: $M=0.5$; right one: $M=1$.

Figure 2: Wormholes in quadratic $F(R)$ gravity with $R_{1}\neq R_{2}$ at the
sides of the throat. Solid lines represent stable solutions while dashed ones
unstable solutions. Meshed zones represent matter satisfying WEC. Grey zones
have no physical meaning. For all plots: $M=1$, $R_{A}=(R_{1}+R_{2})/2$. Upper
row: $R_{1}=0.9R_{A}$ and $R_{2}=1.1R_{A}$; lower row: $R_{1}=0.8R_{A}$ and
$R_{2}=1.2R_{A}$. Left column: $\alpha R_{A}=-0.1$ for $Z>0$ and $\alpha
R_{A}=-0.9$ for $Z<0$; right column: $\alpha R_{A}=-0.2$ for $Z>0$ and $\alpha
R_{A}=-0.8$ for $Z<0$.

###  4.1 Wormholes in general $F(R)$ with $[R]=0$

In the first example, we adopt a symmetric construction across the throat,
i.e. $R_{1}=R_{2}=R_{0}$, $M_{1}=M_{2}=M$, and $Q_{1}=Q_{2}=Q$ so that
$A_{1}(r)=A_{2}(r)=A(r)$. Then, from Eq. (8), by using that
$S_{\hat{i}\hat{j}}=diag(\sigma_{0},p_{0})$, we get the surface energy density
and the pressure for the static solutions

$\sigma_{0}=\frac{F^{\prime}(R_{0})}{\kappa}\frac{A^{\prime}(a_{0})}{\sqrt{A(a_{0})}},$
(23) $p_{0}=-2\frac{F^{\prime}(R_{0})}{a_{0}\kappa}\sqrt{A(a_{0})}.$ (24)

If $\sigma_{0}<0$ or $\sigma_{0}+p_{0}<0$, the matter at the throat is exotic,
because it does not satisfy the weak energy condition (WEC). The continuity of
the trace of the extrinsic curvature imposes that the matter should satisfy
the equation of state $\sigma_{0}-2p_{0}=0$.

From Fig. 1, we see that there are two solutions constituted of normal matter
(satisfying the energy conditions) for a short range of $\sqrt{|R_{0}|}Z$ with
$Z<0$. Hence it requires the presence of ghost fields, i.e.
$F^{\prime}(R_{0})<0$. The solution with the largest radius is stable while
the other one is unstable under radial perturbations. The qualitative behavior
of the solutions is not affected by different values of $M$ which only
represent a change of scale.

###  4.2 Wormholes in quadratic $F(R)$ with $[R]\neq 0$

In the second example, we adopt $M_{1}=M_{2}=M$, and $Q_{1}=Q_{2}=Q$ but
$R_{1}\neq R_{2}$, from Eq. (10) the energy density and pressure at the throat
in the static configuration are

$\sigma_{0}=\frac{1+2\alpha
R_{2}}{\kappa}\left(\frac{A_{2}^{\prime}(a_{0})}{2\sqrt{A_{2}(a_{0})}}\right)+\frac{1+2\alpha
R_{1}}{\kappa}\left(\frac{A_{1}^{\prime}(a_{0})}{2\sqrt{A_{1}(a_{0})}}\right)\,$
(25) $p_{0}=-\frac{1+2\alpha
R_{2}}{\kappa}\left(\frac{\sqrt{A_{2}(a_{0})}}{a_{0}}\right)-\frac{1+2\alpha
R_{1}}{\kappa}\left(\frac{\sqrt{A_{1}(a_{0})}}{a_{0}}\right).$ (26)

In this quadratic case, there exists an external scalar pressure/tension
$\mathcal{T}_{0}$ that satisfies the equation of state
$\sigma_{0}-p_{0}=\mathcal{T}_{0}$ which results from the condition
$[K^{\mu}{}_{\mu}]=0$. The other extra contributions are:
$\mathcal{T}_{\mu}^{(0)}=0$ and
$\mathcal{T}_{\mu\nu}^{(0)}\propto(2\alpha[R]/\kappa)h_{\mu\nu}$.

From Fig. 2, we see that there are two solutions made of normal matter for a
short range of $\sqrt{|R_{A}|}Z$ with $R_{A}=(R_{1}+R_{2})/2$. Since $Z<0$,
there are ghost fields. The solution with the largest radius is stable; the
other is unstable. The larger the difference between $R_{1}$ and $R_{2}$, the
smaller the range of $\sqrt{|R_{A}|}Z$ where solutions exist. The same
qualitative behaviour is obtained when $\alpha$ decreases.

##  5 Summary

We have constructed a family of circularly symmetric wormholes in
(2+1)-dimensional $F(R)$ gravity with constant scalar curvature $R$ at the
sides of the throat. We have studied their stability under perturbations
preserving the symmetry. The junction conditions in $F(R)$ gravity are more
restrictive than in GR. In fact, the extra condition $[K^{\mu}{}_{\mu}]=0$
forces the type of equation of state: for general $F(R)$,
$\sigma_{0}-2p_{0}=0$, and for quadratic $F(R)$,
$\sigma_{0}-2p_{0}=\mathcal{T}_{0}$. We have analyzed two examples: (i)
symmetric wormholes ($R_{1}=R_{2}=R_{0}$; $M_{1}=M_{2}$; $Q_{1}=Q_{2}$) across
the throat in general $F(R)$ and (ii) wormholes asymmetric in the scalar
curvature ($R_{1}\neq R_{2}$; $M_{1}=M_{2}$ and $Q_{1}=Q_{2}$) only in
quadratic $F(R)$. In both cases we have found two solutions made of normal
matter (i.e., satisfying WEC) for a range of values with $Z<0$. However, for
these solution the presence of ghost fields is always required since $Z<0$
implies $F^{\prime}(R_{0})<0$. We found that the solution with the larger
radius is stable while the other is unstable. In particular, for quadratic
$F(R)$ as long as the difference between the constant (negative) scalar
curvatures $R_{1}$ and $R_{2}$ grows, the range of the squared charge
$\mathcal{Q}^{2}$ where stable solutions exist becomes smaller; similar
behavior can be observed when $\alpha$ decreases.

## Acknowledgments

This work has been supported by CONICET and Universidad de Buenos Aires. C. B.
thanks the partial support of the John Templeton Foundation.

## References

  * [1] T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82, 451 (2010).
  * [2] A. De Felice and S. Tsujikawa, Living Rev. Relativ. 13, 3 (2010).
  * [3] S. Nojiri, S.D. Odintsov, and V.K. Oikonomou, Phys. Rep. 692, 1 (2017).
  * [4] G. Darmois, Mémorial des Sciences Mathématiques, Fascicule XXV, Chap. V (Gauthier-Villars, Paris, 1927)
  * [5] W. Israel, II Nuovo Cimento B 44,1 (1966); 48, 463(E) (1967).
  * [6] N. Deruelle, M. Sasaki, and Y. Sendouda, Prog. Theor. Phys. 119, 237 (2008).
  * [7] J. M. M. Senovilla, Phys. Rev. D. 88, 064015 (2013).
  * [8] E. F. Eiroa, G. Figueroa-Aguirre, and J.M.M. Senovilla, Phys. Rev. D 95, 124021 (2017).
  * [9] E. F. Eiroa and G. Figueroa-Aguirre, Eur. Phys. J. C 78, 54 (2018)
  * [10] E. F. Eiroa and G. Figueroa-Aguirre, Eur. Phys. J. C 79, 171 (2019)
  * [11] E. F. Eiroa and G. Figueroa-Aguirre, Eur. Phys. J. Plus 135, 774 (2020).
  * [12] E. F. Eiroa and G. Figueroa-Aguirre, Eur. Phys. J. C 76, 132 (2016).
  * [13] E. F. Eiroa and G. Figueroa-Aguirre, Phys. Rev. D 94, 044016 (2016).
  * [14] J. M. M. Senovilla, Class. Quant. Grav. 31, 072002 (2014).
  * [15] B. Reina, J. M .M. Senovilla, and R. Vera, Class. Quant. Grav. 33, 105008 (2016).
  * [16] S. H. Hendi, B. Eslam Panah, and R. Saffari, Int. J. Mod. Phys. D 23, 1450088 (2014).
  * [17] C. Bejarano, E. F. Eiroa and G. Figueroa-Aguirre, Eur. Phys. J. C 81, 668 (2021).


# Anharmonicity and structural phase transition in the Mott insulator Cu2P2O7

Svitlana Pastukh<EMAIL_ADDRESS>Institute of Nuclear Physics,
Polish Academy of Sciences, PL-31342 Kraków, Poland  Paweł T. Jochym
Institute of Nuclear Physics, Polish Academy of Sciences, PL-31342 Kraków,
Poland  Oleksandr Pastukh  Institute of Nuclear Physics, Polish Academy of
Sciences, PL-31342 Kraków, Poland  Jan Łażewski  Institute of Nuclear Physics,
Polish Academy of Sciences, PL-31342 Kraków, Poland  Dominik Legut
IT4Innovations, VSB - Technical University of Ostrava, 70800 Ostrava, Czech
Republic  Department of Condensed Matter Physics, Faculty of Mathematics and
Physics, Charles University, 121 16 Prague, Czech Republic  Przemysław Piekarz
Institute of Nuclear Physics, Polish Academy of Sciences, PL-31342 Kraków,
Poland

###### Abstract

Ab initio investigations of structural, electronic, and dynamical properties
of the high-temperature $\beta$ phase of copper pyrophosphate were performed
using density functional theory. The electronic band structure shows the Mott
insulating state due to electron correlations in copper ions. By calculating
phonon dispersion relations, the soft mode at the A point of the Brillouin
zone was revealed, showing the dynamical instability of the $\beta$ phase at
low temperatures. The double-well potential connected with the soft mode is
derived and the mechanism of the structural phase transition to the $\alpha$
phase is discussed. The self-consistent phonon calculations based on the
temperature-dependent effective potential show the stabilization of the
$\beta$ phase at high temperatures, due to the anharmonic effects. The
pronounced temperature dependence and the large line width of the soft mode
indicate an essential role of anharmonicity in the structural phase
transition.

##  I Introduction

Transition metal pyrophosphates gained significant interest as promising
materials for energy storage devices [1], solid electrolytes [2],
supercapacitors [3], and photonics [4]. Copper pyrophosphate (Cu2P2O7) has
been extensively studied due to its magnetic, photoluminescence, dielectric,
electrochemical, and catalytic properties [5, 6, 7, 8, 9, 10, 11]. The X-ray
diffraction measurements on Cu2P2O7 revealed the structural phase transition
[12] and the strong negative thermal expansion (NTE) in a broad temperature
range 0–375 K [13]. Both these phenomena are related to the same vibrations of
oxygen atoms, but the details of their mechanisms are not completely
understood.

Copper pyrophosphate crystallizes in two different polymorphic forms: low-
temperature $\alpha$ and high-temperature $\beta$ phase, and the reversible
phase transition occurs around $T=363$ K [14]. As it was found by diffraction
studies, $\alpha$-Cu2P2O7 has the monoclinic $C2/c$ space group with Cu2+
cations in irregular octahedral coordination and the [PO4]3- tetrahedra
sharing one corner with the bridging oxygen atom O1 (see Fig. 1) [12]. At
higher temperatures, the enhanced thermal motion of O1 perpendicular to the
P–O–P bonding makes it positionally disordered and induces transformation to
the $\beta$ phase, which exhibits the monoclinic $C2/m$ structure [12, 15].
Diffuse scattering measurements detected an intermediate state in the
temperature range of 347–363 K [14]. Furthermore, the $\alpha$–$\beta$
structural transformation was classified as a second-order phase transition of
an inhomogeneous type by Robertson and Calvo [12].

Cu2P2O7 exhibits the long-range antiferromagnetic (AFM) order below the Néel
temperature $T_{\text{N}}=27$ K [5, 6]. A two-dimensional model with strongly
coupled magnetic dimers was investigated on the basis of the density
functional theory (DFT) calculations [7]. In our previous studies on
$\alpha$-Cu2P2O7, we demonstrated that local Coulomb interactions in the
Cu($3d$) states are responsible for the Mott insulating state with the AFM
order [16]. We analyzed the lattice’s dynamic properties by calculating the
phonon spectrum for the optimized structure. Recently, our findings were
confirmed by similar results obtained using the hybrid-DFT calculations [17].

The aim of the present first-principles study is to obtain detailed
information on the structural, electronic, and dynamical properties of
$\beta$-Cu2P2O7, and to explain the mechanism of the $\alpha$–$\beta$ phase
transition. The calculations of the phonon dispersion relations revealed the
phonon soft mode, which is responsible for the structural transformation. We
also demonstrate a crucial role of anharmonicity in the stabilization of the
$\beta$ phase at higher temperatures.

Figure 1: Crystal structure of Cu2P2O7. The unit cells of the $\alpha$ and
$\beta$ phases are plotted by the solid black and dashed green lines,
respectively. The following color scheme is used: Cu – pink, P – yellow, and O
– blue. The black arrows indicate oxygen displacements in the phonon soft mode
connected with the $\alpha-\beta$ phase transition. The image was rendered
using VESTA software [18].

##  II Calculation Details

The calculations were performed using the projector augmented-wave method [19]
and the generalized gradient approximation (GGA) [20] implemented in the
Vienna Ab initio Simulation Package [21]. The integration in the $k$-point
space was performed over the $4\times 4\times 8$ Monkhorst-Pack mesh [22] with
the cut-off energy equal to 500 eV. The ground state energy for the AFM order
was found to be consistent with the experimental studies [5, 6] and previous
calculations for the $\alpha$ phase [7, 16, 17]. During our investigation of
the $\alpha$ phase, we utilized the GGA method with $U=J=0$, to reveal a weak
metallic state due to persistent electronic states at the Fermi energy. To
improve this, we applied the GGA+U method [23], using specific parameters:
$U=4$ eV, $J=1$ eV for Cu(3d) states and $U=2$ eV, $J=0.5$ eV for O(2p)
states. This led us to identify an insulating state with the electronic band
gap of $\Delta E=1.66$ eV. However, it’s noteworthy that this gap is too small
in comparison to values reported for copper pyrophosphate materials [24, 25,
17]. To refine our analysis, we applied larger values of the Coulomb
parameter: $U=9$ eV for copper and $U=4$ eV for oxygen. Such values are
commonly used for superconducting cuprates. Additionally, these $U$ and $J$
values yielded good agreement with experimental lattice parameters, as
observed in the $\alpha$ phase analysis [16]. We also included the van der
Waals corrections during optimization of the crystal structure [26, 27].

The interatomic force constants and phonon energies were obtained with the
Alamode software [28]. First, we used the direct method [29] to calculate the
phonon dispersion relations and phonon density of states (PDOS) within the
harmonic approximation. The Hellmann-Feynman forces were obtained by making
small displacements of single atoms in the $2\times 1\times 2$ supercell with
88 atoms using the $2\times 2\times 2$ $k$-point grid for DFT. From this
information, the force constants and dynamical matrices were constructed. The
phonon frequencies and polarization vectors were calculated by the direct
diagonalization of the dynamical matrix. The LO-TO splitting at the $\Gamma$
point is determined by calculating the static dielectric tensor and Born
effective charges within the density functional perturbation theory [30]. This
method was used previously to study lattice dynamical properties of
$\alpha$-Cu2P2O7 [16].

In the present study, we utilize the temperature-dependent effective potential
(TDEP) approach [31, 32] to incorporate anharmonic effects. This method
determines the atomic potential from forces induced by thermal atom
displacements at finite temperatures.

Figure 2: Total and partial electronic density of states for $\beta$ phase of
Cu2P2O7.

The configurations of atomic displacements are generated by the high-
efficiency configuration space sampling (HECSS) method [33], which was
successfully applied to study the structural phase transition in the $A$V3Sb5
($A$ = K, Rb, Cs) superconductors [34] and the negative thermal expansion of
$\alpha$-Sn [35]. In order to investigate the effect of anharmonicity and the
temperature dependence of phonon frequencies, we apply the self-consistent
phonon (SCPH) method [36]. Additionally, we calculate the frequency shift and
line width of phonons using the real and imaginary parts of the phonon self-
energy resulting from the quartic and cubic anharmonicity [37].

##  III CRYSTAL AND ELECTRONIC STRUCTURE

The Cu2P2O7 crystal structure consists of pyrophosphate P2O7 layers parallel
to the (001) crystallographic plane and Cu cations coordinated to four oxygen
atoms (Fig. 1). The diphosphate anion [P2O7]4- is made up of two [PO4]3-
tetrahedra connected by a shared oxygen atom O1. This central oxygen atom in
the P–O–P bridge exhibits the large mean-square displacement in the $\alpha$
phase [13] and positional disorder in $\beta$-Cu2P2O7 [12, 15]. After the
optimization, we obtained the lattice parameters: $a=6.86$ Å, $b=8.09$ Å,
$c=4.32$ Å, and $\beta=109.6$°, in a good agreement with the experimental data
[14]. The size of the unit cell along the $c$ direction in the $\beta$ phase
is approximately two times smaller than in the $\alpha$ phase, and two other
lattice constants $a$ and $b$ have similar values.

Figure 3: (a) Phonon dispersion relations and (b) phonon density of states
obtained for $\beta$-Cu2P2O7 within the harmonic approximation (dotted blue
line) and using the TDEP+SCPH approach (rainbow color bar). (c) The energy
potential of the soft mode at the A point plotted as the function of the
amplitude of the O1 oxygen.

The total and partial element-specific electronic densities of states (EDOS)
calculated for the AFM state in $\beta$-Cu2P2O7 is presented in Fig. 2. Close
to the Fermi energy ($E_{\text{F}}$) the electronic states are dominated by
the O($2p$) orbitals. The lower Hubbard band corresponds to the peak around
$-8$ eV, and it contains mainly the majority-spin Cu($3d$) states with a small
contribution from O($2p$) orbitals. The energy gap appears between the
occupied states with both spin components and the upper Hubbard band with the
empty minority-spin Cu($3d$) states. The value of the energy gap
$E_{\text{g}}=2.94$ eV is slightly reduced compared to the $\alpha$ phase
[16], and it corresponds well to the experimental value for the hydrated
copper pyrophosphate (2.34 eV) [25].

##  IV LATTICE DYNAMICS AND STRUCTURAL PHASE TRANSITION

The phonon dispersion relations of $\beta$-Cu2P2O7 obtained within the direct
and TDEP+SCPH methods are presented in Fig. 3(a). In contrast to the $\alpha$
phase, which is stable at low temperatures [16], the calculation within the
harmonic approximation reveals the soft mode at the A point with the wave
vector $\textbf{q}=(0,0,0.5)$. In Fig. 1, we present the displacements of the
central oxygen atom O1, which show the largest amplitude in the soft mode.
This oxygen atom vibrates along the $b$ direction, perpendicular to the P–O–P
bridge, and the sign of its displacement changes in the following unit cells
along the $c$ direction. Therefore, when this soft mode condensates in the
$\alpha$ phase, the unit cell is doubled along the $c$ direction. As it was
found previously [13], the O1 atoms show very large thermal displacements
along the $b$ direction and their shifts from the central positions are
responsible for the $\alpha$-$\beta$ transition [12, 14]. Here, we confirmed
that vibrations of these atoms are connected with the soft mode, and
therefore, when the temperature is decreased, these displacements along the
$b$ direction determine the structure of the $\alpha$ phase. Indeed, by the
group theory analysis, we verified that the soft mode properly describes the
change of symmetry, $C2/m(\beta)\rightarrow C2/c(\alpha)$.

In Fig. 3(b), we show the total PDOS and partial contributions of Cu, O, and P
vibrations obtained within the direct method. In the low-frequency range, the
phonon states are dominated by Cu and O contributions. At intermediate
frequencies ($\sim$300–600 cm-1), the PDOS is dominated by vibrations of
oxygen atoms, with a smaller contribution of P and Cu atoms. In the narrow
band around 700 cm-1 mainly the vibrations of Cu atoms are present. At the
highest frequencies ($\sim$900–1200 cm-1), the vibrations of all atoms are
present. In the soft mode, we observe the dominant contribution from oxygen
atoms.

We investigate the stability of the $\beta$ phase at higher temperatures by
calculating the phonon dispersion relations using the SCPH procedure. The
potential energy was obtained within the TDEP approach up to the quartic terms
using the atomic displacements generated by the HECSS method. In calculations,
we used $N=100$ configurations of atomic displacements in the supercell with
88 atoms at $T=500$ K. The SCPH equation was solved using the $4\times 4\times
4$ $q$-point grid. Due to the thermal effects, the soft mode at the A point is
stabilized, and its frequency strongly depends on temperature. Its value goes
to zero when the temperature is decreased. This clear dependence of phonon
frequency on temperature indicates a strong anharmonicity of the soft mode.
Also, the lowest phonon modes at the V, L, and M points show a strong shift to
higher frequencies when temperature is increased. This demonstrates the
crucial role of anharmonicity in the stabilization of the high-temperature
$\beta$ phase of copper pyrophosphate.

In the final step, we calculate the phonon spectral function $S(\omega)$ using
the phonon line width ($\Gamma$) and frequency shift ($\Delta$) obtained from
the phonon self-energy including the bubble diagrams [37]. The calculation of
$S(\omega)$ is based on the SCPH results, therefore, the phonon frequencies
are equal to $\omega_{\text{SCPH}}+\Delta$. In order to study the anharmonic
behavior of the soft mode, in Fig. 4 we present the phonon profiles obtained
for the A point with $\bm{q}=(0,0,0.5)$. Most of the phonon peaks depicted in
Fig. 4 show very weak frequency shifts compared to the SCPH values, and their
line widths are rather narrow. A radically different behavior is observed for
the lowest mode, which exhibits strong damping and softening due to the
anharmonic effects. The shape of the soft-mode profile is very irregular and
asymmetric (see the inset in Fig. 4).

Figure 4: The phonon spectral function obtained from the anharmonic self-
energy at the A point with $\bm{q}=(0,0,0.5)$. The positions of phonon peaks
obtained within the SCPH method and with the self-energy correction ($\Delta$)
are shown by the dashed orange and dotted green lines, respectively.

In Fig. 5, we plot the frequency of the soft mode as a function of
temperature. The frequency goes to zero when we approach $T\sim 320$ K, and we
can treat it as the temperature of the structural phase transition. It is only
the approximate estimation of the critical temperature, but it corresponds
quite well to the experimental value.

In order to better understand the structural transformation, we analyze how
the lattice distortion induced by the soft mode modifies the total energy and
position of the O1 atom. Fig. 3(c) presents the change of total energy $\Delta
E$ calculated per one primitive cell of the $\beta$ phase as a function of the
displacement $u$ of the O1 atom, obtained using the polarization vector of the
soft mode. The central position of O1 corresponds to $u=0$. The energy
decreases, reaching the minimal value for the amplitude $u=0.23$ Å, which
agrees quite well with the position of O1 found by the diffraction
measurements for the $\beta$ phase $u_{\text{exp}}(\beta)=0.27$ Å. It
demonstrates that the calculated double-well potential corresponding to the
soft mode determines the shifts of O1 atoms in the high-temperature disordered
phase.

Figure 5: Temperature dependence of the soft-mode frequency.

Relaxation of the distorted structure leads to lower energy, which is marked
in Fig. 3(c) by the dashed violet line. However, this value is higher than the
ground state energy obtained by the full optimization starting from the
experimental values of atomic positions and lattice parameters within the
$\alpha$ phase (dashed blue line). The displacement of O1 from the central
position found for the $\alpha$ phase equals to $u_{\text{th}}(\alpha)=0.34$
Å, and it corresponds well to the experimental value,
$u_{\text{exp}}(\alpha)=0.32$ Å.

Below the phase transition, the soft mode is stabilized and changes into the
optical mode at the $\Gamma$ point. It is the lowest Raman mode of the $A_{g}$
symmetry and frequency 87.06 cm-1. This particular mode was previously
connected with the strong NTE observed in Cu2P2O7 [13]. The transverse motion
of O1 atoms triggers the twist and rotation of the CuO5 and PO4 polyhedra,
which were recognized as the inherent factors for the NTE. The diffraction
studies revealed that the lattice parameters $a$ and $c$ decrease with
increasing temperature, up to the $\alpha$-$\beta$ phase transition.

##  V CONCLUSIONS

In summary, the electronic and lattice properties of the $\beta$-Cu2P2O7 were
investigated by ab initio calculations. The obtained results were compared
with the previous studies performed for the $\alpha$ phase [16] and the
available experimental data. For both phases, the lattice parameters show a
good agreement with the experimental data if local Coulomb interactions and
van der Waals corrections are included in the calculations. The Mott
insulating band gap of $E_{g}=2.94$ eV is slightly smaller than in the
$\alpha$ phase (3.17 eV) and corresponds well to the experimental data [24,
25]. The AFM configurations found in both phases, with the magnetic moments
equal to 0.86 $\mu_{\text{B}}$, are very similar and agree with the long-range
order observed experimentally below $T_{\text{N}}=27$ K. Based on phonon
dispersion relations and PDOS, the dynamical properties of the $\beta$-Cu2P2O7
crystal were analyzed. The soft mode discovered at the A point indicates the
dynamical instability of the $\beta$ phase at low temperatures. It is
responsible for the structural transition, in which the symmetry changes from
the monoclinic $\alpha$ phase ($C2/c$) to the monoclinic $\beta$ phase
($C2/m$). The soft mode induces primarily vibrations of O1 atoms, resulting in
static displacements along the $b$ direction and doubling of the lattice
parameter $c$ in the $\alpha$ phase. The obtained double-well potential
allowed us to estimate the positions of O1 atoms in a good correspondence to
the diffraction data. The calculations performed within the TDEP+SCPH method
demonstrated the strong temperature dependence of the soft mode, which leads
to the stabilization of the $\beta$ phase at higher temperatures. The
transition temperature was estimated to 320 K, close to the experimental value
of approximately 360 K [13, 12, 15, 14]. Strong anharmonicity of the soft mode
is additionally confirmed by the large phonon linewidth obtained from the
phonon self-energy. As it was previously observed, the NTE is related to the
same vibrations of the O1 atoms along the direction perpendicular to the P–O–P
bridge. Our analysis indicates that the lowest optical $A_{g}$ mode, which is
responsible for the NTE, corresponds to the soft mode found in the $\beta$
phase. These observations may be important for further investigations of
structural phase transitions and NTE in other pyrophosphate materials.

The presented results demonstrate that the anharmonic and thermodynamic
properties of complex materials can be successfully studied within the TDEP
and SCPH methods in connection with the recently developed HECSS approach
[33]. It opens another route to investigate even more complicated materials,
molecular crystals, or heterostructures.

## ACKNOWLEDGMENTS

This work was supported by the Ministry of Education, Youth and Sports of the
Czech Republic through the e-INFRA CZ (ID:90254) project and project No.
LUASK22099. P.P. acknowledges the support of the National Science Centre (NCN,
Poland) under Project No. 2021/43/B/ST3/02166. S.P. acknowledges financial
support provided by the Polish National Agency for Academic Exchange NAWA
under the Programme STER– Internationalisation of doctoral schools, Project
no. PPI/STE/2020/1/00020.

## References

  * Pang  _et al._ [2015] H. Pang, Y.-Z. Zhang, Z. Run, W.-Y. Lai, and W. Huang, Amorphous nickel pyrophosphate microstructures for high-performance flexible solid-state electrochemical energy storage devices, Nano Energy 17, 339 (2015).
  * Venckutė  _et al._ [2015] V. Venckutė, P. Dobrovolskis, T. Šalkus, A. Kežionis, A. Dindune, Z. Kanepe, J. Ronis, K.-Z. Fung, and A. F. Orliukas, Preparation and characterization of solid electrolytes based on TiP2O7 pyrophosphate, Ferroelectrics 479, 101 (2015).
  * Patil and Patil [2023] S. S. Patil and P. S. Patil, Transition metal pyrophosphate (MxP2O7): A new arrival in hybrid supercapacitors, Chem. Eng. J. 455, 140639 (2023).
  * Laskowska  _et al._ [2020] M. Laskowska, I. Kityk, O. Pastukh, M. Dulski, M. Zubko, J. Jedryka, K. Cpałka, P. M. Zieliński, and Ł. Laskowski, Nanocomposite for photonics—nickel pyrophosphate nanocrystals synthesised in silica nanoreactors, Micro. Meso. Mater. 306, 110435 (2020).
  * Stiles and Stager [1972] J. A. R. Stiles and C. V. Stager, Magnetic structure of manganese pyrophosphate and copper pyrophosphate, Can. J. Phys. 50, 3079 (1972).
  * Stiles and Stager [1973] J. A. R. Stiles and C. V. Stager, Nuclear magnetic resonance in antiferromagnetic Cu2P2O7, Can. J. Phys. 51, 87 (1973).
  * Janson  _et al._ [2011] O. Janson, A. A. Tsirlin, J. Sichelschmidt, Y. Skourski, F. Weickert, and H. Rosner, Long-range superexchange in Cu2A2O7 (A = P, As, V) as a key element of the microscopic magnetic model, Phys. Rev. B 83, 094435 (2011).
  * Agarwal and Sankapal [2021] A. Agarwal and B. R. Sankapal, Ultrathin Cu2P2O7 nanoflakes on stainless steel substrate for flexible symmetric all-solid-state supercapacitors, Chem. Eng. J. 422, 130131 (2021).
  * In-noi  _et al._ [2021] O. In-noi, P. Daorattanachai, C. Rungnim, K. Prasitnok, B. Rungtaweevoranit, K. Faungnawakij, and P. Khemthong, Insight into fructose dehydration over lewis acid $\alpha$-Cu2P2O7 catalyst, Chem. Nano. Mat. 7, 292 (2021).
  * Sang  _et al._ [2022] J. Sang, P. Wei, T. Liu, H. Lv, X. Ni, D. Gao, J. Zhang, H. Li, Y. Zang, F. Yang, Z. Liu, G. Wang, and X. Bao, A reconstructed Cu2P2O7 catalyst for selective CO2 electroreduction to multicarbon products, Angew. Chem. Inter. Ed. 61, e202114238 (2022).
  * Akiyama  _et al._ [2022] T. Akiyama, M. Shimakawa, and S. Takenaka, Superior performance of copper phosphate $\alpha$-Cu2P2O7 catalysts for partial oxidation of methane into formaldehyde, Chem. Lett. 51, 511 (2022).
  * Robertson and Calvo [1967] B. E. Robertson and C. Calvo, The crystal structure and phase transformation of $\alpha$-Cu2P2O7, Acta Crystallographica 22, 665 (1967).
  * Shi  _et al._ [2020] N. Shi, A. Sanson, Q. Gao, Q. Sun, Y. Ren, Q. Huang, D. O. de Souza, X. Xing, and J. Chen, Strong negative thermal expansion in a low-cost and facile oxide of Cu2P2O7, J. Am. Chem. Soc. 142, 3088 (2020).
  * Pogorzelec-Glaser  _et al._ [2006] K. Pogorzelec-Glaser, A. Pietraszko, H. B., and M. Połomska, Structure and phase transitions in Cu2P2O7, Phase Transitions 79, 535 (2006).
  * Robertson and Calvo [1968] B. E. Robertson and C. Calvo, Crystal structure of $\beta$-Cu2P2O7, Can. J. Chem. 46, 605 (1968).
  * Pastukh  _et al._ [2021] S. Pastukh, M. Laskowska, M. Dulski, T. Krzykawski, K. Parlinski, and P. Piekarz, Ab initio studies for characterization and identification of nanocrystalline copper pyrophosphate confined in mesoporous silica, Nanotechnology 32, 415701 (2021).
  * Yang  _et al._ [2023] X. Yang, P. Zhang, and P. Korzhavyi, Hybrid-density functional calculations of structural, electronic, magnetic, and thermodynamic properties of $\alpha$-Cu2P2O7, Applied Sciences 13 (2023).
  * Momma and Izumi [2011] K. Momma and F. Izumi, Vesta 3 for three-dimensional visualization of crystal, volumetric and morphology data, J. Appl. Cryst. 44, 1272 (2011).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B 13, 5188 (1976).
  * Anisimov  _et al._ [1991] V. I. Anisimov, J. Zaanen, and O. K. Andersen, Band theory and mott insulators: Hubbard U instead of stoner I, Phys. Rev. B 44, 943 (1991).
  * Karaphun  _et al._ [2018] A. Karaphun, P. Chirawatkul, S. Maensiri, and E. Swatsitang, Influence of calcination temperature on the structural, morphological, optical, magnetic and electrochemical properties of Cu2P2O7 nanocrystals, Journal of Sol-Gel Science and Technology 88, 407 (2018).
  * Majtyka-Piłat  _et al._ [2022] A. Majtyka-Piłat, M. Wojtyniak, L. Laskowski, and D. Chrobak, Structure and properties of copper pyrophosphate by first-principle calculations, Materials 15 (2022).
  * Grimme [2006] S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comp. Chem. 27, 1787 (2006).
  * Bučko  _et al._ [2010] T. Bučko, J. Hafner, S. Lebègue, and J. G. Ángyán, Improved description of the structure of molecular and layered crystals: Ab initio dft calculations with van der waals corrections, J. Chem. Phys. A 114, 11814 (2010).
  * Tadano  _et al._ [2014] T. Tadano, Y. Gohda, and S. Tsuneyuki, Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations, J. Phys.: Condens. Matter 26, 225402 (2014).
  * Parlinski  _et al._ [1997] K. Parlinski, Z. Q. Li, and Y. Kawazoe, First-principles determination of the soft mode in cubic ZrO2, Phys. Rev. Lett. 78, 4063 (1997).
  * Gajdoš  _et al._ [2006] M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller, and F. Bechstedt, Linear optical properties in the projector-augmented wave methodology, Phys. Rev. B 73, 045112 (2006).
  * Hellman  _et al._ [2011] O. Hellman, I. A. Abrikosov, and S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B 84, 180301 (2011).
  * Hellman  _et al._ [2013] O. Hellman, P. Steneteg, I. A. Abrikosov, and S. I. Simak, Temperature dependent effective potential method for accurate free energy calculations of solids, Phys. Rev. B 87, 104111 (2013).
  * Jochym and Łażewski [2021] P. T. Jochym and J. Łażewski, High Efficiency Configuration Space Sampling – probing the distribution of available states, SciPost Phys. 10, 129 (2021).
  * Ptok  _et al._ [2022] A. Ptok, A. Kobiałka, M. Sternik, J. Łażewski, P. T. Jochym, A. M. Oleś, and P. Piekarz, Dynamical study of the origin of the charge density wave in $A{\mathrm{V}}_{3}{\mathrm{Sb}}_{5}$ ($A=\mathrm{K}$, Rb, Cs) compounds, Phys. Rev. B 105, 235134 (2022).
  * Jochym  _et al._ [2022] P. T. Jochym, R. K. Kremer, J. Łażewski, A. Ptok, P. Piekarz, E. Brücher, and A. M. Oleś, Influence of anharmonicity on the negative thermal expansion of $\alpha\text{$-$}\mathrm{Sn}$, Phys. Rev. Mater. 6, 113601 (2022).
  * Tadano and Tsuneyuki [2015] T. Tadano and S. Tsuneyuki, Self-consistent phonon calculations of lattice dynamical properties in cubic ${\mathrm{SrTiO}}_{3}$ with first-principles anharmonic force constants, Phys. Rev. B 92, 054301 (2015).
  * Tadano and Tsuneyuki [2018] T. Tadano and S. Tsuneyuki, First-principles lattice dynamics method for strongly anharmonic crystals, Journal of the Physical Society of Japan 87, 041015 (2018).