Novel Dirac Electron in Single-Component Molecular Conductor [Pd(dddt)$_2$]

TL;DR

This paper investigates Dirac electrons in a single-component molecular conductor [Pd(dddt)$_2$] under pressure, revealing a Dirac loop in the 3D Brillouin zone explained by symmetry and parity considerations.

Contribution

It introduces a tight-binding model for [Pd(dddt)$_2$] that explains the Dirac cone formation and the Dirac loop using symmetry and parity analysis.

Findings

Dirac electrons form a loop in the 3D Brillouin zone.

The Dirac cone arises from HOMO-LUMO symmetry properties.

Parity at TRIM points explains the Dirac electron behavior.

Abstract

Dirac electrons in a single-component molecular conductor [Pd(dddt)$_2$] under pressure have been examined using a tight-binding model which consists of HOMO (Highest Occupied Molecular Orbital) and LUMO (Lowest Unoccupied Molecular Orbital) functions in four molecules per unit cell. The Dirac cone between the conduction and valence bands comes from the property that HOMO has ungerade symmetry and LUMO has gerade symmetry. The Dirac point forms a loop in the three-dimensional Brillouin zone, which is symmetric with respect to the plane of $k_y$ = 0 where $k_y$ is the intralayer momentum along the molecular stacking direction, i.e., with the largest (HOMO-HOMO, LUMO-LUMO) transfer energy. The parity at TRIM (time reversal invariant momentum) is calculated using the inversion symmetry around the lattice point of the crystal. It is shown that such an exotic Dirac electron is understood from the parity of the wave function at the TRIM and also from an effective Hamiltonian.

Tables

Table 1. Table 1: Parities P E ​ j ​ ( 𝑮 / 2 ) ( = ± ) annotated subscript 𝑃 𝐸 𝑗 𝑮 2 absent plus-or-minus P_{Ej}(\bm{G}/2)(=\pm) at P 𝑃 P = 8 GPa, where the TRIMs are 𝑮 / 2 = ( 0 , 0 , 0 ) 𝑮 2 0 0 0 \bm{G}/2=(0,0,0) ( Γ Γ \Gamma ), ( π , 0 , 0 ) 𝜋 0 0 (\pi,0,0) (X), ( 0 , π , 0 ) 0 𝜋 0 (0,\pi,0) (Y), ( π , π , 0 ) 𝜋 𝜋 0 (\pi,\pi,0) (M), ( 0 , 0 , π ) 0 0 𝜋 (0,0,\pi) (Z), ( π , 0 , π ) 𝜋 0 𝜋 (\pi,0,\pi) (D), ( 0 , π , π ) 0 𝜋 𝜋 (0,\pi,\pi) (C), and ( π , π , π ) 𝜋 𝜋 𝜋 (\pi,\pi,\pi) (E).

	$E_1$	$E_2$	$E_3$	$E_4$	$E_5$	$E_6$	$E_7$	$E_8$
$P_{Ej}\left(\Gamma \right)$	$+$	$+$	$+$	$-$	$+$	$-$	$-$	$-$
$P_{Ej}\left(X\right)$	$+$	$-$	$-$	$+$	$-$	$+$	$+$	$-$
$P_{Ej}\left(Y\right)$	$-$	$+$	$+$	$-$	$+$	$-$	$-$	$+$
$P_{Ej}\left(M\right)$	$+$	$+$	$+$	$+$	$-$	$-$	$-$	$-$
$P_{Ej}\left(Z\right)$	$+$	$+$	$-$	$-$	$-$	$+$	$-$	$+$
$P_{Ej}\left(D\right)$	$+$	$-$	$+$	$-$	$-$	$+$	$-$	$+$
$P_{Ej}\left(C\right)$	$+$	$-$	$+$	$-$	$-$	$+$	$-$	$+$
$P_{Ej}\left(E\right)$	$-$	$-$	$+$	$+$	$+$	$+$	$-$	$-$

Table 2. Table 2: The parity of P E ​ j ​ ( Γ ) subscript 𝑃 𝐸 𝑗 Γ P_{Ej}(\Gamma) and P E ​ j ​ ( Z ) subscript 𝑃 𝐸 𝑗 Z P_{Ej}({\rm Z}) for P 𝑃 P =0 where those of P E ​ j ​ ( 𝑮 / 2 ) subscript 𝑃 𝐸 𝑗 𝑮 2 P_{Ej}(\bm{G}/2) with X ,Y, M ,D, C and E are the same as Table 1 .

$P=0$	$E_1$	$E_2$	$E_3$	$E_4$	$E_5$	$E_6$	$E_7$	$E_8$
$P_{Ej}\left(\Gamma \right)$	$+$	$+$	$+$	$+$	$-$	$-$	$-$	$-$
$P_{Ej}\left(Z\right)$	$+$	$+$	$-$	$-$	$-$	$-$	$+$	$+$

Equations

