\section{Introduction} 

Electrides are a novel class of crystals with vacant cavities accommodating intrinsic excess electrons, which are not tightly constrained by certain nuclei, but on the contrary, form ionic bonds with the cationic framework, playing the role of anions \cite{dye1990663, dye2003electrons}. In contrast to random defects or color centers, sites of interstitial electrons in electrides adhere to the translational periodicity. Since the experimental synthesis of the first room-temperature stable inorganic electride [Ca$_{24}$Al$_{28}$O$_{64}$]$^{4+}(4e^-)$ (C12A7:2e$^-$) \cite{matsuishi2003high}, electrides have attracted great attention among the communities of materials science, chemistry and condensed matter physics, owing to the intriguing properties exhibited by non-nucleus-bound electrons. For instance, in magnetic electrides the spin-unpaired anionic electrons are well suited for spin injection in a spintronic device \cite{sui2019prediction, wolf2001spintronics}, on account of their low work function. In addition, the interstitial electronic states are found to favor band inversion near the Fermi level, giving rise to the high tendency for electrides towards topological phases \cite{hirayama2018electrides, zhu2019computational}. One could therefore implement electrides as an excellent testbed to study topological physics \cite{PhysRevLett.120.026401, huang2018topological, PhysRevLett.123.206402, zhang2019topological, nie2020application}.

Hitherto, known electrides especially magnetic electrides, however, are extremely scarce\cite{liu2020electrides}. Consecutive efforts have been devoted by a number of researchers to the computational discovery of new electrides. Approaches of the mainstream involve: structure screening from the materials database \cite{zhu2019computational, PhysRevX.4.031023, Tada2014High-Throughput, Burton2018High-Throughput}, crystal structure prediction (CSP) for unknown prototype structures \cite{tsuji2016structural, ming2016first, zhang2017computer, wang2017exploration}, and functionality-tailored electrides design \cite{PhysRevLett.120.026401, PhysRevMaterials.3.024205}. Here we present an alternative computational strategy to design new electrides. It is widely noticed that every natural electride corresponds to a potential oxidized form, such as Ca$_2$N to Ca$_2$NCl, Y$_2$C to Y$_2$H$_2$C. One is therefore motivated to ask whether the reverse direction holds \cite{liu2020electrides}. In other words, if the anions with large electronegativity is taken away from a selected compound (e.g. remove Cl from Ca$_2$NCl), the remaining structure (Ca$_2$N) is suspected to be a possible electride.
To ascertain this question, we mainly focus on ternary halides containing trivalent rare-earth elements, i.e. scandium, yttrium, and lanthanides. Optimizing the reduced structure using a CSP method, a stable or metastable crystal is expected to be the candidate of a binary electride. Chemical substitution of magnetic lanthanides is followed to induce spin-polarized interstitial electronic states, by which magnetic electrides can be constructed. Starting from 11 halides serving as structure hosts, our method predicts 30 thermodynamically stable and metastable nonmagnetic electrides, 29 of which are newly designed. While retaining the thermodynamic (meta)stability, we replace cationic elements with the lanthanide Gd to induce spin polarization for these crystals. This chemical substitution procedure successfully generates 28 candidate structures of magnetic electrides, with 27 of them proposed for the first time. Electronic properties of two representative systems are analyzed in depth. In addition to electride and magnetism features, topological nodal lines in the electronic bands of these two systems are unveiled as well, further confirming the intimate relationship between electrides and topological materials. 

\begin{figure}[htp]
	\centering
	\includegraphics[width=0.95\linewidth]{./figures/flowchart.pdf}
	\caption{Flowchart of designing halide-derived electrides.}
	\label{flowchart}
\end{figure}

\section{Methods}
\subsection{Design scheme}
\label{subsec_design_scheme}

The procedure of our design scheme is presented in Figure \ref{flowchart}, beginning with the selection of halides as electride hosts. For the sake of concreteness, we pick 11 promising ternary halides as host crystals (see Appendix \ref{sec_allcrystals}), which can be classified into 4 stoichiometric families, $A_2$C$X_p$, $A$C$X_p$, $A_2$Ge$X_p$ and $A$Ge$X_p$, where $A$ and $X$ stand for rare-earth and halogen elements, respectively; $p$ be the ratio of mole fraction of the halogen. The halogen atoms in these materials are surrounded by cations, so that the exclusion of halogens gives rise to an intrinsic excess of electrons as well as interstitial sites to accommodate those excess electrons, both of which are necessary for electride formation \cite{zhu2019computational}. We have restricted our scope to rare-earth based halides, considering that rare-earth elements are favorable electron donors in electrides because of their low eletronegativity and standard valence (usually +3). Moreover, close radii among rare-earth atoms form a good condition for chemical substitution (e.g. Y $\rightarrow$ Gd). Owing to intermediate electronegativities (2.55 and 2.01, respectively), C and Ge from group \rom{4}A are chosen as anionic elements, so that electrons could deviate from cationic atoms, and reside in interstitial orbitals at the same time \cite{wang2017exploration}. Si with the electronegativity 1.90 is not considered here, as it is less electronegative compared to C and Ge. In fact, electrons are barely attracted by rare-earth atoms even in the case of Ge, see Sec. \ref{subsec_overview}. Thus, the consideration of Si, which is less likely to attract electrons, is obsolete. We point out that the materials screening here is merely based on rules of thumb, which cannot rule out the possibility of potential electride hosts that disobey these conditions, but would always be an efficient guide of searching electrides \cite{zhang2017computer}.

