# Fate of symmetry protected coherence in open quantum system

Tian-Shu Deng  Institute for Advanced Study, Tsinghua University,
Beijing,100084, China  Lei Pan<EMAIL_ADDRESS>Institute for
Advanced Study, Tsinghua University, Beijing,100084, China

###### Abstract

We investigate the fate of coherence in the dynamical evolution of a symmetry
protected quantum system. Under the formalism of system-plus-bath for open
quantum system, the anti-unitary symmetry exhibits significant difference from
the unitary one in protecting initial coherence. Specifically, taking
advantage of Lindblad master equation, we find that a pure state in the
symmetry protected degenerate subspace will decohere even though both the
system Hamiltonian and system-environment interaction respect the same anti-
unitary symmetry. In contrast, the coherence will persist when the protecting
symmetry is unitary. We provide an elaborate classification table to
illustrate what kinds of symmetry combinations are able to preserve the
coherence of initial state, which is confirmed by several concrete models in
spin-$3/2$ system. Our results could help to explore the possible experimental
realization of stable time-reversal symmetric topological states.

##  I Introduction

Symmetry is one of the greatest unifying themes in modern physics and plays a
fundamental role in classifying quantum phase of matter. In condensed matter
physics, the existence of periodical table in topological free fermions is a
great example in which the anomalous quantum Hall effect, topological
insulators and topological superconductors, and many other interesting
topological phenomena are unified into a single framework of classification
theory. TopoClass1 ; TopoClass2 ; TopoClass3  Despite high degree of
universality on the classification theory, the previous researches about
symmetry analysis and topology in condensed matters mainly focused on the
isolated systemsSymClass1 ; SymClass2 ; Topo1 ; Topo2 ; Topo3 ; Topo4  or
those systems coupled to non-Hermitian potentialsNHTopo1 ; NHTopo2 ; NHTopo3 ;
NHTopo4 ; NHTopo5 ; NHTopo6 ; NHTopo7 ; NHTopo8 .

However, realistic quantum systems inevitably couple to external degrees of
freedom, which, in general, should be described by open quantum systems.
Recently, efforts have been made to study the classification problem on the
steady states in Lindblad equationsCooper2 ; Altland . The key point for
experimentally searching stable symmetry-protected topological states is
maintaining the coherence of quantum state in presence of surroundings. Thus,
a natural question is whether these symmetry protected topological states can
survive from the decoherence induced by environment. In particular, one of the
central issues for open quantum system is the decoherence dynamics concerning
how the quantum coherence evolves and vanishes which is particularly important
for the quantum information and quantum computationNielsen . A crucial
question is how to avoid decoherence which is unavoidable roadblock to quantum
information processing.

Therefore, diagnosing the stability of coherence of symmetry-protected quantum
state has fundamental importance both in theory and practice. Symmetry
analysis plays an important role in this issue. According to Wigner’s
theoremWigner , a symmetry transformation is either unitary or anti-unitary.
Recent theoretical research revealed that the coherence of states underlying
many symmetry protected features, could be fragile when the symmetry is anti-
unitary, even though both the environment and system-environment interaction
respect the same symmetry as the Hamiltonian of systemCooper . But these
features are always robust as long as the protected symmetry is unitary.
Strikingly, when applying this to topological systems protected by time-
reversal symmetry, one concludes that the topological phases could be unstable
to the perturbation of environmentCooper . Understanding these nonequilibrium
dynamics is of fundamental interest and has potential application in quantum
information.

In this work, we investigate the fate of coherence in degenerate subspace
protected by the unitary or anti-unitary symmetry. When the system is weakly
coupled to the environment, the dissipative dynamics under the Born-Markovian
approximation in this open system is governed by Lindblad master equation.
Indeed, the anomalous decoherence or degeneracy breaking in anti-unitary
symmetry has already arisen from the non-Hermitian linear response theoryNHLRT
. At this level, dynamical evolution of density matrix in degenerate subspace
is able to distinguish the difference between anti-unitary symmetry and
unitary symmetry in maintaining coherenceNHLRT2 . Besides, the decoherence is
also related to the number of coupling channels between system and environment
when the all coupling operators are Hermitian onesCooper . We will give a
classification about the maintenance of coherence with different symmetry
combinations respected by system and system-environment interaction.
Furthermore, we take several representative examples in spin-$3/2$ systems to
illustrate our classification. We also provide a general proof about the
origin of decoherence by taking advantage of non-Hermitian Linear response
theory.

The rest of this paper is organized as follows. Sec. II introduces the general
formalism. In Sec. III, we carry out research into the roles of symmetry in
dynamical evolution, and give the classification about coherence of symmetry
protected states. And then in Sec. IV, taking advantage of Lindblad master
equation, we examplify the classification table with spin-$3/2$ models by
calculating the density matrix and von Neumann entropy in degenerate subspace.
Sec. V provides a brief summary and outlook. The analytic proof is provided in
appendixes.

##  II General formalism

In this section, we derive general formalism from microscopical Hamiltonian to
study the dynamical evolution in the system protected by unitary or anti-
unitary symmetry. In order to investigate the robustness of coherence in
degenerate space, we consider a quantum system coupled to a Markov bath whose
Hamiltonian reads

$\displaystyle\hat{H}_{T}$ $\displaystyle=$
$\displaystyle\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{SB},$ (1)

where $\hat{H}_{S}$, $\hat{H}_{B}$ are Hamiltonians belonging to the system
and bath respectively, and $\hat{H}_{SB}$ denotes the interaction between
them. The coupling part $\hat{H}_{SB}$ can be decomposed as

$\displaystyle\hat{H}_{SB}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{M}\hat{\mathcal{A}}_{j}^{\dagger}\otimes\hat{\mathcal{B}}_{j}+\hat{\mathcal{A}}_{j}\otimes\hat{\mathcal{B}}_{j}^{\dagger},$
(2)

where $\hat{\mathcal{A}}_{j}$, $\hat{\mathcal{B}}_{j}$ are operators belonging
to system and bath respectively, and $M$ denotes the number of coupling
channels.

Let us focus on the situation that $\hat{H}_{S}$ has related symmetry in
consideration, and then a natural question is whether the symmetry protected
feature such as degeneracy is maintained or not. It is generally expected that
the symmetry protected feature would be destroyed if the system-bath coupling
$\hat{H}_{SB}$ breaks the related symmetry and it survives as long as
$\hat{H}_{SB}$ respects the same symmetry. This intuitive principle seems to
be always true. However, recently, in their seminal workCooper , McGinley and
Cooper found that it can fail when the respecting symmetry is anti-unitary
even each part ($\hat{\mathcal{A}}_{j}$, $\hat{\mathcal{B}}_{j}$ and
$\hat{H}_{B}$) obeys the same symmetry. This unexpected discovery is deeply
rooted in the Schur’s Lemma for anti-unitary groupanti-Schur ; IrePre  and
immediately indicates the fragility of time-reversal symmetry protected
topological edge statesNHLRT2 . Here we will investigate the fate of coherence
in degenerate subspace of the system. For concreteness, we consider the
following total Hamiltonian

$\displaystyle\hat{H}_{T}=\hat{H}_{S}+\sum_{j=1}^{M}\sum_{\alpha}g_{j,\alpha}\hat{\mathcal{O}}_{j}^{\dagger}\hat{b}_{\alpha}+g_{j,\alpha}^{*}\hat{\mathcal{O}}_{j}\hat{b}_{\alpha}^{\dagger}+\sum_{\alpha}\omega_{\alpha}\hat{b}_{\alpha}^{\dagger}\hat{b}_{\alpha},$
(3)

where the bath is considered as the reservoir of harmonic oscillators in
thermal equilibrium and $\hat{b}_{\alpha}$ ($\hat{b}_{\alpha}^{\dagger}$) is
annihilation (creation) operator of the bath for $\alpha$ mode with bosonic
commutation relation
$\left[\hat{b}_{\alpha},\hat{b}_{\beta}^{\dagger}\right]=\delta_{\alpha\beta}$
and the bath is considered as the reservoir of harmonic oscillators
$\hat{H}_{B}=\sum_{\alpha}\omega_{\alpha}\hat{b}_{\alpha}^{\dagger}\hat{b}_{\alpha}$.
Here we focus on the scenario in which the system Hamiltonian $\hat{H}_{S}$
possesses the symmetry in question exhibiting degeneracy and the coupling
$\hat{H}_{SB}$ and bath $\hat{H}_{B}$ respect the same symmetry. That is to
say the symmetry considered here is represented by a group
$\mathbf{\mathcal{G}}$ and $\hat{\mathcal{O}}$, $\hat{b}_{\alpha}$ are all
invariant under the symmetry operations in $\mathbf{\mathcal{G}}$.

Under this circumstance, a significant practical question is that whether
symmetry protected properties such as degeneracy in system are fundamentally
stable against to perturbations of the environment. To explore this, we
consider the one-channel case $M=1$ for simplicity. Following the standard
procedure to integrate out the bath under Markovian approximation (low
temperature bath and constant noise spectrum) and Born approximation (up to
$g^{2}$ order), one can derive the Lindblad master equation as

$\displaystyle\frac{d\hat{\rho}(t)}{dt}$
$\displaystyle\equiv\mathscr{L}\hat{\rho}(t)$
$\displaystyle=-i\left[\hat{H}_{S},\hat{\rho}(t)\right]-\gamma\left\\{{\hat{\rho}}(t),\hat{\mathcal{O}}^{\dagger}\hat{\mathcal{O}}\right\\}+2\gamma\hat{\mathcal{O}}{\hat{\rho}}(t)\hat{\mathcal{O}}^{\dagger},$
(4)

where $\mathscr{L}$ is Liouvillian superoperator and
$\\{\hat{A},\hat{B}\\}\equiv\hat{A}\hat{B}+\hat{B}\hat{A}$ denotes the anti-
commutator. This equation dominates the dynamics of density matrix. In this
way, the symmetry protected properties in degenerate space are connected with
the dynamics governed by Lindblad master equation. Here we focus on the von
Neumann entropy in degenerate subspace to characterize decoherent process and
the corresponding response of entropy is given by

$\displaystyle\delta S(t)$ $\displaystyle=$ $\displaystyle S(t)-S_{0}(t)$ (5)

where $S_{0}(t)$ is the unperturbed entropy with coherent evolution determined
by $\hat{H}_{S}$. Specifically, we will study von Neumann entropy $S_{\rm
v}(t)=-{\rm Tr}[\hat{\rho}(t)\log\hat{\rho}(t)]$ which can be obtained once
the master equation is solved. The decoherence dynamics in degenerate subspace
and entropy growth can be derived by non-Hermitian linear response theory.
Moreover, the break of degeneracy can be reflected in matrix representation of
Liouvillian $\mathscr{L}$. Using the Choi-Jamiołkowski isomorphismChoi ;
Jamiolkowski , Eq.(4) can be mapped to the following equation

$\displaystyle\frac{d}{dt}\ket{\rho}$
$\displaystyle=\hat{\mathscr{L}}\ket{\rho},$ (6)

with vectorized density matrix $\ket{\rho}$

$\displaystyle\ket{\rho}=\sum_{i,j}\rho_{i,j}|i\rangle\otimes|j\rangle$ (7)

where $\rho_{i,j}=\bra{i}\hat{\rho}\ket{j}$ is matrix element of $\hat{\rho}$.
And $\hat{\mathscr{L}}$ denotes the matrix representation of Liouvillian which
is written as

$\displaystyle\hat{\mathscr{L}}$ $\displaystyle=$
$\displaystyle-i\left(\hat{H}_{S}\otimes\hat{I}-\hat{I}\otimes\hat{H}^{\mathrm{T}}_{S}\right)$
(8) $\displaystyle+$
$\displaystyle\gamma\left[2\hat{\mathcal{O}}\otimes\hat{\mathcal{O}}^{*}-\hat{\mathcal{O}}^{\dagger}\hat{\mathcal{O}}\otimes\hat{I}-\hat{I}\otimes\left(\hat{\mathcal{O}}^{\dagger}\hat{\mathcal{O}}\right)^{\mathrm{T}}\right].$

The Liouvillian superoperator shares the symmetry associated with
$\hat{\mathcal{O}}$ which determines the symmetry of Eq.(4). The fate of
coherence in degenerate subspace depends on whether $\hat{\mathscr{L}}$ is
proportional to identity matrix in degenerate subspace.

