Lattice field theory simulations of Dirac semimetals

TL;DR

This paper uses lattice field theory simulations to analyze Dirac semimetals Na$_3$Bi and Cd$_3$As$_2$, confirming they are in the insulator phase under Coulomb interactions and exploring their phase diagram.

Contribution

It formulates a lattice field theory approach for Dirac semimetals that reproduces low-energy physics and studies their phase diagram with volume-independent results.

Findings

Na$_3$Bi and Cd$_3$As$_2$ are deep in the insulator phase.

The lattice theory reproduces effective low-energy Dirac semimetal behavior.

Results are volume independent, confirming previous findings.

Abstract

In this paper the observed Dirac semimetals Na$_3$Bi and Cd$_3$As$_2$ are studied within lattice simulation. We formulate lattice field theory with rooted staggered fermions on anisotropic lattice. It is shown that in the limit of zero temporal lattice spacing this theory reproduces low energy effective theory of Dirac semimetals. Using this lattice theory we study the phase diagram of Dirac semimetals in the plane effective coupling constant--Fermi velocity anisotropy. Within the formulated theory the results are practically volume independent in contrast with our previous study. Our results confirm our previous finding that within the Dirac model with bare Coulomb interaction both Na$_3$Bi and Cd$_3$As$_2$ lie deep in the insulator phase.

Tables

Table 1. Table 1: The critical value of β ~ ~ 𝛽 \tilde{\beta} and β 𝛽 \beta (Eq. ( 23 )) for different values of ξ = a t a s 𝜉 subscript 𝑎 𝑡 subscript 𝑎 𝑠 \xi=\frac{a_{t}}{a_{s}} (Eq. ( 13 )) and lattice sizes.

$\xi$	$L_s$	$L_t$	${\tilde{\beta }}_c$	${\beta }_c$
1/2	12	12	0.1133(7)	0.0566(4)
1/2	16	16	0.1127(5)	0.0564(2)
1/2	16	32	0.1126(4)	0.0563(2)
1/2	20	12	0.1106(3)	0.0553(2)
1/2	20	24	0.1115(2)	0.05573(14)
1/3	16	48	0.1946(6)	0.0648(2)
1/4	16	64	0.2791(7)	0.06977(16)
1/5	10	30	0.3438(27)	0.0688(5)
1/5	14	48	0.3505(11)	0.0701(2)
1/5	16	80	0.3474(9)	0.0695(2)
1/5	20	30	0.3420(9)	0.0684(2)
1/5	20	48	0.3468(9)	0.0697(2)
1/5	20	100	0.3487(5)	0.06975(11)

Table 2. Table 2: The critical value of β ~ c subscript ~ 𝛽 𝑐 \tilde{\beta}_{c} and β c subscript 𝛽 𝑐 \beta_{c} (Eq. ( 23 )) for different anisotropy of Fermi velocity ζ = v ∥ v ⟂ 𝜁 subscript 𝑣 parallel-to subscript 𝑣 perpendicular-to \zeta=\frac{v_{\parallel}}{v_{\perp}} and lattice sizes.

$\zeta$	$L_s$	$L_t$	${\tilde{\beta }}_c$	${\beta }_c$
1/2	10	40	0.408(3)	0.0816(6)
1/2	10	100	0.433(3)	0.0865(6)
1/2	10	160	0.435(2)	0.0869(5)
1/2	14	60	0.418(2)	0.0836(4)
1/2	16	80	0.4159(12)	0.0832(2)
1/2	16	100	0.4208(13)	0.0842(3)
1/2	16	120	0.4181(11)	0.0836(2)
1/4	12	60	0.476(3)	0.0952(6)
1/4	12	120	0.493(3)	0.0986(5)
1/4	16	80	0.4802(16)	0.0960(3)
1/4	16	100	0.4812(14)	0.0962(3)
1/4	16	140	0.4784(14)	0.0957(3)
1/10	16	100	0.5349(22)	0.1070(4)
1/10	16	160	0.5388(18)	0.1078(4)
1/10	16	200	0.5425(19)	0.1085(4)
1/10	24	140	0.5248(11)	0.10496(22)

Equations

$$E^2=v_{\parallel }^2(k_x^2+k_y^2)+v_{\perp }^2{k}_z^2,$$

$$\displaystyle Z=\int D\psi~{}D\bar{\psi}~{}DA_{4}~{}\exp{\bigl{(}-S_{E}\bigr{)}},$$

$$S_E=S_f+S_g.$$

$$S_g=\frac{1}{8\pi \alpha }\int d^{3}xdt({\partial }_i{A}_4{)}^2,$$

$$S_f$$

$$t\to$$

$$A_4\to$$

$$\begin{array}{r}{S}_E=\sum _{a=1}^{N_f=2}\int d^{3}xdt{\bar{\psi }}_a\left({\gamma }_4({\partial }_4+i{A}_4)+{\zeta }_i{\gamma }_i{\partial }_i\right){\psi }_a+\\ +\frac{1}{8\pi {\alpha }_{eff}}\int d^{3}xdt({\partial }_i{A}_4{)}^2.\end{array}$$

$$\begin{split}Z=\int DA_{4}~{}\exp{\bigl{(}-S_{g}\bigr{)}}\cdot\bigl{(}{\det}D[A_{4}]\bigr{)}^{N_{f}},\end{split}$$