Given the selected halides, we generate initial structures for the CSP. Taken away the halogens, the corresponding 4 categories of binary compounds are $A_2$C, $A$C, $A_2$Ge and $A$Ge, or collectively designated by $A_nB$ ($n=1,2; B=$ C, Ge). Further, the original cation element $A$ is replaced with Y and La separately to develop 2 groups of compositions for each category. For clarity and consistency, throughout the paper we will insist on a classification in four stoichiometric categories, i.e. the member of one category is independent of the chemical nature of $A$. Structure optimization and comparison would be conducted on the Y and La groups of materials individually for each category of $A_nB$. Since neither of Y and La has partially filled 4\textit{f} orbitals, the CSP is entirely freed from spin-polarized calculations to reduce the computational costs. Another reason for this replacement is related to the chemical substitution of magnetic lanthanides, which will be elaborated below.

CSP then follows to eliminate the structure instability caused by the exclusion of halogen atoms. From the viewpoint of CSP, the reduced structure of a known material acts as a reasonable initial guess, assisting CSP to reach the ground state more efficiently, in contrast to random starting structures \cite{zhang2017computer}. Notice that although the size of interstitial regions might be compressed during the structure optimization, proper vacant spaces adopting anionic electrons are often found among energy-low-lying structures in practice. 

Subsequent to the CSP, convex hull diagrams are constructed in order to estimate the stability and potential synthesizability of predicted structures relative to competing phases with different stoichiometries. The most stable structures whose energies above the convex hull (denoted as $E_\text{hull}$) are smaller than 0.20 eV/atom \cite{oganov2019structure, wang2017exploration} are processed for electride identification, in order not to rule out the synthesis possibility of metastable polymorphs \cite{PhysRevLett.116.075503, sun1600225thermodynamic}. Under each category of $A_nB$, the first five promising electrides with the lowest energies are kept as final candidates for both $A=$ Y and La. We note that $E_\text{hull}$ of these candidates mostly turn out to be smaller than 0.10 eV/atom (see Sec. \ref{subsec_overview}). 

Next, for every stoichiometric category of potential electrides $A_nB$, the cations of 10 candidate structures (5 from Y-based systems and 5 from the La-based) are substituted by magnetic lanthanides (we focus on Gd in this work) so that magnetism could be induced. As the ionic radius of Gd is between the ones of Y and La, it is expedient to assume stable or metastable Gd-based phases could be found among structures inherited from those of Y and La, especially in the case that Y$_nB$ and La$_nB$ share an identical crystal structure. This is the other reason why Y and La are simultaneously introduced at the stage of generating initial structures. In addition to the close radii of substitution elements, the stability is also assessed by constructing convex hulls. For simplicity, only ferromagnetic (FM) configurations of each Gd-based electride candidates are computed. At the same time, anionic electrons easily become spin-unpaired in the case of FM configurations of Gd. Applied to all other lanthanides with in-between ionic radii, this element substitution strategy could massively explore stable or metastable phases of lanthanide-based compounds without an expensive CSP for every lanthanide. This is a key point of our work.

Finally, having all the predicted electrides, we concentrate on the structures of interest for detailed analysis, including the dynamical stability, further confirmation of electride features, etc. The magnetic ground state of selected Gd systems is computed (see Sec. \ref{subsec_comp_methods}), estimating the relative stability of the FM state compared to other possible magnetic states. Also, we particularly characterize the topological nature encoded in electronic bands of the designed electrides. Notice that these procedures are applied to only a subset of predicted crystals, as given in the last box of Figure \ref{flowchart}.

\subsection{Computational methods} 
\label{subsec_comp_methods}

All the electronic structure calculations are performed with the density functional theory (DFT) based on the projector augmented wave scheme \cite{PhysRevB.59.1758} as implemented in the Vienna \textit{ab initio} simulation package ({\sc vasp}) \cite{PhysRevB.54.11169}. The exchange-correlation functional is given by the generalized gradient approximation (GGA) of the Perdew-Burke-Ernzerhof (PBE) type \cite{PhysRevLett.77.3865}. Strong intra-atomic interactions in 4\textit{f} orbitals of Gd are approximated by the GGA+U \cite{PhysRevB.52.R5467}. The nonlocal Heyd-Scuseria-Ernzerhof (HSE06) hybrid functional is supplemented to check the band crossing \cite{heyd2003hybrid}. The energy cutoff of the plane wave is chosen as $30\%$ above of the maximum value recommended in the pseudopoetential library. The Brillouin zone is sampled by a $\Gamma$-centered grid with the resolution of $2\pi \times 0.03 \text{ \AA}^{-1}$. All atoms are fully relaxed with a force criterion of 0.001~eV/{\AA}. Pre/Post-processing tools and utilities for solids computation \cite{wang2019vaspkit, herath2020107080, setyawan2010299, togo2018texttt, momma2008vesta} are employed.