##  III Coherence analysis for different symmetry combinations

In this section, we classify the coherent dynamics in degenerate subspace
regarding different combinations of symmetries respected by $\hat{H}_{S}$ and
the $\hat{\mathcal{O}}$. Supposing the system respects unitary or anti-unitary
symmetry characterized by group $\mathbf{\mathcal{G}}$. This means the matrix
representation of $\hat{H}_{S}$ satisfies
$[\hat{H}_{S},\hat{U}_{\mathbf{\mathcal{G}}}]=0$ for any group element
$\hat{U}_{\mathbf{\mathcal{G}}}\in\mathbf{\mathcal{G}}$. If $|\psi\rangle$ is
one of eigenstates of $\hat{H}_{S}$, i.e.,
$\hat{H}_{S}|\psi\rangle=E|\psi\rangle$, then
$\hat{U}_{\mathbf{\mathcal{G}}}|\psi\rangle$ would also be the eigenstate of
$\hat{H}_{S}$. Thus $\\{\hat{U}_{\mathbf{\mathcal{G}}}|\psi\rangle\\}$ spans a
degenerate subspace when
$\hat{U}_{\mathbf{\mathcal{G}}}\in\mathbf{\mathcal{G}}$ satisfying
$\hat{U}_{\mathbf{\mathcal{G}}}|\psi\rangle\neq|\psi\rangle$ exists, which is
also irreducible representation subspace of $\mathbf{\mathcal{G}}$. We call
this symmetry protected degeneracy. A notable example in half-odd integer spin
system is the Kramers’ degeneracy when $\mathbf{\mathcal{G}}$ represents time
reversal symmetry group.

In order to describe the coherent dynamics in degenerate space, we assume that
the $\\{U_{G}|\psi\rangle\\}$ contains two orthogonal basis which are denoted
by $|\phi_{+}\rangle$ and $|\phi_{-}\rangle$ corresponding to two-fold
degenerate subspace. The initial state is prepared to be a pure state
$\hat{\rho}(0)=|\psi(0)\rangle\langle\psi(0)|$ with
$|\psi(0)\rangle=\alpha|\phi_{+}\rangle+\beta|\phi_{-}\rangle$ and then the
dynamics of density matrix will be dominated by Lindblad master equation (4).
We focus on the density matrix in the subspace
$\hat{\rho}_{G}(t)=\hat{\Pi}_{G}\hat{\rho}(t)\hat{\Pi}_{G}$, where the ground
state projective operator is defined by
$\hat{\Pi}_{G}=|\phi_{+}\rangle\langle\phi_{+}|+|\phi_{-}\rangle\langle\phi_{-}|$.

For a general total Hamiltonian, the bath and system-bath coupling don’t
respect the related symmetry of $\hat{H}_{S}$, in which case
$\hat{\rho}_{G}(t)$ evolves into a mixed state where the decoherence happens.
Nevertheless, what we concern here is whether the complete information of
initial state, or at least the coherence, can be maintained when both the
operators {$\hat{a}_{\alpha}$} and {$\hat{\mathcal{O}}_{j}$} respect the same
symmetry as the system Hamiltonian $\hat{H}_{S}$. As stated above, the fate of
coherence is determine by unitarity or anti-unitarity of the symmetry. This
conclusion had been pointed out in Ref.Cooper ; NHLRT2  and here we provide a
brief explanation. For the unitary symmetry,
$\hat{\Pi}_{G}\hat{\mathcal{O}}^{\dagger}\hat{\mathcal{O}}\hat{\Pi}_{G}\propto\hat{\Pi}_{G}$,
$\hat{\Pi}_{G}\hat{\mathcal{O}}^{\dagger}\hat{\Pi}_{G}\propto\hat{\Pi}_{G}$,
and $\hat{\Pi}_{G}\hat{\mathcal{O}}\hat{\Pi}_{G}\propto\hat{\Pi}_{G}$ are
guaranteed by Schur’s lemmaSchur  as long as $\hat{\mathcal{O}}^{\dagger}$,
$\hat{\mathcal{O}}$ respect the same symmetry as $\hat{H}_{S}$. In this case,
Lindblad superoperator $\mathscr{L}$ only acts on the subspace density matrix
as $\hat{\Pi}_{G}\mathscr{L}[\hat{\rho}(t)]\hat{\Pi}_{G}\propto\rho_{G}(0)$.
Accordingly, if we define

$\rho_{G}(t)=\left(\begin{array}[]{cc}\rho_{++}&\rho_{+-}\\\
\rho_{-+}&\rho_{--}\end{array}\right)=\left(\begin{array}[]{cc}\langle\phi_{+}|\hat{\rho}|\phi_{+}\rangle&\langle\phi_{+}|\hat{\rho}|\phi_{-}\rangle\\\
\langle\phi_{-}|\hat{\rho}|\phi_{+}\rangle&\langle\phi_{-}|\hat{\rho}|\phi_{-}\rangle\end{array}\right),$
(9)

then all of the matrix elements synchronously decay with the same rate. That
is to say, renormalized density matrix in subspace
$\tilde{\rho}_{G}=\rho_{G}(t)/{\rm tr}(\rho_{G}(t))$ is always equal to
$\rho_{G}(0)$ under the time-evolution of Lindblad equation which means the
initial coherence is maintained. In addition, the matrix representation of
$\hat{\mathscr{L}}$ in subspace is proportional to identity, namely,

$\displaystyle\hat{\mathscr{L}}_{G}\propto\left(\begin{array}[]{cccc}1&0&0&0\\\
0&1&0&0\\\ 0&0&1&0\\\ 0&0&0&1\end{array}\right),$ (14)

which acts trivially on the state of subspace
$\hat{\mathscr{L}}_{G}\big{[}\rho_{+,+}(t),\rho_{+,-}(t),\rho_{-,+}(t),\rho_{-,-}(t)\big{]}^{\rm
T}\propto\big{[}\rho_{+,+}(0),\rho_{+,-}(0),\rho_{-,+}(0),\rho_{-,-}(0)\big{]}^{\rm
T}$. For the case of anti-unitary symmetry, according to the Schur’s
Lemmaanti-Schur , if $\hat{\mathcal{O}}$ is Hermitian operator obeying this
symmetry, it will be proportional to identity matrix in degenerate subspace,
i.e.,
$\hat{\Pi}_{G}\hat{\mathcal{O}}^{\dagger}\hat{\Pi}_{G}=\hat{\Pi}_{G}\hat{\mathcal{O}}\hat{\Pi}_{G}\propto\hat{I}_{G}$
and the condition (14) is also true which means the density matrix will
maintain its initial coherence. But if $\hat{\mathcal{O}}$ is non-Hermitian
operator, then its matrix presentation is proportional to identity, i.e.,
$\hat{\Pi}_{G}\hat{\mathcal{O}}\hat{\Pi}_{G}\not\propto\hat{I}_{G}$. In this
case, one can find immediately $\tilde{\rho}_{G}=\rho_{G}(t)/{\rm
tr}(\rho_{G}(t))\neq\rho_{G}(0)$ from master equation (4) and decoherence
occurs. Meanwhile $\hat{\mathscr{L}}_{G}$ is no longer proportional to
identity

$\displaystyle\hat{\mathscr{L}}_{G}\not\propto\left(\begin{array}[]{cccc}1&0&0&0\\\
0&1&0&0\\\ 0&0&1&0\\\ 0&0&0&1\end{array}\right).$ (19)

The above statements for unitary or anti-unitary symmetric systems can be
naturally generalized to those with various combinations of symmetries. We
summarize different kinds of symmetry combinations and corresponding fate of
coherence in Table 1, we take time-reversal symmetric group and quaternion
group as candidates of antiuntary and unitary group, respectively. Some
representative examples will be illustrated in next section.

Table 1: Classification table for the fate of initial coherence protected by there kinds of symmetries in $\hat{H}_{s}$ and coupled to operator $\mathcal{O}$ with different types of symmetry combinations. Here $[\hat{H}_{S},\mathcal{Q}]=0$ ($[\hat{H}_{S},\mathcal{T}]=0$) means the system respects ($\mathcal{Q}$-symmetry) time-reversal symmetry. $\hat{I}_{G}$ represents the identity matrix in ground-state subspace. Symmetry of $\hat{H}_{S}$ | Hermiticity of $\hat{\mathcal{O}}$ | Symmetry of $\hat{\mathcal{O}}$ | Coherence/Decoherence |  $\hat{\mathscr{L}}_{G}\propto\hat{I}_{G}\otimes\hat{I}_{G}$ (Yes/No)
---|---|---|---|---
$[\hat{H}_{S},\mathcal{Q}]=0$ | Hermitian | $[\hat{\mathcal{O}},\mathcal{Q}]=0$ | Coherence | Yes
$[\hat{\mathcal{O}},\mathcal{Q}]\neq 0$ | Decoherence | No
Non-Hermitian | $[\hat{\mathcal{O}},\mathcal{Q}]=0$ | Coherence | Yes
$[\hat{\mathcal{O}},\mathcal{Q}]\neq 0$ | Decoherence | No
$[\hat{H}_{S},\mathcal{T}]=0$ | Hermitian | $[\hat{\mathcal{O}},\mathcal{T}]=0$ | Coherence | Yes
$[\hat{\mathcal{O}},\mathcal{T}]\neq 0$ | Decoherence | No
Non-Hermitian | $[\hat{\mathcal{O}},\mathcal{T}]=0$ | Decoherence | No
$[\hat{\mathcal{O}},\mathcal{T}]\neq 0$ | Decoherence | No
$\begin{array}[]{c}[\hat{H}_{S},\mathcal{T}]=0\\\ {}[\hat{H}_{S},\mathcal{Q}]=0\end{array}$ | Hermitian | $[\hat{\mathcal{O}},\mathcal{T}]=0,\quad[\hat{\mathcal{O}},\mathcal{Q}]=0$ | Coherence | Yes
$[\hat{\mathcal{O}},\mathcal{T}]=0,\quad[\hat{\mathcal{O}},\mathcal{Q}]\neq 0$ | Coherence | Yes
$[\hat{\mathcal{O}},\mathcal{T}]\neq 0,\quad[\hat{\mathcal{O}},\mathcal{Q}]=0$ | Coherence | Yes
$[\hat{\mathcal{O}},\mathcal{T}]\neq 0,\quad[\hat{\mathcal{O}},\mathcal{Q}]\neq 0$ | Decoherence | No
Non-Hermitian | $[\hat{\mathcal{O}},\mathcal{T}]=0,\quad[\hat{\mathcal{O}},\mathcal{Q}]=0$ | Coherence | Yes
$[\hat{\mathcal{O}},\mathcal{T}]=0,\quad[\hat{\mathcal{O}},\mathcal{Q}]\neq 0$ | Decoherence | No
$[\hat{\mathcal{O}},\mathcal{T}]\neq 0,\quad[\hat{\mathcal{O}},\mathcal{Q}]=0$ | Coherence | Yes
$[\hat{\mathcal{O}},\mathcal{T}]\neq 0,\quad[\hat{\mathcal{O}},\mathcal{Q}]\neq 0$ | Decoherence | No

##  IV Numerical results with spin-$3/2$ models

In this section, we take spin-$3/2$ models as a concrete example to elucidate
and verify the symmetry classification shown in Table 1. This kind of high-
spin systems coupled to environment was extensively studied in nitrogen-
vacancy (NV) centersNVcenter1 ; NVcenter2 . The matrix representation of
spin-3/2 angular momentum is chosen as

$\displaystyle S_{x}=\frac{1}{2}\left(\begin{array}[]{cccc}0&\sqrt{3}&0&0\\\
\sqrt{3}&0&2&0\\\ 0&2&0&\sqrt{3}\\\ 0&0&\sqrt{3}&0\end{array}\right),$ (24)
$\displaystyle S_{y}=\frac{i}{2}\left(\begin{array}[]{cccc}0&-\sqrt{3}&0&0\\\
\sqrt{3}&0&-2&0\\\ 0&2&0&-\sqrt{3}\\\ 0&0&\sqrt{3}&0\end{array}\right),$ (29)
$\displaystyle S_{z}=\left(\begin{array}[]{cccc}\frac{3}{2}&0&0&0\\\
0&\frac{1}{2}&0&0\\\ 0&0&-\frac{1}{2}&0\\\
0&0&0&-\frac{3}{2}\end{array}\right).$ (34)