$$\displaystyle P[A_{4}]=\exp{\bigl{(}-S_{g}\bigr{)}}\cdot\bigl{(}{\det}D[A_{4}]\bigr{)}^{N_{f}}$$

$$S_m=\sum _{a=1}^{N_f}\int d^{3}xdtm{\bar{\psi }}_a{\psi }_a.$$

$$\sigma =\frac{T}{V}\frac{\partial \ln Z}{\partial m}=\frac{T}{V}⟨\mathrm{Tr}{D}^{-1}[A_4]⟩$$

$$\displaystyle S_{f}=\sum_{x}a_{s}^{3}\biggl{\{}m\bar{\psi}_{x}\psi_{x}+$$

$$\displaystyle+\frac{1}{2}\sum_{\mu=1}^{4}\xi_{\mu}\eta_{\mu}(x)\bigl{[}(\bar{\psi}_{x}\psi_{x+\mu})-(\bar{\psi}_{x+\mu}\psi_{x})\bigr{]}\biggr{\}},$$

$$\begin{array}{rl}{\xi }_{\mu =1,2,3}& =\xi =a_t/a_s\\ {\xi }_{\mu =4}& =1.\end{array}$$

$$x_{\mu }=2y_{\mu }+{\eta }_{\mu },\mu =1,2,3,4.$$

$${\bar{q}}_y^{a\alpha }=\frac{1}{8}\sum _{\eta }{\Gamma }_{\eta ,a\alpha }^{+}{\bar{\psi }}_{2y+\eta },q^{\alpha a}=\frac{1}{8}\sum _{\eta }{\Gamma }_{\eta ,\alpha a}{\psi }_{2y+\eta },$$

$$S_f$$

$${\Delta }_{\mu }f$$

$${\delta }_{\mu }f$$

$$S_f$$

$${\gamma }_5{\gamma }_{\mu =1,2,3}=\left(\begin{array}{cc}i{\sigma }_{\mu }^T& 0\\ 0& -i{\sigma }_{\mu }^T\end{array}\right).$$

$$\det \hat{D}{|}_{a_t\to 0}=$$

$$=\left(\begin{array}{cc}A+{\sum }_{\mu =1,2,3}{B}_{\mu }i{\sigma }_{\mu }^T& 0\\ 0& A-{\sum }_{\mu =1,2,3}{B}_{\mu }i{\sigma }_{\mu }^T\end{array}\right)=$$

$$\displaystyle\quad=\det\biggl{(}A+\sum_{\mu=1,2,3}B_{\mu}i\sigma_{\mu}^{T}\biggr{)}\times$$

$$\displaystyle\hskip 86.72377pt\times\det\biggl{(}A-\sum_{\mu=1,2,3}B_{\mu}i\sigma_{\mu}^{T}\biggr{)}=$$

$$\displaystyle\quad=\biggl{(}\det\bigl{(}AA^{+}+\sum_{\mu=1,2,3}B_{\mu}^{2}\bigr{)}\biggr{)}^{2},$$

$$A$$

$$B_{\mu }$$

$$S_g=\frac{\tilde{\beta }}{2}\sum _{x,i}({\theta }_4(x)-{\theta }_4(x+i){)}^2.$$

$$\begin{array}{r}\tilde{\beta }=\beta \frac{a_s}{a_t}=\frac{v_{\perp }}{4\pi \alpha \xi },\\ \beta =\frac{1}{4\pi {\alpha }_{eff}},\end{array}$$

$${\theta }_4(x)\to {\theta }_4(x)+g(x+\hat{4})-g(x)$$

$$\begin{array}{r}{S}_f={\bar{\psi }}_x{D}_{x,y}{\psi }_y=\sum _x(m{\bar{\psi }}_x{\psi }_x+\\ +\frac{1}{2}{\eta }_4(x)[{\bar{\psi }}_x{e}^{i{\theta }_4(x)}{\psi }_{x+\hat{4}}-{\bar{\psi }}_{x+4}{e}^{-i{\theta }_4(x)}{\psi }_x]+\\ +\frac{1}{2}\sum _{i=1}^3{\xi }_i{\eta }_i(x)[{\bar{\psi }}_x{\psi }_{x+\hat{\text{i}}}-{\bar{\psi }}_{x+\text{i}}{\psi }_x]),\end{array}$$

$${\psi }_x\to e^{i\theta (x)}{\psi }_x,{\bar{\psi }}_x\to {\bar{\psi }}_x{e}^{-i\theta (x)}.$$

Full text

Semimetal, insulator, Coulomb interaction, Monte-Carlo simulations

I Introduction

Recent years are marked with the discovery of numerous materials with interesting properties. Considerable interest is attracted by the materials which low-energy fermionic excitations can be described by massless fermions. Among well-known examples is grapheneNovoselov666 ; Geim2007 - two-dimensional material with two effective massless Dirac fermionsPhysRev.71.622 ; PhysRev.104.666 ; Semenoff:1984dq ; Novoselov2005 ; Zhang2005 . Its three-dimensional analogues include Dirac(Na3BiLiu864 , Cd3As2Neupane2014 ; PhysRevLett.113.027603 ) and Weyl semimetalsXu613 ; Xue1501092 . Such materials provide a perfect opportunity for detailed study of quantum field theory phenomena, which were previously related to high-energy physics.