Crystal structure searches for various compositions of electrides are performed with the minima hopping (MH) \cite{goedecker2004minima} algorithm with a generalized version  to periodic systems employing variable cell shape molecular dynamics \cite{amsler2010crystal,flores2020crystal}. This is a well-tested methodology that aims at finding the global minimum of a complex condensed matter system \cite{flores2020perspective}. For crystals with a small number ($<10$) of atoms per primitive cell, the supercells composed of two even three primitive cells are included as well to extend the searching space. As aforementioned, not only the lowest energy configuration, but also neighboring local minima are evaluated for their potential as electrides. In order to build convex hulls for final candidates, stable phases recorded in Materials Project \cite{Jain2013Materials, Ong2008Li} are adopted. The formation energy is estimated with respect to the corresponding stable elemental polymorph, where the {\sc vasp} pseudopotential C\_$\,$h is employed for all carbon-based systems, as some $A_n$C compounds contain $C_2$ dimers; the plane wave energy cutoffs are 875.0 eV and 320.6 eV for C- and Ge- based materials, respectively. To assess the lattice dynamics, we use the supercell approach, on which the total energy and pattern of displacements (within the harmonic approximation) are evaluated with {\sc vasp}. The harmonic interatomic force constants are built and diagonalized using the {\sc phonopy} software \cite{togo2015first}.

We utilize a recently developed scheme termed CMP+SDFT \cite{huebsch2020benchmark}, that is based on a combination of the cluster multipole (CMP) theory \cite{PhysRevB.95.094406, PhysRevB.99.174407} and calculations in the framework of spin-density functional theory (SDFT), in order to explore the possible stable magnetic states for Gd-based candidates of interest. The CMP basis is an orthogonal basis set of magnetic multipole configurations that can be used to expand an arbitrary magnetic structure. By an analysis of experimental data, it has been shown that the most stable magnetic configurations in nature are linear combinations of only few CMPs. Hence, an exhaustive list of candidate magnetic configurations can be generated with the initial magnetic moments generated by CMPs and their combinations. Subsequently, SDFT calculations are exploited to seek the stable configurations. In the case that the FM state is a local minimum instead of a global one, the energy difference higher than the global minimum is employed to evaluate the energetic stability of the FM structure. We point out again that the magnetic ground state is determined only for a subset of designed magnetic electrides.

\begin{figure*}[htp]
	\centering
	\includegraphics[width=0.98\linewidth]{./figures/overview.pdf}
	\caption{A statistical overview of predicted materials. (a)-(f) Convex hull diagrams of binary compounds $A_nB$ ($A=$ Y, La, Gd; $B=$ C, Ge). Electride candidates are colored by red, and distinguished by the attached numbers that represent their energy rankings. Notice that the number 1 means the corresponding phase has the lowest formation energy merely within the predicted electride structures; there could be a non-electride structure that is more stable. Grey dots stand for the other (meta)stable structures predicted by the CSP, which showing no evidence of electrides or have relatively high formation energy. No grey dots in diagrams of the Gd systems, since the Gd-based structures are obtained by chemical substitution instead of the CSP. (g) Stability (measured by $E_\text{hull}$) distribution of 58 designed electrides. (h) New stable polymorphs ($E_\text{hull} = 0$) of each stoichiometric category discovered by the CSP. Non-electride systems are also presented.}
	\label{overvew}
\end{figure*}

Along the high-symmetry line, crystalline symmetries of energy bands are determined via \texttt{irvsp}, a software package that can obtain irreducible representations of the electronic Bloch states \cite{gao2020irvsp}. Tight-binding models in Wannier representation \cite{wu2018wanniertools} are constructed through Wannier functions without the iterative maximal-localization procedure \cite{PhysRevB.56.12847, PhysRevB.65.035109, pizzi2020wannier90}. Considering the non-nucleus-bound electronic states are mostly \textit{s}-like \cite{matsuishi2003high}, we set the \textit{s} orbital as the trial projection at interstitial sites.  

\subsection{Electrides identification}
The electron localization function (ELF) \cite{becke1990simple, savin2005electron} is one of the mature theoretical descriptors of electrides \cite{liu2020electrides, zhu2019computational, zhao2016electride, dale2018theoretical, wang2017exploration}. The emergence of ELF local maxima (irregardless of its magnitude \cite{zhu2019computational}) at the vacant region indicates a potential electride. For the further identification of candidates of interest, we do the fatband analysis, namely, projecting eigenstates onto the \textit{s} orbital of empty spheres (pseudoatoms) located at the vacancy positions with ELF maxima. The wavefunction on the interstitial-state dominating bands as well as partial charge density close to the Fermi energy represent additional proof.
\vspace{1.5cm}

\section{Results}

\subsection{Overview of predicted electrides}
\label{subsec_overview}

\begin{figure*}[tbp]
	\centering
	\includegraphics[width=0.75\linewidth]{./figures/LaC_bulk.pdf}
	\caption{Electronic properties from the bulk calculation of LaC. (a) Primitive cell of LaC. Yellow bubbles represent the ELF. For a clear visualization, the ELF around atoms is not displayed. (b) The Brillouin zone of LaC. The high-symmetry points and the path chosen are highlighted in red. The $y$-direction projected Brillouin zone is indicated too. (c) Projected band structure of LaC along the path denoted in (b). The Fermi energy is set 0. As the color changes from red to blue, the contribution of interlayer states diminishes. (d) Nodal line of LaC in the reciprocal space. The left and right panels are viewed from \textit{x} and \textit{y} directions respectively. (e) The Zak phase integrated along the \textit{y} direction. Green and white separately stand for $\pi$ and 0 values. High-symmetry points in this reduced reciprocal space are also shown. (f)(g) Real branches of wavefunctions with the corresponding positions indicated in (c). The yellow and cyan mean positive and negative parts, respectively; red and blue stand for their cross-sections.}
	\label{LaC_bulk}
\end{figure*}

