# First-principles estimation of the superconducting transition temperature of
a metallic hydrogen liquid

Haoran Chen  International Center for Quantum Materials, Peking University,
Beijing 100871, China  Xiao-Wei Zhang  International Center for Quantum
Materials, Peking University, Beijing 100871, China  Xin-Zheng Li  State Key
Laboratory for Artificial Microstructure and Mesoscopic Physics, and School of
Physics, Peking University, Beijing 100871, China  Junren Shi
<EMAIL_ADDRESS>International Center for Quantum Materials, Peking
University, Beijing 100871, China  Collaborative Innovation Center of Quantum
Matter, Beijing 100871, China

###### Abstract

We present a full density-functional-theory-based implementation of the
stochastic path-integral approach proposed by Liu et al .[1] for estimating
the superconducting transition temperature ($T_{c}$) of a liquid. The
implementation is based on the all-electron projector-augmented-wave (PAW)
method. We generalize Liu et al.’s formalism to accommodate the pseudo-
description of electron states in the PAW method. A formula for constructing
the overlap operator of the PAW method is proposed to eliminate errors due to
the incompleteness of a pseudo-basis set. We apply the implementation to
estimate $T_{c}$’s of metallic hydrogen liquids. It confirms Liu et al.’s
prediction that metallic hydrogen could form a superconducting liquid.

##  I Introduction

Ever since the discovery of the superconductivity by K. Onnes in 1911, the
search for materials with high superconducting transition temperatures
($T_{c}$’s) had been a constant pursuit. The recent discovery of a new family
of high-$T_{c}$ superconductors in hydrides under high pressure shows the
great promise of first-principles calculations in guiding the search [2, 3].
Based on the celebrated Bardeen-Copper-Schrieffer (BCS) theory of
superconductivity, it was long predicted that metallic hydrogen or hydrogen-
rich materials could exhibit high $T_{c}$’s as light hydrogen atoms provide
both a high Debye frequency of phonons and strong electron-phonon coupling
(EPC) [4, 5, 6, 7]. By using state-of-art first-principles approaches based on
the density-functional theory (DFT), it is now possible to conduct an
extensive search in various hydrides, and predict their structures and
$T_{c}$’s [2, 3, 8, 9, 10].

In spite of the tremendous success, the first-principles prediction of $T_{c}$
for compounds like hydrides is tricky. The widely-adopted algorithm is based
on the Eliashberg theory which relates $T_{c}$ to quantities calculable by DFT
such as the band structure of electrons and the matrix elements of EPC [11,
12]. The theory is built upon the assumption that the vibration amplitudes of
atoms in a solid are small. It thus adopts the harmonic approximation for
describing atom vibrations and a perturbative approach for determining EPC
effects. In hydrides, however, the assumption likely breaks down. Light
hydrogen atoms tend to have large vibration amplitudes, and therefore their
vibrations are subject to anharmonicity [13]. Moreover, hydrogen atoms can
tunnel between lattice sites and could even form a superionic phase in which
they become delocalized and diffuse in a whole lattice [14]. In the mother
compound of metallic hydrides, i.e., the metallic hydrogen, quantum
fluctuations induced by the tunneling is so strong that its melting point is
predicted to be suppressed below the room temperature [15, 16]. The
applicability of the Eliashberg theory in these systems is questionable.

To date, most of the efforts addressing the quantum tunneling and anharmonic
effects in superconductors are based on the approach of the stochastic self-
consistent harmonic approximation (SSCHA) [17] which treats an anharmonic
system as an effective harmonic one with self-consistently determined phonon
modes. $T_{c}$ can then be estimated by applying the standard algorithm for a
harmonic system. SSCHA calculations show that the anharmonicity and quantum
fluctuations could yield different predictions of crystal structures and
$T_{c}$’s [13, 18]. While the approach is effective and efficient, it relies
on a cascade of approximations which could become uncontrolled. In particular,
when a system enters a liquid or superionic phase, the underlying assumption
that the system can be approximated as a harmonic solid obviously breaks down.

Recently, Liu et al. propose a stochastic path-integral approach of
determining $T_{c}$ for general systems including liquids [1]. The core of the
approach is a set of rigorous relations through which an effective pairing
interaction between electrons can be inferred from the fluctuations of
electron-ion scattering $\mathcal{T}$-matrices. The fluctuations of the
$\mathcal{T}$-matrices can be estimated from a stochastic sampling of ion
positions generated by a path integral molecular dynamics (PIMD) simulation
[19]. Furthermore, it is shown that the equations for determining $T_{c}$ in
the Eliashberg theory, i.e., the linearized Eliashberg equations, can be
reestablished rigorously in the new context. As a result, with the effective
pairing interaction, $T_{c}$ can be determined in the exactly same manner as
in the conventional Eliashberg approach. The approach is based on a rigorous
theory and makes no assumption on the nature of atom motion. It should
therefore apply equally well for solids or liquids, with or without
anharmonicity and quantum fluctuations.

Liu et al. present an implementation of the approach for metallic hydrogen
liquids in Ref. 1. They predict that metallic hydrogen would form a
superconducting liquid with $T_{c}$ well above its melting temperature. In the
implementation, hydrogen atom positions are sampled by using the first-
principles PIMD method. However, ionic fields, i.e., the interaction
potentials between hydrogen ions and electrons, are approximated as a bare
Coulomb potential screened by a dielectric function. Such a linear screening
approximation (LSA) should be adequate for metallic hydrogen because hydrogen
ions, i.e., protons, are point particles without inner-core electrons. For
more general systems involving other atom species as in the case of hydrides,
however, the approximation is obviously inadequate. The deficiency limits the
applicability of the implementation to metallic hydrogen only.

In this paper, we develop a full DFT-based implementation of the stochastic
path-integral approach for determining $T_{c}$. We base our implementation on
the Vienna Ab-initio Simulation Package (VASP) [20, 21], and determine ionic
fields by employing the projector augmented-wave (PAW) method [22, 23]. In the
method, the states of valance electrons are described by pseudo (PS)
wavefunctions, and the all-electron (AE) wavefunctions, i.e., the true
wavefunctions of electrons, are related to PS-wavefunctions through a
transformation which depends on the positions of ions. It thus introduces
complexities in directly applying Liu et al.’s formalism. We solve the issue
by generalizing the formalism and accommodating it to the PS description. The
generalized formalism enables our full DFT-based implementation. As a test, we
apply the implementation to metallic hydrogen liquids. It yields results close
to Liu et al.’s results, confirming the prediction that metallic hydrogen
could form a superconducting liquid.

The remainder of the paper is organized as follows. In Sec. II, we outline the
formalism and calculation procedure of the stochastic path-integral approach
proposed by Liu et al., and identify main approximations in their
implementation. In Sec. III, we generalize and implement Liu et al.’s
formalism to accommodate the PAW method. A number of issues relevant to the
implementation are also addressed. In particular, a formula for constructing
the overlap operator in the PAW method is proposed to eliminate errors due to
the incompleteness of PS basis set. In Sec. IV.2, we apply the implementation
to metallic hydrogen liquids, and compare results with Liu et al.’s results.
The convergency of the calculation is also tested. Finally, we summarize our
results in Sec. V.

##  II Stochastic path-integral approach of determining $T_{c}$

In this section, we outline the formalism of the stochastic path-integral
approach of determining $T_{c}$ proposed by Liu et al. and discuss
approximations involved in implementing the approach.

###  II.1 Formalism

In the conventional Eliashberg theory, $T_{c}$ is determined by solving a set
of linearized Eliashberg equations. Liu et al. show rigorously that the same
set of the equations still hold for liquids and general systems, provided that
relevant parameters are properly determined. The set of the linearized
Eliashberg equations reads

$\displaystyle\rho\Delta_{n}$ $\displaystyle=$
$\displaystyle\sum_{n^{\prime}}\left[\lambda(n^{\prime}-n)-\mu^{\star}-\frac{\hbar\beta}{\pi}|\tilde{\omega}(n)|\delta_{nn^{\prime}}\right]\Delta_{n^{\prime}},\quad$
(1) $\displaystyle\tilde{\omega}(n)$ $\displaystyle=$
$\displaystyle\frac{\pi}{\hbar\beta}\left(2n+1+\lambda(0)+2\sum_{m=1}^{n}\lambda(m)\right),n\geq
0\quad$ (2)

and $|\tilde{\omega}(-n)|=|\tilde{\omega}(n-1)|$, where $n$ and $n^{\prime}$
are integers indexing the Matsubara frequencies of Fermions
$\omega_{n}=(2n+1)\pi/\hbar\beta$, $\beta\equiv 1/k_{B}T$ and $T$ is the
temperature. $\\{\lambda(m),\,m\in Z\\}$ is a set of interaction parameters
characterizing the effective pairing interaction, and $\mu^{*}$ is the
empirical Morel-Anderson pseudopotential introduced to take account of the
exchange-correlation effect of electrons [24]. Solving Eq. (1) yields a set of
eigenvalues $\rho$ and eigenvectors $\Delta_{n}$. For normal states, all
eigenvalues are negative. The emergence of a non-negative eigenvalue indicates
that the system enters into the superconducting state. One can thus estimate
$T_{c}$ by solving the equations for a number of temperatures to see when the
maximal eigenvalue changes sign. We note that $\Delta_{n}$ here is not related
to the physical gap.

Two comments are in order. First, in calculations for solids using the
standard approach (see, e.g., Ref. 7), it is common to use the empirical
McMillan formula and estimate $T_{c}$ from a small number of parameters
evaluated at the zero temperature including the mass enhancement factor
$\lambda\equiv\lambda(0)$, an average phonon frequency and $\mu^{*}$. It
simplifies the calculations. Unfortunately, such an approach does not work
well for liquids because (a) unlike harmonic solids, $\lambda$ and the average
phonon frequency of a liquid are temperature-dependent; (b) the empirical
formulas are tuned for solids, and do not fit well with liquids which have
qualitatively different density spectral functions due to the presence of low-
frequency diffusive modes [25, 26]. Second, the choice of parameter $\mu^{*}$
is empirical. For metallic hydrogen systems, we use $\mu^{*}\approx 0.089$ as
suggested in Ref. 7. The empirical parameter introduces an uncertainty in
predicting $T_{c}$. Ideally, one should base calculations on the DFT for
superconductors [27, 28, 29], which provides a rigorous framework for treating
exchange-correlation effects in superconductors. This is not yet implemented
in our current calculation.

To apply the Eliashberg equations, the interaction parameters $\lambda(m)$ are
the central quantities to be determined. They are related to the effective
pairing interaction $W$ by

$\displaystyle\lambda(n-n^{\prime})=-\sum_{\bm{k}^{\prime}}W(\omega_{n}-\omega_{n^{\prime}},\bm{k}-\bm{k}^{\prime})\delta(\tilde{\epsilon}_{\bm{k}^{\prime}}-\epsilon_{F}),$
(3)

where $W(\omega_{n}-\omega_{n^{\prime}},\bm{k}-\bm{k}^{\prime})$ denotes the
matrix element of the effective pairing interaction, which is a function of
the quasi-wave-vector transfer $\bm{k}-\bm{k}^{\prime}$ and the frequency
transfer $\omega_{n}-\omega_{n^{\prime}}$,
$\tilde{\epsilon}_{\bm{k}^{\prime}}$ is the renormalized quasi-electron
dispersion, and $\epsilon_{F}$ is the Fermi energy. To apply the formula, it
is necessary to determine the effective pairing interaction as well as the
renormalized quasi-electron dispersion.

To determine these quantities in a liquid, Liu et al. propose a stochastic
path integral approach. It is based on the classical isomorphism [30], i.e.,
mapping quantum ions to classical ring polymers by applying the imaginary-time
path integral formalism [31]. The spatial configuration of the classical
ensemble of ring polymers is specified by a set of $\tau$-dependent
coordinates $\\{\bm{R}_{a}(\tau)\\}$, where $\tau\in[0,\hbar\beta)$ is the
imaginary time or interpreted as a variable to parameterize the rings. In a
PIMD simulation, the isomorphic ensemble is simulated by using the classical
molecular dynamics. It provides a statistical sampling of the spatial
configurations of the ring polymers.

Quantities related to electrons can be evaluated in the isomorphic ensemble.
Electrons in the ensemble experience a $\tau$-dependent Kohn Sham potential
$V_{\mathrm{KS}}(\tau)$ which is self-consistently determined by the ion
configuration $\\{\bm{R}_{a}(\tau)\\}$ at the specific imaginary-time $\tau$.
As a result, the single-particle Green’s function of electrons with respect to
a given configuration of ring polymers, in the context of the Kohn-Sham
theory, can be determined by

$\left[-\frac{\partial}{\partial\tau}-\frac{\hat{H}(\tau)-\epsilon_{F}\hat{\mathbb{I}}}{\hbar}\right]\hat{\mathcal{G}}(\tau,\tau^{\prime})=\delta(\tau-\tau^{\prime})\hat{\mathbb{I}},$
(4)

where we write the equation in an operator form,
$\hat{H}(\tau)\equiv\hat{K}+\hat{V}_{\mathrm{KS}}(\tau)$ with $\hat{K}$ being
the operators of the kinetic energy. The Green’s function is not a physical
Green’s function for electrons in a liquid. To get the physical one, we need
to further conduct an ensemble average:

$\displaystyle\mathcal{G}_{\mathrm{phy.}}\equiv\bar{\mathcal{G}}(\tau-\tau^{\prime})=\langle\mathcal{G}(\tau,\tau^{\prime})\rangle_{C},$
(5)

