Effective interactions between local hopping modulations on the square lattice

TL;DR

This paper investigates how local hopping modulations influence free fermions on a square lattice with flux and Ising degrees of freedom, revealing a first-order ferromagnetic tendency and long-range Ising interactions in the effective spin model.

Contribution

It introduces a novel effective spin model for fermions coupled to local hopping modulations, highlighting a first-order ferromagnetic term absent in traditional RKKY results.

Findings

First-order ferromagnetic interaction dominates the effective spin model.

Long-range Ising interactions can control the ground state energetics.

Perturbation theory reveals differences from standard RKKY results.

Abstract

We address the problem of free fermions interacting with frozen gauge fields. In particular, we consider a tight-binding model of fermions on the square lattice in which (i) flux 0 or $\pi$ is threaded through each plaquette and (ii) each nearest-neighbor link is decorated with an Ising degree of freedom that describes the local modulation of the hopping amplitude. Following the standard Ruderman--Kittel--Kasuya--Yosida (RKKY) approach, we compute an effective spin model in the coupling strength order by order. Unlike the original RRKY result for site-centered SU(2) spins in which the leading contribution is an effective exchange term at the second-order, perturbation theory in link-centered Z$_2$ case produces a first-order term that favors a collective ferromagnetic moment. If, by some means, an antiferromagnetic configuration can be stabilized, the energetics of ground state is controlled by an effective Ising interaction acting pairwise at the long range across the system.

Tables

Table 1. Table 1: Summarizing coefficients of the effective energy are calculated in both semi-analytic and exact diagonalization methods of the 0- and π 𝜋 \pi -flux models:

Models	0-flux		$\pi$-flux
Methods	$1^{\text{st}}$ order	$2^{\text{nd}}$ order	$1^{\text{st}}$ order	$2^{\text{nd}}$ order
Semi	$-0.811$	$-0.890$	$-0.958$	$-0.48$
Exact	$-0.811$	$-1.02$	$-0.958$	$-0.484$
% error	0.0%	12.7%	0.0%	1.0%

Equations

$$\text{centering}\chi (\mathbf{q})\text{@add@centering}$$

$$\chi (\mathbf{R})$$

$$\begin{array}{rl}\hat{H}& =\sum _{⟨m,n⟩}-t_{mn}(1+\xi {\sigma }_{mn}^z)(c_m^{†}{c}_n+c_n^{†}{c}_m)\\ & =-\sum _{⟨m,n⟩}{t}_{mn}(c_m^{†}{c}_n+c_n^{†}{c}_m)\\ & -\sum _{⟨m,n⟩}{t}_{mn}\xi {\sigma }_{mn}^z(c_m^{†}{c}_n+c_n^{†}{c}_m)\\ & \equiv {\hat{H}}_0+{\hat{H}}_1.\end{array}$$

$${\hat{H}}_{\text{K}}$$

$$E_{0-\text{flux}}^0(\mathbf{k})$$

$$E^{1}_{0-\text{flux}}(\mathbf{q})=\frac{-\xi}{N}\sum_{\mathbf{k}}2n_{F}(E_{\mathbf{k}})\\ \times\Big{[}\cos{k_{x}}\sigma^{z}_{\mathbf{q},1}+\cos{k_{y}}\sigma^{z}_{\mathbf{q},2}\Big{]}\delta_{\mathbf{q},\mathbf{0}},$$

$$\text{centering}{E}_{0-\text{flux}}^2(\mathbf{q})\text{@add@centering}$$

$$J_{\alpha \beta }^{\text{RKKY}}(\mathbf{q},{\Omega }_n)$$

$${\chi }^{\text{F}}(\mathbf{k},\mathbf{q},{\Omega }_n)$$

$$E_{(2,1),\pi-\text{flux}}^{0}(\mathbf{k})=\pm\big{[}2t(\cos^{2}{k_{x}}+\cos^{2}{k_{y}})\big{]}^{1/2}.$$

$$\begin{split}E^{1}_{\pi-\text{flux}}(\mathbf{q})&=\frac{-\xi}{N}\sum_{\mathbf{k}}n_{F}(E_{1,\mathbf{k}})\Bigg{[}\frac{2\cos{k_{x}}u^{2}(\mathbf{k})}{\varv^{2}(\mathbf{k})}\sigma^{z}_{\mathbf{q},1}\\ &\quad\frac{[1+\cos{2k_{y}}]u(\mathbf{k})}{\varv^{2}(\mathbf{k})}(\sigma^{z}_{\mathbf{q},2}+\sigma^{z}_{\mathbf{q},4})\\ &\quad\frac{\cos{k_{x}}[1+\cos{2k_{y}}]}{\varv^{2}(\mathbf{k})}\sigma^{z}_{\mathbf{q},3}\Bigg{]}\delta_{\mathbf{q},\mathbf{0}}.\end{split}$$

$$E_{\pi -\text{flux}}^2(\mathbf{q})$$

$$J_{\alpha \beta }^{\text{RKKY}}(\mathbf{q},{\omega }_n)$$

$${\chi }^{\text{D}}(\mathbf{k},\mathbf{q},{\omega }_n)$$

$$J_{\alpha \beta }^{\text{RKKY}}(\mathbf{R})={\int }_{\mathbf{q}\in \text{BZ}}\frac{d{q}_{x}d{q}_y}{(2\pi {)}^2}{e}^{-i\mathbf{q}\cdot \mathbf{R}}{J}_{\alpha \beta }(\mathbf{q},{\omega }_n),$$

$$J_{11}(\mathbf{q}\to 0,{\omega }_n=0)$$