In this paper we are going to concentrate on studying Dirac semimetals. Low energy spectrum of the observed Dirac semimetals is determined by two Fermi points. In the vicinity of each Fermi point the fermion excitations reveal the properties of massless 3D Dirac fermions with the dispersal relation

[TABLE]

where $v_{\parallel},v_{\perp}$ are Fermi velocities in the $(x,y)$ plane and $z$ direction correspondingly. For the Na3Bi: $v_{\parallel}/c\simeq 0.001,v_{\perp}/v_{\parallel}\simeq 0.1$Liu864 and for the Cd3As2: $v_{\parallel}/c\simeq 0.004,v_{\perp}/v_{\parallel}\simeq 0.25$Liu2014 , where $c$ is the light velocity (further we will express all velocities in the units of $c$).

Electromagnetic interaction between quasiparticles might significantly modify the properties of these systems and thus give rise to a plenty of interesting phenomenaRevModPhys.84.1067 . The reason for that is the smallness of Fermi velocity of quasiparticles $v_{F}\ll 1$. As the result the interaction is reduced to instantaneous Coulomb with effective coupling constant $\alpha_{eff}=\alpha\cdot 1/v_{F}\sim 1$, where $\alpha\approx 1/137$ is the coupling constant of electromagnetic interaction. For instance, the effective coupling constants can be estimated as $\alpha_{eff}\sim 7$ for the Na3Bi and $\alpha_{eff}\sim 2$ for the Cd3As2. So, the coupling constants of both materials are large and they are strongly correlated. It is known that large coupling constant can cause dynamical generation of the fermion mass gap thus leading to the transition from semimetal to the insulator phaseKhveshchenko:2001zz ; Gamayun:2009em ; Sabio:2010yf . The study of this transition is interesting and important from theoretical and experimental points of view.

Since Dirac semimetals are strongly correlated materials investigation of their properties requires nonperturbative consideration. We believe that today the only way to carry out model independent study of strongly coupled systems is lattice Quantum Monte Carlo (QMC) simulation Montvay:1994cy . This approach fully accounts many body effects in strongly coupled systems. Lattice simulation was intensively used to study the properties of graphene Drut:2008rg ; Hands:2008id ; Armour:2009vj ; Drut:2009aj ; Ulybyshev:2013swa ; Boyda:2016emg . Lattice simulation was also used to study topological phase transitions in Weyl semimetalsYamamoto:2016rfr ; Yamamoto:2016zpx . In this paper we are going to use lattice simulation to study semimetals/insulator transition in Dirac semimetals.

The first lattice study of the semimetals/insulator transition in Dirac semimetals was carried out in paper Braguta:2016vhm . In particular, in this paper dynamical chiral symmetry breaking caused by strong interaction was addressed. The result of this study allowed to find critical coupling constant for the semimetals/insulator transition and draw the phase diagram in the plane $(\alpha_{eff},v_{\perp}/v_{\parallel})$. The calculation done in Braguta:2016vhm was carried out at two lattice volumes and exhaustive study of finite volume effects was not done.

In this paper we continue the study of dynamical chiral symmetry breaking caused by strong interaction in Dirac semimetals by means of lattice simulation. We develop an approach that allows to carry out lattice simulation of Dirac semimetals keeping finite volume uncertainties under control. Using this approach we calculate critical coupling for the semimetals/insulator transition as a function of the Fermi velocity asymmetry $v_{\perp}/v_{\parallel}$ and draw the phase diagram of Dirac semimetals.

This paper is organized as follows. In the next section we formulate low energy effective field theory of Dirac semimetals in continuum. Section III is devoted to the description of lattice field theory of Dirac semimetals. Using this lattice theory we study the phase diagram of Dirac semimetals in section IV. In last section we discuss our results.

II Effective field theory

Low-energy fermionic excitations of the discovered semimetals can be described as $N_{f}=2$ massless 3D Dirac fermions. Their dispersion relation is given by formula (1). Important features of the fermionic spectrum are the smallness of the Fermi velocity $v_{F}\ll 1$ and the anisotropy of the Fermi velocity in different directions $v_{\perp}/v_{\parallel}<1$. The smallness of the Fermi velocity allows us to state that magnetic interactions and retardation effects are considerably suppressed in Dirac semimetals. This implies that the electromagnetic interaction between quasiparticles is reduced to the instantaneous Coulomb potential.

The properties described above allow one to build low energy effective quantum field theory of fermionic excitations in Dirac semimetals. In particular, the partition function of theory can be written as a path integral

[TABLE]

where $\bar{\psi},\psi$ are fermion fields, $A_{4}$ is temporal component of the vector potential of the electromagnetic field. The Euclidean action $S_{E}$ can be written as a sum of the contributions of the gauge fields $S_{g}$ and the fermions $S_{f}$

[TABLE]

The gauge field contribution can be written as

[TABLE]

whereas the fermion contribution is

[TABLE]