where $\langle\dots\rangle_{C}$ denotes the average over the ring-polymer
configurations. The averaged (physical) Green’s function
$\bar{\mathcal{G}}(\tau-\tau^{\prime})$ defines the dispersion and lifetime of
a quasi-electron renormalized by the fluctuating ionic field. The quasi-
electron dispersion $\tilde{\epsilon}_{\bm{k}^{\prime}}$ required by Eq. (3)
can then be inferred from it.

For a given spatial configuration of ring polymers, we can determine the
electron-ion scattering $\mathcal{T}$-matrix. The $\mathcal{T}$-matrix is
related to the Green’s function according to the Lippmann-Schwinger equation

$\displaystyle\hat{\mathcal{T}}=\hat{\mathcal{V}}+\frac{1}{\hbar}\hat{\mathcal{V}}\hat{\bar{\mathcal{G}}}\hat{\mathcal{T}},$
(6)

where $\hat{\mathcal{V}}=\hat{V}_{\mathrm{KS}}(\tau)-\hat{\bar{\Sigma}}$ is
the effective scattering potential for electrons, and $\hat{\bar{\Sigma}}$ is
the self-energy inferred from $\hat{\bar{\mathcal{G}}}$.

Liu et al. establish a rigorous relation between the effective pairing
interaction and the $\mathcal{T}$-matrices. One first determines the pair
scattering amplitude $\Gamma$, i.e., the fluctuation of the
$\mathcal{T}$-matrices:

$\displaystyle\Gamma_{11^{\prime}}=-\beta\langle|\mathcal{T}_{11^{\prime}}|^{2}\rangle_{C},$
(7)

where $\mathcal{T}_{11^{\prime}}$ denotes the matrix element of the
$\mathcal{T}$-matrix, and we adopt a short-hand notation by using numbers to
represent the states of electrons: $1\equiv(\bm{k},\omega_{n})$,
$\bar{1}\equiv(-\bm{k},-\omega_{n})$, and similarly for $1^{\prime}$ and
$\bar{1}^{\prime}$. Then, the effective pairing interaction can be obtained by
solving the Bethe-Salpeter equation:

$\displaystyle
W_{11^{\prime}}=\Gamma_{11^{\prime}}+\frac{1}{\hbar^{2}\beta}\sum_{2}\Gamma_{12}|\bar{\mathcal{G}}_{2}|^{2}W_{21^{\prime}}.$
(8)

The obtained $W_{11^{\prime}}$ then enters into Eq. (3), in which we express
it explicitly as $W_{11^{\prime}}\equiv
W(\omega_{n}-\omega_{n^{\prime}},\bm{k}-\bm{k}^{\prime})$. We note that in
general $W_{11^{\prime}}$ depends on $1$ and $1^{\prime}$ independently. For
EPC, it is a good approximation to regard $W_{11^{\prime}}$ as a function of
$\omega_{n}-\omega_{n^{\prime}}$ as long as $\hbar|\omega_{n(n^{\prime})}|$ is
much smaller than the Fermi energy. Moreover, since liquids are isotropic and
only matrix elements with both $\bm{k}$ and $\bm{k}^{\prime}$ residing on the
Fermi surface are needed for evaluating Eq. (3), $W_{11^{\prime}}$ can be
written as a function of $\bm{k}-\bm{k}^{\prime}$.

###  II.2 Approximations

To develop an algorithm for estimating $T_{c}$ based on the the rigorous
formalism shown in the last subsection, Liu et al. employ a number of
approximations. In this subsection, we will discuss these approximations, and
identify the one to be improved.

The most important approximation which makes an implementation possible is the
quasi-static approximation. It is necessary because solving Eqs. (4) and (22)
exactly is deemed to be impossible in practice as it requires a time
resolution $\sim\hbar/\epsilon_{F}$ or a corresponding large cutoff frequency
$\sim\epsilon_{F}/\hbar$. Fortunately, one can exploit the fact that ions move
much slower than electrons. It means that the ionic field only has significant
components for frequencies $\nu_{m}\equiv 2\pi
m/\hbar\beta\lesssim\omega_{\mathrm{ph}}$, where $\omega_{\mathrm{ph}}$ is a
characteristic frequency of phonons. In this case, the equation can be solved
by applying the quasi-static approximation. This is to choose a large $n$ such
that $\omega_{n}\gg\omega_{\mathrm{ph}}$, and solve Eq. (4) as if
$\hat{H}(\tau)$ is time-independent:

$\displaystyle\hat{\mathcal{G}}_{n}(\tau)=\hbar\left[\left(\mathrm{i}\hbar\omega_{n}+\epsilon_{F}\right)\hat{\mathbb{I}}-\hat{H}(\tau)\right]^{-1},$
(9)

The ensemble average
$\hat{\bar{\mathcal{G}}}_{n}=\langle\hat{\mathcal{G}}_{n}(\tau)\rangle_{C}$
will be independent of $\tau$ and a good approximation for the Fourier
transform of $\bar{\mathcal{G}}(\tau-\tau^{\prime})$ at the frequency
$\omega_{n}$ [1]. The $\mathcal{T}$-matrix can then be calculated in the time-
domain in the same manner, and Fourier transformed back to the frequency
domain. See Sec. III.2 for details. For determining the effective pairing
potential, we only need to solve for a specific $n$ with
$\omega_{\mathrm{ph}}\ll\omega_{n}\ll\epsilon_{F}/\hbar$. This is because
$W_{11^{\prime}}\equiv
W(\omega_{n}-\omega_{n^{\prime}},\bm{k}-\bm{k}^{\prime})$ is assumed to be a
function of $\omega_{n}-\omega_{n}^{\prime}\equiv\nu_{m}$.

More approximations are involved in PIMD simulations. It is necessary to
simulate a system in a finite size supercell, and discretize the imaginary
time to a finite number of beads. In Ref. 1, it is found that the
discretization truncates the high-frequency components of the density
correlation function of ions and the effective pairing interaction. It results
in an incorrect asymptotic behavior of $\lambda(m)$ at large $m$. To address
the issue, Liu et al. develop an oversampling scheme to interpolate the
density correlation function in the temporal domain, and obtain an over-
sampled density correlation function. See Ref. 1 for details. To get the
effective pairing interaction, they make use of the empirical relation:

$W(\nu_{m},\bm{q})=|M(\nu_{m},\bm{q})|^{2}\chi(\nu_{m},\bm{q}),$ (10)

where $\chi(\nu_{m},\bm{q})$ is the density correlation function of ions, and
$M(\nu_{m},\bm{q})$ could be regarded as the EPC matrix element as defined in
the Eliashberg theory. It is found numerically that the relation holds well
for metallic hydrogen liquids, and $|M(\nu_{m},\bm{q})|^{2}$ can be obtained
from a linear regression. By substituting the density correlation function in
Eq. (10) with the over-sampled one, one can obtain an over-sampled effective
pairing interaction. It is found that $\lambda(m)$ determined from the over-
sampled effective pairing interaction has a correct asymptotic behavior at
large $m$.

Lastly, Liu et al. adopt the LSA for evaluating the ionic field. This is to
approximate the Fourier transform of the ionic potential as

$V_{\mathrm{ion}}(\bm{q},\tau)=\frac{e^{2}}{\varepsilon_{\mathrm{et}}(q)q^{2}}\rho_{\mathrm{i}}(\bm{q},\tau),$
(11)

where
$\rho_{\mathrm{i}}(\bm{q},\tau)=\sum_{a}\exp[-\mathrm{i}\bm{q}\cdot\bm{R}_{a}(\tau)]$
is the Fourier transform of the density distribution of ions, and
$\varepsilon_{\mathrm{et}}(q)$ is the static electron-test charge dielectric
function [32] which has a form similar to that of the random phase
approximation but with a local field correction factor determined by Ichimaru
et al. [33]. It is derived from the jellium model which models ions as a
uniform positive charge background, an approximation which obviously breaks
down in core regions close to ions. While it is reasonable to expect that such
an approximation could work well for metallic hydrogen in which ions (protons)
have small core volumes and electrons are shared by all ions, it becomes
problematic when other atoms are involved and bonding between atoms is
covalent in nature.

Among the three approximations, the LSA is the one that severely limits the
applicability of Liu et al.’s implementation for systems other than metallic
hydrogen. This will be the focus of this paper aiming to improve the
implementation. The effect of the finite simulation supercell and the finite
number of beads should be assessed to make sure that calculations converge.
This will also be tested. Finally, it is the quasi-static approximation which
makes an implementation feasible for EPC. The approximation is valid because
phonons have an energy scale much smaller than that of electrons, a feature
unique to conventional EPC superconductors.

##  III Formalism for the PAW method

In this section, we develop formalism for a full DFT-based implementation. To
balance speed and accuracy, we base our implementation on the PAW method [22,
23]. As a pseudo-potential method, the PAW method is efficient. Meanwhile, one
can obtain AE wavefunctions from PS wavefunctions by applying a
transformation. The latter is particularly useful for us since the formalism
presented in Sec. II is established in an AE basis. To facilitate an
implementation based on the PAW method, we need formulae of determining the AE
Green’s functions and the $\mathcal{T}$-matrices by using quantities defined
in the PS basis.

In Sec. III.1, we first briefly review the PAW method and outline its
formalism. With the preparation, we derive formulae for determining the AE
Green’s functions and the $\mathcal{T}$-matrices by using the PAW method in
Sec. III.2. In Sec. III.3, we discuss the evaluations of quantities involved
in the formulae, including the transformation matrix and the overlap matrix.

###  III.1 PAW method

In the PAW method, electrons are described by PS wavefunctions. A
transformation operator is introduced to map a PS wavefunction $\tilde{\Psi}$
to an AE wave function $\Psi$:

$\displaystyle\Psi=\hat{T}\tilde{\Psi}\equiv\tilde{\Psi}+\sum_{a}\hat{T}^{a}|\tilde{\Psi}\rangle,$
(12)

where $\hat{T}^{a}$ is an operator transforming the core state of the ion $a$
from the PS space to the AE space, and is nonzero only within an augmentation
sphere surrounding the ion. The radius of the sphere is chosen such that there
is no overlap between spheres associated with different ions. Inside a sphere,
the wave functions are expanded in the basis of partial waves:

$\displaystyle\hat{T}^{a}|\tilde{\Psi}\rangle=|\Psi^{a}\rangle-|\tilde{\Psi}^{a}\rangle=\sum_{i}\left(|\phi_{i}^{a}\rangle-|\tilde{\phi}_{i}^{a}\rangle\right)\langle\tilde{p}_{i}^{a}|\tilde{\Psi}\rangle,\quad$
(13)

where $i$ indexes the partial wave states with angular momentum quantum
numbers $(l_{i}m_{i})$ and at a given reference energy $\epsilon_{i}$. The AE
partial wavefunctions
$\langle\bm{r}|\phi_{i}^{a}\rangle\equiv\phi_{l_{i}\epsilon_{i}}^{a}(r)Y_{l_{i}m_{i}}(\hat{\bm{r}})$
at the radial coordinate $r$ and along the direction $\hat{\bm{r}}$ are
obtained by solving a spherical AE Schrodinger equation at the reference
energy in the presence of an ion at the center. The corresponding PS partial
wavefunctions $\langle\bm{r}|\tilde{\phi}_{i}^{a}\rangle$ are chosen such that
they are smooth and coincide with the AE partial wavefunctions outside the
sphere. $\langle\tilde{p}_{i}^{a}|$ are a set of projectors which are
biorthogonal to the PS partial wave states:

$\displaystyle\langle\tilde{p}_{j}^{a}|\tilde{\phi}_{i}^{a}\rangle=\delta_{i,j}.$
(14)

With a complete set of the basis, the closure relation in each sphere can be
written as:

$\displaystyle\sum_{i}|\tilde{\phi}_{i}^{a}\rangle\langle\tilde{p}_{i}^{a}|=\hat{\mathbb{I}}^{a}.$
(15)

In the PS basis, the Kohn-Sham (KS) equation is transformed to

$\displaystyle\hat{\tilde{H}}\tilde{\Psi}_{n}$ $\displaystyle=$
$\displaystyle\epsilon_{n}\hat{S}\tilde{\Psi}_{n},$ (16a)
$\displaystyle\hat{\tilde{H}}$ $\displaystyle=$
$\displaystyle\hat{T}^{\dagger}\hat{H}\hat{T},$ (16b) $\displaystyle\hat{S}$
$\displaystyle=$ $\displaystyle\hat{T}^{\dagger}\hat{T},$ (16c)

where $\hat{H}$ and $\hat{\tilde{H}}$ denote the Hamiltonians in the AE and PS
basis, respectively. This is a generalized eigenvalue problem, and the
orthogonality condition between PS eigen-wavefunctions becomes

$\displaystyle\langle\Psi_{m}|\Psi_{n}\rangle=\langle\tilde{\Psi}_{m}|\hat{S}|\tilde{\Psi}_{n}\rangle=\delta_{m,n}.$
(17)

We note that in DFT calculations based on the PAW method, the PS Hamiltonian
is directly derived from a total energy functional expressed in terms of the
PS wavefunctions. It is not necessary to construct the AE Hamiltonian first
and then transform it to the PS Hamiltonian by using Eq. (16b).

###  III.2 AE Green’s function and $\mathcal{T}$-matrix

With the preparation, we now proceed to derive formulae for determining the AE
Green’s function and the $\mathcal{T}$-matrix by using the PAW method.

