Magnetic Ordering, Anomalous Lifshitz Transition and Topological Grain Boundaries in Two-Dimensional Biphenylene Network

Graphitic materials have several notable electronic structures such as the Dirac fermionic states in a wide range of energy-momentum space Wallace (1947); Charlier et al. (2007); Castro Neto et al. (2009). These properties are realizable partly owing to the large hopping energy, hexagonal crystal symmetry and negligible spin-orbit interaction (SOI) Charlier et al. (2007); Castro Neto et al. (2009). Recently, a new network of carbon atoms with nonhexagonal rings has been synthesized Fan et al. (2021). The new planar carbon allotrope or biphenylene network (BPN) was grown into the ribbon shape on the Au(111) surface of which electronic properties are dependent of its width Fan et al. (2021). Although they are still one-dimensional, its extension to two-dimensional (2D) form like graphene is surely expected. Considering that carbon materials have been fruitful playgrounds for many fundamental physical concepts Charlier et al. (2007); Castro Neto et al. (2009), we anticipate that BPN could also have potentials to demonstrate new physical phenomena.

In this work, we show that, unlike any other carbon material, a pristine 2D BPN has a saddle shaped energy band at the Brillouin zone (BZ) center and peculiar type-II Dirac fermion. The type-II state here has a nearly flat branch and is shown to support topologically protected grain boundary state. With moderate uniaxial strains, two Dirac cones are driven to merge with the saddle point, thus invoking simultaneous Lifshitz transitions for the van Hove singularity (vHS) and type-II Dirac fermions. Moreover, including effects of Coulomb interactions, the ground state turns out to be ferrimagnetic. Considering control techniques developed for 2D materials Novoselov et al. (2016), we expect the new carbon allotrope Fan et al. (2021) could be one of quintessential materials to study interplays between topological and correlated properties.

We studied electronic properties of BPN based on first-principles methods Giannozzi et al. (2009) and effective models. The generalized gradient approximation (GGA) Perdew et al. (1996) was used for the exchange-correlation functional. A plane-wave basis set with a cutoff energy of 100 Ry is adopted and the ultrasoft pseudopotential Garrity et al. (2014) is used. To consider Coulomb interaction beyond the GGA, we used a newly developed DFT+$U$+$V$ method exploiting the self-consistent evaluation Agapito et al. (2015); Lee and Son (2020); Tancogne-Dejean and Rubio (2020) of the on- and inter-site Hubbard interactions ($U$ and $V$) Campo Jr and Cococcioni (2010). We consider $V$ between the nearest neighbors (n.n.), that was confirmed to achieve desired convergence Lee and Son (2020). The method enable us to include Coulomb interaction in a similar level of $GW$ approximations with a significantly reduced computational cost Lee and Son (2020); Tancogne-Dejean and Rubio (2020); Campo Jr and Cococcioni (2010); Timrov et al. (2021). We note that this method captures the effects of Coulomb interaction in 2D crystals such as black phosphorous Lee and Son (2020) and graphene Tancogne-Dejean and Rubio (2020) very well.

The 2D BPN lattice in Fig. 1(a) has the point group $D_{2h}$ with 4-, 6- and 8-membered rings. It can be regarded as an effective lattice where a square shaped cluster consisting of atoms indexed by $1,\cdots,4$ in Fig. 1 (a) forms a rectangular lattice and a dimer of atoms indexed by 5 and 6 is its basis. By assigning an effective molecular orbital with specific angular momentum for each unit as a basis set, we can describe the low energy electronic states of the BPN lattice. This may offer a unique chance to support interesting correlated states Hlubina et al. (1997); Furukawa et al. (1998); Honerkamp and Salmhofer (2001); Irkhin et al. (2001) as well as symmetry protected energy bands Young and Kane (2015); Kim et al. (2015); Soluyanov et al. (2015); Yan et al. (2017); Armitage et al. (2018).

All structures are fully relaxed and dynamical stability of free standing BPN is confirmed by checking its stable phonon dispersion using density functional perturbation theory Baroni et al. (2001); Giannozzi et al. (2009), agreeing with a recent study Veeravenkata and Jain (2021). Detailed information on crystal structure is in the caption of Fig. 1. Although not all the local bonds are ideal hexagonal $sp^{2}$ bonds, the states are completely decomposed into the $\pi$ and $\sigma$ bonds, as shown by projected local density of states (DOS) in Fig. 1 (b). The non-ideal local $\sigma$ bonds manifest themselves in different charge states depending on the atomic positions. Our calculation shows that the Löwdin charge of the square unit in Fig. 1 (a) is smaller by 0.13$|e|$ ($e$ is an electron charge) per atom than one in the dimer unit, implying the higher site energy of the former than the latter.