In last two equations the following designations are introduced: the $\alpha\approx 1/137$ is the electromagnetic coupling constant, $v_{F,i}$ are the components of Fermi velocity: $v_{F,x}=v_{F,y}=v_{\parallel}$, $v_{F,z}=v_{\perp}$, the $\gamma_{\mu},~{}\mu=1,2,3,4$ are the standard Dirac matrices. The $\psi_{a}$ is four component Dirac spinor, the index $a$ counts the contributions of different Fermi points to the low energy effective action. For further use we also introduced the designation of the Dirac operator $D_{a}(A_{4})$. The fermions in the effective theory (5) are degenerate what allows one to omit the Fermi point index $a$ of the Dirac operator.

Note that the action (3) includes only time-like component $A_{4}$ of the electromagnetic field. As was explained above this is because the contribution of space-like components is suppressed by the factors of $v_{F}$ and can be neglected due to the smallness of Fermi velocity. The action $S_{E}$ is quadratic in the $A_{4}$ so this field can be integrated out what leads to the Coulomb interaction between fermions $\frac{\rho(\vec{r})\rho(\vec{r^{\prime}})}{|\vec{r}-\vec{r^{\prime}}|}$, where $\rho(\vec{r})=\bar{\psi}_{a}(\vec{r})\gamma_{4}\psi_{a}(\vec{r})$ is the charge density. Thus the partition function (2) describes the system of $N_{f}$ Dirac fermions which interact through instantaneous Coulomb law.

Let us now rescale the time coordinate and timelike component of the electromagnetic field:

[TABLE]

what brings us to the action

[TABLE]

In the above equation we introduced Fermi velocity anisotropy factors $\zeta_{i}=\frac{v_{F,i}}{v_{F}}$. It is seen from equation (7) that the strength of the interaction is determined by the effective coupling constant $\alpha_{eff}=\frac{\alpha}{v_{F}}$. For the observed Dirac semimetals the effective coupling constants can be estimated as $\alpha_{eff}\approx 7$ for the Na3Bi and $\alpha_{eff}\approx 2$ for the Cd3As2. So, the semimetals Na3Bi and Cd3As2 belong to the class of strongly coupled systems. The study of such systems is very complicated problem and today there is no analytical approach that can satisfactorily deal with it. For this reason in this paper we are going to apply lattice simulation of Dirac semimetals which exactly regards many body effects in strongly coupled systems.

Similarly to the gauge fields one can integrate out fermion fields from the partition function (3), what leads to the following expression

[TABLE]

Due to the $\gamma_{5}$-hermiticity: $\gamma_{5}D^{{\dagger}}[A_{4}]\gamma_{5}=D[A_{4}]$, the eigenvalues of the Dirac operator go in pairs $\pm i\lambda$. This allows us to state that the ${\det}D[A_{4}]$ is positively defined. The product

[TABLE]

can be interpreted as a probability measure and Monte Carlo methods can be used to calculate partition function.

The chiral symmetry of the action (3) is $SU_{R}(N_{f})\times SU_{L}(N_{f})$111Notice that there is also symmetry $U_{V}(1)\times U_{A}(1)$ of the action (3).. It is known that sufficiently large interaction strength dynamically breaks this symmetry according to the pattern $SU_{R}(N_{f})\times SU_{L}(N_{f})\to SU(N_{f})$Scherer:2002tk ; Pich:1998xt ; Bijnens:1995ww ; Klevansky:1992qe . This leads to the formation of excitonic condensate $\bar{\psi}\psi$ and appearance of dynamical fermionic mass (the mass gap). The system with dynamically generated mass gap is an insulator.

In this paper we are going to study the semimetal/insulator phase transition with dynamical chiral symmetry breaking within lattice simulation of Dirac semimetals. To do this we add the mass term to the action (3)

[TABLE]

Notice that this term violates $SU_{R}(N_{f})\times SU_{L}(N_{f})$ but conserves $SU(N_{f})$, thus playing a role of the seed of the insulator phase. Monte Carlo lattice simulation of Dirac semimetals will be carried out with finite mass $m$ which will be extrapolated to the chiral limit at the end of the calculation. If in the limit $m\to 0$ the chiral condensate is nonzero then the system is in the insulator phase. Otherwise the system is in the semimetal phase.

The order parameter of the semimetal/insulator phase transition with dynamical chiral symmetry breaking in the EFT is the chiral condensate $\sigma=\langle\bar{\psi}\psi\rangle$. If one introduces explicitly non-zero mass term to the action $S_{E}\to S_{E}+m\int d^{3}xdt\bar{\psi_{a}}\psi_{a}$, then the value of the chiral condensate is given by:

[TABLE]

III Lattice field theory for Dirac semimetals

III.1 Staggered fermions without interactions

In this section we are going to formulate low energy effective field theory of Dirac semimetals (7) on discrete lattice. It is known that naive discretization of the continuum fermion action (5) leads to appearance of 16 fermion copies per one lattice fermion. This is well known the fermion doubling problem Montvay:1994cy . To get rid of this problem we are going to use Kogut-Susskind staggered fermions Susskind:1976jm .