$$H=\sum _{i,j}{t}_{i,j;\alpha ,\beta }|i,\alpha ><j,\beta |,$$

$$H=\sum _{\mathbf{k}}|\Phi (\mathbf{k})>\hat{H}(\mathbf{k})<\Phi (\mathbf{k})|,$$

$$t_{\alpha ,\beta }$$

$$\hat{H}(\mathbf{k}){\Psi }_j(\mathbf{k})=E_j(\mathbf{k}){\Psi }_j(\mathbf{k}),$$

$${\Psi }_j(\mathbf{k})=\sum _{\alpha }{d}_{j,\alpha }(\mathbf{k})|\alpha >,$$

$$E_g(\mathbf{k})=\min (E_4(\mathbf{k})-E_5(\mathbf{k})),$$

$$\hat{P}(\mathbf{k})$$

$$\hat{P}(\mathbf{k})u_{\alpha }=p_{\alpha }(\mathbf{k})u_{\alpha },$$

$$\hat{P}(\mathbf{G}/2){\Psi }_j(\mathbf{G}/2)=P_{Ej}(\mathbf{G}/2){\Psi }_j(\mathbf{G}/2)$$

$$P_{Ej}(\mathbf{G}/2)={\Psi }_j(\mathbf{G}/2{)}^{†}\hat{P}(\mathbf{G}/2){\Psi }_j(\mathbf{G}/2)=\sum _{\alpha =H1}^{L4}{p}_{\alpha }(\mathbf{G}/2)|d_{j\alpha }{|}^2.$$

$$P(k_z=0)$$

$$P(k_y=0)$$

$$P(k_z=\pi )$$

$$P_{\delta }=-1,+1.$$

$$D(\omega )$$

$$P(k_z=0)$$

$$P(k_y=0)$$

$$P(k_z=\pi )$$

$$(|d_{j,H_1}{|}^2,|d_{j,H_2}{|}^2,|d_{j,H_3}{|}^2,|d_{j,H_4}{|}^2,|d_{j,L,1}{|}^2,|d_{j,L,2}{|}^2,|d_{j,L,3}{|}^2,|d_{j,L,4}{|}^2)$$

$$\simeq$$

$$=$$

$$H_{\mathrm{eff}}(\mathbf{k})$$

$$t_{H1,H1}$$

$$t_{H1,H2}$$

$$t_{H1,H3}$$

$$t_{H1,H4}$$

$$t_{H1,L1}$$

$$t_{H1,L2}$$

$$t_{H1,L3}$$

$$t_{H1,L4}$$

$$t_{H2,H2}$$

$$t_{H2,H3}$$

$$t_{H2,H4}$$

$$t_{H2,L1}$$

$$t_{H2,L2}$$

$$t_{H2,L3}$$

$$t_{H2,L4}$$

$$t_{H3,H3}$$

$$t_{H3,H4}$$

$$t_{H3,L1}$$

$$t_{H3,L2}$$

Full text

1 Introduction

Since the discovery of the quantum Hall effect in graphene,[1] two-dimensional (2D) massless Dirac fermions have been a fascinating topic. In addition to the graphene with monolayer, a Dirac electron was found in organic conductor $\alpha$-(BEDT-TTF)2I3  as a bulk system, [2] and the properties of molecular Dirac fermion systems have been studied extensively. [3]

Recently, a Dirac electron was found in the single-component molecular conductor [Pd(dddt)2] (dddt=5,6-dihydro-1,4-dithiin-2,3-dithiolate), which shows a constant resistivity with decreasing temperature under pressure. [4, 5] Based on first-principles calculation, which shows the existence of a Dirac cone, [6] a tight-binding model of [Pd(dddt)2] consisting of highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) functions in four molecules per unit cell was proposed.[7] In the crystal, there are two crystallographically independent layers given by layers 1 and 2 (Fig. 1), and the Dirac cone originates from the HOMO-based band in layer 1 and the LUMO-based band in layer 2. [4] The interplay of the interlayer and intralayer transfer integrals is crucial to obtain the Dirac point.[4] The existence of such a Dirac point is clarified as follows for a simple case of 2D momentum with a fixed interlayer momentum. Since the difference in the energy level between the HOMO and LUMO is small, the HOMO band is located higher than the LUMO band around the $\Gamma$ point, and then the overlap between them results in a Fermi line in the absence of the HOMO-LUMO transfer energies. However the Fermi line disappears in the presence of the HOMO-LUMO transfer energies owing to the opening of a gap, which comes from the combined effect of the intralayer HOMO-LUMO transfer energies and the interlayer HOMO-HOMO/LUMO-LUMO transfer energies. Furthermore, there is a line (nodal line) passing through the $\Gamma$ point, on which the HOMO-LUMO coupling vanishes. Thus, the Dirac point is obtained at the intersection of the Fermi line and the nodal line owing to the closing of the gap. [4] Such a nodal line originates from a property of the HOMO-LUMO transfer energy. The [Pd(dddt)2] molecule has an inversion center at the Pd atom. For the HOMO with a bonding $\pi$ character, the inversion produces a phase change for the molecular orbital. This means that the HOMO has ungerade (odd) symmetry. On the other hand, the LUMO with an antibonding $\pi$ character has gerade (even) symmetry, where the inversion results in the same phase for the orbital. Due to these properties, the HOMO-LUMO transfer energies interrelated by the center of symmetry have opposite signs in the [Pd(dddt)2] crystal.