The unitary and anti-unitary symmetries will be discussed respectively.

###  IV.1 Unitary symmetry

We first study the system with unitary symmetry. In order illustrate the
symmetry protected coherence, for simplicity, we choose the quaternion group
$\mathcal{Q}$ as candidate of the unitary group which is a non-Abelian unitary
group whose irreducible representation can be two dimensions, spanning a two-
fold degenerate subspace. The matrix representation of $\mathcal{Q}$ is
displayed explicitly in Appendix B. We then construct a
$\mathcal{Q}$-symmetric Hamiltonian
$\hat{H}_{S}=E_{g}(\hat{S}_{x}\hat{S}_{y}\hat{S}_{z}+\hat{S}_{z}\hat{S}_{y}\hat{S}_{x})$
with twofold degenerate ground states which forms the two-dimensional
irreducible representation space of group $\mathcal{Q}$. For Hermitian
operator $\hat{\mathcal{O}}$ coupled to system respecting
$\mathcal{Q}$-symmetry such as $\hat{\mathcal{O}}=\hat{S}_{y}^{2}$, the
coherence in ground-state subspace maintains all the time but decoherence
occurs in the case of breaking this symmetry, say,
$\hat{\mathcal{O}}=\hat{S}_{x}\hat{S}_{y}+\hat{S}_{y}\hat{S}_{x}$ as shown in
Fig.1(a). The same situation happens to the non-Hermitian coupling operators
that the symmetric and asymmetric $\hat{\mathcal{O}}$ are chosen as
$\hat{\mathcal{O}}=\hat{S}_{x}\hat{S}_{y}\hat{S}_{z}$ and
$\hat{\mathcal{O}}=\hat{S}_{y}\hat{S}_{z}$ respectively (see Fig.1(b)). The
increasing entropy reflects the fact subspace density matrix evolves into a
mixed state. The unchanged von Neumann entropy indicates that the density
matrix in subspace is always a pure state, which means the coherence is
maintained. This is reasonable consequence that the decoherence appears only
when the system is perturbed by symmetry-broken coupling. As mentioned above,
the Schur’s lemma allows us to infer that
$\hat{\Pi}_{G}\hat{\mathcal{O}}\hat{\Pi}_{G}\propto\hat{I}_{G}$ if
$\hat{\mathcal{O}}$ respects $\mathcal{Q}$-symmetry which will lead to
$\hat{\mathscr{L}}_{G}\propto\hat{I}_{G}\otimes\hat{I}_{G}$ according to
expression (8). In contrast,
$\hat{\mathscr{L}}_{G}\not\propto\hat{I}_{G}\otimes\hat{I}_{G}$ if
$\hat{\Pi}_{G}\hat{\mathcal{O}}\hat{\Pi}_{G}\not\propto\hat{I}_{G}$ since
$\hat{\mathcal{O}}$ breaks $\mathcal{Q}$-symmetry which destroys the coherence
and the corresponding entropy of decoherent dynamics reaches steady-state
value $S_{\rm v}(\infty)=\ln 2\approx 0.693$ which is nothing but the maximum
entropy in doublet space.

Figure 1: Time evolution of von Neumann entropy in degenerate subspace
protected by unitary-symmetry ($\mathcal{Q}$-symmetry). (a) Von Neumann
entropy as function of time for hermitian coupling operator: with
$\mathcal{Q}$-symmetry $\hat{Q}=\hat{S}_{y}^{2}$ (red solid line) and without
$\mathcal{Q}$-symmetry $\hat{Q}=\hat{S}_{x}\hat{S}_{y}+\hat{S}_{y}\hat{S}_{x}$
(blue dashed line). (b) Von Neumann entropy as function of time for non-
hermitian coupling operator: with $\mathcal{Q}$-symmetry
$\hat{Q}=\hat{S}_{x}\hat{S}_{y}\hat{S}_{z}$ (red solid line) and without
$\mathcal{Q}$-symmetry $\hat{Q}=\hat{S}_{x}\hat{S}_{y}\hat{S}_{z}$ (blue
dashed line).

###  IV.2 Anti-unitary symmetry

From now on, we focus on the systems respecting time-reversal symmetry. Since
time-reversal operation inverts the angular momentum operators, i.e.,
$\hat{T}\hat{S}_{x,y,z}\hat{T}^{-1}=-\hat{S}_{x,y,z}$, one can construct time-
reversal invariant Hamiltonian as
$\hat{H}_{S}=E_{g}\\{\hat{S}_{x},\hat{S}_{z}\\}$. The ground-state subspace of
$\hat{H}_{S}$ is two-fold degenerate due to Kramers’ theorem and we denote two
degenerate states as $|\phi_{\pm}\rangle$. Similar to the discussion of
unitary symmetric system, here the hermitian and non-hermitian coupling
$\hat{\mathcal{O}}$ will be investigated respectively. In order to verify the
case of $[\hat{H}_{S},\mathcal{T}]=0$ in Table 1, we prepare the initial state
on $|\psi(0)\rangle=\alpha|\phi_{+}\rangle+\beta|\phi_{-}\rangle$, and then
calculate the time evolution of von Neumann entropy in ground-state subspace
by means of Lindblad equation shown in Fig. 2.

For Hermitian case, as shown in Fig.2(a), the initial coherence will maintain
(vanish) if the operator $\hat{\mathcal{O}}$ coupled by system respects
(breaks) time-reversal symmetry. However, there will be a dramatic difference
when this time-reversal symmetric system couples to non-Hermitian operator
$\hat{\mathcal{O}}$. Fig.2(c) plots the von Neumann entropy from which one can
see clearly that the coherence is always destroyed despite whether or not
$\hat{\mathcal{O}}$ respect time-reversal symmetry. Hence, we demonstrate the
statement that coherence could survive only if the time-reversal symmetric
coupling operator is also Hermitian. This is consistent with Schur’s Lemma for
anti-unitary group that an anti-unitary symmetric and hermitian operator is
proportional to identity in degenerate subspace, i.e.,
$\hat{\Pi}_{G}\hat{\mathcal{O}}\hat{\Pi}_{G}\propto\hat{I}_{G}$ which makes
the relation $\hat{\mathscr{L}}_{G}\propto\hat{I}_{G}\otimes\hat{I}_{G}$ true.
When the situation that any of symmetry and hermiticity is not satisfied
happens, it results in
$\hat{\Pi}_{G}\hat{\mathcal{O}}\hat{\Pi}_{G}\not\propto\hat{I}_{G}$ which
causes the decoherence and meanwhile
$\hat{\mathscr{L}}_{G}\not\propto\hat{I}_{G}\otimes\hat{I}_{G}$ .

With regard to density matrix $\hat{\rho}_{G}(t)$ in ground-state space, for
coherent evolution, all matrix elements decay with equal-rate during whole
region of time as shown in Fig.2(b). While for decoherent process, the
nondiagonal elements could decay to zero but its diagonal elements keep finite
which means all coherence is lost and the system reaches maximum entanglement
in groundstate subspace (see Fig.2(d)).

Figure 2: Time evolution of von Neumann entropy and density matrix in
degenerate subspace protected by time-reversal symmetry. (a) Von Neumann
entropy as function of time for hermitian coupling operator: with time-
reversal symmetry $\hat{Q}=S_{x}^{2}$ (red solid line) and without time-
reversal symmetry $\hat{Q}=S_{z}$ (blue dashed line). (b) Matrix elements of
density matrix associated with red solid line in (a) evolves as time. (c) Von
Neumann entropy as function of time for non-Hermitian coupling operator:
without time-reversal symmetry $\hat{Q}=S_{x}S_{y}S_{z}$ (red solid line) and
with time-reversal symmetry $\hat{Q}=iS_{z}$ (blue dashed line). (d) Matrix
elements of density matrix associated with blue dashed line in (c) evolves as
time. The initial value of density matrix is chosen as
$\ket{\psi(0)}\bra{\psi(0)}$ with
$\ket{\psi(0)}=\frac{1}{\sqrt{2}}\big{(}\ket{\phi_{+}}+\ket{\phi_{-}}\big{)}$.

###  IV.3 Both unitary and anti-unitary symmetry

There is another case in which the Hamiltonian respects both time-reversal
symmetry and $\mathcal{Q}$-symmetry, such as
$\hat{H}_{S}=E_{g}\hat{S}_{z}^{2}$ whose degeneracy of groundstates is
protected by both of symmetries. This highly symmetric system allows for more
types of symmetry combinations which contain four circumstances for both
Hermitian and non-Hermitian couplings as shown in the third row of Table 1.
The corresponding entropy dynamics associated with different symmetries is
shown in Fig. (3). For hermitian case, the initial coherence in ground-state
subspace is always maintained whichever symmetry the coupling operator
$\hat{\mathcal{O}}$ has. This is consistent with the conclusions in subsection
IV.1 and IV.2 that
$\hat{\Pi}_{G}\hat{\mathcal{O}}\hat{\Pi}_{G}\propto\hat{I}_{G}$ regardless of
unitarity or antiunitarity which signifies the decoherence occurs only when
the coupling operator breaks both symmetries, as shown in Fig.3.(a). For non-
Hermitian case, the coherence will be destroyed once the
$\mathcal{Q}$-symmetry is broken by coupling operator $\hat{\mathcal{O}}$ as
shown in Fig.3.(b). This reflects the fact that coherence protected only by
antiunitary symmetry is fragile under non-Hermitian perturbation.

Figure 3: The von Neumann entropy as a function of time with the system
Hamiltonian respecting both $\mathcal{Q}$-symmetry and time-reversal symmetry.
The Hermitian and non-Hermitian coupling $\hat{\mathcal{O}}$’s are shown in
(a) and (b) respectively. For hermitian case, the four kinds of coupling
$\hat{\mathcal{O}}$ are chosen as $\hat{S}_{x}^{2}$ (red solid line),
$\hat{S}_{x}\hat{S}_{y}+\hat{S}_{y}\hat{S}_{x}$ (blue dashed line),
$\hat{S}_{x}\hat{S}_{y}\hat{S}_{z}+\hat{S}_{z}\hat{S}_{y}\hat{S}_{x}$ (green
dotted line), and $\hat{S}_{x}$ (dot-dashed line). For non-Hermitian case, the
corresponding choice is
$i(\hat{S}_{x}\hat{S}_{y}\hat{S}_{z}-\hat{S}_{z}\hat{S}_{y}\hat{S}_{x})$,
$\hat{S}_{x}\hat{S}_{y}$, $\hat{S}_{x}\hat{S}_{y}\hat{S}_{z}$, and
$\hat{S}_{x}^{2}\hat{S}_{z}$. The initial value of density matrix is chosen as
$\ket{\psi(0)}\bra{\psi(0)}$ with
$\ket{\psi(0)}=\frac{1}{\sqrt{2}}\big{(}\ket{\phi_{+}}+\ket{\phi_{-}}\big{)}$.

##  V Summary and Outlook

We study the fate of coherence in degenerate subspace which is protected by
unitary symmetry or anti-unitary symmetry. With symmetries lying in system
Hamiltonian, and the interaction between system and environment, we have
demonstrated that the coherence could be fragile when the symmetry is anti-
unitary. We elaborate on classification of various symmetry combinations and
analyze the stability of coherence detailedly which is confirmed by several
spin-3/2 models. The results of this investigation are extensions of
Ref.Cooper  and could be applied to the investigation of robustness in time-
reversal invariant topological edge states. Recent experimental progresses in
controlling and manipulating dissipation in ultracold atoms provide
unprecedented opportunity for understanding the dynamics of open quantum
system which is driven far from of equilibriumBouganne ; Tomita1 ; Tomita2 ;
LuoL ; DissHubbard ; Takahashi1 ; Takahashi2 . We expect our work will guide
the preparation for stable time-reversal symmetry protected topological states
in laboratory.

There is another intriguing issue about the open quantum systems beyond the
Born-Markov approximation. In this case, the coherence protected by anti-
unitary symmetry could also be fragile even the system couples hermitian
operators with the same symmetryfoot2 . What’s more, the current study focuses
on non-interacting systems. An important issue in future study is exploring to
the stability of quantum degeneracy protected by anti-unitary in interacting
systems, especially to interacting topological phases in many-body system such
as Haldane phase in AKLT model chainCaiZi .