We seek for an expression of the Green’s function appropriate for the PAW
method. Under the quasi-static approximation, the AE Green’s function is
determined by Eq. (9). However, it depends on the AE Hamiltonian
$\hat{H}(\tau)$. We can transform it to a form that only depends on quantities
in the PS basis. To see that, we start from the spectral representation of a
general Green’s function
$\hat{\mathcal{G}}(\epsilon)=\hbar(\epsilon\hat{\mathbb{I}}-\hat{H})^{-1}$:

$\displaystyle\hat{\mathcal{G}}(\epsilon)$ $\displaystyle=$
$\displaystyle\hbar\sum_{n}\frac{|\Psi_{n}\rangle\langle\Psi_{n}|}{\epsilon-\epsilon_{n}}$
(18) $\displaystyle=$
$\displaystyle\hat{T}\left[\hbar\sum_{n}\frac{|\tilde{\Psi}_{n}\rangle\langle\tilde{\Psi}_{n}|}{\epsilon-\epsilon_{n}}\right]\hat{T}^{\dagger},$

where $|\Psi_{n}\rangle$ and $|\tilde{\Psi}_{n}\rangle$ are the AE and PS
eigen-wavefunctions, respectively, and $\epsilon_{n}$ is the corresponding
eigen-energy.

We can convert the expression to a more convenient form. By using the
orthonormality condition Eq. (17), it is easy to verify the closure relation

$\displaystyle\sum_{n}\hat{S}|\tilde{\Psi}_{n}\rangle\langle\tilde{\Psi}_{n}|=\hat{\mathbb{I}}.$
(19)

We then have

$\displaystyle\left(\epsilon\hat{S}-\hat{\tilde{H}}\right)^{-1}=$
$\displaystyle\sum_{n}\left(\epsilon\hat{S}-\hat{\tilde{H}}\right)^{-1}\hat{S}|\tilde{\Psi}_{n}\rangle\langle\tilde{\Psi}_{n}|$
$\displaystyle=$
$\displaystyle\sum_{n}\frac{|\tilde{\Psi}_{n}\rangle\langle\tilde{\Psi}_{n}|}{\epsilon-\epsilon_{n}},$
(20)

where we make use of the eigen-equation Eq. (16a). By substituting the
equality into Eq. (18) and setting
$\epsilon=\mathrm{i}\hbar\omega_{n}+\epsilon_{F}$, we obtain the formula of
determining the AE Green’s function in the PAW method:

$\hat{\mathcal{G}}_{n}(\tau)=\hbar\hat{T}(\tau)\left[(\mathrm{i}\hbar\omega_{n}+\epsilon_{F})\hat{S}(\tau)-\hat{\tilde{H}}(\tau)\right]^{-1}\hat{T}^{\dagger}(\tau).\quad$
(21)

We note that $\hat{\tilde{H}}$, $\hat{S}$ and $\hat{T}$ all depend on the
configuration of ions, therefore are time-dependent.

With the Green’s function, the scattering $\mathcal{T}$-matrix can be
determined. To avoid solving the Lippmann-Schwinger equation in the AE space,
we apply an identity

$\hat{\mathcal{T}}_{n}(\tau)=\hbar\hat{\bar{\mathcal{G}}}_{n}^{-1}\hat{\mathcal{G}}(\tau)\hat{\bar{\mathcal{G}}}_{n}^{-1}-\hbar\hat{\bar{\mathcal{G}}}_{n}^{-1}$
(22)

to calculate the $\mathcal{T}$-matirx in time domain. It is easy to prove that
the $\mathcal{T}$-matirx calculated in this way satisfies Eq. (6). We note
that in a liquid, because of the spatial translation symmetry, the average
Green’s function is diagonal in the plane-wave basis. As a result, the matrix
inverse and multiplication can be carried out in a truncated Hilbert space
spanned by the basis without introducing errors.

The matrix elements of $\mathcal{T}$-matrix in frequency domain can then be
determined by applying the Fourier transform:

$\displaystyle\hat{\mathcal{T}}_{\omega_{n}+\nu_{m},\omega_{n}}\approx\frac{1}{\hbar\beta}\int_{0}^{\hbar\beta}d\tau\hat{\mathcal{T}}_{n}(\tau)\mathrm{e}^{\mathrm{i}\nu_{m}\tau}.$
(23)

The matrix elements are all that are needed for evaluating the pair scattering
amplitude and the effective pairing interaction (see Sec. II.1). We note that
for the purpose only the matrix elements with wave-vectors close to the Fermi
surface are needed. Therefore, a large energy cutoff, which is required for
recovering the AE wavefunctions from the PS wavefunctions, is not needed for
our calculation.

###  III.3 Transformation and overlap matrices

To apply the formulae derived in the last subsection, it is necessary to
determine the transformation operator $\hat{T}$ and the overlap matrix
$\hat{S}$.

The transformation operator transforms a PS state to an AE state. From Eq.
(12) and Eq. (13), the operator is

$\hat{T}=\hat{\mathbb{I}}+{\sum_{j,a}}^{\prime}\left(|\phi_{j}^{a}\rangle-|\tilde{\phi}_{j}^{a}\rangle\right)\langle\tilde{p}_{j}^{a}|,$
(24)

where the summation is over ions and partial waves, with the prime indicating
explicitly that the summation is always truncated to a small number of partial
waves in real calculations. For liquids, the basis states in both the AE space
and the PS space are the plane waves. A matrix element in the plane wave basis
is

$\langle\psi_{\bm{k}}|\hat{T}|\tilde{\psi}_{\bm{k}^{\prime}}\rangle=\delta_{\bm{k}\bm{k}^{\prime}}\\\
+{\sum_{i,a}}^{\prime}\left(\langle\psi_{\bm{k}}|\phi_{i}^{a}\rangle-\langle\psi_{\bm{k}}|\tilde{\phi}_{i}^{a}\rangle\right)\langle\tilde{p}_{i}^{a}|\tilde{\psi}_{\bm{k}^{\prime}}\rangle,$
(25)

where $|\psi_{\bm{k}}\rangle$ and $|\tilde{\psi}_{\bm{k}^{\prime}}\rangle$ are
the plane waves in the AE and PS space, respectively. The overlap matrix
element between a plane wave and a partial wave can be evaluated by

$\langle\psi_{\bm{k}}|\phi_{i}^{a}\rangle=\frac{4\pi}{\sqrt{V}}\mathrm{e}^{-\mathrm{i}\bm{k}\cdot\bm{R}_{a}}Y_{l_{i}m_{i}}(\hat{\bm{k}})\\\
\times\int_{0}^{r_{c}}r^{2}\mathrm{d}r(-\mathrm{i})^{l_{i}}j_{l_{i}}(kr)\phi_{l_{i},m_{i}}^{a}(r),$
(26)

and similarly for other matrix elements, where $Y_{lm}$ denotes the spherical
harmonics, $j_{l}(kr)$ is the spherical Bessel function, $r_{c}$ denotes the
cutoff radius of the augmentation sphere in defining $\hat{T}^{a}$, and
$\bm{R}_{a}$ is the position of the ion $a$.

The overlap matrix $\hat{S}=\hat{T}^{\dagger}\hat{T}$ can also be determined
by using the partial waves. From Eq. (24), we have

$\hat{S}=-\hat{\mathbb{I}}+\hat{T}+\hat{T}^{\dagger}\\\
+{\sum_{ij,a}}^{\prime}\Ket{\tilde{p}_{i}^{a}}\Braket{\phi_{i}^{a}-\tilde{\phi}_{i}^{a}}{\phi_{j}^{a}-\tilde{\phi}_{j}^{a}}\Bra{\tilde{p}_{j}^{a}},$
(27)

where we assume that augmentation spheres associating with different ions do
not overlap. We can rewrite it as

$\hat{S}=\hat{S}^{\prime}+\left(\hat{T}-\hat{T}^{\prime}+\mathrm{h.c.}\right),$
(28)

where

$\displaystyle\hat{S}^{\prime}$ $\displaystyle=$
$\displaystyle\hat{\mathbb{I}}+{\sum_{ij,a}}^{\prime}|\tilde{p}_{i}^{a}\rangle\left(\langle\phi_{i}^{a}|\phi_{j}^{a}\rangle-\langle\tilde{\phi}_{i}^{a}|\tilde{\phi}_{j}^{a}\rangle\right)\langle\tilde{p}_{j}^{a}|,$
(29) $\displaystyle\hat{T}^{\prime}$ $\displaystyle=$
$\displaystyle\hat{\mathbb{I}}+{\sum_{ij,a}}^{\prime}|\tilde{p}_{j}^{a}\rangle\left(\langle\tilde{\phi}_{j}^{a}|\phi_{i}^{a}\rangle-\langle\tilde{\phi}_{j}^{a}|\tilde{\phi}_{i}^{a}\rangle\right)\langle\tilde{p}_{i}^{a}|.$
(30)

We note that $\hat{S}^{\prime}$ is the form commonly adopted in literatures.
Actually, by using the closure relation Eq. (15), it is easy to see that
$\hat{T}^{\prime}$ is equal to $\hat{T}$, and therefore
$\hat{S}=\hat{S}^{\prime}$ if the summation is over a complete set of partial
waves. However, since real calculations always use just a small set of partial
waves instead of a complete set, the correction due to
$\hat{T}-\hat{T}^{\prime}$ could be non-negligible. For the determination of
the AE Green’s function and the $\mathcal{T}$-matrix, it is important to
maintain the consistency between $\hat{T}$ and $\hat{S}$ since they both
appear in Eq. (21). For the reason, we use Eq. (28) instead of Eq. (29) for
calculating the overlap matrix.

##  IV Implementation and results

###  IV.1 Implementation

Figure 1:  Flowchart for calculating Green’s function and
$\mathcal{T}$-matrix. It substitutes the corresponding part of the workflow
shown in Fig. 1 of Ref. [1].

We employ the GPU version of VASP [34, 35] for calculating electron
structures. The radius of the augmentation spheres is set to $0.52$ Å, which
may result in overlaps of the spheres in some ionic configurations but saves
computation time. We test a smaller radius of $0.43$ Å and find that
differences in calculation results are negligible. The Perdew-Burke-Ernzerhof
(PBE) exchange-correlation functional [36] is used in all calculations. A
Monkhorst-Pack k-point mesh of spacing not larger than $2\pi\times 0.05$ Å-1
is used to sample the Brillouin zone, and an energy cutoff of $600$ eV is
employed to expand the wavefunctions, yielding a convergence in energy better
than $\sim 2\,\mathrm{meV/atom}$.

To determine the Green’s function and the scattering $\mathcal{T}$-matrix by
using Eqs. (21, 22), we need to obtain the pseudo Hamiltonian
$\hat{\tilde{H}}$, the transformation matrix $\hat{T}$ and the overlap matrix
$\hat{S}$ from the electron structure calculations. Since VASP does not
calculate or output these quantities directly, we modify it to output
intermediate quantities including the local part of the effective KS potential
$\tilde{v}_{\mathrm{loc}}$, the coefficients of the non-local part of the
effective KS potential $(\tilde{v}^{\prime})_{ij}^{a}$ as well as the
projector matrix elements
$\langle\tilde{\psi}_{\bm{k}}|\tilde{p}_{i}^{a}\rangle$. With them, the matrix
elements
$\langle\tilde{\psi}_{\bm{k}}|\hat{\tilde{H}}|\tilde{\psi}_{\bm{k}^{\prime}}\rangle$
of the pseudo Hamiltonian can be determined by using the definition

$\displaystyle\hat{\tilde{H}}=\hat{K}+\tilde{v}_{\mathrm{loc}}+\sum_{ij,a}|\tilde{p}_{i}^{a}\rangle(\tilde{v}^{\prime})_{ij}^{a}\langle\tilde{p}_{j}^{a}|.$
(31)

The matrix elements of the transformation operators $\hat{T}$,
$\hat{T}^{\prime}$ and the overlap matrix $\hat{S}^{\prime}$ can also be
determined by applying Eq. (24), Eq. (30), and Eq. (29), respectively. It is
then straightforward to determine $\hat{S}$ by applying Eq. (28).

We can then proceed to carry out the analysis outlined in Sec. II by following
the procedure detailed in Ref. 1. Figure 1 shows the flowchart of our
implementation.

###  IV.2 Results

Figure 2:  Comparison between the full DFT-based implementation (PAW method)
and the LSA in (a) pair scattering amplitude $\Gamma_{m}(q)$, (b) effective
pair interaction $W_{m}(q)$ and (c) interaction parameters $\lambda(n)$. The
results are for $T=450$ K and $P=500$ GPa. $\Gamma_{m}(q)$ and $W_{m}(q)$ are
obtained by recasting the matrix elements of $\Gamma$ and $W$ as a function of
$q=|\bm{k}_{1}-\bm{k}_{2}|$ for $\bm{k}_{1}$, $\bm{k}_{2}$ in
$0.8k_{F}\leq|\bm{k}_{1}|,|\bm{k}_{2}|\leq 1.14k_{F}$.

To test our first-principles implementation based on PAW method, we apply it
to metallic hydrogen liquids. The ion configurations from the constant-volume
PIMD simulations of 200 hydrogen atoms in Ref. [1] are re-used in the current
study. In this subsection, we summarize our calculation results.

Estimated $T_{c}$’s as well as relevant parameters of metallic hydrogen
liquids under pressures from $0.5$ to $1.5$ TPa are summarized in Table 1. For
comparison, the LSA results are also shown. We find that the results coincide
well with the LSA results. It confirms Liu et al.’s prediction that the
superconducting transition temperature of metallic hydrogen is higher than its
predicted melting temperature [1]. We note that our prediction of $T_{c}$ is
only valid for the liquid phase. $T_{c}$ for the solid phase of metallic
hydrogen had been thoroughly studied in Ref. [7].