In addition to the Dirac cone on the $k_{x}$-$k_{y}$ plane of the intralayer 2D momentum, a novel feature of [Pd(dddt)2] is the formation of a loop of the Dirac point with varying the interlayer momentum $k_{z}$ in the three-dimensional (3D) Brillouin zone. [8] For the case of $\alpha$-(BEDT-TTF)2I3 , the negligible interlayer HOMO-HOMO transfer integral suggests a Dirac line that is open at the boundary of the first Brillouin zone.[9] The Dirac electron in [Pd(dddt)2] is exotic since the Dirac point originates from the interplay of the HOMO and LUMO functions. The purpose of the present paper is to clarify such Dirac electron by examining the Dirac point in the 3D Brillouin zone. The result is analyzed by calculating the parity of the wave function at the TRIM (time reversal invariant momentum), which comes from the inversion symmetry around the molecular site. [10] In Sect. 2, the model and formulation are given. In Sect. 3, Dirac points under a pressure of 8 GPa corresponding to a previous experiment are calculated as a function of $\bm{k}=(k_{x},k_{y},k_{z})$, where $k_{y}$ is the intralayer momentum along the molecular stacking direction, i.e., with the largest transfer energy. It is shown that a pair of Dirac points forms a loop in the 3D Brillouin zone. In Sect. 4, the mechanism for the formation of the loop of the Dirac point is analyzed in terms of the parity at the TRIM, which is calculated for both the pressure of 8 GPa and the ambient pressure to comprehend the emergence of the Dirac point. A summary and discussion in terms of the effective Hamiltonian are given in Sect. 5.

2 Model and Formulation

2.1 Tight-binding model

The crystal structure of [Pd(dddt)2] is shown in Fig. 1, which consists of four molecules (1, 2, 3, and 4) with HOMO and LUMO functions in the unit cell. The crystal structure was determined by the single-crystal X-ray diffraction method at ambient pressure ($P$=0) and was estimated by first-principles density functional theory (DFT) calculations at $P$= 8 GPa.[4] Transfer energies are given by pairs between nearest-neighbor molecules, where those between the layers are expressed as $a$ (molecules 1 and 2 and molecules 3 and 4), and $c$ (molecules 1 and 4 and molecules 2 and 3), those in the same layer are expressed as $p$ (molecules 1 and 3) and $q$ (molecules 2 and 4), and those along the stacking axis are given by $b$ (Fig. 1). Based on the crystal structure, we examine the tight-binding model Hamiltonian given by

[TABLE]

where $i$ and $j$ are the sites of the unit cell with total number $N$, and $\alpha$ and $\beta$ denote the eight molecular orbitals given by the HOMO $(H1,H2,H3,H4)$ and LUMO $(L1,L2,L3,L4)$. The lattice constant is taken as unity. The transfer energies, $t_{i,j;\alpha,\beta}$, are classified as HOMO-HOMO (HH), LUMO-LUMO (LL), and HOMO-LUMO (HL) transfer energies. Taking eV as the unit of energy, the transfer energies $t_{i,j;\alpha,\beta}$ under a pressure of $P$ = 8 GPa[4] (0 GPa) are given by $a_{H}=-0.0345\;(-0.0136)$, $a_{L}=-0.0\;(-0.0049)$, $a_{HL}=0.0260\;(0.0104)$, $b_{1H}=0.2040\;(0.112)$, $b_{1L}=0.0648\;(0.0198)$, $b_{1HL}=0.0219\;(0.0214)$, $b_{2H}=0.0762\;(0.0647)$, $b_{2L}=-0.0413\;(0.0)$, $b_{2HL}=-0.0531\;(-0.0219)$, $c_{H}=0.0118\;(0.0)$, $c_{L}=-0.0167\;(-0.0031)$, $c_{HL}=0.0218\;(0.0040)$, $p_{H}=0.0398\;(0.0102)$, $p_{L}=0.0205\;(0.0049)$, $p_{1HL}=-0.0275\;(-0.0067)$, $p_{2HL}=-0.0293\;(-0.0074)$, $q_{H}=0.0247\;(0.0067)$, $q_{L}=0.0148\;(0.0037)$, $q_{1HL}=-0.0186\;(-0.0048)$, and $q_{2HL}=-0.0191\;(-0.0051)$. The transfer energies were calculated using the extended Hückel method. The gap between the energy of the HOMO and that of the LUMO is taken as $\Delta E=$ 0.696 eV to reproduce the energy band of the first-principles calculation.