## Acknowledgements

We would like to thank Yu Chen for helpful discussions. This work is supported
by Beijing Outstanding Young Scientist Program hold by Hui Zhai. L. P
acknowledges support from the project funded by the China Postdoctoral Science
Foundation (Grant No. 2020M680496).

##  Appendix A Response of density matrix and von Neumann entropy growth

In this appendix, we derive the response of density matrix and von Neumann
entropy growth by means of NHLRT. As discussed in Ref.NHLRT , the whole system
coupling to a bath with white noise can be described by a non-Hermitian
effective Hamiltonian

$\displaystyle\hat{H}_{\rm eff}=\hat{H}_{S}+\hat{H}_{\mathrm{diss}},$ (35)

where $H_{\mathrm{diss}}=\left(-i\gamma\hat{\cal O}^{\dagger}\hat{\cal
O}+\hat{\cal O}^{\dagger}\hat{\xi}+\hat{\xi}^{\dagger}\hat{\cal O}\right)$ and
$\gamma=\pi|g|^{2}\rho$ is dissipation strength where $\rho$ is spectrum
density of bath. $\hat{\xi}(t)$, $\hat{\xi}^{\dagger}(t)$ present the Langevin
noise operators which obey the following relations

$\displaystyle\langle\hat{\xi}(t)\hat{\xi}^{\dagger}(t_{1})\rangle_{\rm
noise}=2\gamma\delta(t-t_{1}),$
$\displaystyle\langle\hat{\xi}(t)\hat{\xi}(t_{1})\rangle_{\rm
noise}=\langle\hat{\xi}^{\dagger}(t)\hat{\xi}(t_{1})\rangle_{\rm
noise}=\langle\hat{\xi}^{\dagger}(t)\hat{\xi}^{\dagger}(t_{1})\rangle_{\rm
noise}=0,$ (36)

where $\langle\cdots\rangle_{\rm noise}$ denotes the noise average foot1 .
This formalism is equivalent to the total Hamiltonian with a white noise bath.
In the interaction picture, the time evolution of density matrix can be
expressed by

$\displaystyle\hat{\rho}(t)$ $\displaystyle=$ $\displaystyle\hat{U}_{\rm
eff}(t)\hat{\rho}_{0}(t)\hat{U}^{\dagger}_{\rm eff}(t)$

where $\hat{U}_{\rm
eff}(t)={\rm\widetilde{T}}\exp{\left(-i\int_{0}^{t}\hat{H}_{\rm
diss}(t^{\prime})dt^{\prime}\right)}$ with anti-time-ordered operator
${\rm\widetilde{T}}$ and
$\hat{\rho}_{0}(t)=e^{-i\hat{H}_{S}t}\hat{\rho}(0)e^{i\hat{H}_{S}t}$ denotes
the time-evolution of density matrix determined by $\hat{H}_{S}$. Taking
$\hat{H}_{\mathrm{diss}}$ as perturbation and then averaging the noise, one
can obtain the density matrix with the first-order correction of $\gamma$

$\displaystyle\rho(t)$ $\displaystyle\equiv$
$\displaystyle\left\langle\hat{\rho}(t)\right\rangle_{\rm noise}=\left\langle
U^{\prime}_{\rm eff}(t)\hat{\rho}_{0}(t)U^{\prime{\dagger}}_{\rm
eff}(t)\right\rangle_{\rm noise}$ (38) $\displaystyle=$
$\displaystyle\Bigg{<}\left(1+\sum_{n=1}^{\infty}(-i)^{n}\int_{t_{1}<\cdots<t_{n}}\hat{H}_{\rm
diss}(t_{1})\cdots\hat{H}_{\rm diss}(t_{n})dt_{1}\cdots
dt_{n}\right)\hat{\rho}_{0}(t)$ $\displaystyle\times$
$\displaystyle\left(1+\sum_{n=1}^{\infty}(i)^{n}\int_{t_{1}<\cdots<t_{n}}\hat{H}^{\dagger}_{\rm
diss}(t_{n})\cdots\hat{H}^{\dagger}_{\rm diss}(t_{1})dt_{1}\cdots
dt_{n}\right)\Bigg{>}_{\rm noise}$ $\displaystyle\approx$
$\displaystyle\rho_{0}(t)-\int_{0}^{t}dt^{\prime}\gamma\Big{\\{}\hat{\cal
O}^{\dagger}(t^{\prime})\hat{\cal
O}(t^{\prime}),\rho_{0}(t)\Big{\\}}+2\int_{0}^{t}dt^{\prime}\gamma\hat{\cal{O}}(t^{\prime})\rho_{0}(t)\hat{\cal
O}^{\dagger}(t^{\prime}),$

where $\langle\cdots\rangle_{\rm noise}$ is the noise average and the
correlation function Eq.(A) has been applied. The linear response of density
matrix is defined by

$\displaystyle\delta\rho(t)\equiv\rho(t)-\rho_{0}(t)$ $\displaystyle=$
$\displaystyle-\int_{0}^{t}dt^{\prime}\gamma\Big{\\{}\hat{\cal
O}^{\dagger}(t^{\prime})\hat{\cal O}(t^{\prime}),\rho_{0}(t)\Big{\\}}$ (39)
$\displaystyle+$ $\displaystyle
2\int_{0}^{t}dt^{\prime}\gamma\hat{\cal{O}}(t^{\prime})\rho_{0}(t)\hat{\cal
O}^{\dagger}(t^{\prime}).$

From the response of density matrix, one can calculate the von Neumann entropy
which characterizes the loss of coherence. Here we focus on the density matrix
in Krammers degenerate subspace $\rho_{\rm K}(t)=\Pi_{\rm K}\rho(t)\Pi_{\rm
K}$ ($\rho_{0,\rm K}(t)=\Pi_{\rm K}\rho_{0}(t)\Pi_{\rm K}$) and the
corresponding response of von Neumann entropy is given by

$\displaystyle\delta S_{\rm v}(t)$ $\displaystyle=$ $\displaystyle S_{\rm
v}(t)-S_{0,\rm R}(t)$ (40) $\displaystyle=$
$\displaystyle-\operatorname{Tr}\left[\frac{\rho_{\rm 0,K}(t)+\delta\rho_{\rm
K}(t)}{\operatorname{Tr}\left(\rho_{\rm 0,K}(t)+\delta\rho_{\rm
K}(t)\right)}\log\left(\frac{\rho_{\rm 0,K}(t)+\delta\rho_{\rm
K}(t)}{\operatorname{Tr}\left(\rho_{\rm 0,K}(t)+\delta\rho_{\rm
K}(t)\right)}\right)\right]+\operatorname{Tr}\left[\rho_{\rm
0,K}(t)\log\rho_{\rm 0,K}(t)\right]$ $\displaystyle=$
$\displaystyle\operatorname{Tr}\Big{[}\rho_{\rm 0,K}(t)\log\rho_{\rm
0,K}(t)-\left(\delta\rho_{\rm 0,K}(t)/\operatorname{Tr}\delta\rho_{\rm
0,K}(t)\right)\log\rho_{\rm 0,K}(t)\Big{]}\operatorname{Tr}\delta\rho_{\rm
0,K}(t).$

We initialize the density matrix in the Kramers’ degenerate space, i.e. ${\rm
Tr}\rho_{0,K}(0)=1$ which means ${\rm Tr}\rho_{0,K}(t)=1$. It is clearly that,
when $\delta\rho_{\rm K}(t)\propto\rho_{\rm 0,K}(t)$ (namely, $\delta\rho_{\rm
K}(t)=\rho_{\rm 0,K}(t)\operatorname{tr}\delta\rho_{\rm K}(t))$ we have
$\delta S_{\rm v}(t)=0$, $\delta S_{\rm R}(t)=0$ but $\delta S_{\rm v}(t)\neq
0$ if $\delta\rho_{\rm K}(t)\neq\rho_{\rm 0,K}(t)$. In other words, the time-
reversal symmetry breaking discussed in main text leads to the growth of
entropy.

##  Appendix B Matrix representation of quaternion group

This appendix provides matrix representation of quaternion group $\mathcal{Q}$
discussed in maintext. The $\mathcal{Q}$-group $\\{Q_{j},j=1,\cdots,8\\}$ is a
non-Abelian group which is isomorphic to subset $\\{1,i,j,k,-1,-i,-j,-k\\}$
whose multiplication table is displayed in Table 2.

Table 2: The multiplication table (Cayley table) of quaternion group. Element | $\mathrm{e}$ | $\overline{\mathrm{e}}$ | $i$ | $\overline{\mathrm{i}}$ | $j$ | $\overline{\mathrm{j}}$ | $k$ | $\overline{\mathrm{k}}$
---|---|---|---|---|---|---|---|---
$\mathrm{e}$ | $\mathrm{e}$ | $\overline{\mathrm{e}}$ | $i$ | $\overline{\mathrm{i}}$ | $j$ | $\overline{\mathrm{j}}$ | $k$ | $\overline{\mathrm{k}}$
-1 | -1 | $\mathrm{e}$ | $\overline{\mathrm{i}}$ | $i$ | $\overline{\mathrm{j}}$ | $j$ | $\overline{\mathrm{k}}$ | $k$
$i$ | $i$ | $\overline{\mathrm{i}}$ | $\overline{\mathrm{e}}$ | $\mathrm{e}$ | $k$ | $\overline{\mathrm{k}}$ | $\overline{\mathrm{j}}$ | $j$
$\overline{\mathrm{i}}$ | $\overline{\mathrm{i}}$ | $i$ | $\mathrm{e}$ | $\overline{\mathrm{e}}$ | $\overline{\mathrm{k}}$ | $k$ | $j$ | $\overline{\mathrm{j}}$
$j$ | $j$ | $\overline{\mathrm{j}}$ | $\overline{\mathrm{k}}$ | $k$ | $\overline{\mathrm{e}}$ | $\mathrm{e}$ | $i$ | $\overline{\mathrm{i}}$
$\overline{\mathrm{j}}$ | $\overline{\mathrm{j}}$ | $j$ | $k$ | $\overline{\mathrm{k}}$ | $\mathrm{e}$ | $\overline{\mathrm{e}}$ | $\overline{\mathrm{i}}$ | $i$
$k$ | $k$ | $\overline{\mathrm{k}}$ | $j$ | $\overline{\mathrm{j}}$ | $\overline{\mathrm{i}}$ | $i$ | $\overline{\mathrm{e}}$ | $\mathrm{e}$
$\overline{\mathrm{k}}$ | $\overline{\mathrm{k}}$ | $k$ | $\overline{\mathrm{j}}$ | $j$ | $i$ | $\overline{\mathrm{i}}$ | $\mathrm{e}$ | $\overline{\mathrm{e}}$

The $\mathcal{Q}$-group contains a two-dimensional irreducible representation
(see Table.3) which can be described as a subgroup of the special linear group
$\mathrm{SL}_{2}(\mathbb{C})$. We can construct the following 4-dimensional
reducible representation

$\displaystyle Q_{1}$ $\displaystyle=I\otimes
I,~{}Q_{2}=-Q_{1},~{}Q_{3}=-I\otimes i\sigma_{z},~{}Q_{4}=-Q_{3},$
$\displaystyle Q_{5}$
$\displaystyle=-i\sigma_{x}\otimes\sigma_{y},~{}Q_{6}=-Q_{5},~{}Q_{7}=-i\sigma_{x}\otimes\sigma_{x},~{}Q_{8}=-Q_{7},$
(41)

where $\sigma_{x,y,z}$ denote Pauli matrices and $I$ is identity matrix. One
can easily find that the matrix representation (41) obeys the multiplication
table (2) and contains two-dimensional irreducible representation.

Table 3: The character table of quaternion group Representation/Conjugacy class | $\\{\mathrm{e}\\}$ | $\\{\overline{\mathrm{e}}\\}$ | $\\{\mathrm{i},\overline{\mathrm{i}}\\}$ | $\\{\mathrm{j},\overline{\mathrm{j}}\\}$ | $\\{\mathrm{k},\overline{\mathrm{k}}\\}$
---|---|---|---|---|---
Trivial representation | 1 | 1 | 1 | 1 | 1
i-kernel | 1 | 1 | 1 | -1 | -1
j-kernel | 1 | 1 | -1 | 1 | -1
k-kernel | 1 | 1 | -1 | -1 | 1
2-dimensional representation | 2 | -2 | 0 | 0 | 0