Within the 4 stoichiometric categories of binary compounds, $A_2$C, $A$C, $A_2$Ge and $A$Ge ($A=$ Y, La, Gd), 58 stable or metastable phases of potential electrides are identified. Among the former three categories, 56 of them are identified for the first time and hence, have not been recorded previously in the Inorganic Crystal Structure Database (ICSD) \cite{bergerhoff1983inorganic}. In $A$Ge no electride appears within the range of $E_\text{hull} < 0.20$ eV/atom, implying that Ge is less favorable for electride phases than C in the rare-earth based systems, which might be due to the smaller electronegativity. Figure \ref{overvew} (a)-(f) illustrate the distribution of predicted stable or metastable structures in convex hull diagrams. The potential electrides are highlighted in red, ranked and labeled according to their formation energies from low to high. The (meta)stable phases that are disregarded because of either showing no evidence to be electrides or yielding a high formation energy compared to other structures are marked in grey. As given in Figure \ref{overvew}(g), all the 58 electride structures have the $E_\text{hull}$ smaller than 0.14 eV/atom, and 49 ($84.5\%$) of them are below 0.10 eV/atom, revealing the possible synthesizability of these candidate structures \cite{oganov2019structure}.

Figure \ref{overvew}(h) lists the the newly discovered materials corresponding to the most stable phase ($E_\text{hull} = 0$). Two non-electride systems are attached too for reference. In the category $A_2$C, being the ground state, La$_2$C with the anti-CdCl$_2$ structure is predicted to be a promising electride. Note that the identical structure is shared by the stable phases of the Y$_2$C and the FM Gd$_2$C, while the latter two systems have been experimentally confirmed as stable 2D electrides \cite{zhang2014two, lee2020ferromagnetic}. Additionally, the FM configuration is verified to be the energetically favorable magnetic state of Gd$_2$C according to CMP+SDFT, in agreement with experiments \cite{lee2020ferromagnetic}. Therefore, the stable phases of $A_2$C represent a benchmark of our design scheme, indicating the repeatability of our CSP strategy. Moreover, the presumption of chemical substitution in Sec. \ref{subsec_design_scheme} is demonstrated to hold, i.e. stable or metastable Gd-based structures could be derived from those of Y and La.

We proceed by introducing the other stable structures in Figure \ref{overvew}(h). In the stoichiometric family $A$C, the possible electrides LaC with the space group $I4_1/amd$ (No. 141) and the FM GdC with the space group $C2/m$ (No. 12) are revealed to be the ground states, neither of which was reported previously. For the category $A_2$Ge, we find an unreported global minimum structure La$_2$Ge that belongs to the space group $Amm2$ (No. 38), though it is unlikely to be an electride. Finally, in the category $A$Ge, despite the absence of electride possibility amongst the most stable phases, a new stable polymorph of YGe is found, sharing the same structure with the known orthorhombic LaGe. Crystal structures of these stable materials are presented in Appendix \ref{sec_stable_structures}.

\begin{figure}[tbp]
	\centering
	\includegraphics[width=0.95\linewidth]{./figures/LaC_surf.pdf}
	\caption{Electronic structure of a 15-layer $(010)$ LaC slab. (a) Electronic band structure of the LaC slab. High-symmetry points in the reduced Brillouin zone are given by Figure \ref{LaC_bulk}(b) or \ref{LaC_bulk}(e). Energy is measured from the Fermi level. The red indicates a large proportion of floating electronic states that localize on the top and bottom of the slab. (b) The squared wavefunction of a surface state localized around the top layer. The corresponding point in the band structure is marked by a green circle in (a).}
	\label{LaC_surf}
\end{figure}

Relative to the corresponding stable phases, metastable phases can exhibit superior physical and chemical characteristics \cite{sun1600225thermodynamic}. In the subsequent discussion, we shall focus on two metastable systems to exemplify the crystal structures and electronic properties of the new electride candidates. One system (denoted as the $A$C system) is represented by the monoclinic LaC with $E_\text{hull} = 4.0$ meV/atom, being a promising 2D electride. The other (denoted as the $A_2$Ge system) is based on an orthorhombic structure, found in Y$_2$Ge and La$_2$Ge concurrently. These two electride systems will be thoroughly discussed in Sec. \ref{subsec_AC} and Sec. \ref{subsec_A2Ge}, respectively.

\subsection{The \textit{A}C system} 
\label{subsec_AC}

\begin{figure*}[htp]
	\centering
	\includegraphics[width=0.7\linewidth]{./figures/GdC_bulk.pdf}
	\caption{Electronic properties of GdC. (a) Primitive cell of GdC. Partial charge densities around the Fermi level ($-0.05 \sim 0.05$ eV) are shown for spin up (left) and spin down (right). Yellow bubbles represent the partial charge density, the green and blue are the cross sections. (b) Electronic band structure of FM GdC for spin up (top) and spin down (bottom). The high-symmetry points are given in (c). Arrows are indicative of the degeneracies between the highest VB and the lowest CB. (c) Nodal lines of FM GdC for the spin-up (left) and spin-down (right) channels. The Brillouin zone together with high-symmetry points are depicted too. (d) The Zak phase for spin up (left) and spin down (right), integrated along the \textit{y} direction, i.e. $\Gamma$-$X$ direction. Green and white mean $\pi$ and 0 values, respectively.}
	\label{GdC_bulk}
\end{figure*}