Using the Fourier transform $|\alpha(\bm{k})>=\sum_{j}\exp[-i\bm{k}\bm{r}_{j}]\;|j,\alpha>$ with wave vector $\bm{k}=(k_{x},k_{y},k_{z})$, Eq. (1) is calculated as[4]

[TABLE]

where $\hat{H}(\bm{k})$ is the Hermite matrix Hamiltonian with the matrix elements $t_{\alpha,\beta}$ defined by

[TABLE]

and the base is given by $<\Phi(\bm{k})|=(<H1|,<H2|,<H3|,<H4|,<L1|,<L2|,<L3|,<L4|).$ The matrix elements of $t_{\alpha,\beta}$ are given in the Appendix. Since the symmetry of the HOMO (LUMO) is odd (even) with respect to the Pd atom, the matrix element of H-L (H-H and L-L) is the odd (even) function with respect to $\bm{k}$. The energy band $E_{j}(\bm{k})$ and the wave function $\Psi_{j}(\bm{k})$, $(j=1,2,\cdots,8)$ are calculated from

[TABLE]

where $E_{1}>E_{2}>\cdots>E_{8}$ and

[TABLE]

with $\alpha=$ H1, H2, H3, H4, L1, L2, L3, and L4. Noting that the band is half-filled owing to the HOMO and LUMO functions, we examine the gap defined by

[TABLE]

for all $\bm{k}$ in the Brillouin zone. The Dirac point $\bm{k}_{D}$ is obtained from $E_{g}(\bm{k}_{\rm D})=0$.

2.2 Parity at TRIM

In order to analyze the Dirac point, we calculate the parity at the TRIM given by $\bm{G}$/2 with $\bm{G}$ being the reciprocal lattice vector, where $\bm{G}/2=(0,0,0)$, $(\pi,0,0)$, $(0,\pi,0)$, and $(\pi,\pi,0)$ correspond to the $\Gamma$, X, Y, and M points, and $\bm{G}/2=(0,0,\pi)$, $(\pi,0,\pi)$, $(0,\pi,\pi)$, and $(\pi,\pi,\pi)$ correspond to the Z, D, C, and E points, respectively. Applying the case of $\alpha$-(BEDT-TTF)2I3  with the 4 $\times$ 4 matrix Hamiltonian[11] to the present case of the 8 $\times$ 8 matrix Hamiltonian, the inversion with respect to molecular site 1 in Fig. 1 gives the matrix for the translation of the base, $\hat{P}(\bm{G}/2)$, expressed as

[TABLE]

The relation $(\hat{P}(\bm{k}))_{{\rm Hj},{\rm Hj}}=-(\hat{P}(\bm{k}))_{{\rm Lj},{\rm Lj}}$ for $j$=1, 2, 3,and 4 originates from the fact that the HOMO has ungerade symmetry and the LUMO has gerade symmetry. The eigenvalue and eigenfunction ($\alpha$ = H1, H2, $\cdots$, L4) are obtained from

[TABLE]

where $p_{\alpha}=(\hat{P}(\bm{k}))_{\alpha,\alpha}$, $u_{\rm H1}(\bm{k})^{t}=(1,0,0,0,0,0,0,0)$, $u_{\rm H2}^{t}=(0,1,0,0,0,0,0,0)$, $\cdots$, and $u_{\rm L4}^{t}=(0,0,0,0,0,0,0,1)$. At the TRIM, one obtains $p_{\alpha}(\bm{G}/2)=+(-)$, which gives the parity. For example, the odd parity $p_{\alpha}(\bm{G}/2)=-$ at the TRIM is given by $\alpha$ = H1, H2, H3, H4 for the $\Gamma$ point, $\alpha$ = H1, H4, L2, L3 for the X point, $\alpha$ = H1, H4, L2, L3 for the Y point, $\alpha$ = H1, H2, H3, H4 for the M point, $\alpha$ = H1, H3, L2, L4 for the Z point, $\alpha$ = H1, H2, L3, L4 for the D point, $\alpha$ = H1, H2, L3, L4 for the C point, and $\alpha$ = H1, H3, L2, L4 for the E point. Note that the wave function at the $\Gamma$ point with odd parity is given only by the HOMO and that at the Z point is given by both the LUMO and HOMO.

Since $[\hat{P}(\bm{G}/2),\hat{H}(\bm{G}/2)]=0$, $\Psi_{j}(\bm{G}/2)$ is also an eigenfunction of $\hat{P}(\bm{G}/2)$. Then, at the TRIM, we obtain the parity from

[TABLE]

with $P_{Ej}(\bm{G}/2)=+(-)$, which denotes the even (odd) parity. Note that $d_{j\alpha}$ in Eq. (5) vanishes for $\alpha$ when the parity $p_{\alpha}(\bm{G}/2)$ is opposite to that of $P_{Ej}(\bm{G}/2)$. Using Eq. (5), Eq. (9) is rewritten as

[TABLE]