The low energy effective action for the staggered fermions (see below) is reduced to the action with four degenerate Dirac fermions. The quantum number which enumerates these degenerate fermions is called taste. Low energy theory of the observed Dirac semimetals contains two Dirac fermions. To get low energy effective action for two Dirac fermions from the action for the staggered fermions one can take square root from the fermion determinant. If the system has the taste symmetry which is the case for the low energy effective action of staggered fermions this procedure is well defined. However, below it will be seen that there are corrections to the low energy effective action of staggered fermions which violate taste symmetry of the fermion action. So, the procedure of taking square root of the staggered fermion determinant is well defined only in the continuum limit. We believe that this also leads to noticeable volume dependence of finite results which were observed in Braguta:2016vhm .

To formulate the theory of the observed Dirac semimetals based on the staggered fermions we are going to use lattice with different spatial and temporal lattice spacings. Below it will be shown that taking continuum limit in temporal direction at arbitrary fixed spatial lattice spacing the procedure of square rooting of the fermion determinant is well defined. The low energy effective action obtained in this case corresponds to two Dirac fermions and finite volume effects are small. Now let us proceed to the details of this approach.

To write the lattice fermion action for staggered fermions let us introduce a regular cubic lattice in four-dimensional space with spatial lattice spacing $a_{s}$ and temporal lattice spacing $a_{t}$. The number of lattice sites is $L_{s}$ in each spatial direction and $L_{t}$ in temporal direction. The sites of the lattice have coordinates $x_{\mu=1,2,3}=0,1,...,L_{s}-1,~{}x_{4}=0,1,...,L_{t}-1$. With this notations the Euclidean action for staggered fermions can be written as

[TABLE]

where $\bar{\psi}_{x},\psi_{x}$ are one component Grassmann fields and the sum is taken over all sites of the lattice, $\eta_{\mu}(x)=(-1)^{\sum_{i<\mu}x_{i}}$, the $m$ is the mass in lattice units, and

[TABLE]

Notice that action (12) is the action for free fermions. Interactions between fermions will be introduced below. In order to reduce the size of the formulas we will take $a_{s}=1$, restoring explicit spatial lattice spacing when necessary.

Action (12) can be considered as some tight-binding model on four-dimensional lattice. Let us determine low energy spectrum of this theory222Here we follow Montvay:1994cy .. To do this we assume that the lattice sizes are even: $L_{s}=2L_{s}^{\prime}$, $L_{t}=2L_{t}^{\prime}$ and divide our lattice to hypercube blocks of size $2^{4}$. The coordinates of the original lattice can be written as

[TABLE]

So, the coordinates $y_{\mu}=0,...,L^{\prime}_{\mu}-1$ label hypercube blocks and the coordinates $\eta_{\mu}=0,1$ label the sites within one hypercube block.

To proceed further we introduce the fields for four tastes of Dirac fermions: $\bar{q}^{a\alpha},q^{\alpha a}$ where the $\alpha=1,2,3,4$ is the Dirac index and the $a=1,2,3,4$ is the taste index. Now one can build the fields $\bar{q}^{\alpha a},q^{\alpha a}$ as a linear combinations of the fields $\bar{\psi}_{x},\psi_{x}$ living in one hypercube block

[TABLE]

where the matrix $\Gamma_{\eta}=\gamma_{1}^{\eta_{1}}\gamma_{2}^{\eta_{2}}\gamma_{3}^{\eta_{3}}\gamma_{4}^{\eta_{4}}$ and $\gamma_{\mu}$ are the Dirac matrices.

With these notations the fermion action can be written in the following form Montvay:1994cy

[TABLE]

In the last equation we used the following designations. $A\otimes B$ means that the matrix $A$ acts on the Dirac indices $\alpha$ and the matrix $B$ acts on the taste indices $a$ of the fermion fields. For further use we also introduced the Dirac matrix $\hat{D}$ which acts in the spaces of the coordinates, Dirac indices and taste indices. The operators $\Delta_{\mu}$ and $\delta_{\mu}$ have the form

[TABLE]

In continuum limit the operators $\Delta_{\mu}$ and $\delta_{\mu}$ are reduced to first and second derivative correspondingly: $\Delta_{\mu}f\to\partial_{\mu}f$, $\delta_{\mu}f\to\partial^{2}_{\mu}f$

Low energy effective action is given by the terms with the smallest order of spatial and temporal derivatives. Thus it is seen from formula (16) that low energy effective action is the sum of the mass term and the terms with the operators $\Delta_{\mu=1,2,3,4}$. So, low energy effective action for the staggered fermions is reduced to the action with four degenerate tastes of Dirac fermions.

In the observed Dirac semimetals low energy effective action contains two Dirac fermions instead of four in the staggered fermions. However, if the action is degenerate in tastes one can take square root of the fermion determinant $\det\hat{D}$ and thus obtain two Dirac fermions from four tastes of staggered fermions.

In addition to the terms degenerate in taste index in formula (16) one can see the terms which are proportional to the operator of the second derivative $\delta_{\mu}$. Although these terms only give corrections to the leading order effective theory, they are accompanied by the matrices $\gamma_{5}\gamma_{\mu}^{T}$ in the taste space. These matrices violate taste degeneracy in action (16). For this reason the procedure of taking square root of the staggered fermion determinant is well defined only in the continuum limit.

In order to get rid of the problem with rooting we are going to use an asymmetric lattice $a_{s}\neq a_{t}$. In particular, on the lattice with finite size of $a_{s}$ we will take limit $a_{t}\to 0$ or equivalently $\xi=a_{t}/a_{s}\to 0$.