Figure 3:  Dependence of $T_{c}$ on (a) different numbers of atoms, and (b)
different numbers of beads for $P=700$ GPa. The error bars indicate the
numerical uncertainties of the estimations of $T_{c}$’s, i.e., the values
shown in the parentheses in Table 1.

In Fig. 2, we show the comparison of the pair scattering amplitude
$\Gamma_{m}(q)$ and the effective pairing interaction $W_{m}(q,\nu)$ between
the two implementations for $P=700$ GPa and $T=450$ K. We find that the
results yielded by the two implementations are very close. We also show the
interaction parameters $\lambda(n)$ in Fig. 2(c). To determine $\lambda(n)$,
we recast the matrix element $W_{11^{\prime}}$ as a function of
$q=|\bm{k}-\bm{k}^{\prime}|$ for $0.8k_{F}\leq|\bm{k}|,|\bm{k}^{\prime}|\leq
1.14k_{F}$, and evaluate Eq. (3) by integrating over $q$ [5]. It is evident
that the results also coincide well. It is not surprising that the LSA can
work well in metallic hydrogen systems because hydrogen has no inner-core
electron and the effective KS potential is dominated by its local part which
can be well approximated by the LSA with an appropriate dielectric function.
For more general systems like hydrides, the LSA is obviously inadequate, and
our implementation based on the standard first-principles approach will be
useful.

We test the convergence of our calculation by varying the size of the
simulation supercell and the number of beads discretizing the imaginary time.
In Ref. 1 and the calculation shown above, we use a supercell with 200
hydrogen atoms and discretize the imaginary time to 24 beads. In Fig. 3(a), we
show the dependence of estimated $T_{c}$ on the number of atoms. We find that
$T_{c}$ slightly decreases with the increasing simulation size, but the
difference are within the error bars. In Fig. 3(b), we show the dependence of
estimated $T_{c}$ on the number of beads. We find that only a relatively small
number of beads ($>16$) is needed to get a converged result. The oversampling
procedure discussed in Sec. II.2 is effective in eliminating the
discretization errors.

| $r_{s}$ | $1.226$ | $1.197$ | $1.17$ | $1.149$ | $1.113$ | $1.049$
---|---|---|---|---|---|---|---
| P | $0.5$ | $0.6$ | $0.7$ | $0.8$ | $1.0$ | $1.5$
PAW method |  350K | $\lambda$ | $10.2(15)$ | $9.0(10)$ | $8.5(10)$ | $7.5(9)$ | $6.3(8)$ | $5.2(5)$
$\overline{\omega}_{2}$ | $102.1(10)$ | $110.9(10)$ | $112.2(15)$ | $124.6(7)$ | $134.6(3)$ | $164.1(30)$
$\rho_{m}$ | $0.37(12)$ | $0.36(11)$ | $0.32(10)$ | $0.31(9)$ | $0.32(11)$ | $0.26(9)$
450K | $\lambda$ | $8.2(16)$ | $7.7(10)$ | $7.1(8)$ | $6.3(10)$ | $5.5(6)$ | $4.6(4)$
$\overline{\omega}_{2}$ | $115.0(1)$ | $114.9(13)$ | $132.0(7)$ | $141.9(4)$ | $147.3(11)$ | $173.1(12)$
$\rho_{m}$ | $-0.10(12)$ | $-0.13(9)$ | $-0.13(8)$ | $-0.10(10)$ | $-0.14(7)$ | $-0.14(6)$
$T_{c}$ | $428(25)$ | $423(19)$ | $421(19)$ | $425(24)$ | $420(19)$ | $415(17)$
Linear screening |  350K | $\lambda$ | $9.5(14)$ | $8.6(11)$ | $8.1(10)$ | $7.2(10)$ | $6.0(8)$ | $5.0(5)$
$\overline{\omega}_{2}$ | $106.7(10)$ | $113.8(1)$ | $115.8(7)$ | $127.5(16)$ | $138.0(1)$ | $167.5(30)$
$\rho_{m}$ | $0.38(12)$ | $0.36(10)$ | $0.32(10)$ | $0.30(8)$ | $0.29(11)$ | $0.24(9)$
450K | $\lambda$ | $7.6(13)$ | $7.4(10)$ | $6.6(8)$ | $6.0(10)$ | $5.2(6)$ | $4.4(4)$
$\overline{\omega}_{2}$ | $119.9(11)$ | $117.8(4)$ | $135.8(3)$ | $144.6(19)$ | $157.1(6)$ | $176.3(12)$
$\rho_{m}$ | $-0.11(11)$ | $-0.13(9)$ | $-0.14(8)$ | $-0.11(9)$ | $-0.14(7)$ | $-0.15(6)$
$T_{c}$ | $428(24)$ | $423(19)$ | $420(19)$ | $424(22)$ | $418(19)$ | $411(17)$
Table 1: Mass enhancement factor $\lambda\equiv\lambda(0)$, average phonon
frequency $\bar{\omega}_{2}$ (in meV) and the maximal eigenvalue $\rho_{m}$ of
the linearized Eliashberg equations for different pressures $P$ (in TPa) and
temperatures $T$ (in K), determined by the standard first-principles approach
(PAW method) and the LSA 222We recalculate the LSA results. They are slightly
different from those presented in Ref. [1] but within the error bars. The
numerical differences are due to different settings in analyzing the PIMD
data.. The transition temperature $T_{c}$, i.e., the temperature at which
$\rho_{m}$ becomes zero, is estimated based on a linear interpolation of
$\rho_{m}$ between the two temperatures. Numerical uncertainties of the
quantities are indicated in parentheses.

##  V Summary

In summary, we present a full DFT-based implementation of the stochastic path-
integral approach for estimating superconducting transition temperatures. We
derive the formulae for determining the Green’s function and the scattering
$\mathcal{T}$-matrix in the pseudo-basis of the PAW method. Properly
calculating the $S$ matrix elements to eliminate the error due to the
incompleteness of the basis is also discussed. We test our implementation in
metallic hydrogen liquid systems and get results close to the those obtained
in Ref. 1. The calculation is in agreement with Liu et al’s prediction, and
suggests that metallic hydrogen could be a superconducting liquid as its
superconducting $T_{c}$ is higher than the predicted melting temperature of
metallic hydrogen. With the full DFT-based implementation, it becomes possible
to apply the stochastic path-integral approach to a wider class of systems
such as hydrides.

This work is supported by the National Basic Research Program of China (973
Program) Grants No. 2018YFA0305603 and No. 2015CB921101 and the National
Science Foundation of China Grant No. 11325416.

## References

  * Liu  _et al._ [2020] H. Liu, Y. Yuan, D. Liu, X.-Z. Li, and J. Shi, Phys. Rev. Research 2, 013340 (2020).
  * Duan  _et al._ [2014] D. Duan, Y. Liu, F. Tian, D. Li, X. Huang, Z. Zhao, H. Yu, B. Liu, W. Tian, and T. Cui, Sci. Rep. 4, 1 (2014), publisher: Nature Publishing Group.
  * Liu  _et al._ [2017] H. Liu, I. I. Naumov, R. Hoffmann, N. Ashcroft, and R. J. Hemley, Proceedings of the National Academy of Sciences 114, 6990 (2017), publisher: National Acad Sciences.
  * Ashcroft [1968] N. W. Ashcroft, Phys. Rev. Lett. 21, 1748 (1968), publisher: American Physical Society.
  * Jaffe and Ashcroft [1981] J. E. Jaffe and N. W. Ashcroft, Phys. Rev. B 23, 6176 (1981).
  * Richardson and Ashcroft [1997] C. F. Richardson and N. W. Ashcroft, Phys. Rev. Lett. 78, 118 (1997).
  * McMahon and Ceperley [2011] J. M. McMahon and D. M. Ceperley, Phys. Rev. B 84, 144515 (2011).
  * Kim  _et al._ [2011] D. Y. Kim, R. H. Scheicher, C. J. Pickard, R. J. Needs, and R. Ahuja, Phys. Rev. Lett. 107, 117002 (2011), publisher: American Physical Society.
  * Flores-Livas  _et al._ [2016] J. A. Flores-Livas, M. Amsler, C. Heil, A. Sanna, L. Boeri, G. Profeta, C. Wolverton, S. Goedecker, and E. K. U. Gross, Phys. Rev. B 93, 020508 (2016), publisher: American Physical Society.
  * Sun  _et al._ [2019] Y. Sun, J. Lv, Y. Xie, H. Liu, and Y. Ma, Phys. Rev. Lett. 123, 097001 (2019).
  * Grimvall [1981] G. Grimvall,  _The electron-phonon interaction in metals_ (North Holland, 1981).
  * Giustino [2017] F. Giustino, Rev. Mod. Phys. 89, 015003 (2017).
  * Errea  _et al._ [2015] I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Phys. Rev. Lett. 114, 157004 (2015).
  * Liu  _et al._ [2018] H. Liu, I. I. Naumov, Z. M. Geballe, M. Somayazulu, J. S. Tse, and R. J. Hemley, Phys. Rev. B 98, 100102 (2018).
  * Chen  _et al._ [2013] J. Chen, X.-Z. Li, Q. Zhang, M. I. J. Probert, C. J. Pickard, R. J. Needs, A. Michaelides, and E. Wang, Nature Communications 4, 2064 (2013).
  * Geng  _et al._ [2015] H. Y. Geng, R. Hoffmann, and Q. Wu, Phys. Rev. B 92, 104103 (2015).
  * Errea  _et al._ [2014] I. Errea, M. Calandra, and F. Mauri, Phys. Rev. B 89, 064302 (2014).
  * Errea  _et al._ [2020] I. Errea, F. Belli, L. Monacelli, A. Sanna, T. Koretsune, T. Tadano, R. Bianco, M. Calandra, R. Arita, F. Mauri, and J. A. Flores-Livas, Nature 578, 66 (2020).
  * Marx and Parrinello [1996] D. Marx and M. Parrinello, Journal of Chemical Physics 104, 4077 (1996).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Comput. Mat. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  * Blöchl [1994] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Morel and Anderson [1962] P. Morel and P. W. Anderson, Phys. Rev. 125, 1263 (1962).
  * Knorr and Barth [1970] K. Knorr and N. Barth, Solid State Communications 8, 1085 (1970).
  * Bergmann [1971] G. Bergmann, Phys. Rev. B 3, 3797 (1971).
  * Gross and Kurth [1991] E. K. U. Gross and S. Kurth, International Journal of Quantum Chemistry 40, 289 (1991).
  * Lüders  _et al._ [2005] M. Lüders, M. A. L. Marques, N. N. Lathiotakis, A. Floris, G. Profeta, L. Fast, A. Continenza, S. Massidda, and E. K. U. Gross, Phys. Rev. B 72, 024545 (2005).
  * Marques  _et al._ [2005] M. A. L. Marques, M. Lüders, N. N. Lathiotakis, G. Profeta, A. Floris, L. Fast, A. Continenza, E. K. U. Gross, and S. Massidda, Phys. Rev. B 72, 024546 (2005).
  * Chandler and Wolynes [1981] D. Chandler and P. G. Wolynes, The Journal of Chemical Physics 74, 4078 (1981).
  * Negele and Orland [1988] J. W. Negele and H. Orland,  _Quantum many-particle systems_ (Addison-Wesley, 1988).
  * Giuliani and Vignale [2005] G. Giuliani and G. Vignale,  _Quantum theory of the electron liquid_ (Cambridge University Press, 2005).
  * Ichimaru and Utsumi [1981] S. Ichimaru and K. Utsumi, Phys. Rev. B 24, 7385 (1981).
  * Hacene  _et al._ [2012] M. Hacene, A. Anciaux-Sedrakian, X. Rozanska, D. Klahr, T. Guignon, and P. Fleurat-Lessard, Journal of Computational Chemistry 33, 2581 (2012).
  * Hutchinson and Widom [2012] M. Hutchinson and M. Widom, Computer Physics Communications 183, 1422 (2012).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996), publisher: American Physical Society.


# Spin-orbit coupling tuned crossover of gaped and gapless topological phases
in the chalcopyrite HgSnX2 (X=N/P): An ab-initio investigation

Surasree Sadhukhan  Corresponding author<EMAIL_ADDRESS><EMAIL_ADDRESS>Sudipta Kanungo<EMAIL_ADDRESS>School of
Physical Sciences, Indian Institute of Technology Goa, Goa-403401, India

###### Abstract

The coupling between electron orbital momentum and spin momentum, known as
spin-orbit coupling (SOC), is a fundamental origin of a multitude of
fascinating physical phenomena, especially it holds paramount significance in
the realm of topological materials. In our work, we have predicted the
topological phase in Hg-based chalcopyrite compounds using the first
principles density functional theory. The initial focus was on HgSnN2,
revealing it to be a nonmagnetic Weyl semimetal, while HgSnP2 displayed
characteristics of a strong topological insulator. What makes our work truly
unique is that despite both compounds having the same SOC strength, arises
from Hg, they exhibit distinct topological phases due to the distinct
hybridization effect of the Hg-$5d$ and X-$p$ bands. This finding can address
a significant factor, i.e., the effect of the band hybridization in deriving
distinct topological phases, keeping the symmetry aspect intact. Our results
indicate that due to the presence of band hybridization between the dominant
X-$np$ orbitals n=2 and 3 for X=N and P respectively and a minor contribution
from Hg-${5d}$, we can tune the topological phase by manipulating SOC
strength, which equivalently achievable by chemical substitutions. This
investigation stands as a remarkable illustration of the unique roles that
hybridization plays in sculpting the topological properties of these compounds
while simultaneously preserving their underlying symmetries.

##  I Introduction