In Eq. (9), we used the notation $P_{Ej}(\Gamma)$, $P_{Ej}(X)$, $P_{Ej}(Y)$, $P_{Ej}(M)$, $P_{Ej}(Z)$, $P_{E_{j}}(D)$, $P_{Ej}(C)$, and $P_{Ej}(E)$, for $\bm{k}$ corresponding to the $\Gamma$, X, Y, M, Z, D, C, and E points, respectively. We define $P_{\delta}$ [ = $P(k_{z}=0)$, $P(k_{y}=0)$, and $P(k_{z}=\pi)$] as

[TABLE]

where each $P_{\delta}$ denotes a quantity assigned on a plane including the four respective TRIMs. The condition for the Dirac point between $E_{4}$ and $E_{5}$ is given by[12, 11]

[TABLE]

When $P_{\delta}=-1(+1)$, the number of pairs of Dirac points between $E_{4}$ and $E_{5}$ is odd (zero or even). This fact can be understood from the idea of the $\pi$ jump for the $Z_{2}$ Berry phase. [13] Note that Eqs.  (11a), (11b), and (11c) describe the condition of the Dirac point on the planes of $k_{c}=0$ (TRIM with the $\Gamma$, X, Y, and M points) and $k_{b}=0$ (TRIM with the Z, $\Gamma$, X, and D points), and $k_{c}=\pi$ (TRIM with the Z, D, C, and E points), respectively.

3 Loop of Dirac Point in Three Dimensions

First we examine the energy bands $E_{4}(\bm{k})$ and $E_{5}(\bm{k})$ with the Dirac point for two typical cases of $k_{z}$ = 0 and $k_{y}$ = 0, which are shown on the $k_{x}$-$k_{y}$ plane in Fig. 2(a) and on the $k_{x}$-$k_{z}$ plane in Fig. 2(b), respectively. A Dirac point with the tilted Dirac cone exists between $E_{4}(\bm{k})$ and $E_{5}(\bm{k})$, where the upper (lower) band denotes $E_{4}(\bm{k})$ ($E_{5}(\bm{k})$) corresponding to the conduction (valence) band. In Fig. 2(a), the center denotes the $\Gamma$ point ($\bm{k}$ = (0, 0, 0)) and the Dirac points are located on the line of $k_{x}$ = 0. The convex upward region in $E_{4}(\bm{k})$ is mainly determined by the HOMO function in layer 1, while the convex downward region in $E_{5}(\bm{k})$ is mainly determined by the LUMO function in layer 2. The energy of $E_{4}(\bm{k})$ ($E_{5}(\bm{k})$) in the remaining region originates from the LUMO (HOMO) function. This means that the HOMO-based band in layer 1 and the LUMO-based band in layer 2 cross around the $\Gamma$ point.[4] In Fig. 2(b), the center denotes the Z point ($\bm{k}$ = (0, 0, $\pi$)) and the Dirac points are disposed symmetrically with respect to the Z point. The saddle point of $E_{4}(\bm{k})$ can be seen at the Z point. The tilted Dirac cone is elongated along the $k_{z}$ axis in Fig. 2(a) and along the $k_{y}$ axis in Fig. 2(b). This suggests that the Dirac cone is overturned with increasing $k_{z}$.

Next we examine the Dirac point in the 3D wave vector space. Since significant interlayer (along the $k_{z}$ direction) transfer integrals exist, the system has a 3D character. In particular, the HOMO-LUMO couplings that play a crucial role in the Dirac cone formation exhibit $k_{z}$ dependence.[4] Figure 3 demonstrates that the Dirac point moves in a manner depending on $k_{z}$ and describes a loop in the extended Brillouin zone. The loop is symmetrical with respect to the plane of $k_{y}$ = 0, while the symmetry of the Dirac point with respect to the $\Gamma$ point is due to the time-reversal symmetry of the Hamiltonian. With increasing $|k_{z}|$, a crossover from the Dirac cone on the $k_{x}$-$k_{y}$ plane to that on the $k_{z}$-$k_{x}$ plane occurs, for example, the Dirac cone is already elongated along the $k_{y}$ axis for the Dirac points $\bm{k}_{\rm D}/(2\pi)$ = $\pm(-0.0425,\pm 0.0265,0.25)$. The loop is almost parallel to the $k_{y}$-$k_{z}$ plane around $k_{z}/2\pi$ = 0, while the loop is turned up from the $k_{y}$-$k_{z}$ plane with increasing $|k_{z}/2\pi|$ up to 0.545 (Fig. 3). Regarding the number of Dirac points on the $k_{x}$-$k_{y}$ plane with a fixed $k_{z}$ value in the reduced Brillouin zone, there are two Dirac points for $0<k_{z}/\pi<0.91$ and $1.09<k_{z}/\pi<2.0$, and four Dirac points for $0.91<k_{z}/\pi<1.09$. In the next section, we show another aspect of two pairs of Dirac points at $k_{z}=\pi$.