Starting from Y$_2$Br$_2$C$_2$ \cite{simon1996supraleitung}, we predict a metastable nonmagnetic crystal LaC , with the formation energy $-0.18$ eV/atom, $4.0$ meV/atom above the convex hull, showing no imaginary phonon modes (see Appendix \ref{sec_phonon}). This is an original material not reported in the database, whose crystal structure is shown in Figure \ref{LaC_bulk}(a). LaC crystallizes in a monoclinic lattice with the space group $P2_1/c$ (No.14), and is composed of four La atoms and four C atoms in a primitive cell. Aside from the time reversal ($\mathcal{T}$) symmetry, the following two symmetries are essential: an inversion $\mathcal{P}$ whose center is at the corner of the unit cell, and an \textit{a}-axis twofold screw operation $\mathcal{S}_{2[100]}$. Indicated by the ELF, LaC is a promising 2D electride, where excess electrons are trapped between cationic layers [Figure \ref{LaC_bulk} (a)]. 

We calculate the electronic band structure along the high-symmetry path given in Figure \ref{LaC_bulk}(b), in the absence of spin-orbit coupling (SOC). Entangled bands around the Fermi level reveal the metallic behavior of LaC. Placing pseudoatoms with the Wigner-Seitz radius 1.50 \AA$\ $at positions with local ELF maxima, i.e. fractional coordinates (0, 0.5, 0) and (0.5, 0.5, 0.5), yields the contribution of anionic electrons to the Bloch wavefunctions. One can notice that interstitial states are closely related to bands around the Fermi energy $E_\text{F}$ and show conducting behavior. The number of valence electrons in this calculation is 60. Hence, the 30 lower bands are counted as valence bands (VBs), while the other higher bands are regarded as conduction bands (CBs). In Figure \ref{LaC_bulk} (f)(g), the real parts of two eigenstates affiliated to the highest VB and the lowest CB are depicted respectively, of which the corresponding positions are marked in the band structure of Figure \ref{LaC_bulk}(c), close to a gapless point. Non-nuclear orbitals are observed for the VB, further confirming the electride characteristics. Nevertheless, this VB claimed by the projection analysis to have little portion of interstitial states, also contains non-nuclear orbitals according to the wavefunction visualization [Figure \ref{LaC_bulk}(g)]. Such contradiction originates from the fact that those non-nuclear orbitals are accidentally distributed in the uncovered region of pseudoatoms. The projection with empty spheres therefore underestimates the contribution from anionic electrons.

Next we move to the topic of band topology. The gapless point between a VB and a CB, mentioned in the last paragraph, is located along the $A$-$X$ line that is invariant under the $\mathcal{S}_{2[100]}$ operation. This degeneracy is protected by the opposite $\mathcal{S}_{2[100]}$ eigenvalues of the two crossing bands, indicated by \texttt{irvsp}. Additionally, the VB and CB remain gapless under the HSE06 hybrid functional. Owing to the $\mathcal{PT}$ symmetry in this spinless system \cite{chan20163, fang2016topological}, this degenerate point extends to general \textit{k}-points and forms a snake-shape nodal line [Figure \ref{LaC_bulk}(d)]. In the manifold of the Brillouin zone, the disconnected nodal line segments should constitute a whole loop. It is noteworthy that this nodal line partly comes from non-nucleus-bound electronic states (e.g. the band crossing in the $A$-$X$ direction), which are insensitive to the relativistic effect of the nuclei field, therefore might exhibit a negligibly small gap even if SOC is present \cite{zhang2019topological}. 

The nodal line here is protected by the $\pi$ Berry phase, defined as
\begin{equation}
	\label{Zak_phase}
	\gamma(\boldsymbol{k}_{\parallel}) = \sum_{n}^{\text{occ.}} \int_{0}^{|\boldsymbol{K}|}\mathrm{d}k_{\perp}\ i\langle u_{n\boldsymbol{k}} | \nabla_{k_{\perp}} | u_{n\boldsymbol{k}} \rangle
\end{equation}
for an arbitrary wave vector $\boldsymbol{k}$, where $\boldsymbol{K}$ stands for the reciprocal lattice vector, $k_{\perp}$ is the magnitude of $\boldsymbol{k}$ component along $\boldsymbol{K}$ and $\boldsymbol{k}_{\parallel}$ is the other component perpendicular to $\boldsymbol{K}$; $|u_{n\boldsymbol{k}}\rangle$ represents the cell-periodic Bloch function, with the gauge $|u_{n,\boldsymbol{k}+\boldsymbol{K}}\rangle = e^{-i\boldsymbol{K}\cdot\boldsymbol{r}}|u_{n\boldsymbol{k}}\rangle$. Note that the Berry phase specified in Eq. (\ref{Zak_phase}), with the integral path along a reciprocal lattice vector, is also known as the Zak phase \cite{PhysRevLett.62.2747}. Making $\boldsymbol{K}$ parallel with the \textit{y} direction of the reciprocal space, yields the Zak phase distribution displayed in Figure \ref{LaC_bulk}(e), which is quantized to be 0 or $\pi$ (mod $2\pi$) as a result of the $\mathcal{P}\mathcal{T}$ symmetry. The nontriviality of this nodal line is demonstrated by the fact that the profile of $\pi$ value region of the Zak phase just matches the nodal line projection along the \textit{y} direction [Figure \ref{LaC_bulk}(e) and the right panel of Figure \ref{LaC_bulk}(d)]. LaC is therefore a TNLE \cite{hirayama2018electrides}.