$$\begin{array}{rl}{E}_{0-\text{flux}}^{\text{RKKY}}& \approx J_{\text{NN}}\sum _{⟨\mu ,\nu ⟩}{\sigma }_{\mu }^z{\sigma }_{\nu }^z+J_{\text{NNN}}\sum _{⟨\zeta ,\eta ⟩}{\sigma }_{\zeta }^z{\sigma }_{\eta }^z.\end{array}$$

$$E^{\text{total}}=E_0+\xi A\sum _{\alpha }{\sigma }_{\alpha }^z+{\xi }^{2}B\sum _{\alpha <\beta }{\sigma }_{\alpha }^z{\sigma }_{\beta }^z,$$

$$\text{centering}{\hat{H}}_0=-t\sum _{\mathbf{R}}[(c_{\mathbf{R}+{\mathbf{a}}_1}^{†}{c}_{\mathbf{R}}+c_{\mathbf{R}}^{†}{c}_{\mathbf{R}+{\mathbf{a}}_1})+(c_{\mathbf{R}+{\mathbf{a}}_2}^{†}{c}_{\mathbf{R}}+c_{\mathbf{R}}^{†}{c}_{\mathbf{R}+{\mathbf{a}}_2})],\text{@add@centering}$$

$$\text{centering}{\hat{H}}_1=-\xi \sum _{\mathbf{R}}[{\sigma }_{\mathbf{R},1}^z(c_{\mathbf{R}+{\mathbf{a}}_1}^{†}{c}_{\mathbf{R}}+c_{\mathbf{R}}^{†}{c}_{\mathbf{R}+{\mathbf{a}}_1})+{\sigma }_{\mathbf{R},2}^z(c_{\mathbf{R}+{\mathbf{a}}_2}^{†}{c}_{\mathbf{R}}+c_{\mathbf{R}}^{†}{c}_{\mathbf{R}+{\mathbf{a}}_2})].\text{@add@centering}$$

$$c_{\mathbf{R}}^{†}$$

$${\sigma }_{\mathbf{R},\alpha }^z$$

$$\hat{H}$$

$$\hat{H_{0}}=-t\sum_{\mathbf{R}}\big{[}(c_{\mathbf{R}+{\mathbf{a}_{1}},\text{A}}^{\dagger}c_{\mathbf{R},\text{A}}+c_{\mathbf{R}-{\mathbf{a}_{1}},\text{A}}^{\dagger}c_{\mathbf{R},\text{A}})+(c_{\mathbf{R},\text{A}}^{\dagger}c_{\mathbf{R},\text{B}}+c_{\mathbf{R},\text{B}}^{\dagger}c_{\mathbf{R},\text{A}})-(c_{\mathbf{R}+\mathbf{a}_{1},\text{B}}^{\dagger}c_{\mathbf{R},\text{B}}+c_{\mathbf{R}-\mathbf{a}_{1},\text{B}}^{\dagger}c_{\mathbf{R},\text{B}})+(c_{\mathbf{R}-{\mathbf{a}_{2}},\text{B}}^{\dagger}c_{\mathbf{R},\text{A}}+c_{\mathbf{R}+{\mathbf{a}_{2}},\text{A}}^{\dagger}c_{\mathbf{R},\text{B}})\big{]},$$

$${\hat{H}}_1$$

$$\begin{array}{rl}{J}_{11}(\mathbf{k},\mathbf{q})& =2+e^{i(2k_x+q_x)}+e^{-i(2k_x+q_x)}=4{\cos }^2(k_x+\frac{q_x}{2}),\\ J_{22}(\mathbf{k},\mathbf{q})& =2+e^{i(2k_y+q_y)}+e^{-i(2k_y+q_y)}=4{\cos }^2(k_y+\frac{q_y}{2}),\\ J_{12}(\mathbf{k},\mathbf{q})& =e^{-i(k_x+k_y+q_y)}+e^{i(k_x+k_y+q_x)}+e^{-i(k_x-k_y)}+e^{i(k_x-k_y+q_x-q_y)},\\ J_{21}(\mathbf{k},\mathbf{q})& =e^{i(k_x+k_y+q_y)}+e^{-i(k_x+k_y+q_x)}+e^{i(k_x-k_y)}+e^{-i(k_x-k_y+q_x-q_y)}.\end{array}$$

$$u(\mathbf{k})$$

$$v(\mathbf{k})$$

$$u(\mathbf{k},\mathbf{q})$$

$$v(\mathbf{k},\mathbf{q})$$

$$\begin{array}{rlr}{M}_{11}(\mathbf{k},\mathbf{q})& =& 4{\cos }^2(k_x+\frac{q_x}{2})\frac{{\cos }^2(k_y){\cos }^2(k_y+q_y)}{v^2(\mathbf{k})v^2(\mathbf{k}+\mathbf{q})},\end{array}$$