To show that this trick helps to solve the problem with rooting, let us concentrate on equation (16). In the limit $a_{t}\to 0$ the term which is proportional to the operator $\delta_{4}$ can be omitted and action (16) becomes

[TABLE]

To proceed let us recall that the matrices $\gamma_{5}\gamma^{T}_{\mu=1,2,3}$ have the form

[TABLE]

Thus the fermion determinant in the taste space can be written in the following form

[TABLE]

where $A$ and $B$ are the matrices acting in the coordinate and Dirac index spaces

[TABLE]

Formula (20) is derived in Appendix.

Notice that the matrix $A$ represents Dirac operator for one Dirac fermion and the matrices $B_{\mu}$ parametrize corrections to the low energy effective action. From equation (20) one sees that the fermion determinant of staggered fermions in the limit $a_{t}\to 0$ represents two copies of the theory with the determinant $AA^{+}+\sum_{\mu=1,2,3}B_{\mu}^{2}$. Thus the square rooting of the staggered fermion determinant in the limit $a_{t}\to 0$ leaves us with one copy of the determinant $\det(AA^{+}+\sum_{\mu=1,2,3}B_{\mu}^{2})$. The terms $\sim B_{\mu}^{2}$ give only corrections to the low energy effective action which is determined by the operator $A$. Thus we have shown that in the limit $a_{t}\to 0$ the rooting procedure gives the low energy effective theory with two Dirac fermions.

III.2 Interacting staggered fermions

First let us consider the action for the electromagnetic field. For discretization of the electromagnetic field we use noncompact action:

[TABLE]

If one replaces $\theta_{4}$ with $a_{t}A_{4}$ and $\tilde{\beta}$ with

[TABLE]

then in the continuum limit one recovers gauge part of the action (7). Here $\beta$ is the inverse effective coupling constant and $\tilde{\beta}$ also includes the ratio of time- and space-like lattice steps $\xi$. In the limit of zero temporal lattice spacing $a_{t}\to 0$ $\tilde{\beta}\to\infty$ and $\beta$ is kept constant.

Note that the action (22) is gauge invariant under the gauge transformations which are described by arbitrary function $g(x)$, where $x$ is lattice coordinate:

[TABLE]

In previous section we considered free fermions. It causes no difficulties to write fermion action for interacting staggered fermions

[TABLE]

where $\theta_{4}(x)$ is the gauge field and $\xi_{i}$ is anisotropy factor. If the Fermi velocity is isotropic in spatial directions, then $\xi_{i}=\xi=a_{t}/a_{s}$. If there is an anisotropy in the Fermi velocity $v_{\perp}\neq v_{\parallel}$ then $\xi_{i}=\xi=a_{t}/a_{s},i=1,2$ and $\xi_{3}=\xi v_{\parallel}/v_{\perp}=(a_{t}v_{\parallel})/(a_{s}v_{\perp})$. Notice that under gauge transformation (24) the fermion fields $\bar{\psi}_{x},\psi_{x}$ transform as

[TABLE]

The requirement of gauge invariance of fermion action unambiguously determines the form of action (25).

Integration over fermionic degrees of freedom is performed explicitly and one gets the following effective action:

[TABLE]

This action corresponds to four degenerate fermion flavorsMontvay:1994cy instead of $N_{f}=2$ observed in Na3Bi and Cd3As2. The reduction $N_{f}=4\to 2$ is performed by the standard rooting procedure. The effective action used in the simulation is

[TABLE]

The generation of the electromagnetic fields $\theta_{4}(x)$ was performed by means of the standard Hybrid Monte-Carlo MethodMontvay:1994cy .

It should be noted that $a_{s}$ plays a role of ultraviolet cutoff in the studied effective field theory. Although the exact value of $a_{s}$ is not known, one expects that it should be of order of the distance between atoms in the crystal structure of the materialBraguta:2013rna . Lattice simulation with action (28) will be performed at few values of the anisotropy coefficient $\xi$ and then the limit $\xi\to 0$ will be taken.

IV Numerical results

In this section we are going to apply lattice field theory developed above to study the semimetal/insulator phase transition in the Dirac semimetals driven by strong Coulomb interaction. The order parameter for this phase transition is the chiral condensate $\sigma$, which is measured directly in our simulations. Its behaviour in the chiral limit $m\to 0$ characterizes the phase of the system. In the semimetal phase $\sigma(m\to 0)\to 0$, while in the insulator phase $\sigma(m\to 0)\neq 0$.

To measure the chiral condensate on the lattice we used stochastic estimator method:

[TABLE]

where $\eta_{x}$ are Gaussian distributed random variables: $\langle\eta_{x}\rangle=0$, $\langle\bar{\eta}_{x}\eta_{y}\rangle=\delta_{x,y}$. Note also that in formula (29) we averaged over configurations of gauge fields $\theta_{4}$. For each set of the parameters we generated $\sim O(200)$ configurations. For the inversion of the Dirac operator $\sum_{y}(D^{-1}[\theta_{4}])_{xy}\eta_{y}$ we used Conjugate Gradient Method.

