Effective Hamiltonian of Topological Nodal Line Semimetal in Single-Component Molecular Conductor [Pd(dddt)$_2$] from First-Principles

Using first-principles density-functional theory calculations, we obtain the non-coplanar nodal loop for a single-component molecular conductor [Pd(dddt)$_2$] consisting of HOMO and LUMO with different parity. Focusing on two typical Dirac points, we present a model of an effective 2 $\times$ 2 matrix Hamiltonian in terms of two kinds of velocities associated with the nodal line. The base of the model is taken as HOMO and LUMO on each Dirac point, where two band energies degenerate and the off diagonal matrix element vanishes. The present model, which reasonably describes the Dirac cone in accordance with the first-principles calculation, provides a new method of analyzing electronic states of a topological nodal line semimetal.

Single-component molecular conductor, [Pd(dddt) 2 ], which is insulating at ambient pressure but becomes conducting and shows almost temperature independent resistivity under a pressure of 12.6 GPa, 1) has been studied using first-principles density-functional theory (DFT) calculations, 2) and linearly dispersed bands with Dirac cones were found at 8 GPa. 1,3) An extended Hückel calculation and a tight-binding analysis of the DFT optimized structure at 8 GPa have shown nodal line semimetal. 4) Relevant physical quantities have been examined to understand the characteristics of such new massless Dirac electron systems. [5][6][7][8] In fact, the multi-orbital nature due to an interplay between HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) on different molecular sites is essential for the emergence of the Dirac electron. 1) Such a state, where the contact points of the Dirac cones form a closed line (loop) called nodal line semimetal, [9][10][11][12][13] is quite different from that of massless Dirac electrons in a two-dimensional molecular conductor. 14) In 1937, Herring 15) claimed accidental degeneracies in the energy bands followed by a closed/open circuit (nodal line) for a system with inversion symmetry. The product of parity eigenvalues at the eight time reversal invariant momentum (TRIM) determines even or odd number of circuits (nodal loop). This formula is quite similar to Fu-Kane's criteria of Z 2 topological insulator, 16) despite assuming the absence of spinorbit coupling (SOC). However, since it should work well for weak SOC materials with light elements such as molecular conductors, the shape of the nodal line and the number of nodal loops would be useful to evaluate the band topology and design new molecular topological materials. The nodal line semimetal is a recent topic related to the new concept associated with the topological property of Dirac electrons. [17][18][19][20] Among them, single nodal loop has been found in pnictides CaAgX (X = P and As), 21) alkaline-earth compounds of AX 2 (A = Ca, Sr, Ba; X = Si, Ge, Sn), 22) CaAs 3 , 23) and Ag 2 S. 24) In * E-mail: tsumu@kumamoto-u.ac.jp † E-mail: reizo@riken.jp ‡ E-mail: suzumura@s.phys.nagoya-u.ac.jp fact, all these systems and [Pd(dddt) 2 ] turn out to be a strong topological insulator 25) when SOC is on. 3,4,26) With the regard to the Dirac electrons in [Pd(dddt) 2 ], it is useful to describe the nodal line semimetal in terms of the effective Hamiltonian with a two-band model. For many nodal line semimetals, first-principles calculations have characterized the shape of the nodal line, and an effective Hamiltonian is analyzed only close to the Γ point due to perturbation. 22,26) However, such a Hamiltonian is insufficient to describe the large nodal line away from TRIM due to the non-coplanar 4) and snake-like loop structures. 22) Further, the model is required to show a reasonable energy band for the Dirac cone along the nodal line in order to examine the physical property, since the chemical potential (i.e., the Fermi energy) crosses the nodal line. Indeed, the model, which provides two kinds of velocity fields associated with the Dirac cone, is useful as shown by the calculation of the topological property of the Berry phase. 8) In single-band molecular conductors, first-principles band structures near the Fermi level are well reproduced by a small number of hopping parameters, 27,28) but the corresponding parameters for a multi-band system of [Pd(dddt) 2 ] are huge and complicated. Therefore, a new method for the direct derivation of the model from the DFT calculation is needed to describe the Dirac cone reasonably for all the Dirac points along the loop.
In this letter, on the basis of first-principles DFT calculations, we calculate a number of Dirac points that form a single non-coplanar nodal loop centered at the Γ-point in the first Brillouin zone (BZ). The characteristics of the two crossing bands near the Dirac points and the velocities of the Dirac cone are examined by focusing on two typical Dirac points on the loop. Then, we demonstrate a new method to obtain the effective Hamiltonian for the nodal line semimetal where the base is not on the TRIM, but on the respective Dirac point. 29,30) A model with analytical matrix elements is derived using the knowledge of the general relation between the nodal line and the velocity fields of the Dirac cone. 8) The crystal structure of [Pd(dddt) 2 ], which is determined by x-ray diffraction at ambient pressure, belongs to a mon-oclinic system with the space group of P2 1 /a. 1) In our previous study, using first-principles calculations, we discovered anisotropic Dirac cones of [Pd(dddt) 2 ] at 8GPa by performing structural optimization under symmetry operations. 1) However, the three-dimensional characteristics of the electronic state and the effective model have not yet been reported from first-principles. The present first-principles DFT calculations are performed based on an all-electron fullpotential linearized plane wave (FLAPW) method as implemented in QMD-FLAPW12. [31][32][33] To determine accurately the anisotropy of the Dirac cone, Fermi velocities, and contact points of Dirac cones along the nodal line, we calculate the eigenvalues around the Dirac cone by using the highly precise FLAPW method. We used a high-dense k-mesh for plotting the 3D band structures. The exchange-correlation functional that is used in the present calculation is the generalized gradient approximation proposed by Perdew, Burke, and Ernzerhof. 34) The k-point sampling grid of 4 × 8 × 4 was used. The BZ integration was performed with an improved tetrahedron method. 35 The shape of the Dirac cone depends on the position of the Dirac point along the nodal line. In the following, we use a wavevector scaled by 2π, lattice constants of the reciprocal lattice taken as unity, and energy in the unit of eV. Figures 2(a) and 2(b) show the Dirac cone with Dirac point (I), k 0 = (0, 0.086, 0), where the axis of the cone is parallel to a vector (-2 −1/2 , 0, 2 −1/2 ). The principal axes of the plane, which are perpendicular to the axis of the Dirac cone, are given by a*+c* (k x+z ) and b* (k y ). Dirac cone (I) shown in Fig. 2(a) is symmetric with respect to k x+z while that shown in Fig. 2 where the axis of the cone is parallel to b* (k y ). The former shows a bird's eye view with tilting along both k x and k z axes. The latter denotes the gap function around Dirac point (II), which is the difference in energy between the valence and conduction bands, E c (k) − E v (k). The principal axes of the Dirac cone (II) is rotated counterclockwise with an angle θ ∼ tan −1 (0.2) and the tilting does not appear. The difference in behavior of the Dirac cone between Dirac points (I) and (II) is discussed in terms of the effective Hamiltonian.   estimation of the effective Hamiltonian, we focus on Dirac points (I) and (II) represented by the numbers 1 and 11, respectively. The trajectory of the Dirac points suggests an accidental degeneracy, which occurs without the symmetry of the crystal. Note that these Dirac points are symmetric with respect to the Γ-point and the plane of k y = 0 (k x -k z is a mirror plane), respectively. Using Fig. 3(a), we examine the effective Hamiltonian of a two-band model, where the base is given by the respective Dirac point k 0 . Note that this method describes the Dirac cone correctly, as shown in Refs. 29 and 30. Two energies consisting of different kinds of parities degenerate on the Dirac point, where arbitrary linear combination of these two wavefunctions becomes possible. The eigenfunctions with k, which is just away from the Dirac point is determined uniquely, and their linear combination rotates around the Dirac point. 36) In the present case, we can choose two eigenfunctions of HOMO and LUMO with different symmetry at the Dirac point. Now we derive the effective Hamiltonian of the two band model explicitly. We start with the full 8 × 8 Hamiltonian that is represented at the Γ-point. Such a Hamiltonina is given by where the 4 × 4 matrix H HH (k) ( H LL (k)) consists of only HOMO (LUMO) and H HL (k) (H HL (k)) denotes the interaction between HOMO and LUMO. The former is the even function of k due to the time reversal symmetry, while the latter is the odd function due to interaction with different parity. Among the energy bands of the four HOMOs and four LU-MOs, we choose two bands given by the top HOMO bands and the bottom LUMO bands, which are assumed to be isolated from the other bands. The Scrödinger equation for these two bands E H (k) and E L (k) is given by where H 0 (k) denotes Eq. (1) without H HL (k) and H LH (k). By using |H(k) > and |L(k) >, the reduced Hamiltonian of the 2 × 2 matrix is written as By defining where k = (k x , k y , k z ) and f 0 (0) is taken as zero. The energy is given by The Dirac point k 0 on the nodal line is obtained from Note that f 0 (k) and f 3 (k) are the even functions of k due to time reversal symmetry. The quantity f 2 (k), which is also real due to the presence of inversion symmetry around the Pd atom, 4) is the odd function of k due to HOMO and LUMO having different parities. The function f 2 (k) is estimated by projecting the nodal line on the k xk z plane in Fig. 4(a), while f 3 (k) is estimated by projecting the nodal line on the k xk y plane in Fig. 4(b). This method reproduces the Dirac points well as shown in Fig. 3(a). We empirically examined the form under the constraint that the Dirac points (0, 0.086, 0) (I) and (-0.1967, 0.000, 0.3924) (II) satisfy Eqs. (6a) and (6b).These functions are obtained as where the velocity of the cone at the Dirac point k 0 is obtained as v 2 = ∇ k 0 f 2 , v 3 = ∇ k 0 f 3 and v 0 = ∇ k 0 f 0 corresponding to Letter Author Name  7a) and (7b)), respectively. The inset shows f 2 (k) around the Γ-point. The arrow, which depicts v 2 at Dirac point (I) with arbitrary scale, corresponds to the horizontal axis a*+c* in Fig. 2(a).
the tilting. Using k 0 = (k 0x , k 0y , k 0z ), v j is calculated as Note that the tangent of the nodal line is parallel to v 2 × v 3 . 8) In the previous work, 30) v j was estimated directly from Eq. (4) using the wavefunction close to k = k 0 . From the energy of the Dirac cone is given by The coefficients of k 3 x , k 5 x and (k x k y ) 2 in Eqs. (7a) and (7b) are determined by fitting the line to the other k 0 . As shown in Figs. 4(a) and 4(b), good coincidence is obtained for both f 2 (k) and f 3 (k) where f 2 (k) reflects the symmetry of k y =0. The arrow in the inset of Fig. 4(a) represents v 2 at Dirac point (I) which is perpendicular to f 2 (k) (solid line). Coefficients C 2 , C 3 , b x , b y , and b z are calculated as follows. First, we examine the Dirac cone at Dirac point (I) where v 2 = C 2 (1, 0, 1), v 3 = C 3 (0, 1, 0), and v 0 = 2b y k 0y (0, 1, 0). Noting v 2 · v 3 = 0, the principal axes are given by k x +k z and k y . By comparing the velocities with those of Figs. 2(a) and 2(b), which give we obtain C 2 ≃ 0.148, b y ≃ −2.6, and C 3 ≃ 0.053. Coefficients C 2 and C 3 can also be obtained from the gap E c (k) − E v (k) at the Γ-point and Z-point, which are estimated as 0.0994 and 0.1117, respectively. The resultant quantities C 2 = 0.137 and C 3 = 0.050 are compatible with those of Eqs. (9a) and (9b), which suggests that the effective Hamiltonian of Eq. (5) is valid not only for k close to the nodal line but also for the extended regions including the Γand Zpoints. Next, we examine the Dirac cone at Dirac point (II) to verify the validity of Eqs. (7a) and (7a) with C 2 = 0.148 and C 3 = 0.053. The Dirac cone is shown on the plane of k x k z , but the principal axes do not coincide with that of v 2 and v 3 since v 2 · v 3 0 from Eqs. (8a) and (8a). We estimate the gap function given by 2∆(k) = E c (k) − E v (k) = 2 (v 2 · δk) 2 + (v 3 · δk) 2 with δk = (x, 0, z), which is calculated as (E c (k) − E v (k)) 2 /4 = 0.45x 2 + 0.117xz + 0.02z 2 (= Ax 2 + 2Cxz + Bz 2 ). This shows the rotation of the principal axes from (x, z) to (X, Z), where x = X cos θ − Z sin θ, and z = X sin θ + Z cos θ. In terms of (X, Z), the gap function is expressed as The cross section of ∆(k) = E 0 shows an ellipse with the radius of the minor [major] axis given . The radius and angle for Dirac point (II) are obtained as a = 1.5E 0 , b = 8.5E 0 , b/a ≃ 5.6 and tan θ = 0.13, which correspond well to the numerical result in Fig. 2(d), i.e., b/a ≃ 7 and tan θ ≃ 0.2. Therefore, our Hamiltonian may be applied to all the Dirac points between points (I) and (II). Further, from Eq. (10), we estimate the area of the ellipse S with the gap 2E 0 for an arbitrary Dirac point, where the axis of the cone is parallel to the tangent of the nodal line ( i.e., ∝ v 2 × v 3 ) and the wavefunction of the Dirac cone is determined by v 2 and v 3 . By taking the axis of the cone as y and the plane of the cone as the xz plane, we obtain the area of the ellipse S with the gap 2E 0 as S (k 0 ) = πab 0 is generally expected except for Dirac point (I), which is located on the symmetric line of Γ -Y.
Finally, we calculate b x and b z in Eq. (7c) from Dirac cone (II), which gives v x = 0.36 ± 0.12 and v z = 0.09 ± 0.06. In a way similar to Eqs. (9a) and (9b), we obtain b x ≃ −0.30 and b z ≃ 0.077 suggesting that the tilting occurs towards the negative (negative) direction for the k x (k z ) axis. This is qualitatively consistent with Fig. 2(c). As shown in Fig. 3(b), our DFT calculation of the energy difference measured from that of Dirac point (I) exhibits nonmonotonical behavior, e.g., -0.0031, -0.0061, -0.0059, -0.0047, and -0.0032 eV for the nodal points No. 3, 5, 7, 9, and 11 in Fig. 3(a), respectively. However, it is quite complicated to reproduce the dispersion quantitatively from Eq. (7c) with a constant C 0 , and such a problem becomes an issue to be solved in the next step.
Here, it should be noted that the SOC was ignored in the present calculation. Since the crystal structure is centrosymmetric, the band structures are in Kramers degeneracy. Then, the loop degeneracies are destroyed by the SOC, leading to an energy gap with ≃ 3 meV at Dirac point (I), 3) which will be examined further.
We comment on the experiment in which the present model is useful for performing the analysis. Since all the directions for the nodal line can be treated on the same footing, the effect of magnetic field is of interest to examine the characteristic behavior of the Landau level. The Hall coefficient is also expected to exhibit anisotropic behavior.
In summary, we performed first-principles calculations for [Pd(dddt) 2 ] at 8GPa, and found a non-coplanar nodal loop within the first BZ. The present effective Hamiltonian obtained from the Dirac points in Fig. 3(a) describes well the electron close to the nodal line. The explicit form of the velocity is obtained by Eqs. (8a)-(8c), with C 2 ≃ 0.148, C 3 ≃ 0.05, b x ≃ −0.31, b y ≃ −2.6, and b z ≃ 0.077, which are useful for future calculations of transport.