The coupling between spin and orbital angular momentum in the relativistic
framework, known as spin-orbit coupling (SOC), is closely intertwined with the
emergence of novel phases in quantum materials[1, 2, 3, 4]. This fascinating
domain of SOC has broadened to encompass a plethora of phenomenon like Kitaev
spin liquid systems[5, 6, 7], unconventional magnetism in heavy elements[8],
superconductivity[9, 10], Rashba splitting[11, 12, 13], chiral magnetic
interaction[14, 15, 16, 17], spin-orbit torque (SOT)[18, 19, 20, 21, 22, 23,
24, 25, 26], spin Hall effect[27, 28], spin-orbit qubits[29], topological
phases of matter[30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41] and many
more. The importance of SOC in relation to topological phases is thoroughly
recognized and meticulously documented, especially the incorporation of SOC
introduces a band gap and leads to band inversion[42, 43, 44, 45, 46] and even
SOC notably plays a significant role in instigating the breakage of band
degeneracy[41]. The strength of SOC is an atomic property that is directly
proportional to the fourth power of atomic number; nonetheless, some possible
approaches, such as doping Pb1-xSnxSe[48] and designing heterostructures[49,
50], have been highly appreciated to engineer the effects of SOC and utilize
in topological-based spintronics devices[3, 51]. Furthermore, both the
magnitude and sign of the SOC strength can be manipulated as experimentally
demonstrated in thallium nitride (TlN)[47]. However, the precise and
definitive connection between topological phases and SOC at a material-
specific level remains a subject yet to be fully explored. Especially in the
context of inter-metallic compounds, the hybridization between bands
originating from distinct chemical elements can effectively account for the
influence of SOC. In this scenario, the identification of the topological
phase is not exclusively dictated by the particular element demonstrating high
SOC. Instead, the compelling impact of hybridization with other elements, even
those with lower SOC strength, comes into play, influencing the determination
of the topological signatures.

To understand the direct impact of SOC in driving the topological phases, we
have selected two sister compounds, HgSnN2 and HgSnP2, having the same
crystalline symmetry, valence electrons, and heavy element (Hg), which
accounts for the dominance of SOC strength. The rationale behind such a choice
of compounds is to ensure everything that determines the topology of the
materials is the same except for the choice of the anion(N/P). Here it is
crucial for meaningful comparisons in our investigation. Moreover, the study
of SOC effects in topological phases of Hg-based heavy elements has found
promising platforms in several cases, including HgTe quantum well, which has
similarity with chalcopyrite[52], MTl4Te3 (M = Cd, Hg)[53], Jacutingaite-
family [M2NX3 (M = Ni, Pt, and Pd; N = Zn, Cd, and Hg; and X = S, Se, and
Te)][54], two dimensional M3C2 (M = Zn, Cd, and Hg)[55], Hg-based materials
HgS, HgSe, and HgTe[56] etc., from both theoretical and experimental
standpoints. Additionally, they provide insights that can aid in predicting
potential Hg-based topological materials in our case. Moreover, it is also
crucial to consider the experimental feasibility of compounds for a sensible
choice of quantum material. The sibling compounds based on N and P in the
chalcopyrite class have already been synthesized successfully, including
ZnSiP2[57], ZnGeP2[58, 59], ZnSnP2[60], CdSiP2[61, 62], CdGeP2[63], CdSnP2[64,
65], MgSiP2[66], ZnSnN2[67], II-Sn-N2 (II= Ca, Mg, Zn)[68], MgSnN2[69],
MgZrN2[70], MnSiN2[71] and many more. While our specific compound hasn’t been
experimentally synthesized yet, there is a significant likelihood of achieving
experimental synthesis for these compounds in the future. It is worth
highlighting the added advantages of small gap semiconductors with broken
inversion symmetry have the potential to serve as promising candidates for
various applications, including photovoltaic cells, optoelectronics (such as
light-emitting diodes and lasers), thermoelectricity, and non-linear optical
properties[72, 73, 74]. Additionally, the ability to engineer a wide range of
band gaps offers opportunities for developing band gap tuning devices driven
by SOC.

In our investigation, we have followed the ab-initio first principles method
to identify two previously unexplored topological materials, namely HgSnN2 and
HgSnP2. To the best of our knowledge, no previous investigations have been
conducted to examine the topological properties of these compounds. Our
findings reveal that HgSnN2 identified as a nonmagnetic Weyl semimetal,
showcasing essential topological features such as band inversion, Fermi arc,
chirality, etc., in the presence of SOC. Similarly, the HgSnP2 shows the
characteristics of a strong topological insulator with a Z2 number of [1,000]
and an odd number of Dirac cones in its surface spectral function when SOC is
incorporated. The distinctive aspect of our research lies in the fascinating
observation that despite the atomic SOC strength being identical for both
compounds, they showcase distinctly different topological phases: one is
gapless, and another is gapped. Such disparity of nontrivial phase arises from
the band hybridization of $N$/$P$ elements with $Hg$, resulting in their
unique behavior of gapless or gaped band structure. To explore the tangible
impact of SOC on hybridized bands, we have artificially modulated the strength
of SOC in the original Hamiltonian. Intriguingly, our investigations revealed
that by escalating the strength of SOC in HgSnN2, we could induce a gap in the
Weyl crossing, leading to a topological phase similar to that of the parent
compound HgSnP2. In a parallel manner, when we decrease the SOC strength in
HgSnP2, we observed the band gap closing, resulting in a Weyl crossing akin to
the parent HgSnN2 compound. This comparative analysis sheds light on the
distinct roles played by hybridization in shaping their topological
properties. The importance of our present work lies in the fact that, as per
the usual notion, the crossover from gap less to gaped topological phases can
be achieved by tuning the governing symmetries; however, in the present case,
we demonstrate how tuning SOC and hybridization can lead to the gap less to
gap crossover in the non-trivial topological phases.

##  II Crystal structure

Figure 1: (Color online)(a) The crystal structure of HgSnX2 [X=N, P] is depicted, where Hg, Sn, and X atoms are represented by blue, green, and red spheres, respectively. We utilized the VESTA package with a ball-and-stick representation to visualize the three-dimensional crystal structure[75].(b) and (c) display the calculated orbital projected density of states (DOS) of HgSnX2, where X=N and P, respectively. The green, blue, and red lines represent the partial DOS for Hg-${5d}$, Sn-${5p}$, and X-${np}$ (n=2 and 3 for X=N and P respectively), respectively. The energy scale is referenced to the Fermi energy (Ef), which is set at zero. Material |  a ($\AA$) |  c ($\AA$) | Ref.
---|---|---|---
HgSnN2 | 5.19 | 9.90 | our calculations
| 6.46 | - | [83]
| 4.96 | 9.65 | [84]
HgSnP2 | 5.98 | 11.77 | our calculations
| 6.17 | - | [83]
| 5.87 | 11.47 | [84]
| 5.95 | 11.83 | [85]
| 5.91 | - | [86]
Table 1: A showcasing of the optimized lattice parameters a and c for both HgSnN2 and HgSnP2. The data provided includes our calculations and the corresponding values found in the literature, with appropriate references. HgSnN2 | Component of | HgSnP2 | Component of
---|---|---|---
(0.016 THz) | displacement vector($\AA$) | (0.010 THz) | displacement vector ($\AA$)
Hg | ( -0.11, 0.15, 0.00) | Hg | ( -0.01, -0.01 0.18 )
Sn | ( -0.09, 0.12, 0.00) | Sn | (-0.01,-0.00, 0.14)
N | (-0.03, 0.04, 0.00 ) | P | (-0.00, -0.00, 0.07)
Table 2: The component of displacement vectors corresponding to Hg, Sn, and X (where X = N or P) are presented in the acoustic mode with the lowest frequency at the $\Gamma$ point. HgSnN2 | Component of | HgSnP2 | Component of
---|---|---|---
(18.97 THz) | displacement vector($\AA$) | (9.71 THz) | displacement vector ($\AA$)
Hg | (0.00, 0.00, 0.00) | Hg | (0.01, 0.01, 0.00)
Sn | ( 0.00, 0.00, 0.00) | Sn | (-0.07, -0.05, -0.00)
N | ( 0.15, 0.00, -0.09) | P | (0.12, -0.05 -0.11)
| (0.05, 0.15, 0.09) |  | (-0.01, 0.13, 0.11)
| (-0.05, -0.15, 0.09) |  | ( -0.01 0.13 -0.11)
| ( -0.15, -0.00 -0.09) |  | ( 0.12, -0.05 0.12)
Table 3: The component of displacement vectors corresponding to Hg, Sn, and X
(where X = N or P) are presented in the optical mode with the highest
frequency at the $\Gamma$ point.

HgSnX2 (X=N, P) is a non-centrosymmetric intermetallic alloy belonging to the
chalcopyrite class with space group I$\bar{4}$2d (space group no. 122). The
unit cell possesses the body-centered tetragonal structure having a
tetrahedral unit with the cation Hg, Sn, and the anion X (X=N, P) as depicted
in Fig. 1(a). The optimized lattice parameters are listed in Table I. The
Wyckoff positions of the three atoms are assigned the position as Hg atom at
(0, 0, 0), Sn atom at (0, 0, 0.5), and X atom at ($x$, 0.25, 0.125), where $x$
referred as the anion displacement parameter. For HgSnN2 and HgSnP2, the
values of $x$ are 0.7921 and 0.7548, respectively. We have observed slight
deviations in the lattice parameters of both compounds in the literature
compared to our optimized values, which we listed in Table I. Therefore, we
have examined the lattice parameters versus energy plot very carefully and
confirmed that the lattice parameters listed in Table I are energetically
lower compared to any other choice of lattice parameters reported so far, as
demonstrated in Fig. S1(a) and Fig. S1(b) for HgSnN2. Similar analyses have
been conducted for HgSnP2, as shown in Fig. S1(d) and S1(e). The calculated
formation energies for HgSnN2 and HgSnP2 are -1.6374 eV/f.u. and -6.91705
eV/f.u., respectively, which show that the HgSnP2 is more stable than HgSnN2.
Since this compound has yet to be synthesized experimentally, the choice of
lattice parameters is critical for the chalcopyrites class, especially in
addressing the topological properties[87] and exploring its dynamical
stability and thermodynamic characteristics is of utmost importance. We have
calculated the phonon density of states and confirmed the absence of any
negative frequency phonon DOS. The negative Helmholtz free energy as a
function of temperature and positive entropy, depicted in the inset of Fig.
S1(c), further supports the thermodynamic stability of HgSnN2. Similarly, Fig.
S1(f) illustrates the phonon density of states without any negative frequency
region and the temperature-dependent positive evolution of Helmholtz free
energy for HgSnP2. These calculations assure the stability of our structure
from a thermodynamic standpoint.

Additionally, we have depicted the phonon dispersion of HgSnN2 and HgSnP2 in
Fig. S2(a) and Fig. S2(b), respectively, with two distinct modes—acoustic and
optical—i.e., respective the lowest and the highest frequencies observed at
the $\Gamma$ point in the phonon dispersion. To enhance comprehension of both
optical and acoustic modes, we have provided a visual representation of
displacement vectors within the primitive cell arrangement. In Fig. S2(a), the
acoustic mode registers a frequency of 0.016 THz, while the optical mode is
observed at 18.97 THz. Interestingly, no phonon bands are within the frequency
range of 6.48 THz to 8.80 THz. Examining the acoustic mode of HgSnN2, we note
that mainly cations (Hg and Sn atoms) exhibit prominent displacement vectors
along the negative X and positive Y directions, with no displacement in the
Z-direction. The specific components of displacement vectors associated with
each atom are detailed in Table 2, and the corresponding displacement vector
is illustrated in Fig. S2(c). Conversely, in the optical mode, cations (Hg and
Sn atoms) display no displacement vector. Within this mode, anions (N)
predominantly contribute to the vibrations. The detailed components of
displacement vectors are listed in Table 3, and the corresponding displacement
vector is depicted in Fig. S2(d).

In the phonon dispersion of HgSnP2, illustrated in Fig. S2(b), the acoustic
mode manifests at a frequency of 0.010 THz, while the optical mode is
discerned at 9.71 THz. Notably, there are no phonon bands within the frequency
range of 4.52 THz to 7.20 THz. Upon closer examination of the acoustic mode in
HgSnP2, it is observed that cations (Hg and Sn atoms) prominently exhibit
displacement vectors along the Z-direction with tiny displacement of anions
(P). The specific components of these displacement vectors are comprehensively
outlined in Table 2, with the corresponding visual representation of the
displacement vector provided in Fig. S2(e). Conversely, in the optical mode,
cations (mainly Sn atoms) showcase relatively less displacement vector
compared to anions (P). Within this mode, anions (P) play a predominant role
in contributing to the vibrations dynamics. The intricate details of the
displacement vectors are cataloged in Table 3, and the corresponding
displacement vector is portrayed in Fig. S2(f). The notable resemblance
between the two sibling compounds lies in the dominance of cations (Hg and Sn)
contributions in the acoustic mode. Also, in the optical mode, the
participation is primarily from anions, with X (where X = N/P) playing a major
role. A significant disparity is observed in the displacement vectors in the
acoustic mode for HgSnN2 Hg, and Sn shows displacement along X and Y
directions, but for HgSnP2, it is in the Z-direction. Another crucial
observation is Sn has taken part in optical mode only for HgSnP2. We have also
determined the elastic stiffness constants by assessing the energy variation
by applying minor strains to the equilibrium lattice configuration to ensure
the elastic stability of HgSnX2 (X=N/P). The count of elastic stiffness
constants increases as the symmetry of the crystal structure is reduced. In
the case of our system, which belongs to the tetragonal group, there are six
independent elastic stiffness constants namely C11, C12, C13, C33, C44,
C66[88, 89]. The value of the elastic constant is tabulated in Table 4. A
material is considered elastically stable when the elastic constants meet the
criteria (1) to (4) for a tetragonal structure with the space group 122. We
have found that both compounds satisfy the four necessary and sufficient
conditions, which ensures the mechanical stability of these two compounds.