To check finite volume effects we performed numerical simulations with different lattice sizes $L_{s}^{3}\times L_{t}$ (see Tab. 1 and  2). We have found, that for all values of lattice size under study, the measured values of the chiral condensate as well as the extracted value of $\beta_{c}$ coincide with good accuracy, thus we expect that finite volume effects are small. For our final predictions we took the results for $L_{s}=16$, while $L_{t}$ was different for different set of parameters.

IV.1 Isotropic Fermi velocity

First we present the results for the case of isotropic Fermi velocity $\xi_{1}=\xi_{2}=\xi_{3}=\xi$. In order to study the limit $a_{t}\to 0$ we performed simulations for four values of $\xi=1/2,1/3,1/4,1/5$. To control finite size effects simulations for the smallest and the largest value of $\xi=1/2,1/5$ were performed at different lattice sizes $L_{s}^{3}\times L_{t}$. We have found, that for lattice sizes used in our simulations the measured values of the chiral condensate coincide within the errorbars. Based on this observation we expect that finite volume effects are very small and can be neglected. Values of lattice sizes used in our simulations are presented in Tab. 1.

The dependence of the chiral condensate $\sigma$ on $\tilde{\beta}$ for different values of mass $m$ is presented in the Fig. 1. On the plots the data points for different lattice sizes almost coincide. All plots exhibit the similar behavior showing the formation of the chiral condensate and phase transition from semimetal to insulator phase. The critical value of $\tilde{\beta}$ can be roughly estimated to be $\tilde{\beta}_{c}(\xi=1/2)\sim 0.12$, $\tilde{\beta}_{c}(\xi=1/3)\sim 0.2$, $\tilde{\beta}_{c}(\xi=1/4)\sim 0.3$, $\tilde{\beta}_{c}(\xi=1/5)\sim 0.35$.

For more precise determination of the critical coupling constant we used the Equation of State(EoS), which was successfully applied to strongly coupled QED Gockeler:1996me and grapheneDrut:2008rg :

[TABLE]

where $X(\tilde{\beta})=X_{0}+X_{1}(1-\tilde{\beta}/\tilde{\beta}_{c})$, $Y(\tilde{\beta})=Y_{1}(1-\tilde{\beta}/\tilde{\beta}_{c})$ are expanded in the vicinity of critical $\tilde{\beta}_{c}$. For the functions $f_{1}$ and $f_{3}$ we allowed nonclassical critical exponents: $f_{1}(\sigma)=\sigma^{b}$, $f_{3}(\sigma)=\sigma^{\delta}$. Classical values of $b=1$ and $\delta=3$ can be easily visualized, if one plots $\sigma^{2}$ as a function of $m/\sigma$(Fisher plot). The resulting dependence $\sigma^{2}(m/\sigma)$ form straight lines, crossing the origin at critical value $\beta_{c}$. The deviations from straight lines might be attributed to non-classical values of $b$ and $\delta$. The Fisher plots for the considered values of $\xi$ are presented in Fig. 2. The lines on the plots correspond to the fit by EoS (30). One sees small deviations from classical values $b=1$, $\beta=3$. The critical values of $\tilde{\beta}_{c}$, obtained with the help of this fit, are also presented in Tab. 1. One can easily obtain critical value of inverse effective coupling by the formula $\beta_{c}=\tilde{\beta}_{c}\xi$.

We plot the dependence of the critical value $\beta_{c}$ on $\xi$ in the Fig. 3. Note, that the limit $\xi\to 0$ should be taken. One sees that for $\xi=1/4$ and $\xi=1/5$ the critical values $\beta_{c}$ are almost the same. As the final estimate we take half of the sum of $\beta_{c}$ obtained at $\xi=1/4$ (lattice size $16^{3}\times 64$) and $\xi=1/5$ (lattice size $16^{3}\times 80$). The obtained critical value is $\beta_{c}=0.06964(13)$, which corresponds to the critical coupling $\alpha_{c}=1.1414(21)$. Note that this critical coupling constant is smaller, than was obtained earlier in Braguta:2016vhm . We believe that this discrepancy can be attributed to finite volume effects which are under control in the present study.

IV.2 Fermi velocity anisotropy

In the discovered materials Na3Bi and Cd3As2 there is an anisotropy in the Fermi velocity $v_{\parallel}\neq v_{\perp}$. For this reason it is interesting to study how this anisotropy affects the phase diagram of the system.

The calculations presented in the previous section show that for $\xi=1/5$ we are close to the limit $a_{t}\to 0$. Accordingly, here we fix two values $\xi_{1}=\xi_{2}=\xi=1/5$ and vary $\xi_{3}=\zeta\xi$, choosing $\zeta\in(0;1)$. One can expect that in this case we are also close to the limit $a_{t}\to 0$. Note that it corresponds to the relation $v_{\parallel}=\zeta v_{\perp}$ between Fermi velocity in different directions.

For checking finite volume effects we also performed simulations at various lattice sizes, which are listed in the Tab. 2.

In Fig. 4 the dependences of the chiral condensate on $\beta$ as well as Fisher plots are shown for three values of $\zeta=\frac{1}{2},\frac{1}{4},\frac{1}{10}$. By fitting data to the same EoS (30), critical $\beta_{c}$ was extracted. Obtained values of critical $\beta_{c}$ are presented in Tab. 2. As in the isotropic case the presented results suggest that finite volume effects are small. Using these values the phase diagram of the system is drawn in Fig. 5. For this plot we took spatial lattice size $L_{s}=16$ and the largest temporal lattice size $L_{t}=120$ for $\zeta=1/2$, $L_{t}=140$ for $\zeta=1/4$, $L_{t}=200$ for $\zeta=1/10$.