While underlying the fingerprint of nodal lines in the bulk, the Zak phase also discloses information on the surface. To elaborate this point, we construct a 15-layer (010) slab with the electron layer terminated, and compute the electronic band structure using DFT, as presented in Figure \ref{LaC_surf}(a), whose high-symmetry points in the surface Brillouin zone are given by Figure \ref{LaC_bulk}(b) or \ref{LaC_bulk}(e). Contribution from excess electrons localized on the surfaces is highlighted by projecting eigenstates onto pseudoatoms on the top and bottom ends of the slab. At the wave vector $\boldsymbol{k}_{\parallel}$ with $\pi$ Zak phase [e.g., $\bar{\Lambda}_1$ in Figure \ref{LaC_bulk}(b)], the bulk polarization is nonzero and quantized, and the topological charge $e/2$ (mod $e$) would be induced at the slab surface by this bulk polarization \cite{PhysRevB.48.4442, hirayama2017topological}. As the two surfaces of the slab are equivalent on account of the $\mathcal{P}$ symmetry, the total charge at $\bar{\Lambda}_1$ is different from that in the region with 0 Zak phase [e.g., $\bar{Z}$ in Figure \ref{LaC_bulk}(b)] by \textit{e}. The charge difference can be evidenced by the slab band structure. Along the $\bar{Z}\rightarrow\bar{\Gamma}$ direction, one unoccupied state traverses the bulk line node and merges into the valence spectrum in the end, as marked by the blue circle in Figure \ref{LaC_surf}(a). As a result, in comparison with $\bar{Z}$, there exists one more VB at $\bar{\Lambda}_1$, being responsible for the charge difference \textit{e}. In addition to the surface charge, midgap states emerging from the line node are observed. On the left side of the line node, the midgap states are dominated by interstitial orbitals. In Figure \ref{LaC_surf}(b), visualization of a midgap state in real space unveils the fact that this boundary state is not tightly attached to any atoms, but floats on the surface instead. Such floating surface states are also predicted in other TNLEs like Y$_2$C \cite{hirayama2018electrides}.

\begin{figure*}[htbp]
	\centering
	\includegraphics[width=0.75\linewidth]{./figures/Y2Ge_bulk.pdf}
	\caption{Electronic properties from the bulk calculation of Y$_2$Ge. (a) Primitive cell of Y$_2$Ge. Yellow bubbles are ELF. The isosurface value in the bottom panel is smaller than those in the top and middle panels. Connected interstitial states are highlighted by a dashed red circle in the bottom panel. (b) Projected band structure of Y$_2$Ge with the high-symmetry points denoted in (c). The Fermi energy is set 0. As the color changes from red to blue, the contribution of interstitial states decays. Some gapless points are indicated by arrows. (c) Nodal lines of Y$_2$Ge accompanied by the Brillouin zone. Two inequivalent groups of nodal lines are labels as NL1 and NL2, respectively. Spinless Weyl points are also displayed, with the orange and magenta representing the monopole charge $+1$ and $-1$ separately. (d)(e) show the real branches of the wavefunctions at the positions indicated in (b); the positive part is in yellow and the corresponding cross-sections are colored by red; the cyan and blue are for the negative part.}
	\label{Y2Ge_bulk}
\end{figure*}

We now discuss the crystal after chemical substitution Y $\rightarrow$ Gd, namely, the isostructural GdC. As described in Sec. \ref{subsec_design_scheme}, we have restricted our discussion to the FM configuration of Gd atoms, in which case the interstitial electrons can be easily magnetized. For the GdC here, although an antiferromagnetic (AFM) configuration is predicted by CMP+SDFT to be the magnetic ground state, it is not highly energetic-favorable compared to the FM state (see Appendix \ref{sec_ground_mag}). This makes our 'FM restriction' more plausible. GdC possesses a formation energy equal to $-0.22$ eV/atom and $E_\text{hull}= 0.032$ eV/atom. To look into the electronic structure, the Hubbard $U$ for Gd 4\textit{f} orbitals is included, with $U - J = 6$ eV and $J = 0.7$ eV \cite{PhysRevX.4.031023}; other values ($U - J = 3$ eV and 9 eV) are also checked, showing that only the dense flat $4f$ bands shift while other bands near the Fermi level remain unchanged. Figure \ref{GdC_bulk}(a) presents the partial charge density of the FM GdC, in the vicinity of $E_{\text{F}}$ ($-0.05 \sim 0.05$ eV) for two spin channels. While only intralayer electrons can be seen for spin up, excess electrons rest in the interlayer gap for spin down. Therefore, spin-polarized interstitial states are achieved in this FM system, although the macroscopic magnetism is dominated by Gd atoms. Figure \ref{GdC_bulk}(b) gives the electronic band structure of the FM GdC without SOC. One can notice that bands of two spin channels resemble those in LaC, as a result of the close electronic configuration between Y and Gd, and the nearly identical crystal structure. Such a similarity has been also observed in the electride systems Y$_2$C and Gd$_2$C \cite{PhysRevX.4.031023, liu2020ferromagnetic}. Hence, we expect GdC exhibits a similar electronic behavior to LaC, especially regarding the band topology. Band degeneracies around the Fermi level are displayed in Figure \ref{GdC_bulk}(b), all of which are preserved at the level of HSE06 functionals. In parallel with LaC, these degenerate points turn out to extend in the reciprocal space, forming nodal rings for both spin up and down [Figure \ref{GdC_bulk}(c)]. All the degeneracies still exist in the HSE06 band structure. The nontriviality is again verified by the Zak phase [Figure \ref{GdC_bulk}(d)]. In spite of the broken $\mathcal{T}$ symmetry, the remaining $\mathcal{P}$ symmetry still protects the $\pi$ Berry phase nodal lines. This time reversal breaking TNLE state was also reported in Gd$_2$C \cite{liu2020ferromagnetic}.