Here we note a slight variation of the energy on the loop, which gives an electron pocket around $k_{z}$ = 0 (Fig. 2(a)) and a hole pocket around $k_{y}$ = 0 (Fig. 2(b)), implying a nodal line semi-metal.[14, 15] In fact, the energy at the Dirac point in Fig. 2(a) (Fig. 2(b)) is $\simeq 0.002$ V ($\simeq 0.001$ eV) below (above) the Fermi energy[16].

In order to understand the electronic properties associated with such a loop of the Dirac point, we examine the density of states (DOS) per spin and per site, which is given by[16]

[TABLE]

where the Fermi energy, $E_{\rm F}(=0.5561)$, is obtained from $\int_{-\infty}^{E_{\rm F}}d\omega D(\omega)$ = 2 due to a half-filled band. The quantity $[\cdots]$ corresponds to the 2D DOS for the fixed $k_{z}$, which is proportional to $|\omega-E_{\rm D}|$ with $E_{\rm D}(=E_{4}(\bm{k}_{\rm D})=E_{5}(\bm{k}_{\rm D}))$ close to the Dirac point. As shown in the inset of Fig. 4, $E_{\rm D}$ exhibits a monotonic increase with increasing $k_{z}/\pi$ from 0 to 1.09 due to the $k_{z}$ dependence of $\bm{k}_{\rm D}$. The total DOS obtained after the summation of $k_{z}$ is shown in Fig. 4, where $\omega$ is measured from $E_{\rm F}$. This variation of $E_{\rm D}$ results in $D(0)\not=0$, i.e., metallic behavior for small $\omega$. The location of the dip due to the Dirac cone is slightly above $E_{\rm F}$ and the width of the dip is smaller than that of $E_{\rm D}$ since the number of Dirac points above $E_{\rm F}$ is larger than that below $E_{\rm F}$ as seen from Fig. 3. This asymmetry of the DOS may give rise to a characteristic temperature ($T$) dependences of the magnetic susceptibility[16] and specific heat compared with those of the 2D Dirac cone, which are proportional to $T$ and $T^{2}$, respectively.

4 Dirac Point vs. Parity at TRIM

Now we examine the Dirac points using the condition Eq. (12). Since there is a degeneracy between $E_{1}$ and $E_{2}$ ($E_{3}$ and $E_{4}$ ) at certain TRIMs, a small amount of site potential is added to the diagonal elements in $\hat{H}(\bm{k})$ to obtain the parity of Eq. (10). In the present case, we add a potential $-0.0001$ eV to $t_{H3,H3}$, $t_{H4,H4}$, $t_{L3,L3}$, and $t_{L4,L4}$ in Eq. (3). Some of the $P_{Ej}(\bm{G}/2)$ depend on the choice of such potential but Eq. (12) remains the same. From Eq. (10), we obtain 64 parities, where half of them correspond to those of the filled band. Note that $\Pi_{j=1}^{8}P_{Ej}(\bm{G}/2)$ = + since half of the eight eigenvalues of Eq. (LABEL:eq:P) are negative for every $\bm{G}/2$. The respective parities $P_{Ej}(\bm{G}/2)$ at $P$ = 8 GPa are listed in Table 1.

Substituting the parities of Table 1 into Eqs. (11a), (11b), and (11c), we obtain

[TABLE]

Equations (14a) and (14b) are consistent with the fact that there is a pair of Dirac points for $k_{z}=0$ and for $k_{y}=0$ as shown in Figs. 2(a) and 2(b), respectively. The loop of the Dirac point in the Brillouin zone can be understood from Eqs. (14a) and (14b) as follows. Equation (14a) shows a pair of Dirac points located on the planes of $k_{z}=0$, while Eq. (14b) displays another pair of Dirac points on the plane of $k_{y}=0$ being perpendicular to that of $k_{z}=0$. Noting that the lines connecting these Dirac points are symmetric with respect to $k_{y}=0$, the existence of the Dirac point on the plane of $k_{y}=0$ suggests the loop of the Dirac point. Equation (14c) is also consistent with the fact that two pairs of Dirac points exist for $k_{z}=\pi$ as mentioned in the previous section.

Here we numerically examine the difference in the parity between $P_{Ej}(\Gamma)$ and $P_{Ej}({\rm Z})$ with $j$= 4 and 5, which gives Eqs. (14a) and (14c). The components of the wave function $d_{j,\alpha}(\bm{k})$ in Eq. (5) are estimated as

[TABLE]