## References

  * (1) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008).
  * (2) A. Kitaev, AIP Conf. Proc. 1134, 22 (2009).
  * (3) S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, New J. Phys. 12, 065010 (2010).
  * (4) M. R. Zirnbauer, J. Math. Phys. (N.Y.) 37, 4986 (1996).
  * (5) A. Altland and M. R. Zirnbauer, Phys. Rev. B 55, 1142 (1997).
  * (6) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (7) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (8) B. Bernevig and T. Hughes, Topological Insulators and Topological Superconductors (Princeton University Press,2013).
  * (9) C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * (10) Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Phys. Rev. X 8, 031079 (2018).
  * (11) K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Phys. Rev. X 9, 041015 (2019).
  * (12) C.-H. Liu and S. Chen, Phys. Rev. B 100, 144106 (2019).
  * (13) H. Zhou and J. Y. Lee, Phys. Rev. B 99, 235112 (2019).
  * (14) K. Kawabata, T. Bessho, and M. Sato, Phys. Rev. Lett. 123, 066405(2019).
  * (15) K. Kawabata, S. Higashikawa, Z. Gong, Y. Ashida and M. Ueda, Nat Commun 10, 297 (2019).
  * (16) C.-H. Liu, H. Jiang, and S. Chen, Phys. Rev. B 99, 125103 (2019).
  * (17) J. Y. Lee, J. Ahn, H. Zhou, and A. Vishwanath,Phys. Rev. Lett. 123, 206404 (2019).
  * (18) S. Lieu, M. McGinley, and N. R. Cooper, Phys. Rev. Lett. 124, 040401 (2020).
  * (19) A. Altland, M. Fleischhauer, S. Diehl, arXiv:2007.10448 (2020).
  * (20) M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, 2010).
  * (21) E. Wigner, Gruppentheorie (Vieweg, Braunschweig, 1931).
  * (22) M. McGinley and N. R. Cooper, Nat. Phys. 16, 1181 (2020).
  * (23) L. Pan, X. Chen, Yu, Chen and H. Zhai, Nat. Phys. 16, 767 (2020).
  * (24) Tian-Shu Deng, Lei Pan, Yu Chen, Hui Zhai, arXiv:2009.13043 (2020).
  * (25) J. O. Dimmock, J. Math. Phys. 4, 1307 (1963).
  * (26) Y. Jian and Z.-X. Liu, J. Phys. A 51, 025207 (2017).
  * (27) It is worth emphasizing that the Langevin noise is essential since it preserves canonical commutation relation. Hamiltonian without Langevin noise contradicts the principle of quantum mechanics.
  * (28) M.-D. Choi, Lin. Alg. Appl. 10, 285 (1975).
  * (29) A. Jamiołkowski Rep. Math. Phys. 3, 275 (1972).
  * (30) E. P. Wigner, Group Theory and Its Application to the Quantum Mechanics of Atomic Spectra (Academic Press, New York, 1959).
  * (31) For the bath beyond Born-Markov approximation at finite temperature, the anti-unitary protected coherence in system could also be unstable against to hermitian coupling for the multi-channel case. This can be illustrated by a two-channel model with two independent coupling operators which give rise to a non-hermitian operator $c_{1}\mathcal{\hat{O}}_{1}+c_{2}\mathcal{\hat{O}}_{2}$ (with real numbers $c_{1},c_{2}$) acting on the system. As discussed in section IV.2, this non-Hermiticity leads to decoherence despite the possession of anti-symmetry.
  * (32) J. Choi, S. Choi, G. Kucsko, P. C. Maurer, B. J. Shields, H. Sumiya, S. Onoda, J. Isoya, E. Demler, F. Jelezko, N. Y. Yao, and M. D. Lukin, Phys. Rev. Lett. 118, 093601 (2017).
  * (33) G. Kucsko, S. Choi, J. Choi, P. C. Maurer, H. Zhou, R. Landig, H. Sumiya, S. Onoda, J. Isoya, F. Jelezko, E. Demler, N. Y. Yao, and M. D. Lukin, Phys. Rev. Lett. 121, 023601 (2018).
  * (34) T. Tomita, S. Nakajima, I. Danshita, Y. Takasu, and Y. Takahashi, Sci. Adv. 3, e1701513 (2017).
  * (35) T. Tomita, S. Nakajima, Y. Takasu, and Y. Takahashi, Phys.Rev. A 99, 031601(R) (2019).
  * (36) J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar, and L. Luo, Nat. Commun.10, 855 (2019).
  * (37) K. Sponselee, L. Freystatzky, B. Abeln, M. Diem, B. Hundt, A. Kochanke, T. Ponath, B. Santra, L. Mathey, K. Sengstock, and C. Becker, Quantum Science and Technology 4, 014002 (2019).
  * (38) R. Bouganne, M. B. Aguilera, A. Ghermaoui, J. Beugnon, and F. Gerbier, Nat. Phys. 16, 21 (2020).
  * (39) F. Schäfer, T. Fukuhara, S. Sugawa, Y. Takasu, and Y. Takahashi, Nat. Rev. Phys. 2, 411 (2020).
  * (40) Yosuke Takasu, Tomoya Yagami, Yuto Ashida, Ryusuke Hamazaki, Yoshihito Kuno, Yoshiro Takahashi, arXiv:2004. 05734 (2020).
  * (41) Z. Wang, Q. Li, W. Li, Z. Cai, arXiv:2102.06383 (2021).


# A first-principles investigation of the linear thermal expansion
coefficients of BeF2: Giant thermal expansion

Chee Kwan Gan<EMAIL_ADDRESS>Institute of High Performance
Computing, 1 Fusionopolis Way, #16-16 Connexis 138632, Singapore  Abdullah I.
Al-Sharif, Ammar Al-Shorman  Department of Physics, Yarmouk University,
Irbid-21163, Jordan  Abdallah Qteish<EMAIL_ADDRESS>Department of Physics,
Yarmouk University, Irbid-21163, Jordan

###### Abstract

We present the results of a theoretical investigation of the linear thermal
expansion coefficients (TECs) of BeF2, within a direct Grüneisen formalism
where symmetry-preserving deformations are employed. The required physical
quantities such as the optimized crystal structures, elastic constants, mode
Grüneisen parameters, and phonon density of states are calculated from first-
principles. BeF2 shows an extensive polymorphism at low pressures, and the
lowest energy phases [$\alpha$-cristobalite with space group (SG)
$P4_{1}2_{1}2$ and its similar phase with SG $P4_{3}2_{1}2$] are considered in
addition to the experimentally observed $\alpha$-quartz phase. For
benchmarking purposes, similar calculations are performed for the rutile phase
of ZnF2, where the volumetric TEC ($\alpha_{v}$), derived from the calculated
linear TECs along the $a$ ($\alpha_{a}$) and $c$ ($\alpha_{c}$) directions, is
in very good agreement with experimental data and previous theoretical
results. For the considered phases of BeF2, we do not find any negative
thermal expansion (NTE). However we observe diverse thermal properties for the
distinct phases. The linear TECs are very large, especially $\alpha_{c}$ of
the $\alpha$-cristobalite phase and its similar phase, leading to giant
$\alpha_{v}$ ($\sim 175\times 10^{-6}$ K-1 at 300 K). The giant $\alpha_{v}$
arises from large Grüneisen parameters of low-frequency phonon modes, and the
$C_{13}$ elastic constant that is negatively signed and large in magnitude for
the $\alpha$-cristobalite phase. The elastic constants, high-frequency
dielectric constants, Born effective charge tensors, and thermal properties of
the above phases of BeF2 are reported for the first time and hence serve as
predictions.

phonon calculations, lattice dynamical studies, small-displacement methods,
crystal symmetry, negative thermal expansion, Grüneisen parameters

##  I Introduction

Beryllium fluoride (BeF2) is known to exist in glass and crystalline phases
and has a variety of technological applications. BeF2 glass has a large
bandgap of about $13.8$ eV, the lowest refractive index and highest Abbe
number of any inorganic material, and exceptional resistance to damage. These
properties have enabled the manufacturing of special glasses (from BeF2 and
its mixtures with fluorides and other difluorides) that have excellent
transmittance in the UV region[1, 2] and for use in high-power laser
systems[3]. The LiF–BeF2 mixture is a primary coolant and fuel solvent in
molten salt nuclear reactors[4]. In protein crystallography, BeF2 is used to
restrict protein motion to facilitate the crystallography process[5]. Very
recently, crystalline BeF2 is predicted to be a better neutron filter than
MgF2, which has been considered an effective neutron filter candidate[6]. The
main aim of this work is to investigate, for the first time, the linear
thermal expansion coefficients (TECs) of a few low-energy crystalline phases
of BeF2. We also consider a benchmark system ZnF2 that has exceptional
electric and optical properties, and interesting technological applications
ranging from catalysis to spectroscopy and laser applications[7].