${\color[rgb]{0,0,0}C_{11}>|C_{12}|}$ (1)
${\color[rgb]{0,0,0}(C_{11}+C_{12})C_{33}>2(C_{13})^{2}}$ (2)
${\color[rgb]{0,0,0}C_{44}>0.0}$ (3) ${\color[rgb]{0,0,0}C_{66}>0.0}$ (4)

Material | C11 (GPa) | C12 (GPa) | C13 (GPa) | C33 (GPa) | C44 (GPa) | C66 (GPa)
---|---|---|---|---|---|---
HgSnN2 | 206.294 | 74.055 | 64.916 | 223.207 | 101.362 | 48.705
HgSnP2 | 338.194 | -26.859 | -54.115 | 124.782 | 57.694 | 30.692
Table 4: The six independent elastic stiffness constants for HgSnN2 and HgSnP2
shows the elastic stability of the materials.

##  III Calculation methodology

We conducted electronic structure calculations based on Density Functional
Theory (DFT) using the Vienna ab-initio simulation package (VASP) with a
plane-wave basis set and pseudo-potential framework[76, 77]. The calculations
employed the Perdew–Burke–Ernzerhof (PBE) prescription, which is a generalized
gradient exchange-correlation approximation (GGA) functional[78]. The
influence of SOC was considered by applying a relativistic correction to the
original Hamiltonian. A cut-off energy of 600 eV was used for the plane-wave
basis. We relaxed the internal atomic positions to optimize the structure
until the Hellman–Feynman force was less than 0.001 eV/$\AA$. The Brillouin
zone was sampled using an 8 $\times$ 8 $\times$ 6 k-point mesh, and the
electronic convergence criteria were set at 10-9 eV for self-consistent
calculations. We constructed a tight-binding model using the Wannier function
basis set to study the topological properties. This energy-selective method
mainly addresses the low-energy, few-band model Hamiltonian defined in the
effective Wannier function basis. Using a re-normalization method, we
integrated out the irrelevant degrees of freedom to keep the orbitals that are
solely responsible for topological effects. In the down-folding calculations,
we considered only the Hg-${5d}$, Sn-${5p}$, X-${np}$, X-${ns}$ (n=2,3 for
X=N/P) orbitals as active degrees of freedom for HgSnN2 and HgSnP2
respectively, while other degrees of freedom were down-folded as they made
less significant contributions in determining the topological properties. To
analyze the topological aspects, we calculated the topological index (Z2
number), examined the surface states, studied the distribution of normalized
Berry curvature in the 2D momentum plane, and investigated the Fermi arc at
the particular chemical potential of Weyl crossing. These calculations were
performed using Wannier90[79, 80, 81] and WannierTool[82], starting from the
full converged self-consistent DFT outputs.

##  IV Electronic structure

Figure 2: (Color online) (a) and (c) depict the calculated bulk electronic
band structures for HgSnN2 and HgSnP2, respectively, employing the GGA (blue
lines) and GGA+SOC (red lines). The inset images in (a) and (c) provide a zoom
version of the gap-less and gaped states in the GGA+SOC band structures for
HgSnN2 and HgSnP2, respectively. (b) and (d) illustrate the calculated orbital
projected band structures with GGA+SOC for the Hg and X (X=N, P) atoms,
showing the band inversion. The Brillouin zone (BZ) is shown in the lower
inset of (c). The energy scale is referenced to the Fermi energy(Ef), which is
set at zero.

The GGA density of states (DOS) for HgSnX2 (X=N, P) in the nonmagnetic state
has been successfully determined using first-principles density functional
theory, shown in Fig. 1(b)-(c). By analyzing the DOS of HgSnN2 in Fig.1(b), it
becomes evident that it exhibits a semi-metallic nature due to the absence of
a noticeable gap in DOS at the Fermi energy, primarily influenced by the
significant $p$-orbital contribution of N. However, upon replacing N with P,
the contribution from the $p$-orbital of P shifts slightly away from the Fermi
energy, as depicted in Fig. 1(c), opening a small gap at Fermi energy in DOS.
In both cases, the DOS shows that the main contributions near the Fermi energy
originate from Hg-${5d}$, Sn-${5p}$, and X-${np}$ (n=2,3 for X=N/P) orbitals,
with minor contributions from Sn-$5s$, and X-$ns$ (n=2,3 for X=N/P) orbitals.
However, these minor parts are not shown in the DOS plots. Notably, X-$np$
(n=2,3 for X=N/P) orbitals dominate over Hg-${5d}$ orbitals near the Fermi
energy, and the presence or absence of states in DOS near the Fermi energy is
mostly determined by the $p$-orbital contribution from the X atom, which is
further crucial for understanding the electronic band structure.

To gain a comprehensive understanding of the topological aspects, we have
analyzed the electronic band structure of HgSnX2 (X=N, P). We have plotted the
band structures of HgSnX2 (X=N, P) using both the GGA and GGA+SOC methods,
focusing on the few high symmetry points, as shown in Fig. 2. A typical
connection of both bands is the removal of degeneracy in the GGA bands in the
presence of SOC, which leads to the involvement of these bands in generating
nontrivial topological signatures and forming a gapless state (band crossing)
or band gap. Specifically, in Fig. 2(a), we have noticed that the
incorporation of SOC has a fascinating effect on the degeneracy of the GGA
band structure. SOC has broken the two-fold bands in both the valence and
conduction bands near Fermi energy. Subsequently, a noteworthy occurrence
takes place when a particular section of those valence and conduction bands
come close to each other and again form a point degeneracy in the band
structure, which leads to Weyl crossing protected by two-fold degeneracy at a
non-high symmetry point in momentum space. These coordinates are found to be
(-0.05534, 0.07620, -0.04912) along the path from S_0 to $\Gamma$. It is
important to highlight that the significant role played by SOC is the creation
of point degenerate Weyl crossing for HgSnN2, while the gapped band structure
for HgSnP2.

Additionally, by examining the orbital-projected bands, we have found that the
major contribution in the Hg atom (like $\Gamma_{6}$) comes from the Hg-$6s$
orbital, which is lower in energy than the major contribution from the N-$2p$
orbital (like $\Gamma_{8}$). This observation leads to a negative band gap at
the $\Gamma$ symmetry point, supporting a nontrivial topological phase. We
have also observed a contribution from the N-$2s$ orbital in addition to the
Hg-$6s$ orbital, as depicted in Fig. S3(a). The $p$ orbitals of N (2py and
2pz) are also marked in Fig. S3(b) and Fig. S3(c), respectively. Therefore,
conducting a detailed investigation of the band-crossing phenomenon is crucial
in light of the relevant topological signatures. In the case of HgSnP2, an
interesting observation was made when incorporating the SOC strength. It led
to the elimination of the two-fold band degeneracy in both the valence and
conduction bands; consequently, the specific regions of these bands come close
to each other, resembling the behavior observed in HgSnN2. However, unlike
HgSnN2, the bands in HgSnP2 do not exhibit a Weyl crossing. Instead, a notable
distinction was noticed where a global band gap of approximately 14 meV
emerges between the valence and conduction bands at the $\Gamma$ point, as
illustrated in Fig. 2(c). There is still band inversion at the high symmetry
$\Gamma$ point, where the major contribution comes from the s orbitals of Hg
and P ($\Gamma_{6}$), which are lower in energy than the major p orbital
contribution ($\Gamma_{8}$), as shown in Fig. 2(d), Fig. S3(d), Fig. S3(e),
and Fig. S3(f). The point to be noted is that it is essential to acknowledge
that the presence of SOC primarily alters the characteristics of the GGA
bands. The introduction of SOC leads to the breaking of degeneracy in the GGA
bands for both compounds, resulting in a state with either no energy gap or a
noticeable gap in the GGA+SOC bands. This intriguing observation prompts us to
conduct a more in-depth analysis of the topological characteristics exhibited
by both compounds and, simultaneously, to see the impact of SOC on their
respective topological phases. Interestingly, although the strength of SOC is
identical in both compounds, but they still display distinct topological
phases. Therefore, in the next section, our investigation aims to uncover the
underlying reasons behind this disparity, and the anticipation regarding the
role of SOC also demands further exploration.

##  V Topological investigation

In this section, we have extensively examined the topological properties of
both compounds by implementing an ab-initio-derived model Hamiltonian
approach. Our main objective was to discern and understand the distinctive
topological traits displayed by these materials.

###  V.1 Topological signatures of HgSnN2

Figure 3: (Color Online) (a) The calculated Fermi arcs, with a chemical
potential positioned approximately 7 meV above the Fermi energy (Ef),
establish connections between the opposite chiral Weyl points in the momentum
space. These Weyl points are indicated by red (source) and blue (sink) filled
circles in the $K_{y}-K_{z}$ momentum plane. In the same momentum plane, (b)
illustrates the corresponding distribution of normalized Berry curvatures in
the outward (source) and inward (sink) direction of fluxes between the two
opposite chiral Weyl points. Displayed in (c) is the flow of the Wannier
Charge Centers (WCC) that move in opposite directions for source and sink
points.

In HgSnN2, we have successfully identified a band crossing phenomenon in the
band structure along the S_0 to $\Gamma$ path, as illustrated in Fig. 2(a),
which exhibits all the characteristic features of a time-reversal symmetry
protected Weyl semimetallic phase. In a Weyl semimetal, there exists a Fermi
arc connecting two Weyl crossings of opposite chirality in momentum space.
These Weyl crossings occur in pairs and can be viewed as magnetic monopoles in
momentum space. For HgSnN2, we discovered a similar outcome that two Weyl
crossings, WP+ with a momentum space coordinate of (-0.0553, 0.0762, -0.0491)
and WP- with a coordinate of (-0.0550, 0.0759, 0.0493) which are connected by
a Fermi arc in the $K_{y}-K_{z}$ momentum plane, as depicted in Fig. 3(a),
specifically around an energy contour close to the chemical potential,
approximately 7 meV above the Fermi energy level. The opposite chirality of
the Weyl crossings is clearly evident from their normalized Berry curvature
distribution, as displayed in Fig. 3(b) in the same momentum plane. The Berry
curvature near the blue-filled circle is inward, acting as a sink in momentum
space, while the filled red circle represents a source with an outward field
line direction. To validate our analysis further, we examined the evolution of
the Wannier Charge Center (WCC) throughout the Brillouin Zone (BZ). Our
analysis revealed that the average position of the WCC, calculated using the
Wilson-loop method, shows a discontinuous jump marked by a double-headed
arrow. Notably, the average position of the WCC for the sink (upper panel) is
in the opposite direction compared to the source (lower panel). We have
enlisted the coordinates of the Weyl crossing with its chirality and energy
positions with respect to the Fermi level in Table V.

| Coordinates | Chemical | Chirality
---|---|---|---
Compound | in K-space (1/$\AA$) | Potential |
|  | (meV) |
HgSnN2 | -0.0550 | 0.0759 | 0.0493 | 6.8 | -1 (WP-)
| -0.0553 | 0.0762 | -0.0491 | 7.9 | +1 (WP+)
Table 5: In the momentum space, the Weyl node coordinates (in unit of 1/$\AA$)
of HgSnN2 with their chemical potentials relative to the Fermi energy, and
chiralities are listed.

We have also calculated the spectral function using the slab model approach as
implemented in WannierTool[82] and found that two Weyl points are connected by
a Fermi arc- related band in the surface state, as shown in supplementary Fig.
S5(a). Collectively, all these findings in HgSnN2 firmly establish it as a
nonmagnetic Weyl semi-metal, possessing all the key topological
characteristics. Due to the proximity of Weyl crossing to the Fermi energy,
this compound can be regarded as an excellent candidate for an Hg-based ideal
Weyl semi-metal [90], making it highly suitable for experimental
investigations of the Fermi arc. Furthermore, the Weyl crossing observed in
this compound is sufficiently distinct and isolated from the surrounding bulk
bands, facilitating experimental observations and analysis.

###  V.2 Topological signatures of HgSnP2

Figure 4: (Color online) (a) Depicts the flow of average Wannier Charge
Centers (WCC) position for the highest occupied valence band of HgSnP2,
obtained through the Wilson-loop method. Notably, the flow showcases a
continuous trajectory of the WCC. (b) The surface spectral function across the
surface Brillouin zone (BZ) is represented in this figure for HgSnP2. A green
circle highlights the presence of the odd numbers of the Dirac crossing. (c)
shows the Dirac cone at $\Gamma$ point shown along [100] for clear
visualization. The energy scale is referenced to the Fermi energy(Ef), which
is set at zero.

The presence of a narrow band gap of 14 meV and a band inversion between the
$\Gamma_{6}$ and $\Gamma_{8}$ bands, as depicted in Fig. 2(c) and Fig. 2(d),
prompted us to investigate the topological properties of HgSnP2. While we did
not observe any band crossings similar to its sister compound HgSnN2, but the
existence of band inversion with a small band gap between the maxima of
valence and minima of the conduction band is a significant indication that
could be categorized as a topological insulator. Following standard practices
for studying topological behaviors, we examined the evolution of the WCC along
the Ky direction of the Brillouin Zone (BZ) in Kx=0.0 plane, as shown in Fig.
4(a). We found that the WCC exhibits continuity, indicative of its nontrivial
topological nature. Additionally, we plotted the same quantities, i.e., WCC,
along different directions in the BZ plane, as shown in the supplementary Fig.
S5(b) other momentum planes in the Kx=0.5, Ky=0.5, and Kz=0.5 which are
discontinuous. Based on these findings, we conclude that HgSnP2 possesses a
nonzero topological index with a $Z_{2}$ number [1,000] which falls under the
category of strong topological insulator.