We first highlight two novel features in energy bands of a single layer BPN. The saddle shaped band intersects with the almost flat band near $\frac{2}{5}\Gamma\text{X}$, forming a type-II Dirac state Soluyanov et al. (2015); Yan et al. (2017); Armitage et al. (2018). As shown in Fig. 1 (b), the nearly flat band is completely projected on the $p_{z}$-orbitals localized within the square unit. The almost flat dispersion indicates that the state is close to a transition between the type-I and II Dirac states Volovik and Zhang (2017); McCormick et al. (2017). The vHS at the BZ center is also uncommon since it is usually at the zone boundary Hlubina et al. (1997); Furukawa et al. (1998); Honerkamp and Salmhofer (2001); Irkhin et al. (2001); Kim and Son (2021). These topological bands should be maintained under strong Coulomb interactions unless ground state symmetries are broken. Like other graphitic materials Charlier et al. (2007); Castro Neto et al. (2009), however, the Coulomb interactions will alter energies and slopes of DFT-GGA bands.

Including Coulomb interactions beyond the GGA, we obtain a ferrimagnetic ground state. Our self-consistent evaluation of Hubbard interactions are $U=5.93~{}(5.98)$ eV and $V=3.05~{}(3.09)$ eV for $p_{z}$-orbitals in the square (dimer) unit. These values are comparable to $U$ and $V$ of graphene Wehling et al. (2011); Tancogne-Dejean and Rubio (2020). Unlike graphene where all atoms have the same Hubbard interactions, atoms belonging to different units of BPN lattice have different values. Moreover, diverging DOS near $E_{F}$ can invoke large exchange splitting. Our ab initio calculation indeed obtains a net magnetic moment of $0.33~{}\mu_{B}$ per unit cell where a moment per atom belonging to square (dimer) unit is $0.12~{}(-0.07)~{}\mu_{B}$. The total energy for the magnetic state is lower than nonmangetic one by 14.2 meV per atom. As shown in Fig. 1 (c), spin resolved band structures have four crossing points where two of them are quite close to the zone center and $E_{F}$. The energy level of vHS is also closer to the $E_{F}$ than DFT-GGA, as also can be seen in 2D Fermi surfaces for each case shown in Figs. 1 (d) and (e), respectively. We note that our result is consistent with several theoretical studies Hlubina et al. (1997); Furukawa et al. (1998); Honerkamp and Salmhofer (2001); Irkhin et al. (2001); Ziletti et al. (2015) on effects of interaction near vHS. While the strong repulsion within a single vHS alone may not be so effective in realizing spin ordering Hlubina et al. (1997), the half-filled vHS could realize a magnetic instability within the Stoner limit Ziletti et al. (2015). We also note that the computed work function of 4.24 eV with extended Hubbard interactions slightly differ from 4.32 eV with DFT-GGA.

As demonstrated in the other 2D single layer magnets Gong et al. (2017); Kim et al. (2019), our magnetic state may not be so stable in the limit of single layer owing to intrinsic fluctuations with negligible SOI here. The magnetic state may stabilize by an enhanced anisotropy from external fields or substrate. So, hereafter, we include the $U$ and $V$ without magnetism to take into account the effects of Coulomb interactions on band energies and group velocities like previous calculations for nonmagnetic materials Lee and Son (2020); Tancogne-Dejean and Rubio (2020).

Having discussed unique low energy bands, now we turn to their control. Like graphene Pereira and Castro Neto (2009); Choi et al. (2010) or black phosphorous Baik et al. (2015); Doh and Choi (2017), strain can alter electronic and topological properties of a single layer BPN. We consider uniaxial strains of $\varepsilon_{xx}~{}(\varepsilon_{yy})$ along $x(y)$-direction. With increasing $\varepsilon_{xx}(\varepsilon_{yy})$, Dirac points moves to the BZ boundary (center) [Figs. 2 (a)]. The energy difference ($\Delta^{\Gamma}_{\alpha\beta}$ in Fig. 2 (a)) between the nearly flat and saddle shaped bands at the $\Gamma$ point also increases (decreases) as $\varepsilon_{xx}~{}(\varepsilon_{yy})$ increases as shown in Fig. 2 (c). With increasing $\varepsilon_{yy}$, the Dirac point located on $\Gamma$X and its time reversal partner approach each other and eventually merge together at the $\Gamma$-point. Unlike typical merging Murakami et al. (2007); Son et al. (2011) or generation Baik et al. (2015); Doh and Choi (2017) of a pair of Dirac cones, they annihilate each other at the saddle point as shown in Figs. 2 (a) and (c). Thus, its DOS at the transition diverges unlike the usual semi-Dirac cone at the critical point Dietl et al. (2008); Son et al. (2011) as shown in Fig. 2 (d). The sign change of $\Delta^{\Gamma}_{\alpha\beta}$ near $\varepsilon_{yy}\simeq 9\%$ in Fig. 2 (c) signals gap opening at the zone center, occurring the Lifshitz transition. This transition is also accompanied with charge transfer between the square and dimer units. As shown in Fig. 2(e), the charge difference between two units ($\delta\rho_{L}$) changes its sign at the $\varepsilon_{yy}$ for gap opening. Thus, phenomenologically, we correlate the site energy or local charge reversal between two units with the topological transition.