\subsection{The \textit{A}$_2$Ge system}
\label{subsec_A2Ge}

\begin{figure*}[htbp]
	\centering
	\includegraphics[width=0.7\linewidth]{./figures/Gd2Ge_bulk.pdf}
	\caption{Electronic properties of Gd$_2$Ge. (a) Primitive cell of Gd$_2$Ge. Partial charge densities around the Fermi level ($-0.05 \sim 0.05$ eV) are presented for spin up (left) and spin down (right). The yellow bubbles represent the partial charge density, and the green means the cross-section. (b) Electronic band structure of FM GdC for the spin-up (top) and spin-down (bottom) parts. The high-symmetry points are indicated in (c). Degenerate points between the highest VB and lowest CB are marked by arrows. (c) Nodal lines of FM GdC for spin up (left) and spin down (right). The Brillouin zone, together with high-symmetry points, are depicted too. The blue or green distinguishes nodal lines that are distributed on the vertical or horizontal mirror plane. Weyl points are also displayed, with the orange and magenta representing the monopole and antimonopole, respectively.}
	\label{Gd2Ge_bulk}
\end{figure*}

Y$_2$Ge is another new inorganic electride, derived from the halide La$_2$I$_2$Ge\cite{mattausch2005reduced}, with the structure metastability supported by the negative formation energy $-0.62$ eV/atom ($0.13$ eV/atom higher than the hull) and the absence of imaginary phonon frequencies (see Appendix \ref{sec_phonon}). The base-centered orthorhombic lattice of Y$_2$Ge belongs to the space group $Amm2$ (No. 38). The same structure emerges as a metastable phase of La$_2$Ge too. Figure \ref{Y2Ge_bulk}(a) depicts the primitive cell of Y$_2$Ge that contains nine Y atoms and three Ge atoms. The symmetry operations are: a twofold rotation $\mathcal{C}_{2[0\bar{1}1]}$, two mirror reflections $\mathcal{M}_{[011]}$ and $\mathcal{M}_{[100]}$, whose axes are perpendicular to each other. Notice that the inversion symmetry $\mathcal{P}$ is absent. The spatial distribution of the ELF is given in Figure \ref{Y2Ge_bulk}(a). Two ELF maxima related by the $\mathcal{C}_{2[0\bar{1}1]}$ or  $\mathcal{M}_{[011]}$ symmetry exist in the cationic cavities (the top and middle panels). These two interstitial sites are not isolated but connected with each other, which is uncovered in the bottom panel with a smaller isosurface value of ELF. Now that these interstitial states are connected through the one-dimensional (1D) channel, we classify this candidate as a 1D electride. Set pseudoatoms with the Wigner-Seitz radius 0.5 \AA$\ $at the interstitial locations with local maxima of the ELF, yielding the projected band structure of Y$_2$Ge (without SOC) given by Figure \ref{Y2Ge_bulk}(b), in which the high-symmetry points are denoted in Figure \ref{Y2Ge_bulk}(c). Eigenstates around the Fermi level mostly emanate from Y $4d$ and interstitial \textit{s} orbitals. Again the number of VBs is selected to match the half number of valence electrons in the primitive cell.

Along the $T$-$Z$-$\Gamma$-$Y$ path, the gap between the highest VB and the lowest CB disappears at four positions, which have been labeled by arrows NL1 or NL2 in Figure \ref{Y2Ge_bulk}(b). While the crossings marked by NL1 on the $T$-$Z$ and $Z$-$\Gamma$ paths are confirmed by the band structure computed with the HSE06 hybrid functional, the other two crossings NL2 on the $\Gamma$-$Y$ path disappear when employing HSE06 (see Appendix \ref{sec_hse_bands}). For completeness, we take the latter two degenerate points into account too. Since the $T$-$Z$-$\Gamma$-$Y$ line is located on the $\mathcal{M}_{[011]}$-invariant plane $(011)$, each state on this line has a well-defined eigenvalue $\pm1$ of $\mathcal{M}_{[011]}$, as verified by \texttt{irvsp}. It turns out that every pair of crossing bands have opposite mirror eigenvalues, protecting the degeneracy. In Figure \ref{Y2Ge_bulk}(d)(e), two eigenstates---one is reflection antisymmetric and the other is symmetric with respect to the (011) plane---are depicted. Their corresponding positions have been indicated in Figure \ref{Y2Ge_bulk}(b), belonging to a CB and a VB, respectively. Such opposite crystalline symmetries guarantee the band crossing in the $T$-$Z$ direction. Moreover, one of the crossing bands is dominated by anionic electrons [Figure \ref{Y2Ge_bulk}(e)], again suggesting the significant role of interstitial states for band inversion in electrides. Figure \ref{Y2Ge_bulk}(c) reveals the nodal lines stemming from those gapless points are confined within the mirror plane $(011)$. The big ring (labeled as NL1) crosses the boundary of the Brillouin zone; two small rings (labeled as NL2) are the inversion of each other, associated by the $\mathcal{T}$ symmetry. They are related to the band crossings NL1 and NL2 in Figure \ref{Y2Ge_bulk}(b), respectively. As the degenaracies on the $\Gamma$-$Y$ line disappear at the HSE06 level, the small nodal rings NL2 might be absent under an experimental observation. We then implement the $\mathbb{Z}$-type reflection invariant to characterize the mirror-symmetry protected nodal lines \cite{PhysRevB.90.205136, fang2016topological}. Within the mirror plane, we pick two points $p_a$ and $p_b$ on the two sides of the nodal line, at which the VBs and CBs are separated. Then count the number of VBs with the mirror eigenvalue $+1$, denoted by $N_{a,b}$. The reflection invariant reads
\begin{equation}
	\label{Z-invariant}
	\zeta = N_a - N_b.