Substituting Eqs. (15)-(15d) into Eq. (10), we obtain $P_{E4}(\Gamma)=-$, $P_{E5}(\Gamma)=+$, $P_{E4}({\rm Z})=-$, and $P_{E5}({\rm Z})=-$, respectively, where $(p_{\rm H1}(\bm{G}/2),p_{\rm H2}(\bm{G}/2),p_{\rm H3}(\bm{G}/2),p_{\rm H3}(\bm{G}/2),$ $p_{\rm L1}(\bm{G}/2),p_{\rm L2}(\bm{G}/2),p_{\rm L3}(\bm{G}/2),p_{\rm L4}(\bm{G}/2))$ = $(-,-,-,-,+,+,+,+)$ for the $\Gamma$ point and $(-,+,-,+,+,-,+,-)$ for the Z point from Eq. (LABEL:eq:P). Both $E_{4}(\Gamma)$ and $E_{4}({\rm Z})$ are determined by the HOMO function while both $E_{5}(\Gamma)$ and $E_{5}({\rm Z})$ are determined by the LUMO function. The difference in the parity between $E_{4}$ and $E_{5}$ at the $\Gamma$ point is understood from the LUMO and HOMO functions, whose parities are different from each other. However, at the Z point, the parity of $E_{4}$ is the same as that of $E_{5}$, although the former (latter) is described by the HOMO (LUMO) function. Actually, the same parity of $P_{E4}({\rm Z})=P_{E5}({\rm Z})=-$ is obtained from Eqs. (15c) and (15d) with $p_{\rm H1}({\rm Z})=p_{\rm H3}({\rm Z})=p_{\rm L2}({\rm Z})=p_{\rm L4}({\rm Z})=-$. We also see, from Eq. (LABEL:eq:P), that the product of the parities of $\Psi_{4}$(Z), $\Psi_{6}$(Z), $\Psi_{7}$(Z), and $\Psi_{8}$(Z) becomes positive since all of these wave functions are given by the HOMO function. Thus, we obtain $P(k_{z}=\pi)=+1$ owing to $P_{E4}({\rm Z})=P_{E5}({\rm Z})=-$, although $\Psi_{5}$(Z) is given by the LUMO function. Furthermore we note that the loop of the Dirac point is associated with the fact that $\Psi_{5}$(Z) is determined by $L2$ and $L4$, i.e., the LUMO function is determined by layer 2.

Finally, we note the emergence of the Dirac point where a pair of Dirac points appears at one of the TRIMs.[17] In Table 2, the parities at the $\Gamma$ point and the Z point at $P$=0 (ambient pressure) are shown to compare with those at $P$=8 GPa in Table 1. It turns out that the Dirac point is absent at $P$ = 0 owing to $P(k_{z}=0)$= + (i.e., $P_{E4}(\Gamma)$= + and $P_{E5}(\Gamma)=-$ ). Using a linear interpolation for the transfer energy between $P=0$ and 8 GPa, we obtain that the emergence of the loop of the Dirac point occurs at $P$ $\simeq$ 7.6 GPa. For the $\Gamma$ point, the level crossing between $E_{4}(\Gamma)$ and $E_{5}(\Gamma)$ followed by the emergence of a pair of Dirac points occurs owing to the different parity. For the Z point, the level crossing between $E_{4}$(Z) and $E_{5}$(Z) also occurs for $P\simeq$ 7.8 GPa owing to the different functions of the HOMO and LUMO but with the same parity, resulting in a loop within the first Brillouin zone for $7.6\;{\rm GPa}<P<7.8\;{\rm GPa}$.

5 Summary and Discussion

We examined the Dirac point in the single-component molecular conductor [Pd(dddt)2] within a tight-binding model that consists of HOMO and LUMO functions in four molecules in the unit cell. It is crucial for the present Dirac electron that the HOMO has the ungerade symmetry and the LUMO has the gerade symmetry. We obtained a loop of the Dirac point in the 3D Brillouin zone as a result of the combined effect of the interlayer HOMO-HOMO/LUMO-LUMO and intralayer HOMO-LUMO couplings. The present Dirac point is exotic since a conventional molecular Dirac electron system with only a single molecular orbital gives a line of the Dirac point extending in the Brillouin zone as shown in $\alpha$-(BEDT-TTF)2I3. From the calculation of the Dirac point for both fixed $k_{z}$ and fixed $k_{y}$, we found that the plane displaying the Dirac cone rotates with increasing $k_{z}$, which comes from the combined effect of the interlayer and intralayer matrix elements.

The Dirac point was analyzed using the parity of the wave function at the TRIM. The parity was calculated from the matrix of Eq. (LABEL:eq:P), which comes from the inversion symmetry with respect to the lattice site of the Pd atom and describes the difference in the symmetry between the HOMO and the LUMO. The behavior of the parity well explains the loop and the emergence of the Dirac point. The conditions of the Dirac point given by Eqs. (14a), (14b), and (14c) support the existence of a loop of the Dirac point. The parities of $E_{4}$ and $E_{5}$ are different for the $\Gamma$ point while they are the same for the Z point. This difference in the parity between the $\Gamma$ and Z points is compatible with the respective behaviors of the Dirac point at the $\Gamma$ and Z points, where, for a fixed $k_{z}$, a pair of Dirac points exists for the former case and two pairs of Dirac points exist for the latter case.