With $U$ and $V$, all the variations by strain discussed above still hold. The band width increases by including inter-site Hubbard interaction of $V$ Lee and Son (2020); Tancogne-Dejean and Rubio (2020); Campo Jr and Cococcioni (2010) so that the Dirac points are closer to $E_{F}$ and $\Gamma$-point than ones from DFT-GGA, resulting in a smaller critical strain for the topological transition. So, as shown in Figs. 2 (b) and (c), the Lifshitz transition achieves with $\varepsilon_{yy}\simeq 2.5\%$ for DFT+$U$+$V$ bands while $\varepsilon_{yy}\simeq 9\%$ for DFT-GGA bands as also shown in sign change of $\Delta^{\Gamma}_{\alpha\beta}$ in Fig. 2 (c). Like $\delta\rho_{L}$ for DFT-GGA bands, the difference of on-site Hubbard repulsion between two units ($\delta U$) also changes its sign near the transition as shown in Fig. 2 (e). Unlike semimetallic nature of strained BPN within DFT-GGA calculations, we obtain indirect band gaps with the Hubbard interactions when $\varepsilon_{yy}>3\%$ [Fig. 2 (f)] owing to enhanced band repulsion by $U$. So, the annihilation of Dirac cones by strain induces metal-insulator transition with the maximum gap of 0.78 eV at $\varepsilon_{yy}=8\%$.

Having described variations in electronic energy bands, we briefly discuss control of magnetic moments under external perturbations. As shown in Fig. 2 (g), the magnetic moment quickly disappears as the energy gap develops with $\varepsilon_{yy}$ while it maintains the value with $\varepsilon_{xx}$. The variation of magnetic moments also can be achieved with external dopings. In Fig. 2 (h), it is shown that the moment decreases gradually with increasing hole doping while it remains with electronc doping since the large difference between hole and electron side density of states shown in Fig. 2 (d).

Hereafter, We discuss origins for the novel bands and their topological consequences with simpler models. All features found in our ab initio bands can be retrieved by tight-binding (TB) approximations with the $p_{z}$-orbitals. Since the fine-tuned hopping energy depending on n.n. distances only changes the minute details, we set them to be a same value of $-t_{0}$. A hopping energy of $-t_{d}$ along diagonal directions in the square unit is included to obtain realistic hopping energies. To account for $\delta\rho_{L}$ or $\delta U$ in Fig. 2 (e), the site energy of $\epsilon_{D}$ for the atoms in the dimer unit is indispensable. With these parameters, a TB Hamiltonian can be written as

$${\mathcal{H}}_{\text{TB}}=-t_{0}\sum_{\langle i,j\rangle}c^{\dagger}_{i}c_{j}-t_{d}\sum_{\langle\langle i,j\rangle\rangle\in{\mathcal{S}}}c^{\dagger}_{i}c_{j}+\epsilon_{D}\sum_{i\in{\mathcal{D}}}c^{\dagger}_{i}c_{i}+(\text{c.c.}),$$

(1)

where $\mathcal{S(D)}$ denotes a set of atoms belonging to the square (dimer) units and $\langle\cdots\rangle$ and $\langle\langle\cdots\rangle\rangle$ indicate the n.n. and next n.n. pairs, respectively. Comparing TB bands with DFT-GGA ones, we obtained $t_{0}=2.59$ eV, $t_{d}=0.45$ eV and $\epsilon_{D}=-0.95$ eV, respectively, and the TB bands agree with our DFT-GGA bands [Fig. 3 (a)]. Even without $t_{d}$, Eq. 1 can reproduce almost same bands and topological characteristics. For such a case, however, the fitted values are unrealistic. We also note that a long-range hopping between rectangular units along ${\bf a}_{1}$ direction makes the flat $\alpha$-band in Fig. 3 (a) be dispersive to match our DFT-GGA bands better.