$$M_{33}(\mathbf{k},\mathbf{q})=\big{[}2+2\cos(2k_{x}+q_{x})\big{]}\frac{u^{2}(\mathbf{k})u^{2}(\mathbf{k}+\mathbf{q})}{v^{2}(\mathbf{k})v^{2}(\mathbf{k}+\mathbf{q}},$$

$$M_{13}(\mathbf{k},\mathbf{q})$$

$$M_{31}(\mathbf{k},\mathbf{q})$$

$$M_{14}(\mathbf{k},\mathbf{q})$$

Full text

I Introduction

Novel properties of two dimensional (2D) lattices such as the integer quantum Hall effect in graphene [1] and superconducting phenomenon in the $\mathrm{CuO_{2}}$ plane (high-$T_{\text{c}}$ superconductors) are attractive to investigate in detail. Recently, the observation of practical single-layer ferromagnet with the finite Curie temperature $T_{\text{C}}$ can be applied for fabricating spintronic and magneto-electric devices [2]. The tight-binding models in the square lattice with adding some generic interacting terms are used to describe the emergent properties of the superconducting order in cuprate- and iron-based superconductors [3].

Observing gapless Fermi surface in the cuprate compounds by “resonating valence bond” theory was the primary reason that Affleck et al. proposed the $\pi$-flux model on the square lattice [4]. By applying perpendicularly the constant magnetic field through the lattice, the hopping integral between fermions at the adjacent sites is coupled with the gauge-invariant magnetic flux [5]. When the total flux passing through each plaquette is exact $\pi$, the gauge-inequivalent state is formed, and the unit cell tunes from single basis into bipartite one. The Fermi surface with finite volume shrinks into four Dirac points in the energy dispersion [6; 7].

The models of 0-flux [3] and the $\pi$-flux [8; 7; 9] square lattices reemerging in the context of coupling fermions with the transverse field Ising spins show various exotic ground-state phases and quantum phase transitions. In those models, an Ising spin variable on the link modifies the hopping amplitude of fermions at the neighboring sites. The first example of Ising-nematic quantum phase in the square lattice has been proposed to understand the formation of phase diagram of the Fe-based superconductors [3]. Second, Assaad et al. has presented the existence of different ground-state phases and quantum phase transitions by changing the number of fermions per site [8]. The rich variety of quantum phase transitions such as the first-order transition of two different topological ferromagnetic orders [7] and the transition between BCS (Bardeen–Cooper–Schrieffer model) and BEC states (Bose–Einstein condensation) [9] have also been observed in these models. However, none of above works show the effective interaction of Ising spins at the weak coupling limit between fermion and Ising spin.

The Ruderman–Kittel–Kasuya–Yosida (RKKY) interactions between localized spins are formed indirectly via nonlocal electrons [10; 11; 12]. That is one of the most important prototypes to explain the formation of magnetic order in pure rare-earth elements (e.g Gd, Sm, and Dy), their alloys, heavy fermion materials, diluted magnetic semiconductors, and impurities in graphene [13; 14; 15; 16; 17]. Therefore, the derivation of RKKY interaction is varied from one system to the other with a main spin susceptibility, or static Lindhard function in the momentum space $\mathbf{q}$ and the real space $\mathbf{R}$:

[TABLE]

In the Eq. (1), the momentum space of Lindhard function $\chi(\mathbf{q})$ is calculated by integrating over the first Brillouin zone (BZ), where $n_{F}(E_{\mathbf{k}})$ is the Fermi-Dirac distribution function, and $E(\mathbf{k})$ is the energy dispersion relation in the momentum $\mathbf{k}$. The spin susceptibility consists the singular or maximum point in the momentum space that define the magnetic ordering vector of spin system. The real-space susceptibility $\chi(\mathbf{R})$ is achieved by taking the Fourier transformation of the momentum space one with spatial dimensions $d$ Eq. (2). The real-space $\chi(\mathbf{R})$ function is oscillatory and decaying with the distance $R$ of two localized spins. For example, the conventional real RKKY interaction of spins showed the sign-changing oscillation and decaying rate of $1/R^{d}$ ($d$ is the dimension of crystal lattice) [18]. However, real RKKY interaction between magnetic impurities in the graphene shows no sign-changing oscillation due to vanishing of the Fermi surface, and the $1/R^{3}$ decaying rate instead of $1/R^{2}$ for 2D [17; 15]. The microscopic magnetic interactions of impurities in the graphene are known as FM and AFM couplings for spins on the same and different sublattices, respectively. The exact macroscopic magnetic order of impurity-doped graphene is quite controversial discussion such as AFM [19; 20] and FM orders [21].

This letter, we have constructed the total Hamiltonian, including fermions at the sites and Ising spins at the bonds, of the 0-flux and $\pi$-flux square lattices. We solve problem by two different approaches: semi-analytic integration and exact diagonalization. In the semi-analytic method, we separate the Hamiltonian to unperturbing (pure hopping between fermions) and interacting parts (coupling between fermion and Ising spin). Using the perturbation theory, the interacting Hamiltonian is treated at the weak limit. The first- and second-order coupling terms are derived analytically, then computed numerically in the momentum and real spaces. Since the frustrating effect emerges in our models due to competing effect between nearest and next-nearest couplings, the exact diagonalization method is used to verify semi-analytic calculations. Based on microscopic discussions, we will determine the final magnetic order and the ground-state energy for each case.

II Model and Methods

II.1 Model

The general Hamiltonians $\hat{H}$ of the 0- and $\pi$-flux square lattices are formed:

[TABLE]

Here, spin-1/2 fermion $c^{\dagger}_{m}=(c^{\dagger}_{m\uparrow},c^{\dagger}_{m\downarrow})$ lives at the site $m$ of the square lattice, and $\langle m,n\rangle$ is the nearest-neighbor pair of the orbital fermions. The Ising spin $\sigma_{mn}^{z}$ is positioned at the bond between two fermions $m$ and $n$ – link-centered $\text{Z}_{2}$ model. For the 0-flux lattice, there are two Ising spins in each unit cell $i$: $\sigma^{z}_{i,1}$ and $\sigma^{z}_{i,2}$ (Fig. 1a). We set the hopping integral $t_{mn}=t$ and the coupling parameter $t_{mn}\xi$ ($t$ and $\xi$ are real positive values). For the $\pi$-flux model, there are four Ising spins in each unit cell: $\sigma^{z}_{i,1}$, $\sigma^{z}_{i,2}$, $\sigma^{z}_{i,3}$ and $\sigma^{z}_{i,4}$ (Fig. 1b). Along the solid line in the Fig. 1b, $t_{mn}=t$, and along the dotted line $t_{mn}$ is reversed sign due to the $\pi$-flux passing through each plaquette. (Note: in 0-flux lattice, the number of site equals to number of unit cell, so $i=m$. However, in $\pi$-flux lattice, the number of site is doubled number of unit cell, so two values are different.)

We separate the full Hamiltonian into the tight-binding part $\hat{H}_{0}$ and interacting part $\hat{H}_{1}$ with the small perturbation variable $\xi$. The Hamiltonian is written explicitly in the real space with the lattice translation vector $\mathbf{R}$, then transformed into the momentum space $\mathbf{k}$ (see Appendix A). In our model, the interacting term represents coupling between Ising spins and the itinerant fermions compared with the spin of fermions and localized spins in Kondo lattice – site-centered SU(2) case:

[TABLE]

where $\mathbf{s}_{m}$ is spin of itinerant electron, and $\mathbf{S}_{m}$ is the localized spin at the position $m$. This term is used to derive the second-order conventional spin susceptibility of RKKY interaction [22].

Below, we show the main steps of two calculating methods.

II.2 Semi-analytic integration

II.2.1 Momentum-dependent 0-flux model

The momentum-dependence of the non-interacting Hamiltonian $\hat{H}_{0}$ is diagonalized to give the energy dispersion $E^{0}_{0-\text{flux}}(\mathbf{k})$:

[TABLE]

In the 0-flux model, the expectation energies of the interacting Hamiltonian $\hat{H}_{1}(\mathbf{k},\mathbf{q})$ are calculated using perturbation theory with the Fermi-sea ground state. The first-order energy is found:

[TABLE]

where $N$ is number of unit cells in the lattice, and $n_{F}(E_{\mathbf{k}})$ is Fermi–Dirac function. The second-order perturbation provides a $2\times 2$ interacting matrix $J^{\text{RKKY}}_{\alpha\beta}(\mathbf{q},\omega_{n})$ inside the effective energy $E^{2}_{0-\text{flux}}(\mathbf{q})$:

[TABLE]

Here, $\alpha$ and $\beta$ are labeled of Ising spin in the unit cell of 0-flux lattice (Fig. 1a). $\chi^{\text{F}}(\mathbf{k},\mathbf{q},\omega_{n})$ is the Lindhard function for Fermi metallic band or intra-band interaction. All elements of $J_{\alpha\beta}(\mathbf{k},\mathbf{q})$ matrix are listed in the Sect. 1 of Appendix B.

To compute $\mathbf{q}$-space RKKY interaction, we take integration over the whole square Brillouin zone for momentum $\mathbf{k}$ and diagonalize $J^{\text{RKKY}}_{\alpha\beta}(\mathbf{q},\omega_{n})$ matrix to obtain two eigenvalues Eigen1 and Eigen2. We divide $k_{x}$ and $k_{y}$ from $-\pi$ to $\pi$ into $N$ intervals, and $N=L_{x}=L_{y}=400$ ($L_{x}$ and $L_{y}$ are the lengths of real lattice). For Matsubara frequency $\Omega_{n}=2nT$, we select the value of $\Omega_{n}/T=10^{-5}$ and integer number $n$. Two eigevalues as the functions of $\mathbf{q}$ are plotted along high symmetric points in the Brillouin zone such as $\Gamma=(0,0)$, $M=(\pi,0)$ and $K=(\pi,\pi)$ in Fig. 2a.

II.2.2 Momentum-dependence of the $\pi$-flux model

We transform the unperturbed Hamiltonian $\hat{H}_{0}$ into the momentum space and diagonalize it to obtain the energy dispersion with two bands:

[TABLE]

Similar to the 0-flux model, the interacting Hamiltonian of $\pi$-flux model gives the first- and second-order effective energies. The first-order is:

[TABLE]

Here, $u(\mathbf{k})$ and $\varv(\mathbf{k})$ functions are defined in the Sect. 2 of Appendix B.

The second-order effective energy $E^{2}_{\pi-\text{flux}}(\mathbf{q})$ of the interacting Hamiltonian includes a $4\times 4$ matrix with common Lindhard function:

[TABLE]

Here, $\alpha$ and $\beta$ are labeled of Ising spins in the unit cell of $\pi$-flux lattice (Fig. 1b). $\chi^{\text{D}}(\mathbf{k},\mathbf{q},\omega_{n})$ is the Lindhard function for Dirac semimetallic bands. The values of $s$ and $s^{\prime}$ represent the energy bands. For the half-filling case, there is only one situation with $s=1$ and $s^{\prime}=2$ results in the non-zero value of Lindhard function, or it is interband interaction [20]. (All the terms of matrix $M_{\alpha\beta}(\mathbf{k},\mathbf{q})$ are in the Sect. 3 of Appendix B.)

We find the eigenvalue spectrum of $J^{\text{RKKY}}_{\alpha\beta}(\mathbf{q},\omega_{n})$ by taking integration over the rectangular Brillouin zone and diagonalizing $4\times 4$ functional matrix. Analogous to the 0-flux model, that function shows four different eigenvalues such as Eigen1, Eigen2, Eigen3 and Eigen4 that plot along the high symmetric points such as $\Gamma=(0,0)$, $X=(\pi,0)$, $Y=(0,\pi/2)$, and $K=(\pi,\pi/2)$ in Fig. 2b.

II.2.3 Real-space calculation

Previous formulas are written in the momentum space. We compute numerically the real-space RKKY interaction between Ising spin pair by taking the Fourier transformation of the momentum space of $J_{\alpha\beta}(\mathbf{q},\omega_{n})$ function:

[TABLE]

where $\mathbf{R}=\mathbf{R}_{j}-\mathbf{R}_{i}$ is the distance vector of two $i$ and $j$ unit cells. For the 0-flux lattice, the $\alpha,\beta=1$ and 2 (or Ising spins 1 and 2). For the $\pi$-flux lattice, the $\alpha,\beta=1,2,3$ and 4 (or Ising spins 1, 2, 3 and 4). All interacting functions are plotted via the real distance $R=|\mathbf{R}|$ along the x- or y-directions. We use the lattice size of $L_{x}=L_{y}=160$, and $L_{x}=160$ and $L_{y}=80$ for 0- and $\pi$-flux models, respectively. From those calculations, we have constructed the real-space effective interacting energy of the system.

II.3 Exact diagonalization

Following the illustrated square lattices in the Fig. 1, the real Hamiltonian matrices are constructed based on the hopping amplitude $t$ and coupling parameter $\xi$ with their sizes of $N=L_{x}\times L_{y}=160\times 160$. With the 0-flux model, the hopping terms along the x- and y-directions are $-t-\xi\sigma^{z}_{1}$ and $-t-\xi\sigma^{z}_{2}$, respectively. We set the value of $t=1$, and the total energy $E^{\text{total}}=\langle\Psi_{\text{gs}}|\hat{H}|\Psi_{\text{gs}}\rangle$ depends on the Ising configuration and parameter $\xi$.

Since the lattice size is large $L_{x}\times L_{y}=160\times 160$, the number of Ising spin in this lattice is $N_{\text{Ising}}=51200$ spins. It is impossible to find the minimal energy of the system by optimizing over all Ising spin configurations. Because the definition of the magnetic ordering vector $\mathbf{Q}$ would be at the high symmetric points of Brillouin zone, we have found the period of Ising spin configuration using formula $\mathbf{Q}\cdot\mathbf{r}_{m}=2\pi m$ with integer number $m$ [23; 24]. For example, in the 0-flux model, if we choose the magnetic ordering vector is at the M point in Fig. 2a, $\mathbf{Q}=(\pi,0)$, and set up Ising 1: spin-up (or $\sigma^{z}_{1}=1$) and Ising 2: spin-down (or $\sigma^{z}_{1}=-1$), all the spins are reversed directions at next unit cell on the right. Along the y-direction, they copy the similar configuration. The obtaining result of numerical method is used to compare with the semi-analytic calculation one to determine the correct spin orders for both models.

III Results and Discussions

III.1 Energy dispersions

Figure 2 shows the 2D contour plots of Fermi band and Dirac bands for the 0-flux and $\pi$-flux lattices, respectively. The energy dispersion $E^{0}_{0-\text{flux}}=-2(\cos{k_{x}}+\cos{k_{y}})$ of 0-flux lattice exhibits the continuous metallic band and diamond shape Fermi surface (red line in Fig. 2a). (The 3D surface is also plotted in my MSc dissertation [25]). At the half-filling or chemical potential $\mu=0$, each site of lattice is occupied exactly one electron. The Fermi surface of the 0-flux square lattice provides nesting property. There are nesting vector $\mathbf{Q}_{0}=(\pm\pi,\pm\pi)$ that connect all points on the Fermi surface [22]. The high symmetric points in this Brillouin zone include $\Gamma=(0,0)$, $\text{M}=(\pi,0)$, and $\text{K}=(\pi,\pi)$.

2D contour of the $\pi$-flux lattice is plotted with the energy dispersions E_{(2,1),\pi-\text{flux}}^{0}=\pm\big{[}2(\cos^{2}{k_{x}}+\cos^{2}{k_{y}})\big{]}^{1/2}. The upper $E_{2}^{0}$ and lower $E_{1}^{0}$ bands contact each other at four Dirac points $D=(\pm\pi/2,\pm\pi/2$) (white region in the Fig. 2b). So, the upper band $E_{2}^{0}$ is empty state (electron band), and lower band $E_{1}^{0}$ is completely filled (hole band). That is a typical band structure of semimetal (an example of graphene [26]). Four distinct symmetric points in the rectangular Brillouin zone are $\Gamma=(0,0)$, $\text{X}=(\pi,0)$, $\text{Y}=(0,\pi/2)$, and $\text{K}^{\prime}=(\pi,\pi/2)$.

III.2 Eigenvalue spectra in the momentum space

Figure 3 shows the momentum-dependence of the second-order susceptibility of the 0-flux and $\pi$-flux models. The eigenvalue spectrum of $J^{\text{RKKY}}_{\alpha\beta}(\mathbf{q})$ ($\alpha,\beta=$ 1 and 2) matrix is plotted along the high symmetric points of the square zone such as $\Gamma=(0,0)$, $\text{M}=(\pi,0)$ and $\text{K}=(\pi,\pi)$. Because of symmetric properties between the x- and y- directions, we only calculate it at the M$(\pi,0)$ point. That spectrum shows degenerate maximum value at the K = ($\pi$, $\pi$) point. According to the She’s discussions [24], the ordering vector of magnetic interaction should be defined at the maximum point K of eigenvalue plot. That is due to the maximum value with minus sign of coupling constant $-\xi^{2}/N^{2}$ gives the minimum magnetic energy or stable system (see equations 7 and 13) [24]. However, there is the other hidden singular point which does not show in the way we calculate eigenvalue spectrum. RKKY interaction is called the static function as $\omega_{n}=0$ and $\mathbf{q}\to\mathbf{0}$. The $J_{11}(\mathbf{q})$, for example, is:

[TABLE]

The final integral form of $J_{11}(\mathbf{q}\to\mathbf{0})$ is diverged at those integral limits. It is an approximation we compute the RKKY interaction in the 0-flux square lattice (see the Sect. 1 of Appendix C for detail calculation).

The eigenvalue spectrum of $J^{\text{RKKY}}_{\alpha\beta}(\mathbf{q})$ ($\alpha,\beta=$ 1, 2, 3 and 4) matrix of $\pi$-flux model is showed in Fig. 3b. Four distinct eigenvalue curves includes $Eigen1$, $Eigen2$, $Eigen3$ and $Eigen4$ are plotted along $\Gamma=(0,0)$, $\text{X}=(\pi,0)$, $\text{K}^{\prime}=(\pi,\pi/2)$ and $\text{Y}=(0,\pi/2)$. The $Eigen4$ plot shows the highest value at the $\text{X}(\pi,0)$ point. Since the Lindhard function of $\pi$-flux model is interband transition, there is no appearance of zero value at its numerator and denominator. The maximum point defines the magnetic order $\mathbf{Q}=(\pi,0)$ of the system.

III.3 Real-space semi-analytic calculations

From the first-order and second-order perturbations in the momentum space, we perform their Fourier transformations into the real space. The first-order effective energies of 0- and $\pi$-flux models exist only when $\mathbf{q}=\mathbf{0}$ ($\delta_{\mathbf{q},\mathbf{0}}$ terms in the equations (6) and (11)). So, they are calculated exactly by taking integration over the Brillouin zone to provide a unique value depending linearly on the parameter $\xi$. The second-order energies are computed approximately according to the strength of Ising spin pair interactions.

Figure 4 shows the real interacting pairs in 0-flux lattice such as Ising 1–1 (a pair of two Ising spins 1) and 1–2 couples (Figs. 4a and 4b). Both interactions are decaying rapidly within one or two lattice distance $R$. Comparing with the normal RKKY interaction in the square lattice [27], our interacting functions are anisotropic and distinct along the x- and y-directions. For example, with the distance of one unit cell along the x-direction, the $J_{11}^{\text{RKKY}}(R=1)$ value shows a negative sign with the factor $-(\xi/4\pi^{2})^{2}$, that gives antiferromagnetic coupling between two Ising spins. However, along the y-direction, that coupling shows a positive sign that corresponds to FM coupling with three times larger magnitude (inset of Fig. 4a). However, the interaction between different spins are antiferromangetically along both the x- and y-directions (Fig. 4b). The long-tail interaction of Ising spins are considered carefully by taking $\log-\log$ plot in Fig. 4c. Our results are consistent with the power decaying rate $R^{-1.85}$ that is closed to known value of RKKY interaction in two-dimensional lattice ($R^{-2}$) [27].

Because of fast decaying in the real-space interaction, we take the value of the closest distance with considering of different Ising spin coupling is the nearest-neighbor $J_{\text{NN}}=J_{\alpha\beta}^{\text{RKKY}}=317.57(\xi/4\pi^{2})^{2}$ with $\alpha\neq\beta$, and same Ising spins coupling is the next-nearest-neighbor $J_{\text{NNN}}$ along the x- and y-directions with $J_{\alpha\alpha-x}^{\text{RKKY}}=17.57(\xi/4\pi^{2})^{2}$, and $J_{\alpha\alpha-y}^{\text{RKKY}}=-60.3(\xi/4\pi^{2})^{2}$, respectively. The magnitude of the nearest-neighbor coupling is much larger than the next-nearest one. The second-order effective energy for the 0-flux lattice is calculated by semi-analytic method:

[TABLE]

Where $\langle\mu,\nu\rangle$ and $\langle\zeta,\eta\rangle$ are the sum of all Ising spin pairs with distances of $\sqrt{2}/2$ and 1, respectively. With the coexistence of the nearest and next-nearest couplings between spins, the formation of magnetic order becomes frustrasted. That effect is similar to the observation of spin-nematic model in heavy fermion $\mathrm{LiCuVO_{4}}$ compounds [28]. This compound is a typical example of competing effect in spin chain with the exisence of nearest-neighbor FM and next-nearest-neighbor AFM couplings. We need the exact diagonalization method to solve that problem and define the correct magnetic order.

In Fig. 5, we calculate different coupling terms between four Ising spins of the $J_{\alpha\beta}^{\text{RKKY}}$ matrix (with $\alpha,\beta=$1, 2, 3 and 4). Quite similar to the 0-flux case, $\pi$-flux results show two distinct interacting directions. Along the x-direction, the effective interactions is decaying rapidly within four to five lattice distance (Fig. 5a). However, we observe the strong sign-changing oscillation along the y-direction even though the interacting distance between spins is raised in double (Fig. 5b). That phenomenon is completely opposite to the no sign-changing oscillation of RKKY interaction in graphene explained by collapsing of Fermi surface to Dirac points [17]. We believe that the sign-changing oscillation is also determined by the momentum function $M_{\alpha\beta}(\mathbf{k},\mathbf{q})$ outside of the Lindhard function. Interestingly, this model shows the magnitude of the nearest-neighbor are smaller than the next-nearest-neighbor one. For example, the nearest-neighbor coupling of Ising spin 1–2, 2–3 pairs are nearly vanished (Fig. 5c) whereas the magnitude of Ising spin 3–4 are smaller $J_{34}=6(\xi/2\pi^{2})^{2}$ comparing with $J_{11-y}=45(\xi/2\pi^{2})^{2}$. Along the x-direction, interactions of Ising spin 1–1 and 3–3 couples have similar strength and ferromagnetic, but spin 2–2 and 4–4 couple ones are antiferromagnetic (inset of Fig. 5a). Along the y-direction, both of them behave similarly as ferromagnetic interactions with larger value than x-direction ones (inset of Fig. 5b). Because there are of ten couplings in two different directions, the magnetic order of the system is hard to predict by this method.

Figure 5d shows $\log-\log$ plots of $J_{\alpha\beta}^{\text{RKKY}}$ via distance $R$. We have a trouble to find the consistent decaying rate of that model. Along the y-direction, because of sign-changing oscillation, we take the maximum peaks of plot in Fig. 5b to calculate $\log$ value, so they decay with power of $R^{-1.08}$. The long-tail of nearest coupling 3–4 pair is the fastest decaying in the x-direction with $R^{-4.3}$ rate. Other x-direction couplings are decaying rate with power of $R^{-3}$. Our results show different behavior comparing with the long-distance limit of RKKY interaction in the graphene with decaying rate $R^{-3}$, even though we have similar form of interband Lindhard function and Dirac points. That issue may be interpreted by the effect of the momentum function $M_{\alpha\beta}(\mathbf{k},\mathbf{q})$ outside of Lindhard function.

III.4 Energy configurations and magnetic orders

We take advantage of the numerical calculation to search for the ground-state energy for each model. The total energy is calculated by exact diagonalization, and depended on the initial setup Ising spins at one unit cell, magnetic ordering vector $\mathbf{Q}$ and coupling parameter $\xi$. Figures 6a and 6b show different energy configurations of the 0- and $\pi$-flux lattices calculating both semi-analytic and exact diagonalization methods. We observe that the FM order is the lowest trivial energy (ground state) that comprises in both methods for two lattices. Our results show perfect matching between the semi-analytic method – First Order line and exact calculation – FM-q00 open square in 0-flux model, or First Order line and open circle FM-00-1111 in the $\pi$-flux (see Fig. 6a(2) and Fig. 6b(2) or $0\%$ error in Table. 1). It is meant that the non-interacting energy dispersion of the system is stretched out with $1+\xi$ amplitude because each hopping integral is increased linearly with factor $\xi$. If the ferromagnetic order is set up, the first-order perturbing energy is dominated.

We visualize the FM orders of two lattices in Fig. 7a and Fig. 7b. In the 0-flux lattice, there is a single fermion basis that surrounding by four Ising spin-ups (red arrow in the Fig. 7a). Similarly, two different fermions in the $\pi$-flux couple four Ising spins-ups in four directions (dash violet and green boxes in the Fig. 7b).

The total energy of the system is generalized by:

[TABLE]

where $A$ and $B$ are some constant coefficients that are extracted from the effective calculation in the semi-analytic method. $\sum_{\alpha}$ is the sum of all Ising spins in lattices, and $\sum_{\alpha<\beta}$ is the sum of Ising pairs. So, $A\sum_{\alpha}\sigma_{\alpha}^{z}$ and $B\sum_{\alpha<\beta}\sigma_{\alpha}^{z}\sigma_{\beta}^{z}$ terms are exact first and second coefficients, respectively (listed in the Table 1). For exact diagonalization method, they are interpolated from the plot of $E^{\text{total}}-\xi$. Detail calculations of $A$ and $B$ coefficients are found in the section 2 of Appendix C. The $E^{\text{total}}$ includes the linear or quadratic forms of coupling parameter $\xi$.

If our lattices are set up antiferromagnetically, the first-order perturbation energy is vanished ($\sum_{\alpha}\sigma_{\alpha}^{z}=0$). For semi-analytic method, we follow the equation (17) to calculate the effective energy of the 0-flux model. Unlike to unique FM magnetic order, several AFM configurations appear in the 0-flux and $\pi$-flux models.

In the Fig. 6a(1), we have seen that two AFM orders of 0-flux lattice: AFM-q00 (black star symbol) and AFM-q$\pi\pi$ (blue line and circle) behave similarly. However, when those plots are zoomed in the high enough resolution, the energy curve of AFM-q00 order is lower value than AFM-q$\pi\pi$ one (Fig. 6a(3) ). So, the advantage of exact diagonalization method over semi-analytic one is neglecting the singular point, and we have found the correct magnetic ordering vector $\mathbf{Q}=(0,0)$. Figure 6a(2) shows agreement between semi-analytic and exact calculation results with $\xi\leq 0.2$ known as the fundamental property of perturbation theory – weak coupling limit. The percent error of two methods is quite high as $12.7\%$ (Table 1) because of the existence of singular point and approximation in the semi-analytic computation.

That AFM order with the magnetic ordering vector $\mathbf{Q}=(0,0)$ is visualized in the Fig. 7c, and there are two spin-ups along the x-direction and two spin-downs along the y-direction. If we look at the $45^{\circ}$ rotation of the square lattice, that magnetic order is exact Néel state that is comprised with typical RKKY ground state in the square lattice [22]. The effect of spin-up enhances hopping (or weak bond), and the spin-down retards the hopping (or strong bond) of fermion from one site to the other. If the lattice is applied the electric field, the electric current will flow anisotropic, and be higher conductivity in the x-direction than y-direction.

In Fig. 6b, we plot the energy $E^{\text{total}}$ versus coupling parameter $\xi$ curves of the AFM configurations in $\pi$-flux lattice. Because of four different Ising spins at one unit cell, we have investigated more spin configurations than the 0-flux one. Figure 6b(2) zooms in the AFM energy curves calculated from Second order plot by semi-analytic method (red solid line) and AFM$\pi 0\_$1-111 plot by exact diagonalization (open square symbol) with setting up of three spin-ups positioned at Ising 1, 3, and 4, one spin-down for Ising 2, and the magnetic ordering vector $\mathbf{Q}=(\pi,0)$. Two energetic ($E^{\text{total}}-\xi$) plots of the methods match very well with $1.0\%$ error of the different coefficients (Table 1). We believe thi magnetic order is the ground-state energy of AFM configurations in the Fig. 6b(3). That order preserves the $Z_{2}$ symmetry when I reversed all signs of Ising spins. Analogous to the 0-flux lattice, the semi-analytic result is only fitted with exact diagonalization at the limit of $\xi\leq 0.3$.

The AFM order is visualized in Fig. 7d with two distinct fermions (in the dash violet and green boxes). Each fermion is coupled with four surrounding Ising spins. Fermion A (filled black circle) first coupled with two spin-ups in the East and South directions, and two spin-downs in the West and North directions (green box), so the next-nearest fermion A along the x-direction reverses all signs of Ising spins. Fermion B has a similar behavior (open circle in the violet box). However, along the y-direction, the the sign of Ising spin are unchanged. That follows the magnetic ordering order $\mathbf{Q}=(\pi,0)$, or it is like a spin wave along the x-direction. Therefore, the AFM state of $\pi$-flux is different from 0-flux one.

IV Conclusions

We consider the model of Ising spins on the links of the square lattices, coupled to the fermion charge fluctuations for the 0-flux and $\pi$-flux cases. The unit cell of the $\pi$-flux lattice doubles of 0-flux one. The Brillouin zone changed from symmetrical square of 0-flux to asymmetrical rectangular of $\pi$-flux lattice. At the half-filling, the tight binding parts of Hamiltonian provide continuous metallic band and semimetallic Dirac bands for the 0- and $\pi$-flux lattices, respectively. We observe the Fermi surface with the nesting vectors, and Fermi points in Dirac bands for two cases.

Use the semi-analytic perturbation theory and exact diagonalization method, we achieve different magnetic orders and ground state energies based on spin configurations. With the first-order perturbation, both kinds of models show the trivial FM ground state, and ground state energy that is linear dependent of coupling parameter $\xi$. That is due to spin-up enhances the hopping magnitude with the factor of $1+\xi$. We can use our models or depleted Anderson’s model [29] to construct the FM magnetic order in the square lattice.

When the Ising spin configuration is set up antiferromagnetically, the second-order perturbation is considered. The ground state energy now is following the quadratic function of coupling parameter $\xi$. The strong agreement between the semi-analytic and exact diagonalization methods only exists with the coupling parameter $\xi\ll 1$. We finally construct lowest energy AFM configuration with the Néel state of the $45^{\circ}$ rotation $\mathbf{Q}=(0,0)$, and spin-wave AFM oscillating along the x-direction with the ordering vector $\mathbf{Q}=(\pi,0)$ for the 0-flux and $\pi$-flux models, respectively.

Our achieving results contribute to the rich variety of phenomena in the Ising-nematic square lattices [7; 8; 9]. Particularly, the FM state exists in the 2D experimental observation can be explained by our models [2]. Moreover, we consider a case of coupling $\xi$ between spin and fermion in the weak limit. Our future work will continue with inserting the other interaction like Hubbard term [6; 30], and tuning the filling factor [9] to those lattices.

V Acknowledgement

Numerical calculations have been carried out on machines hosted at the Mississippi Center for Supercomputing Research, University of Mississippi.

Appendix A Explicit real-space Hamiltonian

A.1 Zero-flux model

The Hamiltonian of 0-flux lattice is written explicitly in real space with the lattice translation vector $\mathbf{R}=m\mathbf{a}_{1}+n\mathbf{a}_{2}$ ($m$, $n$ are integers):

[TABLE]

[TABLE]

A.2 Discrete Fourier transformation identities

The discrete Fourier transformation for fermion $c_{\mathbf{R}}^{\dagger}$ and Ising spin $\sigma_{\mathbf{R},\alpha}^{z}$ from the real space to momentum space are:

[TABLE]

A.3 $\pi$-flux model

The general Hamiltonian of the $\pi$-flux lattice is the form of vector $\mathbf{R}$ in the real space:

[TABLE]

We separate into two parts: tight-binding Hamiltonian $\hat{H}_{0}$ and interacting Hamiltonian $\hat{H}_{1}$:

[TABLE]

[TABLE]

The real-space Hamiltonians are then transformed into to the momentum space using the discrete Fourier transformations above. After that, we use the Fermi-sea ground state to calculate the first- and second-order effective interactions and energies. (Detail calculation is described in my master thesis [25], or in the reference [22]).

Appendix B Function definitions

B.1 Second-order RKKY interaction in the 0-flux lattice

All terms of 2$\times$2 $J_{\alpha\beta}(\mathbf{k},\mathbf{q})$ matrix are listed below:

[TABLE]

B.2 Second-order RKKY interaction in the $\pi$-flux lattice

There are some useful functions:

[TABLE]

All terms of 4$\times$4 $M_{\alpha\beta}(\mathbf{k},\mathbf{q})$ matrix are listed below:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix C Example calculation

C.1 Derivation of static function $J_{11}(\mathbf{q},\omega_{n}=0)$

If we consider the $\mathbf{q}$ is very small, we can expand the Fermi-Dirac function $n_{F}(E_{\mathbf{k}+\mathbf{q}})$ and energy dispersion relation $E(\mathbf{k}+\mathbf{q})$ of the Lindhard function in terms of $\mathbf{q}$. So the momentum-dependence of two spin 1 interaction is written:

[TABLE]

C.2 Calculation of second-order effective energy using semi-analytical methods

For 0-flux lattice, the real-space interactions calculating by semi-analytic method include the nearest-neighbor and next-nearest-neighbor pairs. When the lattice size of the system is $L$, the total number of Ising spins living at the bonds is $N_{\text{spin}}=2L^{2}$. We have $N$ number of Ising 1 and $N$ number of Ising spin 2. There are $2N$ numbers of antiferromagnetic Ising 1–2 and 2–1 pairs, and $N/2$ number of Ising 1–1 along the x- and y- directions.

[TABLE]

So, the second-order coefficient in the 0-flux lattice is $B\sum_{\alpha<\beta}\sigma_{\alpha}^{z}\sigma_{\beta}^{z}=-0.89$. The coefficients of the $\pi$-flux are calculated in the similar way.

For exact diagonalization, when we obtain the total energy data, we plot them via the coupling parameter $\xi$ for different Ising spin configurations. Using power fitting, we find the first- and second-order coefficients.