\end{equation}
A nonzero $\zeta$ confirms the nontrivial nodal line, which means that the band inversion cannot be adiabatically eliminated without breaking the mirror symmetry. In our case, we obtain $\zeta_1 = \zeta_2 = 1$ for the big and small nodal rings. As a result, the stability of the line nodes in Y$_2$Ge is protected by the reflection symmetry $\mathcal{M}_{[011]}$.

Aside from the line nodes, there appear six pairs of point nodes (spinless Weyl points) owing to the absence of $\mathcal{P}$ symmetry. We compute the monopole charge (also known as the chirality) of each node using the Berry flux across an enclosed surface accommodating the node \cite{doi:10.7566/JPSJ.87.041002}. Figure \ref{Y2Ge_bulk}(c) depicts the positions and monopole charges of the Weyl points. Each pair of Weyl points are linked by the $\mathcal{T}$ symmetry, possessing the same monopole charge.

Next, we replace the Y with Gd while retaining the crystal structure, obtaining a metastable magnetic electride whose formation energy is $-0.58$ eV/atom and $E_\text{hull}=0.12$ eV/atom. Figure \ref{Gd2Ge_bulk}(a) displays the partial charge density distribution in the primitive cell, within the energy range $-0.5\text{ eV} < E-E_{\text{F}} < 0.5\text{ eV}$. Similar to GdC, interstitial states prefer the minority spin (spin down) in Gd$_2$Ge. That is, the magnetization opposite to Gd is induced among anionic electrons. Utilizing the GGA+U method, we show the electronic band structure in Figure \ref{Gd2Ge_bulk}(b) without SOC. Compared to the bands of Y$_2$Ge, the spin-up bands are lowered and the spin-down bands are lifted as expected. Degenerate points along high-symmetry lines are labeled by arrows. All the degenaracies are confirmed by the computation using the HSE06 functional, except for the right two band crossings on the $T$-$Z$ path of the spin-down band structure, which are replaced by one crossing (see Appendix \ref{sec_hse_bands}). These gapless points are guaranteed by the mirror symmetries $\mathcal{M}_{[011]}$ and $\mathcal{M}_{[100]}$. Figure \ref{Gd2Ge_bulk}(c) presents the nodal lines generated from the band crossings within the mirror planes. In contrast to Y$_2$Ge, the $k_x=0$ (spin up)  and $k_x=\pi$ (spin down) planes also hold nodal lines (highlighted in green), allowed by the $\mathcal{M}_{[100]}$ symmetry. Additionally, the perpendicular mirror planes give rise to two sets of nodal chains in the reciprocal space of spin down. In parallel with Y$_2$Ge, Weyl points are observed for both spin channels [Figure \ref{Gd2Ge_bulk}(c)]. The intriguing properties arose out of the coexistence of Weyl points and nodal chains remain to be explored.

\section{Conclusion} 
In summary, we put forward a new electride design scheme, predicting electrides as the reduced form of halides. To demonstrate the success of this design strategy, three stoichiometric categories of binary compounds, $A_2$C, $A$C, $A_2$Ge ($A$ = Y, La, Gd), with 56 newly discovered stable or metastable phases derived from ternary halides, are identified as promising electrides based on first-principles crystal structure prediction and chemical substitution. Of all the candidate structures, the energies above the convex hull are smaller than $0.14$ eV/atom, indicating potential synthesizability. Benefiting from the similarities of atomic and chemical properties of rare-earth elements, stable or metastable Gd-based systems are generated by element substitution from the Y- and La-based structures so that magnetism can manifest itself. Subsequently, we study the ferromagnetic case, a situation where the anionic electrons can easily be magnetized. Two electride systems, the monoclinic $A$C and the orthorhombic $A_2$Ge, have been discussed in detail. According to our results for $A=$ Gd, their magnetic states are suitable for spin injection. The topological phases investigated show that both of the two designed electride systems manifest topological nodal lines in the reciprocal space, which cements the close relation between topological materials and electrides because of the interstitial states around the Fermi level. 

Our strategy to design magnetic and topological electrides can be easily generalized to other sets of materials.  For example, one could expand the searching space to include hydrides, oxides, nitrides, etc. accessible through material databases in combination with high-throughput techniques \cite{setyawan2010high-throughput}. The replacement of other lanthanides in host systems is likely to open new avenues to synthesize magnetic electrides. One future task would be collaborating with experimentalists to synthesize the designed electrides, and devising incorporation into spintronic devices.

\begin{acknowledgments}
M. H. acknowledges financial support from JSPS KAKENHI Grants No. 20K14390. J. A. F.-L. acknowledges computational resources provided by the Swiss National Supercomputing Center withing the project s970. T. N. acknowledges financial support from JST PRESTO (Grant No. JPMJPR20L7). M.-T. H. and R. A. were supported by a Grant-in-Aid for Scientific Research (No. 19H05825, and No. 16H06345) from the Ministry of Education, Culture, Sports, Science and Technology.	
\end{acknowledgments}

\newpage