The topological transitions also can be described by Eq. 1. Not by varying hoppings according to strain in detail, the transition can be simply controlled by sign of $\epsilon_{D}$. As shown in Fig. 3 (b), the $\epsilon_{D}$ following the sign of $\delta\rho_{L}$ or $\delta U$ in Fig. 2 (e) realises the transition: with $\epsilon_{D}<0$, the type-II Dirac point is on the flat band and with $\epsilon_{D}>0$, gap opens. When $\epsilon_{D}=0$ implying the same site energy for two sublattice units, the two tilted Dirac cones and saddle point merge at the zone center resulting in coexistence of critical states and vHS.

We find that wave functions (WFs) at the $\Gamma$-point for the vHS and type-II Dirac states are linear combinations of characteristic states in the square ($\varphi_{s}$ and $\varphi_{p}$) and dimer units ($\varphi_{s^{\prime}}$). The states are expressed with a linear combination of $p_{z}$ orbitals of atom $i$ ($\phi_{i}$) where $i=1,\cdots,6$ in Fig. 1(a) such that $\varphi_{s(p)}=\frac{1}{2}(\phi_{1}+\phi_{2}\pm\phi_{3}\pm\phi_{4})$ and $\varphi_{s^{\prime}}=\frac{1}{\sqrt{2}}(\phi_{5}+\phi_{6})$. These WFs for the units can be mapped onto effective orbitals as described in Fig. 3 (c) such that $\varphi_{s(s^{\prime})}$ is regarded as an effective $s$-orbital and $\varphi_{p}$ as $p_{y}$-orbital, respectively. Then we immediately notice that the hopping channel between neighboring $\varphi_{p}$ orbitals through $\varphi_{s^{\prime}}$ along $x$ direction is forbidden, realizing the flat band along $k_{x}$ direction and localization of its WFs within the square unit as shown in the ab initio bands [Fig. 1(b)]. We also can explain the zone-center vHS by considering right combination of hopping between $\varphi_{p}$ and $\varphi_{s^{\prime}}$ along $y$-direction and one between $\varphi_{s}$ and $\varphi_{s^{\prime}}$ along $x$-direction because their hopping energies should have opposite signs to each other.

The qualitative discussions above can be explicitly expressed by constructing effective models. For a basis set of $\Psi=(\varphi_{p},\varphi_{s},\varphi_{s^{\prime}})$ using the orbitals in Fig. 3 (c), a Hamiltonian for three effective orbitals can be written as

$$\mathcal{H}_{\text{orb}}({\bf k})=\begin{pmatrix}\varepsilon_{p}&0&-if({\bf{k}})\\ 0&\varepsilon_{s}&-g({\bf{k}})\\ if({\bf{k}})&-g({\bf{k}})&\varepsilon_{s^{\prime}}\end{pmatrix},$$

(2)

where $\hbar=1$, ${\bf k}=(k_{x},k_{y})$, $f({\bf{k}})\equiv 4t^{\prime}\cos\frac{k_{x}a_{1}}{2}\sin\frac{k_{y}a_{2}}{2}$, $g({\bf{k}})\equiv 4t\cos\frac{k_{x}a_{1}}{2}\cos\frac{k_{y}a_{2}}{2}$, and $\varepsilon_{l}$ $(l=s,s^{\prime},p)$ is a site energy for $\varphi_{l}$. $t(t^{\prime})$ is a hopping energy between $\varphi_{s}$ and $\varphi_{s^{\prime}}(\varphi_{p})$. Without loss of generality, we can set $\varepsilon_{p}=0$, $t=t^{\prime}$ and $\varepsilon_{D}\equiv\varepsilon_{s^{\prime}}-4t=\varepsilon_{s}-4t$, respectively. Then, it can be further reduced to a simpler Hamiltonian by projecting Eq. 2 onto a basis of $\Psi_{\text{eff}}=(\varphi_{p},\tilde{\varphi_{s}})$ where $\tilde{\varphi_{s}}=\frac{1}{\sqrt{2}}(\varphi_{s}+\varphi_{s^{\prime}})$. The resulting ${\mathcal{H}}_{\text{eff}}$ near $\Gamma$-point is

$$\mathcal{H}_{\text{eff}}({\bf k})\simeq\begin{pmatrix}-\frac{\delta_{y}}{2m_{y}}k_{y}^{2}&iv_{y}k_{y}\\ -iv_{y}k_{y}&\varepsilon_{D}+\frac{k_{x}^{2}}{2m_{x}}+\frac{k_{y}^{2}}{2m_{y}}\end{pmatrix},$$

(3)