Single crystal BeF2 has been grown and found to have a crystal structure
remarkably similar to that of the $\alpha$-quartz (SiO2) structure[8], which
has a trigonal symmetry with space group (SG) $P3_{1}21$ (#152). A recent
first-principles study has revealed that BeF2 shows extensive polymorphism at
low pressures[9]. Interestingly, three crystal phases [namely, (i) the
$\alpha$-cristobalite phase that has a tetragonal symmetry with SG
$P4_{1}2_{1}2$ (#92), (ii) a similar phase to the $\alpha$-cristobalite phase
(hereafter referred to as the $\alpha^{\prime}$-cristobalite phase) with SG
$P4_{3}2_{1}2$ (#96) , and (iii) the $C2/c$-$2\times 4$ phase with SG $C12/c1$
(#15)] are predicted to be energetically more stable than $\alpha$-quartz.
However, these phases have a very small stability pressure range (less than
$0.4$ GPa), and the $\alpha$-quartz phase transforms to the coesite-I phase SG
C2/c at $3.1$ GPa. The high-pressure phases of BeF2 have been the subject of
other first-principles calculations[10]. Very recently, first-principles
calculations have also been employed to construct the $P$-$T$ phase diagram of
BeF2[11]. The HSE06 optical bandgap of the $\alpha$-quartz structure is found
to be about $10.6$ eV, and increases by increasing the applied pressure[9].
The lattice vibrations, inelastic scattering cross-sections, and neutron
transmission of BeF2 have been thoroughly investigated[6] using first-
principles calculations and compared to those of MgF2.

The benchmark system ZnF2 crystallizes in the tetragonal rutile structure with
SG $P4_{2}/mnm$ (#136). Very recently, Raman scattering measurements with the
use of the diamond anvil cell have been employed to investigate the structural
phase transformations of ZnF2 under high pressures[12]. This experimental work
is supplemented by first-principles calculations. In addition to the
structural stability and pressure variation of the Raman active phonon modes,
the electronic bandgap of the considered phases as a function of pressure has
been investigated at the HSE06 level. Neutron diffraction has been employed to
study the temperature dependence of the lattice parameters and unit cell
volume of ZnF2[13], and NTE has been observed in a small temperature range
(below 75 K). This NTE behavior has been supported by first-principles
calculations[13, 14]. However, only the volumetric TEC has been theoretically
investigated.

In the present work, the linear TECs of BeF2 and ZnF2 are investigated by
employing the recently introduced direct approach[15] in which the symmetry of
the deformed structures could be preserved. Since this approach has not been
applied to systems with tetragonal symmetry, ZnF2 is thus chosen as a suitable
benchmark system because of its tetragonal crystal structure, in addition to
the existence of experimental and previous theoretical results of its
volumetric thermal expansion. The elastic constants and phonon frequencies
required to compute linear TECs are calculated from first-principles. For
BeF2, the $\alpha$-quartz, $\alpha$-cristobalite and
$\alpha^{\prime}$-cristobalite phases are considered. Moreover, the relative
stability of the above three phases of BeF2 will be also investigated using
different levels of approximation of the exchange-correlation potential.

##  II Methodology

The linear TECs of the considered phases of ZnF2 and BeF2 are calculated
within the Grüneisen formalism following the procedure described in Refs. [16,
17, 18, 19, 20, 21, 22]. To compute the mode Grüneisen parameters we
considered two types of symmetry-preserving deformations obtained by changing
the in-plane ($a$) or out-of-plane ($c$) lattice parameters by $\pm 0.5\%$.
These deformations allow for the full utilization of the tetragonal or
trigonal point-group symmetry[23] of the considered systems, which minimizes
the required number of independent atomic displacements (i.e., number of
supercells) to calculate the phonon frequencies within the direct method[24,
25, 26, 27, 28]. The amplitude of atomic displacements, from the corresponding
equilibrium positions, is $0.015$ Å. The supercell sizes are of $2\times
2\times 3$ for ZnF2, and $2\times 2\times 2$ for the $\alpha$-quartz,
$\alpha$-cristobalite, and $\alpha^{\prime}$-cristobalite phases of BeF2. The
adequacy of these supercells have been checked by considering larger ones for
each of these systems, and we found that such actions do not alter appreciably
our main results and conclusions. The determination of the linear TECs also
requires the elastic constants that may be obtained through fittings of energy
versus strain curves[29, 30]. Specifically, these TECs at temperature $T$ in
the $a$ ($\alpha_{a}$) and $c$ ($\alpha_{c}$) directions of the above systems
are given by

$\left[\begin{matrix}\alpha_{a}(T)\\\
\alpha_{c}(T)\end{matrix}\right]=\frac{1}{\Omega
D}\left[\begin{matrix}C_{33}&-C_{13}\\\
-2C_{13}&(C_{11}+C_{12})\end{matrix}\right]\left[\begin{matrix}I_{a}(T)\\\
I_{c}(T)\end{matrix}\right]$ (1)

where $C_{ij}$ are the elastic constants,
$D=(C_{11}+C_{12})C_{33}-2C_{13}^{2}$ and $\Omega$ is the volume of the
primitive cell. The phonon density of states (PDOS) weighted by the Grüneisen
parameters are

$\Gamma_{i}(\nu)=\frac{\Omega}{(2\pi)^{3}}\sum_{\lambda}\int_{\rm
BZ}d{\boldsymbol{q}}\
\delta(\nu-\nu_{\lambda{\boldsymbol{q}}})\gamma_{i,\lambda{\boldsymbol{q}}},$
(2)

with the Grüneisen parameter
$\gamma_{i,\lambda{\boldsymbol{q}}}=-\nu_{\lambda{\boldsymbol{q}}}^{-1}\partial\nu_{\lambda{\boldsymbol{q}}}/\partial\epsilon_{i}$
for the deformation of type $i$ (with $i=a$ for in-plane and $i=c$ for out-of-
plane deformations). The derivative
$\partial\nu_{\lambda{\boldsymbol{q}}}/\partial\epsilon_{i}$ measures the
change of the frequency $\nu_{\lambda{\boldsymbol{q}}}$ with respect to the
strain parameter $\epsilon_{i}$[15]. The summation is over all frequencies
$\nu_{\lambda{\boldsymbol{q}}}$ for the phonon band index $\lambda$ and
${\boldsymbol{q}}$ vector in the Brillouin zone (BZ). The $I_{i}(T)$ are
calculated from

$I_{i}(T)=\int_{\nu_{\rm min}}^{\nu_{\rm max}}d\nu\ \Gamma_{i}(\nu)c(\nu,T),$
(3)

where $c(\nu,T)=k_{B}(r/\sinh r)^{2}$ is the contribution of the phonon modes
with frequency $\nu$ to the specific heat. Here, $r=h\nu/2k_{B}T$, and $h$ and
$k_{B}$ are the Planck and Boltzmann constants, respectively. The maximum
(minimum) frequency is denoted by $\nu_{\rm max}$ ($\nu_{\rm min}$).

The DFT calculations of the optimized structural parameters, phonon
frequencies, and elastic constants are performed by employing the projector
augmented wave (PAW) method, as implemented in the Vienna Ab-Initio Simulation
Package (VASP). A relatively high cutoff energy of $600$ eV is used for the
plane-wave basis. Geometry optimization is stopped when the maximum force on
each atom is less than $10^{-3}$ eV/Å. We find that phonon frequency shifts
are more consistent when we use the local density approximation (LDA) for
BeF2, and the PBE$\\_$sol functional of the generalized gradient approximation
(GGA) for ZnF2. Therefore, for the linear TECs only the results of these
calculations are reported.

Table 1:  Calculated lattice constants and internal parameters of the $\alpha$-quartz and $\alpha$-cristobalite phases of BeF2. Also shown are the relative energy ($\Delta E$) of the $\alpha$-cristobalite phase with respect to that of $\alpha$-quartz, and the available experimental data (measured at 100 K) and other theoretical results. Phase | Approach | Lattice constants |  | Internal parameters |  $\Delta E$ (meV)
---|---|---|---|---|---
|  |  $a$ (Å) |  $c$ (Å) |  | $x_{1}$ | $x_{2}$ | $y_{2}$ | $z_{2}$ |
$\alpha$-quartz | LDA | 4.5958 | 5.0529 |  | 0.4579 | 0.4098 | 0.2867 | 0.2290 |
| PBE$\\_$sol | 4.7301 | 5.1814 |  | 0.4662 | 0.4138 | 0.2737 | 0.2188 |
| PBE | 4.8497 | 5.3070 |  | 0.4756 | 0.4176 | 0.2575 | 0.2053 |
| PBE[6] | 4.8282 | 5.2837 |  | 0.4740 | 0.4171 | 0.2601 | 0.2075 |
| LDA[31] | 4.6663 | 5.1608 |  |  |  |  |  |
| Expt.[8] | 4.7390 | 5.1875 |  | 0.4700 | 0.4164 | 0.2671 | 0.2131 |
$\alpha$-cristobalite | LDA | 4.5967 | 6.1773 |  | 0.3226 | 0.2230 | 0.1454 | 0.2000 | 25
| PBE$\\_$sol | 4.8087 | 6.5984 |  | 0.3044 | 0.2378 | 0.1112 | 0.1825 | $-2$
| PBE | 4.8934 | 6.7428 |  | 0.2988 | 0.2400 | 0.1001 | 0.1769 | $-9$
| LDA[31] | 4.695 | 6.318 |  |  |  |  |  |
| LDA[11] | 4.684 | 6.373 |  |  |  |  |  |
| PBE[11] | 4.960 | 6.910 |  |  |  |  |  |
Figure 1: The crystal structures of (a) rutile ZnF2, (b) $\alpha$-quartz BeF2,
and (c) $\alpha$-cristobalite BeF2. The $\alpha^{\prime}$-cristobalite BeF2
structure which is similar to the $\alpha$-cristobalite BeF2 structure is not
shown. The small gray balls represent the F atoms. The $c$ axis of the three
crystal structures is along the vertical direction.

##  III Results and discussion

###  III.1 Structural properties

We start with the benchmark system, ZnF2, which at ambient conditions
crystallizes in the rutile structure, shown in Fig. 1(a). This crystal
structure has a tetragonal symmetry and six atoms per primitive unit cell. The
two Zn atoms occupy the Wyckoff $2a(0,0,0)$ sites and the four F atoms are
located at the $4f(x,x,0)$ sites. Therefore, this structure is characterized
by three crystallographic parameters: two lattice parameters ($a$ and $c$) and
an internal parameter for the coordinates of the four F atoms ($x$). The
structural parameters obtained using the PBE$\\_$sol functional are
$(a,c,x)=(4.7194~{}{\rm\AA},3.1376~{}{\rm\AA},0.3037)$, while the
corresponding LDA results are (4.6373 Å, 3.0990 Å, 0.3033). These calculated
values are in very good agreement with the experimental data (4.7038 Å, 3.1336
Å, 0.3035)[32] and (4.7034 Å, 3.1335 Å, 0.303)[13]. As expected, the LDA
lattice parameters are underestimated while those of the PBE$\\_$sol are
slightly overestimated.

For BeF2, the three considered crystal structures are $\alpha$-quartz
[Fig.1(b)], $\alpha$-cristobalite [Fig.1(c)] and
$\alpha^{\prime}$-cristobalite. The $\alpha$-quartz phase has trigonal
symmetry and nine atoms in the conventional hexagonal unit cell. The three Be
atoms occupy the Wyckoff $3a(x_{1},x_{1},0)$ sites and the six F atoms occupy
the $6c(x_{2},y_{2},z_{2})$ sites. On the other hand, the
$\alpha$-cristobalite phase has a tetragonal symmetry and twelve atoms per
primitive unit cell. The four Be atoms occupy the Wyckoff $4a(x_{1},x_{1},0)$
sites and the eight F atoms occupy the $8c(x_{2},y_{2},z_{2})$ sites.
Therefore, each of these two structures has six crystallographic parameters:
two lattice parameters ($a$ and $c$), and four internal parameters (denoted as
$x_{1}$, $x_{2}$, $y_{2}$ and $z_{2}$). The atomic coordinates of the
$\alpha^{\prime}$-cristobalite phase can be obtained from those of the
$\alpha$-cristobalite by mirror-image transformation
$(x,y,z)\rightarrow(-y,-x,z)$, and the lattice parameters of the two
structures are identical. Therefore, only the structure parameters of the
first two crystal structures are reported. Our LDA, PBE$\\_$sol, and PBE
results, shown in Table 1, are in good agreement with available experimental
data and other theoretical calculations. The PBE$\\_$sol results lie between
the corresponding LDA and PBE results and show the best agreement with the
experimental data for the $\alpha$-quartz.

Table 2:  Calculated elastic constants of the $\alpha$-quartz and $\alpha$-cristobalite phases of BeF2, and the rutile phase of ZnF2.  System | Phase | Approach | Elastic constants (GPa)
---|---|---|---
|  |  | $C_{11}$ | $C_{12}$ | $C_{13}$ | $C_{14}$ | $C_{33}$ | $C_{44}$ | $C_{66}$
BeF2 |  $\alpha$-quartz | LDA | 46.975 | 14.223 | 12.067 | $-6.401$ | 75.287 | 31.745 |
|  | PBE$\\_$sol | 42.278 | 4.991 | 3.760 | $-9.077$ | 53.009 | 30.548 |
|  $\alpha$-cristobalite | LDA | 33.309 | 7.300 | $-5.087$ |  | 22.487 | 35.810 | 13.731
|  | PBE$\\_$sol | 32.589 | 4.839 | $-5.336$ |  | 24.412 | 37.813 | 16.078
ZnF2 | Rutile | LDA | 139.442 | 121.550 | 109.127 |  | 220.673 | 36.583 | 91.826
|  | PBE$\\_$sol | 128.470 | 98.717 | 94.547 |  | 200.902 | 35.523 | 83.001
|  | Expt.[33] | 125.5 | 91.8 | 83.0 |  | 192.2 | 39.5 | 80.7

###  III.2 Elastic properties and stability

We show in Table 1 the LDA, PBE$\\_$sol and PBE relative energies ($\Delta E$)
of the $\alpha$-cristobalite phase with respect to those of $\alpha$-quartz.
According to the PBE and PBE$\\_$sol calculations the latter phase is slightly
more stable, in accordance with the Nelson et al.[9] GGA calculations.
However, the LDA calculations lead to an opposite conclusion. Similar
conclusions have recently been reached by Masoumi[11], using both LDA and GGA
calculations. These results show that these two phases have extremely close
cohesive energies.

The $\alpha$-quartz phase of BeF2 with a trigonal crystal symmetry has six
independent elastic constants[34]. On the other hand, the rutile phase of ZnF2
and $\alpha$-cristobalite phase of BeF2 (both have a tetragonal crystal
symmetry) also have six independent elastic constants[34]. The elastic
constants of these phases, obtained by using the LDA and PBE$\\_$sol
functionals are shown in Table 2. There are two features to note from this
table. First, our results for the rutile ZnF2 are in very good agreement with
the available experimental values[33]. Secondly, the PBE$\\_$sol values are
systematically smaller than the corresponding LDA values, which is expected
since the PBE$\\_$sol GGA functional leads to softer materials than LDA (see
above). The calculated elastic constants are used in the calculations of the
linear TECs (see Sec. II).

The elastic constants could be used to investigate the mechanical stability of
the crystal structure. For $\alpha$-quartz structure, the Born stability
criteria[34] are

$D=(C_{11}+C_{12})C_{33}-2C_{13}^{2}>0,\\\ $ (4)

and

$(C_{11}-C_{12})C_{44}-2C_{14}^{2}>0.$ (5)

Note that the same expression of $D$ appears in Eq. 1. As for the rutile and
$\alpha$-cristobalite phases (of tetragonal (I) class)[35], the Born stability
criteria are: $C_{11}>|C_{12}|$, $D>0$, $C_{44}>0$ and $C_{66}>0$. The elastic
constants reported in Table 2 show that these criteria are satisfied, and
therefore the considered phases of ZnF2 and BeF2 are mechanically stable. The
mechanical stability of these crystal structures can also be inferred from
phonon dispersion relations, discussed below.

Table 3:  Calculated high-frequency dielectric constants (DCs) and Born effective charges of the $\alpha$-quartz (AQ) and $\alpha$-cristobalite (AC) phases of BeF2, and the rutile phase of ZnF2. System | Phase | Approach | DC |  | Born effective charge
---|---|---|---|---|---
|  |  | xx=yy | zz | atom | xx | xy | xz | yx | yy | yz | zx | zy | zz
BeF2 | AQ | PBE$\\_$sol | 1.902 | 1.912 | Be | 1.728 | 0.000 | 0.000 | 0.000 | 1.914 | 0.081 | 0.000 | -0.079 | 1.866
|  |  |  |  | F | $-0.747$ | 0.221 | $-0.115$ | 0.213 | $-1.074$ | 0.337 | $-0.095$ | 0.345 | $-0.933$
|  | LDA | 1.887 | 1.896 | Be | 1.727 | 0.000 | 0.000 | 0.000 | 1.918 | 0.080 | 0.000 | $-0.0758$ | 1.864
|  |  |  |  | F | $-0.746$ | 0.226 | $-0.124$ | 0.218 | $-1.076$ | 0.343 | $-0.104$ | 0.356 | $-0.932$
BeF2 | AC | PBE$\\_$sol | 1.723 | 1.718 | Be | 1.832 | 0.005 | $-0.045$ | 0.005 | 1.832 | 0.049 | 0.100 | $-0.101$ | 1.810
|  |  |  |  | F | $-1.221$ | $-0.117$ | 0.377 | $-0.101$ | $-0.611$ | 0.061 | 0.380 | 0.100 | $-0.905$
|  | LDA | 1.827 | 1.818 | Be | 1.823 | 0.005 | $-0.036$ | 0.005 | 1.823 | 0.036 | 0.110 | $-0.110$ | 1.790
|  |  |  |  | F | $-1.199$ | $-0.135$ | 0.329 | $-0.120$ | $-0.624$ | 0.076 | 0.333 | 0.115 | $-0.896$
ZnF2 | Rutile | PBE$\\_$sol | 2.549 | 2.664 | Zn | 2.222 | $-0.162$ | 0.000 | $-0.162$ | 2.222 | 0.00 | 0.000 | 0.000 | 2.424
|  |  |  |  | F | $-1.111$ | $-0.409$ | 0.000 | -0.409 | $-1.111$ | 0.000 | 0.000 | 0.000 | $-1.200$
|  | LDA | 2.547 | 2.653 | Zn | 2.206 | $-0.1493$ | 0.000 | $-0.1493$ | 2.206 | 0.000 | 0.000 | 0.000 | 2.392
|  |  |  |  | F | $-1.103$ | $-0.395$ | 0.000 | $-0.395$ | -1.103 | 0.000 | 0.000 | 0.000 |  $-1.19$6
|  | Expt.[36] | 2.6 | 2.1 |  |  |  |  |  |  |  |  |  |
Figure 2: Calculated phonon dispersion relations and PDOS of ZnF2, with (black
solid curves) and without (red dashed curves) non-analytic correction (NAC).
The Zn and F projected PDOS, with NAC, and also shown. Symbols: available
experimental data[37, 38].

###  III.3 Phonon dispersion relations

Since the considered crystals are polar in character we perform non-analytical
correction (NAC) to their dynamical matrices. To do that, we have calculated
the high-frequency dielectric constant and Born effective charge tensors, and
the results are listed in Table 3. The features to note from this table are
the following. (i) The calculated results have a very weak dependence on the
used exchange-correlation functional. (ii) The effective charges in the ZnF2
are larger than in BeF2, which shows that the ionicity of Zn-F bond is larger
than that in Be-F. This is consistent with the corrected Allred-Rochow
electronegativity values[39] (larger for Be). (iii) Our calculated $xx$ and
$yy$ components of the dielectric constant are in very good agreement with
available experimental data[36], while that of $zz$ is larger than the
measured one. However, the comparable values of diagonal components of the
dielectric constant in our calculations are consistent with the experimental
and calculated values for other metal fluorides crystallizing in the rutile
structure (such as MgF2 and FeF2)[40].

Fig. 2 shows the calculated phonon dispersion relations and PDOS of the rutile
ZnF2, with and without NAC. Also shown are the calculated Zn and F projected
PDOS, with NAC, and the available experimental data[37, 38]. The frequency
spans across an interval of about 500 cm-1. The features to note from this
figure are the following. (i) As expected, the NAC leads to longitudinal
optical-traverse optical (LO-TO) splitting, near the $\Gamma$ point. The
strongest effects are felt by high-frequency optical modes. However, the
effects of the NAC on the calculated PDOS are quite small. (ii) Experimental
data are available only for infrared[37] and Raman[38, 12] active modes at the
$\Gamma$-points. The reported frequencies of the latter modes are in very good
agreement with each other, and hence only those of Ref. [38] are shown in Fig.
2. For the designation of these phonon modes see Ref.[40]. Fig 2 shows that
these experimental data agree reasonably well with our first-principles
results.

Figure 3: Calculated phonon dispersion relations of the (a) $\alpha$-quartz
and (b) $\alpha$-cristobalite phases of BeF2. (c) PDOS, and the Be and F
projected PDOS of the $\alpha$-cristobalite phase.

Fig. 3 shows the phonon dispersion relations of the $\alpha$-quartz and
$\alpha$-cristobalite phases of BeF2, taking into account the NAC. Also shown
are the PDOS, and Be and F projected PDOS of the $\alpha$-cristobalite phase.
The results of the $\alpha^{\prime}$-cristobalite phase are very similar to
those of $\alpha$-cristobalite and hence are not shown. The features to note
from this figure are the following. (i) The very wide frequency range of the
phonon modes in these systems, compared to that of ZnF2. This can be
understood as a consequence of the rather large mass difference between Be and
Zn atoms. (ii) The frequency range of both phases of BeF2 can be separated,
according to the character of the phonon modes, into three sub-regions. (a)
The lower frequency region between 0 and about 700 cm-1, where the phonon
modes are mainly due to the vibrations of F atoms. The contribution of the Be
atoms becomes appreciable above 300 cm-1. It is worth noting that in the case
of ZnF2 the dominance of the vibrations of the F atoms occurs in the upper
part of the frequency range because the Zn atom is heavier than the F atom.
(b) A narrow intermediate region at about 770 cm-1, where the rather localized
phonon modes originate from vibrations involving both Be and F atoms. (c) The
upper-frequency region, where the phonon modes originate mainly from
vibrations of Be atoms. (iii) The opening of two frequency gaps, between (a)
and (b), and between (b) and (c) sub-regions. These frequency gaps can be
understood as a consequence of the localization of phonon modes in the (b)
sub-region.

Figure 4: Calculated linear and volumetric TECs of ZnF2 using the LDA,
compared to the available experimental data[13] and previous theoretical
results[13, 14].

###  III.4 Thermal expansion

The calculated linear and volumetric TECs of the considered phases of ZnF2 and
BeF2, according to the procedure described in Sec. II, are depicted in Figs. 4
and 5, respectively. It is worth mentioning that, as expected, the NACs to the
dynamical matrices have negligible effects on the calculated linear TECs.

Figure 5: Calculated linear and volumetric TECs with PBE$\\_$sol of the three
considered phases of BeF2. Note the difference in the scales of the two
panels.

We will first look at the TECs of ZnF2. The important features to note from
Fig. 4 are the following. (i) The NTE at low temperatures is mostly due to
$\alpha_{a}$. The negative values of $\alpha_{c}$ are smaller (in magnitude)
than those of $\alpha_{a}$ and lie in a considerably shorter $T$-range. This
is clear from the magnitude and the location of the minimum values:
$\alpha_{a}\sim-1.05\times 10^{-6}$ K-1 at $55$ K, and
$\alpha_{c}\sim-0.5\times 10^{-6}$ K-1 at $40$ K. These results are consistent
with observed $T$-variations of the $a$ and $c$ lattice parameters at low
temperatures (see Fig. 3 of Ref. [13]). (ii) The calculated $\alpha_{v}$ from
the linear TECs (i.e., $\alpha_{v}=2\alpha_{a}+\alpha_{c}$) are in good
agreement with the previous direct theoretical calculations[13, 14], and the
results of all these theoretical calculations are in a qualitative agreement
with experimental data[13]. This finding reflects the accuracy and reliability
of our calculated linear TECs. (iii) $\alpha_{c}$ is systematically and
appreciably larger than $\alpha_{a}$. For example, at $300$ K the calculated
value of $\alpha_{c}$ (of $8.9\times 10^{-6}$ K-1) is about $60$% larger than
that of $\alpha_{a}$ (of $5.6\times 10^{-6}$ K-1).

As for the thermal expansion of the considered phases of BeF2, the features to
note from Fig. 5 are the following. (i) Unlike ZnF2, the calculated values of
both $\alpha_{a}$ and $\alpha_{c}$ are always positive for all of the
considered phases of BeF2. (ii) Both $\alpha_{c}$ and $\alpha_{a}$ of the
$\alpha$-cristobalite structure are very close to those of
$\alpha^{\prime}$-cristobalite, and hence only those of the former phase will
be discussed below. (iii) In the considered $T$-range, $\alpha_{a}(T)$ of the
$\alpha$-quartz structure is slightly larger than $\alpha_{c}(T)$, whereas
$\alpha_{a}(T)$ of the $\alpha$-cristobalite phase is much smaller than
$\alpha_{c}(T)$. (iv) The large $\alpha_{c}$ and $\alpha_{a}$ lead to very
large $\alpha_{v}$ for both phases of BeF2. For example, at $300$ K, the
values of $\alpha_{v}$ are of $77.6$ and $169.9$ $\times 10^{-6}$ K-1
respectively for the above two phases of BeF2, compared to that of $20.0\times
10^{-6}$ K-1 for ZnF2. Our largest calculated linear TEC is of $\sim 95\times
10^{-6}$ K-1 at $300$ K for $\alpha_{c}$ of the $\alpha$-cristobalite phase.
This is indeed large compared to the experimental linear TECs at $300$ K of
four fluorites, i.e., CaF2, SrF2, BaF2, and PbF2[41] that range between $18.1$
and $29~{}\times 10^{-6}$ K-1, but still is somewhat smaller than the measured
linear TEC value of $163.9\times 10^{-6}$ K-1 of an Ti-Nb alloy[42].

Figure 6: The PDOS weighted by the Grüneisen parameters of the rutile ZnF2,
and the $\alpha$-quartz and $\alpha$-cristobalite phases of BeF2.

The key physically insightful quantity for the interpretation of the above
results is the PDOS weighted by the Grüneisen parameters, $\Gamma_{i}(\nu)$,
defined in Eq. 2. Fig. 6 shows $\Gamma_{i}(\nu)$ of ZnF2, and the
$\alpha$-quartz and $\alpha$-cristobalite phases of BeF2. The important
features to note from this figure are the following. (i) For ZnF2, the low-
frequency modes ($\nu<150$ cm-1) have negative Grüneisen parameters, which
lead to negative $\Gamma_{i}(\nu)$ in this $\nu$ range. Since low-frequency
modes are easily thermally excited, this finding explains the observed NTE in
ZnF2. Moreover, by inspecting the differences between $\Gamma_{a}(\nu)$ and
$\Gamma_{c}(\nu)$ one can easily understand why $\alpha_{a}$ is always lower
than $\alpha_{c}$. (ii) The $\Gamma_{i}(\nu)$ of the considered phases of BeF2
are always positive, which reflects the dominance of positive mode Grüneisen
parameters in these phases. This explains the lack of NTE in the considered
phases of BeF2. (iii) The $\Gamma_{a}(\nu)$ and $\Gamma_{c}(\nu)$ of
$\alpha$-quartz BeF2 have comparable magnitudes, with $\Gamma_{c}(\nu)$ being
smaller than $\Gamma_{a}(\nu)$ for $\nu<100$ cm-1, which explain the
comparable magnitudes and ordering of its $\alpha_{c}$ and $\alpha_{a}$. (iv)
The peak in $\Gamma_{c}(\nu)$ of the $\alpha$-cristobalite phase around
$\nu\sim 34$ cm-1 is much higher than that of the $\Gamma_{a}(\nu)$, which
results in a large $\alpha_{c}$ compared to $\alpha_{a}$. This finding means
that, in this $\nu$-range, the positive mode Grüneisen parameters associated
with the out-of-plane deformation are significantly larger than those
associated with the in-plane deformation. The large Grüneisen parameters can
be viewed as a manifestation of strong anharmonic effects in the
$\alpha$-cristobalite and $\alpha^{\prime}$-cristobalite structures of BeF2.
However, it should be noted that large Grüneisen parameters are not the only
factor that is responsible for the giant $\alpha_{v}$ of the
$\alpha$-cristobalite: the elastic property via the negative (and with a large
magnitude) $C_{13}$ elastic constant (see Table 2 and Eq.1) plays also a major
role. The above findings explain the much larger volumetric TEC of the
$\alpha$-cristobalite phase of BeF2, compared to that of the $\alpha$-quartz
phase, which, in turn, is larger than that of ZnF2.

##  IV Summary

First-principles calculations are performed to investigate the structural,
elastic, and vibrational properties of the rutile structure of ZnF2 and three
crystal structures of BeF2 ($\alpha$-quartz, $\alpha$-cristobalite and its
similar phase with space group $P4_{3}2_{1}2$). The so-obtained phonon density
of states, mode Grüneisen parameters, and elastic constants are used to study
the linear thermal expansion coefficients (TECs) of the compounds mentioned
above, within a Grüneisen formalism. We have used deformations that preserve
the symmetry of the crystal to obtain the Grüneisen parameters. The considered
crystal structures of both systems are found to be mechanically stable. The
calculated physical quantities for both systems are in very good agreement
with the available experimental data and previous theoretical results. For
ZnF2, the calculated linear TECs, $\alpha_{a}$ and $\alpha_{c}$, along the $a$
and $c$ directions are consistent with the experimental $T$-variations of the
corresponding lattice parameters, respectively. The volumetric TEC
$\alpha_{v}$ computed from these linear TECs is in qualitative agreement with
experiment at low temperatures, including negative thermal expansion (NTE)
behavior. The considered phases of BeF2 are not NTE materials, and their
linear TECs are much higher than those of ZnF2, especially for the
$\alpha$-cristobalite phase. The elastic constants, high-frequency dielectric
constants, Born effective charge tensors, and TECs of the considered phases of
BeF2 are reported in this work for the first time and could serve as
predictions.

##  V Acknowledgment

We thank the National Supercomputing Center, Singapore (NSCC) and A*STAR
Computational Resource Center, Singapore (ACRC) for computing resources. We
also thank Iyad Al-Qasir for fruitful discussions. This work is supported by
RIE2020 Advanced Manufacturing and Engineering (AME) Programmatic Grant No
A1898b0043.

## References

  * Parker and France [2011] J. Parker and P. France, in  _Glasses and glass-ceramics, edited by M. H. Lewis_ (Chapman and Hall, 2011) p. 156.
  * Gan [1995] F. Gan, Optical properties of fluoride glasses: a review, Journal of Non-Crystalline Solids 184, 9 (1995).
  * Weber  _et al._ [1978] M. J. Weber, C. F. Cline, W. L. Smith, D. Milam, D. Heiman, and R. Hellwarth, Measurements of the electronic and nuclear contributions to the nonlinear refractive index of beryllium fluoride glasses, Appl. Phys. Lett. 32, 403 (1978).
  * Mei  _et al._ [2013] L. Mei, X. Cai, D. Jiang, J. Chen, W. Guo, and W. Xiong, Investigation of thermal neutron scattering data for BeF2 and LiF crystals, J. Nucl. Sci. Tech. 50, 419 (2013).
  * Kagawa  _et al._ [2004] R. Kagawa, M. Montgomery, K. Braig, A. Leslie, and J. Walker, The structure of bovine F1-ATPase inhibited by ADP and beryllium fluoride, EMBO J. 23, 2734 (2004).
  * Al-Qasir and Qteish [2017] I. Al-Qasir and A. Qteish, Neutron filter efficiency of beryllium and magnesium fluorides, J. Appl. Cryst. 50, 441 (2017).
  * Kaawar and Paulus [2015] Z. Kaawar and B. Paulus, A computational study of the structure of zinc fluoride surfaces, AIP Conf. Proc. 1653, 020051 (2015).
  * Ghalsasi and Ghalsasi [2011] P. Ghalsasi and P. S. Ghalsasi, Single crystal X-ray structure of BeF2: $\alpha$-quartz, Inorg. Chem. 50, 86 (2011).
  * Nelson  _et al._ [2017] J. R. Nelson, R. J. Needs, and C. J. Pickard, High-pressure phases of group II difluorides: polymorphism and superionicity, Phys. Rev. B 95, 054118 (2017).
  * Rakitin  _et al._ [2015] M. S. Rakitin, A. R. Oganov, H. Niu, M. M. D. Esfahani, X.-F. Zhou, G.-R. Qian, and V. L. Solozhenko, Novel phase of beryllium fluoride at high pressure, Phys. Chem. Chem. Phys. 17, 26283 (2015).
  * Masoumi [2021] N. Masoumi,  _First-principles DFT study of imide and fluoride analogs of silicon oxide, silicon oxynitride, and their alloys_ , Ph.D. thesis, Arizona State University, Cambridge, MA (2021).
  * Kurzydlowski  _et al._ [2020] D. Kurzydlowski, A. Oleksiak, S. B. Pillai, and P. K. Jha, High-pressure phase transitions of zinc difluoride up to 55 GPa, Inorg. Chem. 59, 2584 (2020).
  * Chatterji  _et al._ [2011] T. Chatterji, M. Zbiri, and T. C. Hansen, Negative thermal expansion in znf2, Appl. Phys. Lett. 98, 181911 (2011).
  * Wang  _et al._ [2012] L. Wang, P.-F. Yuan, F. Wang, Q. Sun, E.-J. Liang, and Y. Jia, Theoretical study of negative thermal expansion mechanism of ZnF2, Materials Research Bulletin 47, 1113 (2012).
  * Gan and Chua [2019] C. K. Gan and K. T. E. Chua, Large thermal anisotropy in monoclinic niobium trisulfide: A thermal expansion tensor study, J. Phys.: Condens. Matter 31, 265401 (2019).
  * Grüneisen [1926] E. Grüneisen, Handb. Phys. 10, 1 (1926).
  * Barron  _et al._ [1980] T. H. K. Barron, J. G. Collins, and G. K. White, Thermal expansion of solids at low temperatures, Adv. Phys. 29, 609 (1980).
  * Schelling and Keblinski [2003] P. K. Schelling and P. Keblinski, Thermal expansion of carbon structures, Phys. Rev. B 68, 035425 (2003).
  * Ding and Xiao [2015] Y. Ding and B. Xiao, Thermal expansion tensors, Grüneisen parameters and phonon velocities of bulk MT2 (M=W and Mo; T=S and Se) from first principles calculations, RSC Adv. 5, 18391 (2015).
  * Gan  _et al._ [2015] C. K. Gan, J. R. Soh, and Y. Liu, Large anharmonic effect and thermal expansion anisotropy of metal chalcogenides: The case of antimony sulfide, Phys. Rev. B 92, 235202 (2015).
  * Liu  _et al._ [2017] G. Liu, H. M. Liu, J. Zhou, and X. G. Wan, Temperature effect on lattice and electronic structures of WTe2 from first-principles study, J. Appl. Phys. 121, 045104 (2017).
  * Liu and Allen [2018] J. Liu and P. B. Allen, Internal and external expansions of wurtzite ZnO from first principles, Comput. Mater. Sci. 154, 251 (2018).
  * Gan  _et al._ [2021] C. K. Gan, Y. Liu, T. C. Sum, and K. Hippalgaonkar, Efficacious symmetry-adapted atomic displacement method for lattice dynamical studies, Comput. Phys. Comm. 259, 107635 (2021).
  * Frank  _et al._ [1995] W. Frank, C. Elsässer, and M. Fähnle, Ab initio force-constant method for phonon dispersions in alkali metals, Phys. Rev. Lett. 74, 1791 (1995).
  * 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).
  * Kresse  _et al._ [1995] G. Kresse, J. Furthmüller, and J. Hafner, Ab initio force constant approach to phonon dispersion relations of diamond and graphite, Europhys. Lett. 32, 729 (1995).
  * Gan  _et al._ [2010] C. K. Gan, X. F. Fan, and J.-L. Kuo, Composition-temperature phase diagram of BexZn1-xO from first principles, Comput. Mater. Sci. 49, S29 (2010).
  * Liu  _et al._ [2014] Y. Liu, K. T. E. Chua, T. C. Sum, and C. K. Gan, First-principles study of the lattice dynamics of Sb2S3, Phys. Chem. Chem. Phys. 16, 345 (2014).
  * Dal Corso [2016] A. Dal Corso, Elastic constants of beryllium: a first-principles investigation, J. Phys.: Condens. Matter 28, 075401 (2016).
  * Gan and Lee [2018] C. K. Gan and C. H. Lee, Anharmonic phonon effects on linear thermal expansion of trigonal bismuth selenide and antimony telluride crystals, Comput. Mater. Sci. 151, 49 (2018).
  * F.Yu  _et al._ [2013] F.Yu, M. Xu, M. Jiang, and J.-X. Sun, The phase transitions and electronic structures of crystalline BeF2 under high-pressure: First-principle calculations, Solid State Commun. 169, 14 (2013).
  * O’Toole and Streltsov [2001] N. O’Toole and V. Streltsov, Synchrotron x-ray analysis of the electron density in CoF2 and ZnF2, Acta Cryst. B 57, 128 (2001).
  * Rimai [1977] D. S. Rimai, Elastic properties of ZnF2 between 4.2 and 300 K, Phys. Rev. B 16, 4069 (1977).
  * Pabst and Gregorová [2013] W. Pabst and E. C. Gregorová, Elastic properties of silica polymorphs - A review, Ceramics-Silikáty 57, 167 (2013).
  * Mouhat and Coudert [2014] F. Mouhat and F.-X. Coudert, Necessary and sufficient elastic stability conditions in various crystal systems, Phys. Rev. B 90, 224104 (2014).
  * Alaoui-Bichri  _et al._ [1980] R. E. Alaoui-Bichri, J. Giordano, R. Almairac, C. Benoit, and P. Nassiri, Properties, J. Physique 41, 543 (1980).
  * A. S. Barker [1964] J. A. S. Barker, Transverse and longitudinal optic mode study in MgF2 and ZnF2, Phys. Rev. 136, A1290 (1964).
  * Porto  _et al._ [1967] S. P. S. Porto, P. A. Fleury, and T. C. Damen, Raman spectra of TiO2, MgF2, ZnF2, FeF2, and MnF2, Phys. Rev. 154, 522 (1967).
  * Qteish [2019] A. Qteish, Electronegativity scales and electronegativity-bond ionicity relations: A comparative study, J. Phys. Chem. Solids 124, 186 (2019).
  * Benoit and Giordano [1988] C. Benoit and J. Giordano, Dynamical properties of crystals of MgF2, ZnF2 and FeF2. II. Lattice dynamics and infrared spectra, J. Phys. C: Solid State Phys 21, 5209 (1988).
  * Roberts and White [1986] R. B. Roberts and G. K. White, Thermal expansion of fluorites at high temperatures, J. Phys. C: Solid State Phys. 19, 7167 (1986).
  * Bönisch  _et al._ [2017] M. Bönisch, A. Panigrahi, M. Stoica, M. Calin, E. Ahrens, M. Zehetbauer, W. Skrotzki, and J. Eckert, Giant thermal expansion and $\alpha$-precipitation pathways in Ti-alloys, Nature Commun. 8, 1429 (2017).