V Conclusions and discussion

This paper is devoted to the study of electronic properties of recently discovered Dirac semimetals Na3Bi and Cd3As2. We formulate the low energy effective theory for these materials, which is described by $N_{f}=2$ massless 3D Dirac fermions with Coulomb interaction. We construct lattice discretized version of this theory and show that in the limit of zero temporal lattice spacing $a_{t}\to 0$ it reproduces low energy effective theory with the correct number of degrees of freedom. Methods for numerical simulation of this lattice theory are briefly summarized.

In the second part of the paper we studied the phase diagram of Dirac semimetals using supercomputer simulations of the formulated lattice effective field theory. In particular, we concentrate on the semimetal-insulator phase transition driven by Coulomb interaction. Systematic uncertainties due to finite volume effects and nonzero temporal lattice spacing $a_{t}\neq 0$ are under control. By means of numerical simulations we determine the value of the effective critical coupling and its dependence on the Fermi velocity anisotropy. Based on these results we draw tentative phase diagram for the Dirac semimetals.

The phase diagram of the Dirac semimetals discussed in this paper was studied earlier in the strong coupling limitSekine:2014yna ; Araki:2015php , by means of Dyson-Schwinger equationsGonzalez:2015iba and within ladder approximationGonzalez:2015tsa . It is interesting to discuss their results and to compare them with ours. In Araki:2015php it was shown that in the limit of infinitely strong coupling $\beta=0$ for all values of Fermi velocity anisotropy $v_{\perp}/v_{\parallel}$ there is a gap, which is in line with our results. Note that the authors of Araki:2015php also discussed the importance of anisotropic lattices $a_{t}/a_{s}\to 0$. The authors of Gonzalez:2015iba studied the phase diagram of Dirac semimetals without Fermi velocity anisotropy within Dyson-Schwinger equations. The critical value of the coupling constant obtained by such approach is much larger $\alpha\approx 14.7$. Later the same authors revised this result within ladder approximation Gonzalez:2015tsa and obtained the critical coupling to be equal $\alpha_{eff}^{c}=1.8660$ which is close to the value $\alpha_{c}=1.1414(21)$ obtained here.

It is also worth mentioning the paper 2017arXiv170203551S , where the authors studied the graphene under uniaxial strain, leading to anisotropic Fermi velocity. The results, obtained in this paper are in qualitative agreement with ours: large Fermi velocity anisotropy $v_{\perp}/v_{\parallel}\ll 1$ leads to smaller effective coupling constant.

Finally we would like to mention that the results obtained in this paper are approximately by a factor of $\sim 1.5$ smaller than that obtained in our previous paper Braguta:2016vhm for all values of the Fermi velocity anisotropy. We believe that the discrepancy can be attributed to the finite volume effects which are under control in the present paper.

Despite this discrepancy our main conclusions remain the same. The first one is that the Fermi velocity anisotropy leads to decrease of critical coupling constant. The second one is that the observed Dirac semimetals Na3Bi and Cd3As2 lie deep in the insulator phase what contradicts to experiment. A possible resolution of this discrepancy is that in real world the interaction potential at small distances is screened by bound electrons what can lead to the shift of the critical coupling to larger values. Another possible explanation is that due to the renormalization effects strong interaction can considerably modify the basic parameters of the theory. Although the study of different explanation of the raised question is very important it is beyond the scope of this paper.

Acknowledgments

The authors would like to thank Z.V. Khaidukov for useful discussions. The work on numerical simulation and the determination of the critical coupling constants at different Fermi velocity asymmetries, and the formulation of lattice field theory on anisotropic lattice were supported by the RSF grant under contract 16-12-10059. The work of M.I.K. was supported by Act 211 of the government of the Russian Federation, Contract No. 02.A03.21.0006. Numerical simulations were carried out on GPU cluster of NRC Kurchatov Institute and at MSU supercomputer ”Lomonosov”.

Appendix A

In Appendix we are going to derive formula (20). We start from expression

[TABLE]

where the following designation is used

[TABLE]

Notice that the Pauli matrices $\sigma_{\mu}$ and the operator $\hat{\delta}$ act in the space of taste indices.

In the calculation we used Euclidean $\gamma$-matrices which are defined $\gamma_{4}=\gamma_{0}^{M},~{}\gamma_{\mu}=i\gamma_{\mu}^{M},~{}\mu=1,2,3$, where $\gamma^{M}_{\mu}$ are $\gamma$-matrices in Minkowski space. Expanding $\gamma$-matrices one can write

[TABLE]

where the Pauli matrices $\tau_{\mu}$ act in the space of spinor indices. Now it is evident that the diagonal operators $m\pm\Delta_{4}$ (which are unity operators in spin and taste indices) commute with the non-diagonal ones. For this reason one can write

[TABLE]

where the unity matrices ${\bf 1}$ act in spinor and taste spaces. Evidently last equation can be rewritten as

[TABLE]

Similarly

[TABLE]