To further support our conclusions, we calculated the spectral function or
surface state in the [001] and [100] direction, scanning the entire BZ at
time-reversal invariant momenta (TRIM) points. The analysis revealed a gap-
less surface state with an odd number of Dirac crossings, as illustrated in
Fig. 4(b)-(c). We analyzed the spectral distribution at the surface of the
semi-infinite slab of HgSnP2. The presence of a Dirac-type gap-less crossing
in the spectral function serves as a crucial indicator for identifying the
topological state. This observation aligns with our earlier findings related
to bulk bands, affirming the consistency of our results. We utilized Green’s
function techniques, employing the WannierTool[82] to extract the tight-
binding model through maximally projected Wannier functions (WFs). In our
visual representation, the red bands correspond to surface bands against a
blue background of bulk bands. For example, no crossing is observed at the M
high symmetry point of the bulk band, approximately -0.15 eV. However, in the
surface bands, we identified a Dirac-like band crossing. These crossings are
denoted as SS1, SS2, and SS3, with SS1 located at the N high symmetry point,
SS2 at the $\Gamma$ high symmetry point, and SS3 at the M high symmetry point
along the [001] and [100] direction. Since SS1 is not very clear along the
[001] direction, we have identified it clearly in the [100] directions. We
have put an image along [100] to show the Dirac cone at $\Gamma$ point for
clear visualization in Fig. 4(c). Consequently, our analysis indicates an odd
number of Dirac cones in the surface state, providing valuable insights into
the topological characteristics of the system.

The $Z_{2}$ number [1,000] and the existence of an odd number of Dirac cones
unequivocally establish HgSnP2 as a strong topological insulator,
characterized by a bulk band gap and a metallic surface state protected by
time-reversal symmetry. We have also observed spin momentum locking, as shown
in supplementary Fig. S5(c). While the chalcopyrite class offers numerous
examples of topological insulators, our current understanding suggests a lack
of comprehensive studies focusing on the topological aspects of Hg-based
chalcopyrites in the existing literature. In a nutshell, the relevant
signature indicates that HgSnP2 is a strong topological insulator.

##  VI Effect of SOC modulation

Figure 5: (Color online) Calculated GGA+SOC band structure for (a) HgSnN2 and
(b) HgSnP2 with different SOC strength ($\lambda$). (c) and (d) represent the
projected band structures which show the band inversion at 3$\lambda$
($\xi$=3) and 0.3 $\lambda$ ($\xi$=0.3) for HgSnN2 and HgSnP2 respectively.
Here $\lambda$ is the original SOC strength, and $\xi$ is the scaling factor
of SOC. The energy scale is referenced to the Fermi energy(Ef), which is set
at zero.

The primary focus revolves around comprehending the fundamental factors that
give rise to the distinct (gap and gap less) topological phases observed in
two closely related compounds, namely, HgSnN2 and HgSnP2, which have the same
chemical composition except for the anion (N/P). These compounds share the
same space group with broken inversion symmetry and preserving time-reversal
symmetry. In this context, it is essential to emphasize that Hg has strong
atomic SOC strength, while atomic N or P does not contribute to SOC strength
significantly. Interestingly, despite both compounds having the same atomic
SOC strength, but they still display different topological phases, which
prompts us to delve deeper into the underlying reasons behind this non-
identical behavior.

To reinforce our analysis further, we have taken into account a $\xi$ term in
the relativistic Hamiltonian, H${}_{\text{SOC}}$ = $\xi$$\lambda$
$\vec{L}$.$\vec{S}$, where $\vec{L}$ and $\vec{S}$ represent orbital and spin
angular momentum, respectively as followed by the prescription[91, 92] and
$\lambda$ is the original SOC strength while $\xi$ serve as a modulating
factor of original SOC strength. Generally, SOC is an inherent atomic property
of the element that remains unaltered with the same nominal valance state.
However, we can artificially tune the SOC by scaling the relativistic term in
the DFT calculations. Indeed, the artificial increase or decrease of SOC is
not factually correct and does not represent the physical picture. However,
such model calculations will be very effective in understanding the impact of
the SOC in comparison to the other energy scales involved in the materials. In
the case of original SOC strength $\xi$=1, but we can arbitrarily vary the
value of $\xi$ to probe its effect. We have explored the consequences by
sequentially scaling up and down the value of $\xi$ in Hg-based materials. The
strength of the SOC is dependent on the fourth power of the atomic number,
meaning that the effect of SOC becomes more dominant as we move down the
periodic table and very strong for Hg. For HgSnN2, we have utilized values of
1.2$\lambda$ (i.e $\xi$=1.2) and 3$\lambda$ (i.e $\xi$=3), while for HgSnP2,
the effect has been reduced to 0.8$\lambda$ (i.e $\xi$=0.8) and 0.3$\lambda$
(i.e $\xi$=0.3). Remarkably, we have observed that we can artificially
calibrate multiple topological phases by scaling up and down $\xi$ parameter
in both materials.

Now, we investigated the effect of varying SOC strength on the Weyl crossing
in HgSnN2. In Fig. 5(a), we increased the strength of SOC ($\xi$), say, up to
three times the original SOC strength (3$\lambda$) and observed that at
1.2$\lambda$, the gapless feature disappears, and a small gap begins to form,
which becomes more prominent at 3$\lambda$. These two artificial SOC-driven
band structures at 1.2$\lambda$ and 3$\lambda$ in HgSnN2 indicate that
relativistic SOC indeed plays a crucial role in maintaining the gapless nature
of the band structure. We performed a similar analysis on HgSnP2 for cross-
validation, as shown in Fig. 5(b). In this case, we reduced the strength of
SOC, say, to 0.8$\lambda$ and 0.3$\lambda$. We found that the band gap at the
same momentum coordinate vanishes. At 0.8$\lambda$, the band gap between the
valence and conduction bands is reduced; at 0.3$\lambda$, it becomes gap-less.
However, it is essential to note that the band structure at 0.3$\lambda$ in
HgSnP2 is already artificially influenced in terms of SOC, so further
verification of its topological characteristics is not meaningful. However,
from first-principles DFT calculations, we can verify the nontrivial
properties by observing the band ordering in terms of orbital contributions,
whether they are normal or inverted. We observed that in HgSnN2 at SOC
strength 3$\lambda$, the major Hg-$6s$ orbital contribution ($\Gamma_{6}$) is
lower than the major N-$p$ orbital contribution ($\Gamma_{8}$) as shown in
Fig. 5(c). Furthermore, the involvement of the N-$2s$ orbital has been
identified alongside the Hg-$6s$ orbital, illustrated in Fig. S4(a). The
N-$2p$ orbitals (2py and 2pz) have also been designated in Fig. S4(b) and Fig.
S4(c) correspondingly. Hence, a negative band gap is present at the high-
symmetry point ($\Gamma$), similar to the original band structure as shown in
Fig. 2(b). Similarly, in HgSnP2 at 0.3$\lambda$, the major Hg-$6s$ orbital
contribution from the Hg atom ($\Gamma_{6}$) is lower than the major P-$3p$
orbital contribution ($\Gamma_{8}$) as shown in Fig. 5(d). We have observed an
additional involvement of the P-$3s$ orbital alongside the Hg-$6s$ orbital,
which is illustrated in Fig. S4(d). The $3p$ orbitals of P, namely 3py and
3pz, have also been indicated in Fig. S4(e) and Fig. S4(f), respectively.
Notably, both band structures at artificially modified SOC strengths also
exhibit band inversion. The study of these sister compounds enables us to
examine the role of SOC in driving the topological phase in quantum materials
from a modeling perspective. Another critical perspective is that we have
noticed a striking resemblance between the gapped band structure of HgSnN2 at
0.3$\lambda$ and the parent band structure of HgSnP2 at its original SOC
strength. Similarly, a comparable equivalence can be established for HgSnP2,
where the gapless band structure at 3$\lambda$ closely resembles the parent
state of HgSnN2 at its original SOC strength. In summary, by modulating the
SOC strength, we observe the crossover from gapped to gapless and vice-versa
topological phases can emerge in these two sister compounds, which can also be
achievable by tuning the chemical composition via N or P.

Figure 6: (Color online)The GGA+SOC orbital projected hybridized band
structure of HgSnN2 with major N-$p$ (grey) and minor Hg-${5d}$(blue)
orbitals. The energy scale is referenced to the Fermi energy(Ef), which is set
at zero (Ef=0.0 eV)

Ideally, both the compound HgSnN2 and HgSnP2 should be the same in terms of
SOC. To correlate such a strong effect of SOC in these compounds, we found a
notable feature that the band crossing or Weyl crossing observed in HgSnN2
mainly arises from the hybridization of the major N-$p$ orbital and minor
Hg-${5d}$ orbital states, as shown in Fig. 6. This hybridization plays a
crucial role in explaining the topological characteristics of the respective
compounds. The SOC strength is particularly effective in the hybridized band
due to the presence of Hg-${5d}$ orbitals. The minor contribution from
Hg-${5d}$ orbitals may modify the band dispersion without substantial changes
in the overall band structure. The impact of SOC on the topological phase is
widely recognized, but this study has elegantly illuminated the crucial role
of hybridization in modifications to the topological phase.

##  VII Summary and Conclusion

In this paper, we have undertaken a comprehensive ab-initio analysis of two
topological materials, namely HgSnN2 and HgSnP2, highlighting their intriguing
topological properties driven by SOC. We confirmed the structural stability of
the compounds by the absence of phonon-negative frequency modes. In our
investigation, we found that HgSnN2 exhibits characteristics of a nonmagnetic
Weyl semi-metal, and the Weyl crossing is protected by a point band
degeneracy. In the case of HgSnP2, we found a similar band inversion even
after the N is replaced by P. This change in chemical composition effects the
topological properties, mainly opening a small band gap. In particular, HgSnP2
is identified as a strong topological insulator. The point to be noted is that
although both the compounds offer the same atomic SOC strength, but the
hybridization of bands between the elements present in our compounds plays a
critical role in deciding whether the topological phase is a topological
insulator or semi-metal. Furthermore, we have observed that by adjusting the
strength of SOC in the model calculation, for the case of HgSnN2, it is
possible to generate equivalent topological signatures to those of the parent
HgSnP2, and vice versa. Such tuning of SOC can be accomplished through doping,
substitution, and the design of heterostructures in the quantum materials[48,
49, 50].

Exploring topological properties in Hg-based chalcopyrite has been relatively
less reported compared to other chalcopyrites. We have observed that
chalcopyrites based on N and P have been successfully synthesized, including
the formation of single crystals[61, 63]. In contrast, both structures lack
inversion symmetry but preserve time-reversal symmetry, making our study
simple to capture the dependence of SOC effectively. This type of work can be
further extended for any suitable choice of such sibling magnetic compounds to
capture the interplay of topology and magnetism. Furthermore, the remarkable
proximity of the Weyl crossings to the Fermi energy in the case of HgSnN2
suggests the potential for experimental detection. These findings pave the way
for exciting advancements in these fields.

##  VIII Acknowledgments

SS expresses gratitude to IIT Goa for providing the research fellowship under
the Government of India.