Here, an effective Hamiltonian for the bands $E_{4}(\bm{k})$ and $E_{5}(\bm{k})$ is briefly discussed to comprehend the behavior of the Dirac point. Noting that $P_{E4}(\Gamma)=-$ and $P_{E5}(\Gamma)=+$, we apply the method used for the case of $\alpha$-(BEDT-TTF)2I3.[18] After replacing $\hat{H}(\bm{k})$ of Eq. (2) by $\tilde{H}(\bm{k})(=\hat{P}(\bm{k})^{1/2}\hat{H}(\bm{k})\hat{P}(\bm{k})^{-1/2})$ to obtain the real matrix elements, $\tilde{H}(\bm{k})$ is rewritten as $\hat{H}(\bm{k})=\tilde{H}^{\rm HH}(\bm{k})+\tilde{H}^{\rm LL}(\bm{k})+\tilde{H}^{\rm HL}(\bm{k}),$ where the matrix elements of $\tilde{H}^{\rm HH}(\bm{k})$, $\tilde{H}^{\rm LL}(\bm{k})$, and $\tilde{H}^{\rm HL}(\bm{k})$ are expressed in terms of the transfer energies of HOMO-HOMO, LUMO-LUMO, and HOMO-LUMO, respectively. Note that $\tilde{H}^{\rm HH}(\bm{k})=\tilde{H}^{\rm HH}(-\bm{k})$ and $\tilde{H}^{\rm LL}(\bm{k})=\tilde{H}^{\rm LL}(-\bm{k})$, and that the relation $\tilde{H}^{\rm HL}(\bm{k})=-\tilde{H}^{\rm HL}(-\bm{k})$ comes from the difference in the symmetry between the LUMO and HOMO. Defining $|A>$ as the wave function for the maximum eigenvalue of $\tilde{H}^{\rm HH}(\bm{k})$ and $|B>$ as that for the minimum eigenvalue of $\tilde{H}^{\rm LL}(\bm{k})$, the 2 $\times$ 2 effective Hamiltonian is given by

[TABLE]

where $h_{11}(\bm{k})=<A|\tilde{H}^{\rm HH}(\bm{k})|A>$, $h_{22}(\bm{k})=<B|\tilde{H}^{\rm LL}(\bm{k})|B>$, and $h_{12}(\bm{k})=h_{21}^{*}(\bm{k})=<A|\tilde{H}^{\rm HL}(\bm{k})|B>$. The quantity $h_{12}(\bm{k})$ is determined by the combined effect of the interlayer and intralayer couplings.[4] Noting that $E_{4}(\bm{k})$ and $E_{5}(\bm{k})$ are eigenvalues of Eq. (16), the Dirac point ($E_{4}(\bm{k})=E_{5}(\bm{k})$) is obtained from $h_{11}(\bm{k})=h_{22}(\bm{k})$ and $h_{12}(\bm{k})=0$, which give the Fermi surface and nodal plane, respectively. The loop of the Dirac point is obtained by the intersection of these two planes, which are shown in Figs. 5(a)-5(d) on a 2D plane with a reduced zone. Figures 5(a) and 5(b) correspond to Figs. 2(a) and 2(b), respectively. For $k_{z}=0$, the Fermi line is an ellipsoid while the nodal line is given by $k_{x}=0$ (Fig. 5(a)). For $k_{y}=0$, the Fermi line is almost parallel to $k_{x}=0$ with a bottleneck for $k_{z}\sim\pi$, while there are two nodal lines, $k_{x}=0$ and $k_{z}-\pi\simeq-0.6\;k_{x}$ (Fig. 5(b)). With increasing $k_{z}$ the Fermi surface remains almost the same but the nodal line varies as follows. For $0.91<k_{z}/\pi<1.09$, there are two kinds of nodal lines as shown in Fig. 5(c), while only the left line intersects with the Fermi line for $0<k_{z}/\pi<0.91$. In Fig. 5(d), the case of $k_{z}/\pi=1$ is shown, where the line $k_{y}\simeq\pm 0.6\;k_{x}$ is consistent with the parity $P_{E4}(Z)=-$ and $P_{E5}(Z)=-$. These behaviors suggest that the nodal plane is periodic with respect to $k_{z}$ and the overlap between two nodal planes occurs in the interval region of $0.91<k_{z}/\pi<1.09$. It is found that the Dirac points obtained from Figs. 5(a)-5(d) well reproduce those of Fig. 3.

Acknowledgements.

The authors are grateful to T. Tsumuraya for useful discussions in the early stage of the present work. One of the authors (Y.S.) thanks T. Kariyado for useful comments and also C. Hotta and T. Osada for helpful comments. This work was supported by JSPS KAKENHI Grant Numbers JP15H02108, JP26400355, and JP16H06346.

Appendix A Matrix elements of Hamiltonian

The matrix elements of Eq. (3) are given by

[TABLE]

where $t_{\beta,\alpha}=t_{\alpha,\beta}^{*}$.