where $\delta_{y}=\frac{2t}{\varepsilon_{D}+8t}$, $v_{y}=\sqrt{2}ta_{2}$, $m_{x}^{-1}=2ta_{1}^{2}$ and $m_{y}^{-1}=2ta_{2}^{2}$. By matching $\varepsilon_{D}=-\Delta^{\Gamma}_{\alpha\beta}\simeq\frac{2}{3}\epsilon_{D}$, and $t=t_{0}$, the energy bands from Eq. 3 reproduce the topological transitions for vHS and Dirac states obtained from Eq. 1 very well as shown in Fig. 3(b). The critical type-II Dirac Hamiltonian can also be obtained by projecting Eq. 2 to $\Psi_{\text{eff}}$ around the Dirac point of ${\bf k}_{D}=(\pm k_{D},0)$ such as

$$\mathcal{H}_{\text{crit}}(q_{x}\pm k_{D},q_{y})\simeq\pm v_{x}q_{x}\sigma_{0}\mp v_{x}q_{x}\sigma_{z}-v_{y}q_{y}\sigma_{y},$$

(4)

where $k_{D}=\frac{2}{a_{1}}\arccos(1+\frac{\varepsilon_{D}}{4t})$, $v_{x}=\frac{a_{1}}{2}[-\varepsilon_{D}(\varepsilon_{D}+8t)]^{\frac{1}{2}}$, $\sigma_{0}$ is an identity and $\sigma_{z(y)}$ Pauli matrix. Eq. 4 indeed satisfies the condition for the transition between types-I and II Dirac states Volovik and Zhang (2017); McCormick et al. (2017) and ensures the pseudospin structures for nontrivial topological properties Soluyanov et al. (2015); Yan et al. (2017).

Having discussed the effective Hamiltonian of BPN, we consider its topological consequences. Topological characteristics regarding the band crossing at the type-II Dirac point is revealed by calculating the Zak phase. Here we consider the Zak phase of each band along the $k_{y}$ direction with fixed $k_{x}$, $\gamma_{n}(k_{x})=i\int_{-\pi}^{\pi}dk_{y}\;\langle u_{n}({\bf{k}})|{\partial_{k_{y}}}|u_{n}({\bf{k}})\rangle$ where $u_{n}({\bf{k}})$ is the cell periodic function of the $n$-th Bloch state of Eq. 1. This quantity represents the Wannier center of the corresponding band along the $y$ direction and should be quantized by either 0 or $\pi$ (mod $2\pi$) due to the mirror symmetry of the system Lau et al. (2016). The Zak phases of the two bands that form the tilted Dirac cone are quantized as shown in Fig. 4 (a). Each band has a different quantized value that are exchanged at the crossing point. The discrete jump of the Zak phase makes the Berry phase along the circular path around the Dirac point be either $\pi$ or $-\pi$. The non-trivial Berry phase indicates pseudo-spin texture winding around the Dirac point as already conceived in Eq. 4. Moreover, once we translate half of the unitcell, the Zak phase of each band jumps from 0 to $\pi$ or vice versa compared to that from the original unitcell, because the Wannier center is translated by the same amount. Therefore, the domain boundaries between half-translated unitcells are the topological domain boundary where two bands with different Zak phases encounter and the topological states associated with the difference should occur. In Figs. 4 (b) and (c), two different topological domain boundaries between half-translated unit cell structures can be constructed. In each domain boundary, topologically protected edge states that connect two opposite Dirac points emerge (Figs. 4(d) and (e)). The protected states on the domain boundaries are expected to show characteristic scanning tunneling microscopy images and spectroscopic signals as demonstrated recently for other 2D topological materials Kim et al. (2020).

In conclusion, we have shown that the new nonbenzenoid 2D carbon crystal has interesting low energy states such as the zone-center vHS and near-critical tilted Dirac cones. These two states undergo topological transitions by merging themselves under strains. We also expect that several anomalies related with the zone-center saddle point such as fluctuating magnetic states Hlubina et al. (1997); Furukawa et al. (1998); Honerkamp and Salmhofer (2001); Irkhin et al. (2001) as well as nontrivial magnetic responses Giura and Wanderlingh (1968); Henrich (1971); Nikolaev (2021) could be investigated in the new carbon allotrope.

Y.-W.S. was supported by NRF of Korea (Grant No. 2017R1A5A1014862, SRC program: vdWMRC center) and KIAS individual Grant No. (CG031509). H.J. was supported by NRF of Korea (Grant Nos. 2021M3H4A1A03054864, 2019R1A2C1010498, and 2017M3D1A1040833). Computations were supported by the CAC of KIAS.