## References

  * [1] E E Krasovskii, J. Condens. Matter Phys. $\bf{27}$, 493001 (2015).
  * [2] X. Fan, H. Celik, J. Wu, C. Ni, K. J. Lee, V. O. Lorenz and J. Q. Xiao, Nat. Commun. $\bf{5}$, 3042 (2014).
  * [3] A. Soumyanarayanan, N. Reyren, A. Fert and C. Panagopoulos, Nature $\bf{539}$, 509 (2016).
  * [4] E. G. Moon, C. Xu, Y. B. Kim and L. Balents, Phys. Rev. Lett. $\bf{111}$, 206401 (2013).
  * [5] B. Huang, Y. B. Kim, and Y. M. Lu, Phys. Rev. B $\bf{95}$, 054404 (2017).
  * [6] B. Huang, Y. B. Kim, and Y. M. Lu, npj Quantum Mater. $\bf{3}$, 23 (2018).
  * [7] K. Kitagawa, T. Takayama, Y. Matsumoto, A. Kato, R. Takano, Y. Kishimoto, S. Bette, R. Dinnebier, G. Jackeli and H. Takagi, $\bf{554}$, 341 (2018).
  * [8] H. S. Kim, V. Shankar V., A. Catuneanu, and H. Y. Kee, Phys. Rev. B $\bf{91}$, 241110 (2015).
  * [9] A. B. Shick and W. E. Pickett, Phys. Rev. B $\bf{100}$, 134502 (2019).
  * [10] S. Chakrabortty, R. Kumar and N. Mohapatra, Phys. Rev. B $\bf{107}$, 024503 (2023).
  * [11] A. Manchon, H. C. Koo, J. Nitta, S. M. Frolov and R. A. Duine , Nat. Mater. $\bf{14}$, 871 (2015).
  * [12] A. D. Caviglia, M. Gabay, S. Gariglio, N. Reyren, C. Cancellieri and J.-M. Triscone, Phys. Rev. Lett. $\bf{104}$, 126803 (2010).
  * [13] M. Governale, Phys. Rev. Lett. $\bf{89}$, 206802 (2002).
  * [14] T. Yu, Z. Luo and G. E.W. Bauer , Phys. Rep. $\bf{1009}$, 1 (2023).
  * [15] K. W. Kim, H. W. Lee, K. J. Lee and M. D. Stiles, Phys. Rev. Lett. $\bf{111}$, 216601 (2013).
  * [16] S. Grytsiuk, J. P. Hanke, M. Hoffmann, J. Bouaziz, O. Gomonay, G. Bihlmayer, S. Lounis, Y. Mokrousov and S. Blügel , Phys. Rev. Lett. $\bf{11}$, 511 (2020).
  * [17] S. Banerjee, J. Rowland, O. Erten, and M. Randeria, Phys. Rev. X $\bf{4}$, 031045 (2014).
  * [18] K. W. Kim, K. J. Lee, J. Sinova, H. W. Lee and M. D. Stiles, Phys. Rev. B $\bf{96}$, 104438 (2017).
  * [19] A. Brataas and K. M. D. Hals, Nat. Nanotechnol. $\bf{9}$, 86 (2014).
  * [20] L. Zhu, D.C. Ralph and R. A. Buhrman, Phys. Rev. Lett. $\bf{112}$, 077201 (2019).
  * [21] S. h. C. Baek, V. P. Amin, Y. W. Oh, G. Go, S. J. Lee, G. H. Lee, K. J. Kim, M. D. Stiles, B. G. Park and K. J. Lee, Nat. Mater. $\bf{17}$, 509 (2018).
  * [22] C. Ciccarelli, L. Anderson, V. Tshitoyan, A. J. Ferguson, F. Gerhard, C. Gould, L. W. Molenkamp, J. Gayles, J. Železný, L. Šmejkal, Z. Yuan, J. Sinova, F. Freimuth and T. Jungwirth , Nat. Phys. $\bf{12}$, 855 (2016).
  * [23] N. Tesařová, P. Němec, E. Rozkotová, J. Zemen, T. Janda, D. Butkovičová, F. Trojánek, K. Olejník, V. Novák, P. Malý and T. Jungwirth , Nat. Photonics $\bf{7}$, 492 (2013).
  * [24] S. Lee, M. G. Kang, D. Go, D. Kim, J. H. Kang, T. Lee, G. H. Lee, J. Kang, N. J. Lee, Y. Mokrousov, S. Kim, K. J. Kim, K. J. Lee and B. G. Park , Commun. Phys. $\bf{4}$, 234 (2021).
  * [25] C. Song, R. Zhang, L. Liao, Y. Zhou, X. Zhou, R. Chen, Y. You, X. Chen and F. Pan, Prog. Mater. Sci. $\bf{118}$, 100761 (2021).
  * [26] Y. Li, K. W. Edmonds, X. Liu, H. Zheng and K. Wang, Adv. Quantum Technol. $\bf{2}$, 1800052 (2019).
  * [27] C. L. Kane and E. J. Mele, Phys. Rev. Lett. $\bf{95}$, 226801 (2005).
  * [28] A. Hoffmann, IEEE Trans. Magn. $\bf{49}$, 5172 (2013).
  * [29] S. N. Perge, S. M. Frolov, E. P. A. M. Bakkers and L. P. Kouwenhoven, Nature $\bf{468}$, 1084 (2010).
  * [30] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. $\bf{82}$, 3045 (2010).
  * [31] X. L. Qi and S. C. Zhang, Rev. Mod. Phys. $\bf{83}$, 1057 (2011).
  * [32] W. Beugeling, N. Goldman and C. M. Smith, Phys. Rev. B $\bf{86}$, 075118 (2012).
  * [33] X. Zhang, J. Liu and F. Liu, Nano Lett. $\bf{22}$, 9000 (2022).
  * [34] M. Kim, Z. Wang, Y. Yang, H. T. Teo, J. Rho and B. Zhang, Nat. Commun. $\bf{13}$, 3499 (2022).
  * [35] S. M. Young, S. Chowdhury, E. J. Walter, E. J. Mele, C. L. Kane and A. M. Rappe, Phys. Rev. B $\bf{84}$, 085106 (2011).
  * [36] J. O. Island, X. Cui, C. Lewandowski, J. Y. Khoo, E. M. Spanton, H. Zhou, D. Rhodes, J. C. Hone, T. Taniguchi, K. Watanabe, L. S. Levitov, M. P. Zaletel and A. F. Young , Nature $\bf{571}$, 85 (2019).
  * [37] M. Kheirkhah, Z. Yan, Y. Nagai and F. Marsiglio, Phys. Rev. Lett. $\bf{125}$, 017001 (2020).
  * [38] S. Fang, L. Ye, M. P. Ghimire, M. Kang, J. Liu, M. Han, L. Fu, M. Richter, J. Brink, E. Kaxiras, R. Comin and J. G. Checkelsky, Phys. Rev. B $\bf{105}$, 035107 (2022).
  * [39] L. Rademaker, Phys. Rev. B $\bf{105}$, 195428 (2022).
  * [40] N. T. Dejean, F. G. Eich and A. Rubio , Npj Comput. Mater. $\bf{8}$, 145 (2022).
  * [41] L. Tian, Y. Liu, W. Meng, X. Zhang, X. Dai and G. Liu, J. Phys. Chem. Lett. $\bf{11}$, 10340 (2020).
  * [42] D. Pesin and A. H. MacDonald, Nat. Mater. $\bf{11}$, 409 (2012).
  * [43] D. Pesin and A. H. MacDonald, J. Phys. Soc. Jpn. $\bf{84}$, 121003 (2015).
  * [44] G. Bian, T. R. Chang, R. Sankar, S. Y. Xu, H. Zheng, T. Neupert, C. K. Chiu, S. M. Huang, G. Chang, I. Belopolski, D. S. Sanchez, M. Neupane, N. Alidoust, C. Liu, B. Wang, C. C. Lee, H. T. Jeng, C. Zhang, Z. Yuan, S. Jia, A. Bansil, F. Chou, H. Lin and M. Z. Hasan, Nat. Commun. $\bf{7}$, 10556 (2016).
  * [45] B. Singh, S. Mardanya, C. Su, H. Lin, A. Agarwal and A. Bansil, Phys. Rev. B $\bf{98}$, 085122 (2018).
  * [46] S. Howlader, S. Saha, R. Kumar, V. Nagpal, S. Patnaik, T. Das and G. Sheet, Phys. Rev. B $\bf{102}$, 104434 (2020).
  * [47] X. L. Sheng, Z. Wang, R. Yu, H. Weng, Z. Fang and X. Dai, Phys. Rev. B $\bf{90}$, 245308 (2014).
  * [48] P. Dziawa, B. J. Kowalski, K. Dybko, R. Buczko, A. Szczerbakow, M. Szot, E. Łusakowska, T. Balasubramanian, B. M. Wojek, M. H. Berntsen, O. Tjernberg and T. Story , Nat. Mater. $\bf{11}$, 1023 (2012).
  * [49] F. Mahfouzi, N. Nagaosa and B. K. Nikolić, Phys. Rev. Lett. $\bf{109}$, 166602 (2012).
  * [50] T. Naimer and J. Fabian, Phys. Rev. B $\bf{107}$, 195144 (2023).
  * [51] A Belayadi and P Vasilopoulos, Nat. Nanotechnol. $\bf{34}$, 365706 (2023).
  * [52] M. Zholudev, F. Teppe, M. Orlita, C. Consejo, J. Torres, N. Dyakonova, M. Czapkiewicz, J. Wróbel, G. Grabecki, N. Mikhailov, S. Dvoretskii, A. Ikonnikov, K. Spirin, V. Aleshkin, V. Gavrilenko and W. Knap, Phys. Rev. B $\bf{86}$, 205420 (2012).
  * [53] Y. Li and G. Xu, Phys. Rev. Mater. $\bf{6}$, 104201 (2022).
  * [54] F. C. D. Lima, R. H. Miwa and A. Fazzio, Phys. Rev. B $\bf{102}$, 235153 (2020).
  * [55] P. F. Liu , L. Zhou , S. Tretiak and L. M. Wu , J. Mater. Chem. C $\bf{5}$, 9181 (2017).
  * [56] F. Virot, R. Hayn, M. Richter and J. Brink, Phys. Rev. Lett. $\bf{111}$, 146803 (2013).
  * [57] S. A. Mughal, A. J. Payne and B. Ray, J. Mater. Sci. $\bf{4}$, 895 (1969).
  * [58] G. A. Verozubova, A.O. Okunev, A.I. Gribenyukov, A.Y. Trofimiv, E.M. Trukhanov and A.V. Kolesnikov, J. Cryst. Growth $\bf{312}$, 1122 (2010).
  * [59] K. Masumoto, S. Isomura and W. Goto, J. Phys. Chem. Solids $\bf{27}$, 1939 (1966).
  * [60] M. Rubenstein and R W Ure Jr., J. Phys. Chem. Solids $\bf{29}$, 551 (1968).
  * [61] K. T. Zawilski, P. G. Schunemann, T. C. Pollak, D. E. Zelmon, N. C. Fernelius and F. K. Hopkins, J. Cryst. Growth $\bf{312}$, 1127 (2010).
  * [62] G. Zhang, X. Tao, H. Ruan, S. Wang and Q. Shi, J. Cryst. Growth $\bf{340}$, 197 (2012).
  * [63] P. Vohl, J. Electron. Mater. $\bf{8}$, 517 (1979).
  * [64] E. Buehler and J.H. Wernick, J. Cryst. Growth $\bf{8}$, 324 (1971).
  * [65] S. Shirakata and S. Isomura, J. Cryst. Growth $\bf{99}$, 781 (1990).
  * [66] R.T. Trykozko, Mater. Res. Bull. $\bf{10}$, 489 (1975).
  * [67] I. S. Khan, K. N. Heinselman and A. Zakutayev, J. Phys. Energy $\bf{2}$, 032007 (2020).
  * [68] F. Kawamura, H. Murata, M. Imura, N. Yamada and T. Taniguchi, Inorg. Chem. $\bf{60}$, 1773 (2021).
  * [69] A. L. Greenaway, A. L. Loutris, K. N. Heinselman, C. L. Melamed, R. R. Schnepf, M. B. Tellekamp, R. W. Robinson, R. Sherbondy, D. Bardgett, S. Bauers, A. Zakutayev, S. T. Christensen, S. Lany and A. C. Tamboli, J. Am. Chem. Soc. $\bf{142}$, 8421 (2020).
  * [70] C. L. Rom, M. J. Fallon, A. Wustrow, A. L. Prieto and J. R. Neilson,Chem. Mater. $\bf{33}$, 5345 (2021).
  * [71] S. Esmaeilzadeh, U. Hålenius and M. Valldor, Chem. Mater. $\bf{18}$, 2713 (2006).
  * [72] F. J.Fan, L. Wu and S. H. Yu , Energy Environ. Sci. $\bf{7}$, 190 (2014).
  * [73] M. C. Ohmer and R. Pandey, MRS Bull. $\bf{23}$, 16 (1998).
  * [74] B. Sadhukhan, Y. Zhang, R. Ray and J. Brink, Phys. Rev. Mater. $\bf{4}$, 064602 (2020).
  * [75] K. Momma and F. Izumi, J. Appl. Crystallogr.$\bf{41}$, 653 (2008).
  * [76] G. Kresse and J. Hafner, Phys. Rev. B $\bf{47}$, 558(R) (1993).
  * [77] G. Kresse and J. Furthmuller, Phys. Rev. B $\bf{54}$, 11169 (1996).
  * [78] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. $\bf{77}$, 3865 (1996).
  * [79] N. Marzari and D. Vanderbilt, Phys. Rev. B $\bf{56}$, 12847 (1997).
  * [80] I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B $\bf{65}$, 035109 (2001).
  * [81] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, et al., J. of Phys.: Cond. Matt. $\bf{32}$, 165902 (2020).
  * [82] Q. Wu, S. Zhang, H.-F. Song, M. Troyer, and A. A. Soluyanov, Computer Physics Communications $\bf{224}$, 405 (2018).
  * [83] C. Suh and K. Rajan, Appl. Surf. Sci. $\bf{223}$, 148 (2004).
  * [84] Y. M. Basalaev and M. Y. Basalaeva, J. Struct. Chem. $\bf{61}$, 1007 (2020).
  * [85] Y. Zhong, H. Mei, D. He, X. Du and N. Cheng, J. Phys. Chem. Solids $\bf{134}$, 157 (2019).
  * [86] J. E. Jaffe and A. Zunger,Phys. Rev. B $\bf{29}$, 1882 (1984).
  * [87] S. Sadhukhan, B. Sadhukhan and S. Kanungo, Phys. Rev. B $\bf{106}$, 125112 (2022).
  * [88] F. Mouhat and F. X. Coudert, Phys. Rev. B $\bf{90}$, 224104 (2014).
  * [89] A. H. Baghdad, H. Bouafia, B. Djebour, B. Sahli, S. Hiadsi, M. Elchikh and M. Attou , J. Phys. Chem. Solids $\bf{167}$, 110756 (2022).
  * [90] J. Ruan, S. K. Jian, D. Zhang, H. Yao, H. Zhang, S. C. Zhang and D. Xing,Phys. Rev. Lett. $\bf{116}$, 226801 (2016).
  * [91] S. Kanungo, K. Mogare, B. Yan, M. Reehuis, A. Hoser, C. Felser and M. Jansen, Phys. Rev. B $\bf{93}$, 245148 (2016).
  * [92] R. Roy and S. Kanungo, Phys. Rev. B $\bf{106}$, 125113 (2022).