$\eta$ -pairing on bipartite and non-bipartite lattices

Yutaro Misu¹, Shun Tamura², Yukio Tanaka² and Shintaro Hoshino¹ ¹Department of Physics, Saitama University, Saitama 338-8570, Japan
²Department of Applied Physics, Nagoya University, Nagoya 464-8603, Japan

Abstract

The $\eta$ -pairing is a type of Cooper pairing state in which the phase of the superconducting order parameter is aligned in a staggered manner, in contrast to the usual BCS superconductors with a spatially uniform phase. In this study, we search for a characteristic $\eta$ -pairing state in a triangular lattice where a simple staggered alignment of the phase is not possible. As an example, we consider the attractive Hubbard model on both the square and triangular lattices under strong external Zeeman field. Using the mean-field approximation, we have identified several $\eta$ -pairing states. Additionally, we have examined the electromagnetic stability of the pairing state by calculating the Meissner kernel. Odd-frequency pairing plays a crucial role in achieving diamagnetic response if the electrons experience a staggered superconducting phase during the propagation of current.

I Introduction

The diversity of superconducting phenomena has been attracting continued attention. The superconducting state of matter is characterized by the properties of Cooper pairs, which can be classified based on their space-time and spin structures. With regard to their space structure, Cooper pairs are typically classified as $s$ -wave, $p$ -wave, or $d$ -wave pairs depending on their relative coordinate structure. As for their center-of-mass coordinate, while it is usually assumed to be zero in most superconductors, it is possible to consider the existence of a finite center-of-mass momentum. One example of this is the Flude-Ferrell-Larkin-Ovchinnikov (FFLO) state [1, 2], in which the Cooper pair has a small but finite center-of-mass momentum under the influence of a magnetic field. More generally, the magnitude of the center-of-mass momentum can be larger and of the order of the reciprocal lattice vector $\sim\pi/a$ , where $a$ is a lattice constant. This type of pairing state is known as $\eta$ -pairing, a concept first proposed by C. N. Yang, which forms a staggered alignment of the superconducting phase on a bipartite lattice [3]. The spatially modulating order parameter is known also as the pair density wave, and has been discussed in relation to cuprate superconductors [4].

The actual realization of the $\eta$ -pairing has been proposed for the correlated electron systems such as the attractive Hubbard (AH) model with the magnetic field [5], the single- and two-channel Kondo lattices [6, 7], the Penson-Kolb model [8], and also the non-equilibrium situation [9, 10, 11, 12, 13, 14]. Since the phase of the superconducting order parameter can be regarded as the XY spin, the $\eta$ -pairing is analogous to an antiferromagnetic state of the XY spin model. Hence, the $\eta$ -pairing state should be strongly dependent on the underlying lattice structure and we naively expect a variety of the $\eta$ -pairing state if we consider the geometrically frustrated lattice such as the triangular lattice since the simple staggered state cannot be realized.

In this paper, we deal with the AH model on the non-bipartite lattice in order to search for possible new superconducting states depending on the feature of the non-bipartite lattice structure in equilibrium. Already in the normal state without superconductivity, it has been pointed out that the non-bipartite lattice generates a non-trivial state of matter. For example in the Kondo lattice, a partial-Kondo-screening, which has a coexisting feature of Kondo spin-singlet and antiferromagnetism, is realized [15]. Also in the AH model at half-filling, charge-density-wave (CDW) is suppressed due to the frustration effect [16]. The $\eta$ -pairing that appears in a photodoped Hubbard model on the triangular lattice has been studied recently [14]. In the equilibrium situation, the properties of the AH model have been studied on bipartite lattices [5], but the model on a non-bipartite lattice has not been explored.

As shown in the rest of this paper, there are several types of $\eta$ -pairings on the triangular lattice of the AH model under the Zeeman field. One of the $\eta$ -pairing states is regarded as a 120^∘-Néel state. Since the relative phase between the nearest neighbor sites is neither parallel nor anti-parallel, the inter-atomic Josephson current is spontaneously generated. This state can also be regarded as a staggered flux state, where the flux is created by the atomic-scale superconducting loop current. While the staggered flux state has been studied so far [17, 18, 19, 20, 21, 22, 23], the staggered flux in this paper is induced by the Josephson effect associated with superconductivity and has a different origin.

For the analysis of the AH model, we employ the mean-field approximation in this paper. It has been suggested that a simple $\eta$ -pairing shows a paramagnetic Meissner state [24]. Hence it is necessary to investigate the electromagnetic stability of the solution for superconductivity. We evaluate the Meissner kernel whose sign corresponds to the diamagnetic (minus) or paramagnetic (plus) response of the whole system, where the physically stable state should show diamagnetism. We confirm that if the mean-field $\eta$ -pairing state has the lowest energy compared to the other ordered states, the calculation of the Meissner kernel shows the diamagnetic response. It is also notable that the odd-frequency pairing amplitude, which has an odd functional form with respect to the frequency [25, 26, 27, 6, 28, 29, 30], can contribute to the diamagnetism in the $\eta$ -pairing state. This is in contrast to the usual superconductivity with the uniform phase where the conventional even-frequency pairing contributes to the diamagnetism. It has been shown that the odd-frequency pairing induced at the edge, interface or junctions [31, 32, 33, 34, 35, 36] shows a paramagnetic response [37, 38, 39, 40, 41]. In this paper, by contrast, we consider the odd-frequency pairing realized in bulk, which shows a qualitatively different behavior.

This paper is organized as follows. We explain the model and method for the AH model in Sec. II, and the Meissner kernel in Sec. III. The numerical results for the AH model are shown in Sec. IV, and we summarize the paper in Sec. V.

II Attractive Hubbard model

II.1 Hamiltonian

We consider the Hamiltonian of the AH model with magnetic field $\bm{h}$ which induce Zeeman effect only (Zeeman field) :

	$\displaystyle{\cal{H}}=-t\displaystyle\sum\limits_{\langle i,j\rangle\sigma}c_{i\sigma}^{\dagger}c_{j\sigma}+\mathrm{H.c.}$	$\displaystyle+U\displaystyle\sum\limits_{i}n_{i\uparrow}n_{i\downarrow}$
		$\displaystyle-\mu\displaystyle\sum\limits_{i}n_{i}-\bm{h}\cdot\displaystyle\sum\limits_{i}\bm{s}_{i},$		(1)

where $c_{i\sigma}^{\dagger}$ and $c_{i\sigma}$ are the creation and annihilation operators of the $i$ -th site with spin $\sigma$ , respectively. The symbol $\langle i,j\rangle$ represents a pair of the nearest-neighbor sites. Here, the parameter $t$ is the nearest-neighbor single-electron hopping integral. $U$ ( $=-|U|$ ) is the onsite attractive interaction. The spin operator is defined as $\bm{s}_{i}=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}c_{i\sigma}^{\dagger}\bm{\tau}_{\sigma\sigma^{\prime}}c_{i\sigma^{\prime}}$ , where $\bm{\tau}$ is the Pauli matrix, and the number operator of electrons is denoted as $n_{i}=n_{i\uparrow}+n_{i\downarrow}=\sum_{\sigma}c_{i\sigma}^{\dagger}c_{i\sigma}$ . The electron concentration is controlled by adjusting the chemical potential $\mu$ .

Refer to caption — Figure 1: Sketches of the phase diagrams for the repulsive Hubbard model [44] (left panel) and AH model (right panel). $n_{c}$ is the electron concentration and $m$ is the magnetization. When the interaction $|U|$ is large, the ground state in the repulsive Hubbard model is ferromagnet (FM), while the ground state in the AH model is $\eta$ -pairing.

The AH model has been successfully used to elucidate several important and fundamental issues in superconductors [42]. The model on a bipartite lattice at half filling is theoretically mapped onto the repulsive Hubbard model by the following partial particle-hole transformation [43]

\displaystyle c_{i\uparrow}^{\dagger}\rightarrow~{}c_{i\uparrow}^{\dagger},~{}c_{i\downarrow}^{\dagger}\rightarrow~{}c_{i\downarrow}\mathrm{e}^{{\rm{i}}\bm{Q}\cdot\bm{R}_{i}}.

(2)

The reciprocal vector $\bm{Q}$ satisfies the condition $\mathrm{e}^{{\rm{i}}\bm{Q}\cdot\bm{R}_{i}}=(-1)^{i}$ that takes $\pm 1$ depending on $\bm{R}_{i}$ belonging to A or B sublattice on the bipartite lattice. Then, the $\eta$ -pairing appears in the region that corresponds to a ferromagnet with transverse magnetization in the repulsive model [5]. In a mean-field theory, the phase diagram for the repulsive Hubbard model without the magnetic field is shown in the left panel of Fig. 1 [44]. From this figure, we find that the ferromagnet is located in the regime where the repulsive interaction $U>0$ is large and the electron concentration is not half-filled. Hence, the $\eta$ -pairing phase is located in the regime where the attractive interaction $U<0$ is large and the magnetization is finite. The phase diagram of the AH model at half filling is shown in the right panel of Fig. 1. In principle, an attractive interaction large enough to realize $\eta$ -pairing could be realized in artificial cold atom systems [45].

The Cooper pair is formed by the two electrons with $(\bm{k}\uparrow,~{}-\bm{k}+\bm{q}\downarrow)$ where $\bm{q}$ is the center-of-mass momentum. The FFLO state and the $\eta$ -pairing are distinguished by the magnitude of $|\bm{q}|$ . In $\eta$ -pairing, the center-of-mass momentum of the Cooper pair is the order of the reciprocal lattice vector, while the momentum of the FFLO state is much smaller and the spatial modulation is slowly-varying compared to the atomic scale. Although the large center-of-mass momentum is usually not energetically favorable, a strong attractive interaction can make it stable.

II.2 Mean-field theory

By applying the mean-field approximation, we obtain the mean-field Hamiltonian

	$\displaystyle{\cal{H}}^{\rm{MF}}$	$\displaystyle=-t\displaystyle\sum\limits_{\langle i,j\rangle\sigma}c_{i\sigma}^{\dagger}c_{j\sigma}+\mathrm{H.c.}-\mu\displaystyle\sum\limits_{i}n_{i}-\bm{h}\cdot\displaystyle\sum\limits_{i}\bm{s}_{i}$
		$\displaystyle~{}~{}-\displaystyle\sum\limits_{i}\left(v_{i}n_{i}+\bm{H}_{i}\cdot\bm{s}_{i}-\Delta_{i}c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}-\Delta^{*}_{i}c_{i\downarrow}c_{i\uparrow}\right).$		(3)

The order parameters are given by the self-consistent equations

	$\displaystyle v_{i}\equiv\dfrac{\|U\|}{2}\langle{n_{i}}\rangle,$		(4)
	$\displaystyle\Delta_{i}\equiv-\|U\|\langle c_{i\downarrow}c_{i\uparrow}\rangle,$		(5)
	$\displaystyle\bm{m}_{i}=\dfrac{1}{2}\displaystyle\sum\limits_{\sigma\sigma^{\prime}}\langle{c_{i\sigma}^{\dagger}\bm{\tau}_{\sigma\sigma^{\prime}}c_{i\sigma^{\prime}}}\rangle,~{}~{}\bm{H}_{i}=~{}-2\|U\|\bm{m}_{i},$		(6)

where $\displaystyle\langle A\rangle=\mathrm{Tr}\left[A\mathrm{e}^{-{\cal{H}}^{{\rm{MF}}}/T}\right]/\mathrm{Tr}\left[\mathrm{e}^{-{\cal{H}}^{{\rm{MF}}}/T}\right]$ is a quantum statistical average with the mean-field Hamiltonian and $T$ is temperature. $\Delta_{i}$ is the order parameter for $s$ -wave singlet superconductivity (pair potential). The phase $\theta_{i}\in[0,2\pi)$ of the pair potential $\Delta_{i}=|\Delta_{i}|\mathrm{e}^{\mathrm{i}\theta_{i}}$ is dependent on the site index and will be represented by the arrow in a two-dimensional space. The mean-fields for the charge and spin are given by $v_{i}$ and $\bm{H}_{i}$ , respectively, at each site. The derivation of the self-consistent equations is summarized in Appendix A. We will consider the AH model both on the two-dimensional square and triangular lattices.

III Meissner kernel for a general tight-binding lattice

III.1 Definition

As we explained in Sec. I, it is necessary to calculate the Meissner kernel to determine whether the mean-field solution for $\eta$ -pairing is electromagnetically stable. In the tight-binding model, the electromagnetic field appears as Peierls phase:

\displaystyle{\cal{H}}_{\rm{kin}}=-t\displaystyle\sum\limits_{\langle i,j\rangle\sigma}\mathrm{e}^{{\rm{i}}A_{ij}}c_{i\sigma}^{\dagger}c_{j\sigma}+\rm{H.c.}.

(7)

The Meissner effect is examined by the weak external orbital magnetic field applied perpendicular to the plane, while the $\eta$ -pairing is stabilized only under a strong Zeeman field. In order to make these compatible, we apply the Zeeman field parallel to the plane $\bm{h}=(h,0,0)$ , which does not create the orbital motion of the tight-binding electrons. Thus, the weak magnetic field that triggers the Meissner effect is applied perpendicular to the plane in addition to the in-plane magnetic field. While the out-of-plane Zeeman effect is also induced by the weak additional field, it is neglected since the dominant Zeeman field already exists by the strong in-plane magnetic field.

Let us formulate the Meissner response kernel on a general tight-binding model. We apply the formulation in Refs. [46, 47, 48] to the present case with sublattice degrees of freedom. The current density operator between two sites is defined as

	$\displaystyle\bm{j}_{ij}$	$\displaystyle=\dfrac{\partial{\cal{H}}_{\rm{kin}}}{\partial A_{ij}}\hat{\bm{\delta}}_{ij}$
		$\displaystyle=-\mathrm{i}t\displaystyle\sum\limits_{\sigma}\left(c_{i\sigma}^{\dagger}c_{j\sigma}\mathrm{e}^{\mathrm{i}A_{ij}}-c_{j\sigma}^{\dagger}c_{i\sigma}\mathrm{e}^{-\mathrm{i}A_{ij}}\right)\hat{\bm{\delta}}_{ij},$		(8)

where $\bm{\delta}_{ij}=\bm{R}_{i}-\bm{R}_{j}$ is the inter-site lattice vector between $i$ -th and $j$ -th sites, and hat ( $\hat{\ }$ ) symbol means a unit vector. In the linear response theory, the current operator which appears as a response to the static magnetic field in equilibrium is written as

$\displaystyle\bm{j}_{ij}$	$\displaystyle\simeq-{\rm{i}}t\displaystyle\sum\limits_{\sigma}(c_{i\sigma}^{\dagger}c_{j\sigma}-c_{j\sigma}^{\dagger}c_{i\sigma})\hat{\bm{\delta}}_{ij}$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}+t\displaystyle\sum\limits_{\sigma}(c_{i\sigma}^{\dagger}c_{j\sigma}+c_{j\sigma}^{\dagger}c_{i\sigma})\hat{\bm{\delta}}_{ij}A_{ij}$
	$\displaystyle\equiv\bm{j}^{\rm{para}}_{ij}+\bm{j}^{\rm{dia}}_{ij}.$	(9)

The first term is called the paramagnetic term and the second term is diamagnetic. The Fourier-transformed paramagnetic and diamagnetic current density operators are written as $\bm{j}^{\rm para}(\bm{q})$ and $\bm{j}^{\rm dia}(\bm{q})$ . The linear response kernel is then defined by $\langle{j_{\nu}(\bm{q})}\rangle=\sum_{\mu}K_{\nu\mu}(\bm{q})A_{\mu}(\bm{q})$ , where $\nu,\mu=x,y$ is the direction. We evaluate the kernel $K_{\nu\mu}(\bm{q}\rightarrow\bm{0})\equiv K_{\nu\mu}$ when investigating the stability of superconductivity. This is called the Meissner kernel, which is proportional to the superfluid density.

The Meissner kernel is separated into paramagnetic and diamagnetic terms as $K_{\nu\mu}=\left(K_{\rm{para}}\right)_{\nu\mu}+\left(K_{\rm{dia}}\right)_{\nu\mu}$ . The paramagnetic kernel is given by

\displaystyle\left(K_{{\rm{para}}}\right)_{\nu\mu}=\dfrac{1}{N}\displaystyle\int_{0}^{1/T}d\tau\langle j_{\nu}^{{\rm{para}}}(\bm{q}=0,\tau)j_{\mu}^{{\rm{para}}}(\bm{q}=0)\rangle,

(10)

where $N=\sum_{i}1$ is the number of sites. The Heisenberg representation with the imaginary time $\tau$ is defined as $A(\tau)=\mathrm{e}^{{\cal{H}}\tau}A\mathrm{e}^{-{\cal{H}}\tau}$ . The form of the diamagnetic kernel is obvious from Eq. (9).

We note that if the sign of the Meissner kernel $K$ is negative, the superconducting state is electromagnetically stable and is also called a diamagnetic Meissner state, which expels magnetic flux. On the other hand, if the sign is positive, the superconducting state is called the paramagnetic Meissner state, which attracts magnetic flux. For a stable thermodynamic superconducting state, the negative value of $K$ is required.

III.2 Method of evaluation

The actual evaluation of the kernels is performed based on the wave-vector representation. Here, the physical quantities are described by the operator $c_{\bm{k}\sigma}^{\alpha}$ where $\alpha$ distinguishes the sublattice. Note that the Brillouin zone is folded by $\sum_{\alpha}1$ times. The diamagnetic kernel is rewritten as

\displaystyle\left(K_{{\rm{dia}}}\right)_{\nu\mu}=\dfrac{1}{N}\displaystyle\sum\limits_{\alpha,\beta}\displaystyle\sum\limits_{\bm{k}\sigma}\left(m^{-1}_{\bm{k}\alpha\beta}\right)_{\nu\mu}\langle{c_{\bm{k}\sigma}^{\alpha\dagger}c_{\bm{k}\sigma}^{\beta}}\rangle.

(11)

The inverse mass tensor $m^{-1}_{\bm{k}\alpha\beta}$ , which reflects the characteristics of the lattice shape, are given by

\displaystyle\left(m^{-1}_{\bm{k}\alpha\beta}\right)_{\nu\mu}\equiv t\displaystyle\sum\limits_{\langle i_{\alpha},j_{\beta}\rangle}\left({\hat{\bm{\delta}}}_{i_{\alpha}j_{\beta}}\right)_{\nu}\left({\hat{\bm{\delta}}}_{i_{\alpha}j_{\beta}}\right)_{\mu}\mathrm{e}^{-i\bm{k}\cdot\bm{R}_{{i}_{\alpha}{j}_{\beta}}},

(12)

where $i_{\alpha}$ is the $i$ -th unit cell with sublattice $\alpha$ . The symbol $\langle i_{\alpha},j_{\beta}\rangle$ represents a pair of the nearest-neighbor sites and $\bm{R}_{{i}_{\alpha}{j}_{\beta}}$ is the vector between the unit lattice with the $i$ -th sublattice $\alpha$ and the unit lattice with the $j$ -th sublattice $\beta$ .

The paramagnetic term has the form of a current-current correlation function. We can calculate this term by using the Green’s function matrix

\displaystyle\check{\cal{G}}_{\bm{k}}(\tau)\equiv-\langle T_{\tau}\bm{\psi}_{\bm{k}}(\tau)\bm{\psi}_{\bm{k}}^{\dagger}\rangle

(13)

where $\bm{\psi}_{\bm{k}}=(c_{\bm{k}\uparrow}^{\alpha},c_{-\bm{k}\downarrow}^{\alpha\dagger},\cdots)^{T}$ is the Nambu-spinor. $T_{\tau}$ is time-ordering operator regrading $\tau$ . Each component of the Green’s function matrix is given by the diagonal and off-diagonal Green’s functions:

	$\displaystyle G_{\sigma\sigma^{\prime}}^{\alpha\beta}(\bm{k},\tau)\equiv-\langle{T_{\tau}c_{\bm{k}\sigma}^{\alpha}(\tau)c_{\bm{k}\sigma^{\prime}}^{\beta\dagger}}\rangle,$		(14)
	$\displaystyle\bar{G}_{\sigma\sigma^{\prime}}^{\alpha\beta}(\bm{k},\tau)\equiv-\langle{T_{\tau}c_{\bm{k}\sigma}^{\alpha\dagger}(\tau)c_{\bm{k}^{\prime}\sigma^{\prime}}^{\beta}}\rangle,$		(15)
	$\displaystyle F_{\sigma\sigma^{\prime}}^{\alpha\beta}(\bm{k},\tau)\equiv-\langle{T_{\tau}c_{\bm{k}\sigma}^{\alpha}(\tau)c_{-\bm{k}\sigma^{\prime}}^{\beta}}\rangle,$		(16)
	$\displaystyle F_{\sigma\sigma^{\prime}}^{\alpha\beta\dagger}(\bm{k},\tau)\equiv-\langle{T_{\tau}c_{-\bm{k}\sigma}^{\alpha\dagger}(\tau)c_{\bm{k}\sigma^{\prime}}^{\beta\dagger}}\rangle.$		(17)

The anomalous part of Green’s function [Eq. (16)] is also called the pair amplitude. The paramagnetic kernel in Eq. (10) can be divided into the normal ( $G$ ) and anomalous ( $F$ ) Green’s function contributions as

$\displaystyle\left(K_{\rm{para}}\right)_{\nu\mu}$	$\displaystyle=-\dfrac{1}{N}\displaystyle\sum\int_{0}^{1/T}\mathrm{d}\tau\left(\bm{v}_{\bm{k}\alpha\beta}\right)_{\nu}\cdot\left(\bm{v}_{\bm{k}\alpha^{\prime}\beta^{\prime}}\right)_{\mu}\times\left(\bar{G}^{\alpha\beta^{\prime}}_{\sigma\sigma^{\prime}}(\bm{k},\tau)G^{\alpha^{\prime}\beta}_{\sigma\sigma^{\prime}}(\bm{k},\tau)+\bar{G}^{\alpha\beta^{\prime}}_{\sigma\sigma^{\prime}}(-\bm{k},\tau)G^{\alpha^{\prime}\beta}_{\sigma\sigma^{\prime}}(-\bm{k},\tau)\right)$
	$\displaystyle\hskip 8.53581pt-\dfrac{1}{N}\displaystyle\sum\int_{0}^{1/T}\mathrm{d}\tau\left(\bm{v}_{\bm{k}\alpha\beta}\right)_{\nu}\cdot\left(\bm{v}_{-\bm{k}\alpha^{\prime}\beta^{\prime}}\right)_{\mu}\times\left(F^{\beta\alpha\dagger}_{\sigma^{\prime}\sigma}(\bm{k},-\tau)F^{\alpha^{\prime}\beta^{\prime}}_{\sigma,\sigma^{\prime}}(\bm{k},\tau)+F^{\beta\alpha\dagger}_{\sigma^{\prime}\sigma}(-\bm{k},-\tau)F^{\alpha^{\prime}\beta^{\prime}}_{\sigma,\sigma^{\prime}}(-\bm{k},\tau)\right)$
	$\displaystyle\equiv K^{G}_{\rm{para}}+K^{F}_{\rm{para}}.$	(18)

The summation $\sum$ is performed over the indices which appears only in the right-hand side. The velocity vector $\bm{v}_{\bm{k}\alpha\beta}$ is defined by

\displaystyle\left(\bm{v}_{\bm{k}\alpha\beta}\right)_{\nu}\equiv t\displaystyle\sum\limits_{\langle i_{\alpha},j_{\beta}\rangle}\left(\hat{\bm{\delta}}_{i_{\alpha}j_{\beta}}\right)_{\nu}\mathrm{e}^{-\mathrm{i}\bm{k}\cdot\bm{R}_{i_{\alpha}j_{\beta}}}.

(19)

In order to perform the integral with respect to $\tau$ in Eq. (18), we define the Fourier-transformed Green’s function as

\displaystyle g_{\bm{k}}({\rm{i}}\omega_{n})\equiv\displaystyle\int_{0}^{1/T}d\tau g_{\bm{k}}(\tau)\mathrm{e}^{{\rm{i}}\omega_{n}\tau},

(20)

where $g_{\bm{k}}$ represents one of Eqs. (14)-(17) and $\omega_{n}=(2n+1)\pi T$ is fermionic Mastubara frequency. Moreover, the Fourier-transformed Green’s function matrix is given by using the matrix representation of mean-field Hamiltonian Eq. (3) as

\displaystyle\check{\cal{G}}_{\bm{k}}({\rm{i}}\omega_{n})

\displaystyle=\left[{\rm{i}}\omega_{n}\check{1}-\check{\cal{H}}^{{\rm{MF}}}_{\bm{k}}\right]^{-1}=\check{U}_{\bm{k}}\left[{\rm{i}}\omega_{n}\check{1}-\check{\Lambda}_{\bm{k}}\right]^{-1}\check{U}^{\dagger}_{\bm{k}},

(21)

where $\check{\Lambda}_{\bm{k}}$ and $\check{U}_{\bm{k}}$ are, respectively, a diagonal eigenvalue matrix and a unitary matrix satisfying $\check{U}^{\dagger}\check{\cal{H}}_{\bm{k}}^{{\rm{MF}}}\check{U}=\check{\Lambda}_{\bm{k}}=\mathrm{diag}(\lambda_{\bm{k}1},\lambda_{\bm{k}2},\ldots)$ . From Eq. (21), $K_{\rm{para}}$ can be calculated as

\displaystyle\left(K_{\rm{para}}\right)_{\nu\mu}

\displaystyle=-\dfrac{1}{N}\sum\left[\left(\bm{v}_{\bm{k}\alpha\beta}\right)_{\nu}\cdot\left(\bm{v}_{\bm{k}\alpha^{\prime}\beta^{\prime}}\right)_{\mu}{\cal{U}}_{\bm{k}p}^{\beta^{\prime}\sigma^{\prime},\alpha\sigma}{\cal{U}}_{\bm{k}p^{\prime}}^{\alpha^{\prime}\sigma,\beta\sigma^{\prime}}+\left(\bm{v}_{\bm{k}\alpha\beta}\right)_{\nu}\cdot\left(\bm{v}_{-\bm{k}\alpha^{\prime}\beta^{\prime}}\right)_{\mu}{\cal{U}}_{\bm{k}p}^{\beta\sigma^{\prime},\alpha\sigma}{\cal{U}}_{\bm{k}p^{\prime}}^{\alpha^{\prime}\sigma,\beta^{\prime}\sigma^{\prime}}\right]\dfrac{f\left(\lambda_{\bm{k}p}\right)-f\left(\lambda_{\bm{k}p^{\prime}}\right)}{\lambda_{\bm{k}p}-\lambda_{\bm{k}p^{\prime}}}+\rm{c.c.}

(22)

where $f(\lambda_{\bm{k}p})=\frac{1}{\mathrm{e}^{\lambda_{\bm{k}p}/T}+1}$ is the Fermi-Dirac distribution function and we have defined the coefficient ${\cal{U}}_{\bm{k}p}^{\alpha\sigma,\beta\sigma^{\prime}}\equiv\left[\check{U}_{\bm{k}}\right]_{\alpha\sigma,p}\left[\check{U}^{\dagger}_{\bm{k}}\right]_{p,\beta\sigma^{\prime}}$ .

The anomalous part of Eq. (18) $K_{\rm{para}}^{F}$ is further decomposed into the contributions $K^{\rm{EFP}}$ and $K^{\rm{OFP}}$ from the even-frequency pair (EFP) and odd-frequency pair (OFP) amplitudes defined by

	$\displaystyle F^{\rm{EFP}}(\bm{k},\mathrm{i}\omega_{n})\equiv\dfrac{F(\bm{k},\mathrm{i}\omega_{n})+F(\bm{k},-\mathrm{i}\omega_{n})}{2},$		(23)
	$\displaystyle F^{\rm{OFP}}(\bm{k},\mathrm{i}\omega_{n})\equiv\dfrac{F(\bm{k},\mathrm{i}\omega_{n})-F(\bm{k},-\mathrm{i}\omega_{n})}{2}.$		(24)

Then, we obtain $K^{\rm{EFP}}$ and $K^{\rm{OFP}}$ by using Eqs. (23) and (24) as

$\displaystyle K^{\rm{EFP,OFP}}_{\nu\mu}$	$\displaystyle=-\dfrac{1}{2N}\displaystyle\sum\limits_{\bm{k}}\displaystyle\sum\limits_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\left(\bm{v}_{\bm{k}\alpha\beta}\right)_{\nu}\cdot\left(\bm{v}_{-\bm{k}\alpha^{\prime}\beta^{\prime}}\right)_{\mu}$
	$\displaystyle\times\displaystyle\sum\limits_{\sigma\sigma^{\prime}}\displaystyle\sum\limits_{pp^{\prime}}{\cal{U}}_{\bm{k}p}^{\beta\sigma^{\prime},\alpha\sigma}{\cal{U}}_{\bm{k}p^{\prime}}^{\alpha^{\prime}\sigma,\beta\sigma^{\prime}}$
	$\displaystyle\times\left[\dfrac{f\left(\lambda_{\bm{k}p}\right)-f\left(\lambda_{\bm{k}p^{\prime}}\right)}{\lambda_{\bm{k}p}-\lambda_{\bm{k}p^{\prime}}}\mp\dfrac{f\left(\lambda_{\bm{k}p}\right)-f\left(-\lambda_{\bm{k}p^{\prime}}\right)}{\lambda_{\bm{k}p}+\lambda_{\bm{k}p^{\prime}}}\right]$
	$\displaystyle+\rm{c.c.},$	(25)

where the minus ( $-$ ) sign in the square bracket is taken for EFP contribution and the plus ( $+$ ) for OFP pairing. These quantities are numerically calculated as shown in the next section. Note that the cross term of the EFP and OFP terms of Green’s functions vanishes after the summation with respect to the Matsubara frequency.

III.3 Paramagnetic Meissner response of a simple $\eta$ -pairing state

Before we show the results of the AH model, let us show that a simple $\eta$ -pairing state leads to the paramagnetic response which would not arise from thermodynamically stable states [24, 49]. We consider the simple bipartite lattice with staggered ordering vector $\bm{Q}$ . The anomalous contribution to the Meissner kernel may be written as [49]

\displaystyle K_{{\rm para},xx}^{F}

\displaystyle=-T\sum_{n\bm{k}\bm{k}^{\prime}\sigma\sigma^{\prime}}v_{\bm{k}}^{x}v_{\bm{k}^{\prime}}^{x}F^{*}_{\sigma^{\prime}\sigma}(\bm{k}^{\prime},\bm{k},\mathrm{i}\omega_{n})F_{\sigma\sigma^{\prime}}(\bm{k},\bm{k}^{\prime},\mathrm{i}\omega_{n}).

(26)

This contribution must be negative (diamagnetic response) in order to dominate over the paramagnetic contribution. For a purely $\eta$ -pairing state, we assume the relation $F_{\sigma\sigma^{\prime}}(\bm{k},\bm{k}^{\prime})=F_{\sigma\sigma^{\prime}}(\bm{k})\delta_{\bm{k}^{\prime},-\bm{k}-\bm{Q}}$ , and obtain

\displaystyle K_{{\rm para},xx}^{F}

\displaystyle=-T\sum_{n\bm{k}\sigma\sigma^{\prime}}(v_{\bm{k}}^{x})^{2}F^{*}_{\sigma^{\prime}\sigma}(\bm{k},\mathrm{i}\omega_{n})F_{\sigma\sigma^{\prime}}(\bm{k},\mathrm{i}\omega_{n}),

(27)

where we have used $v^{x}_{-\bm{k}-\bm{Q}}=v^{x}_{\bm{k}}$ valid for square lattice, which is in contrast to the relation $v^{x}_{-\bm{k}}=-v^{x}_{\bm{k}}$ for the uniform pairing with additional minus sign [24]. We separate the spin-singlet and triplet parts as $F_{\sigma\sigma^{\prime}}=F_{s}\mathrm{i}\tau^{y}_{\sigma\sigma^{\prime}}+\bm{F}_{t}\cdot(\bm{\tau}\mathrm{i}\tau^{y})_{\sigma\sigma^{\prime}}$ , and then obtain

\displaystyle K_{{\rm para},xx}^{F}

\displaystyle=2T\sum_{n\bm{k}}(v_{\bm{k}}^{x})^{2}\Big{[}|F_{s}(\bm{k},\mathrm{i}\omega_{n})|^{2}-|\bm{F}_{t}(\bm{k},\mathrm{i}\omega_{n})|^{2}\Big{]}.

(28)

If we consider the simple $\eta$ -pairing with only spin-singlet part ( $\bm{F}_{t}=\bm{0}$ ), it leads to the paramagnetic response (positive).
Thus, a simple $s$ -wave spin-singlet $\eta$ -pairing is unlikely realized as a stable state. On the other hand, in the AH model with magnetic field, the spin-triplet pair contribution is substantially generated by the Zeeman field, which plays an important role for the diamagnetic response as shown below.

IV Numerical result for AH Model

IV.1 Square lattice

IV.1.1 Prerequisites

Let us begin with the analysis of the AH model on the square lattice. We consider the two-sublattice structure to describe the staggered ordered phase such as a $\eta$ -pairing. While the superconducting states in the attractive model are interpreted in terms of the magnetic phases of the repulsive model by the particle-hole transformation in Eq. (2), the response functions such as the Meissner kernel are specific to the attractive model and have not been explored.

In the following, we choose the band width $W=1$ as the unit of energy. We fix the value of the attractive interaction $U=-1.375$ . The electron concentration is fixed as $n_{c}=1$ , and the temperature is taken to be $T=1.0\times 10^{-3}$ unless otherwise specified. We will investigate the change of the Meissner kernel for $\eta$ -pairing as a function of magnetic field strength $h=|\bm{h}|$ . In this paper, the mean-field solutions are calculated using the $60\times 60$ mesh in $\bm{k}$ -space. The result of the Meissner kernel for $\eta$ -pairings is calculated with the mesh $300\times 300$ . We also checked that the behaviors remain qualitatively unchanged when these numbers are increased. The self-consistent equations in Eqs. (4)-(6) are computed by using an iterative method. In the following subsection IV.1.2, we restrict ourselves to the analysis of two-sublattice mean-field solutions, and in IV.1.3, we examine the solutions when the two-sublattice constraint is relaxed.

IV.1.2 Two-sublattice solution

Before investigating the electromagnetic stability, we clarify the regime where the $\eta$ -pairing becomes the ground state. In this paper, we assume that the internal energy in Eq. (1) is approximately equal to the free energy in the low temperature region. The upper panel of Fig. 2 shows the internal energy of several ordered states measured from the normal-state energy as a function of the Zeeman field $h$ . Here, the $\eta$ -pairing solution is obtained by solving the self-consistent equation with imposing the constraint of the staggered phase of the pair amplitude. A constraint is also used for the calculation of the other types of order parameters. Our calculations have not found any ordered states other than the types shown in Fig. 2 even when a random initial condition is employed.

We determine the thermodynamically stable ground state by comparing the internal energies. In low magnetic fields, BCS and CDW are degenerated ground states. On the other hand, we find that the $\eta$ -pairing becomes the ground state in the magnetic field located in $1.063<h<1.875$ . The $\eta$ -pairing solution itself is found in the wider regime although the internal energy is not the lowest one. It has been known that the attractive Hubbard model under a magnetic field also shows the FFLO state [50], but this possibility cannot be considered when we take the two-sublattice condition. This point will be revisited in the next subsection where the two-sublattice condition is relaxed.

The lower panel of Fig. 2 shows the density of state (DOS) at the Fermi level for each state. The result indicates that there is no energy gap in the $\eta$ -pairing state, in contrast to the conventional BCS pairing state. There exists the regime where the DOS at the Fermi level for $\eta$ -pairing is larger than that of normal metal ( $1.25\lesssim h\lesssim 1.5$ ). This is due to the van-Hove singularity of the square lattice model as shown in FIG. 3. We also perform the calculation for the cubic lattice where the van-Hove singularity is absent at zero energy and confirm in this case that the DOS is smaller than the normal state (see Appendix B).

The stability of the $\eta$ -pairing depends upon the magnitude of the magnetic field as seen in the Meissner response kernel $K$ ( $=K_{xx}=K_{yy}$ ) (green symbol) in Fig. 4(a). The contributions from the paramagnetic ( $K_{\rm{para}}$ , positive) and diamagnetic ( $K_{\rm{dia}}$ , negative) parts are also separately plotted in the figure. In the regime with $h\leq 1.125$ and $1.75\leq h$ , the $\eta$ -pairing is electromagnetically unstable, while it is stable in $1.125<h<1.75$ . In Fig. 4, the yellow shaded rectangle indicates the regime where the $\eta$ -pairing becomes the ground state as seen from Fig. 2. We find a narrow region where $\eta$ -pairing is regarded as the ground state but is not an electromagnetically stable state around $h=1.125$ . From these results, we see that the $\eta$ -pairing is not necessarily electromagnetically stable even if it becomes the ground state in a two-sublattice calculation. As we shall see later, the simple $\eta$ -pairing in this narrow regime does not necessarily exist if we relax the two-sublattice condition of the mean-field solution.

We also show in Fig. 4(a) the contributions from the even- and odd-frequency pairs defined in Eqs. (23) and (24). The negative sign of the kernel, which means the response is diamagnetic, is partly due to the odd-frequency component of the pair amplitude, ( $K^{\rm{OFP}}<0$ ). This is in contrast to the FFLO state whose Meissner kernel is also negative due to the even-frequency component [51]. Hence, it implies that the mechanism of the diamagnetism is different between the FFLO and $\eta$ -pairing states.

In addition to the Meissner kernel, we calculate the local pair amplitudes which are shown in FIG. 4(b). Here the left- and right-panels represent the spin-triplet and spin-singlet components of the local pair amplitude, respectively. The triplet component $\sum_{\sigma\sigma^{\prime}}(\tau^{\mu}\mathrm{i}\tau^{y})_{\sigma\sigma^{\prime}}F_{\sigma\sigma^{\prime}}(\mathrm{i}\omega_{n})$ with $\mu=x$ has a finite imaginary part and zero real part, which represents the odd-frequency pair. The other $\mu=y,z$ components are zero. On the other hand, the singlet component has a finite real part and zero imaginary part and is the even-frequency pair. We can see that the maximum value of the spin-triplet component of the pair amplitude is largest at the magnetic field $h=1.375$ , where the magnitude of $K^{\rm{OFP}}$ is largest. It is also notable that the magnitude of the odd-frequency pair amplitude correlates with the magnitude of DOS at zero energy as seen by comparing Figs. 3 and 4.

We comment on the singular behavior of $K^{\rm{OFP}}$ at the magnetic field $h=1.375$ , although it does not affect the total Meissner kernel $K$ . This anomalous feature is related to the van Hove singularity of the DOS at zero energy as shown in FIG. 3, which shows a sharp peak at the Fermi level.

IV.1.3 Beyond two-sublattice

In order to clarify the stable ordered state where the Meissner kernel is positive (paramagnetic), we investigate mean-field solutions on finite-sized lattice where the two-sublattice condition is not imposed. We have numerically solved the Eqs. (4)-(6) self-consistently by using the mean-field solutions of the $\eta$ -pairing obtained for two-sublattice as an initial condition.

Figure 5 shows the spatial distribution of the phase of the gap function when the number of sites is $8\times 8$ . At $h=0.5$ in (a), where the $\eta$ -pairing is not a ground state, the uniform BCS pairing state is realized as expected. With increasing the magnetic field, the longer-periodicity structures are found as shown in Figs. 5(b), (c) and (d). At $h=1.375$ in (c), where the $\eta$ -pairing solution has the lowest energy and the electromagnetic response is well diamagnetic, we obtain the staggered alignment of the phases. When the parameters are close to the edges of the yellow-highlighted region in Fig. 4, the complex structures are formed as shown in (b) and (d). The behavior in (b) is interpreted as due to the competing effect where the simple uniform and staggered phases are energetically close to each other.

We also investigate the case with the other choice of parameters: $U=-1.25$ and $h=1.25$ . In this case, we find the staggered flux state where the phase of pair potential is characterized by 90^∘-Néel ordering as in Fig. 6(a). This ordered state cannot be described in the mean-field theory with two sublattices. Owing to a non-colinear 90^∘-Néel ordering vector, the spontaneous clockwise or counterclockwise loop currents arise by the inter-atomic Josephson effect. The current density is calculated by

\displaystyle j_{ij}

\displaystyle=-{\rm{i}}t\displaystyle\sum\limits_{\sigma}\langle c_{i\sigma}^{\dagger}c_{j\sigma}-c_{j\sigma}^{\dagger}c_{i\sigma}\rangle

(29)

which is identical to the expression of the paramagnetic current in the linear response theory. We can also evaluate the flux for each plaquette, which is define by

\displaystyle\Phi

\displaystyle=\sum_{(i,j)\in{\rm plaquette}}j_{ij}

(30)

This expression is similar to the flux $\displaystyle\oint_{C}\bm{j}\cdot d\bm{s}=\int_{S}\bm{b}\cdot d\bm{S}$ ( $\bm{j}=\bm{\nabla}\times\bm{b}$ ) defined in a continuum system, where $\bm{b}$ is a flux density. The flux is aligned in a staggered manner on a dual lattice as indicated in Fig. 6(b). The staggered flux originating from the normal part has been studied before [20, 21, 22, 23], while the staggered flux shown in Fig. 6(b) has a different origin: it arises from the superconductivity associated with the off-diagonal part in the Nambu representation.

We also comment on a feedback effect to the electromagnetic field from the supercurrent. Since the characteristic length scale for the magnetic field in layered superconductor becomes long [52], each magnetic flux on the plaquette is smeared out with this length. Hence we expect that the net magnetic field is not created from the staggered superconducting flux.

IV.2 Triangular lattice

IV.2.1 Mean-field solution

Now we search for the $\eta$ -pairing reflecting the characteristics of a geometrically frustrated triangular lattice at the half-filling ( $n_{c}=1.0$ ). We choose the parameters $U=-1.83$ and $T=1.0\times 10^{-3}$ . We consider the cases of two- and three-sublattice structures. For a usual antiferromagnet, the typical ordered state in the two-sublattice case has a stripe pattern, while in the three-sublattice case we expect a 120^∘-Néel state. Below we study the superconducting $\eta$ -pairing phases within the mean-field theory.

We have found the four types of superconducting states reflecting the characteristics of the triangular lattice, which are referred to as the $\eta$ -pairing I, II, III, and IV. The schematic pictures for these four states are shown in Fig. 7(a), where the arrow indicates the phase of the superconducting order parameter at each site. We make a few general remarks: the three-sublattice structure is assumed for I, II, III, while the two sublattice is employed for IV. The type-I has a non-colinear structure, and in the other $\eta$ -pairings the vectors are aligned in a colinear manner. We also note that CDW accompanies the $\eta$ -pairings II and III, where the number of local filling is indicated by the size of the filled circle symbols in Fig. 7(a).

Figure 7(b) shows the internal energy of the ordered states measured from the normal state (Upper panel) and from the $\eta$ -pairing I (Lower panel). From the lower panel of Fig. 7(b), we can identify the ground state. With increasing the magnetic field, the ground state changes as BCS $\rightarrow$ $\eta$ -pairing II $\rightarrow$ $\eta$ -pairing I $\rightarrow$ $\eta$ -pairing IV $\rightarrow$ $\eta$ -pairing I $\rightarrow$ normal. Figure 7(c) shows the particle density and $x$ -direction magnetization $m_{i}^{x}$ of each sublattice for $\eta$ -pairing II (Upper panel) and $\eta$ -pairing IV (Lower panel). The values of $m_{i}^{y}$ and $m_{i}^{z}$ are zero because the Zeeman field $\bm{h}$ is applied along the $x$ -direction. Below, we explain the characteristic features for each $\eta$ -pairing state.

$\eta$ -pairing-I state.— The $\eta$ -pairing I has 120^∘ Néel ordering vector (Green pentagon in Fig. 7(b)). The spontaneous supercurrent appears in this non-colinear state as schematically shown in Fig. 8(a). This superconducting state forms a staggered flux state, where the flux is aligned on a honeycomb dual lattice, which is similar to the $\eta$ -pairing with 90^∘-Néel ordering vector on the square lattice shown in Fig. 6(b). Figure 8(b) displays the values of spontaneous loop current density as a function of the magnetic field.

$\eta$ -pairing-II state.— The $\eta$ -pairing II has the structure with up-up-down colinear phases plus CDW (Red hexagon in Fig. 7(b)). There is the relation $n_{\rm{A}}=n_{\rm{B}}<n_{\rm{C}}$ for the electron filling at each sublattice shown in Fig. 7(c). We note that this site-dependent feature is characteristic for the II (and IV) state. The phases of the pair potential at A and B sublattices are “ferromagnetic”, while the phase at C sublattice is “antiferromagnetic”. The resulting ordered state is regarded as the emergence of the honeycomb lattice formed by equivalent A and B sublattices.

$\eta$ -pairing-III state.— This is the $\eta$ -pairing with a staggered ordering vector and CDW (Magenta square in Fig. 7(b)). The order parameter $\Delta$ at C sublattice is zero, but the others (A,B) are finite. The electron-rich sublattices A and B form a simple bipartite $\eta$ -pairing state on an emergent honeycomb lattice. Since this state does not become a ground state anywhere for the present choice of $U=-1.83$ , we do not further investigate this state in the following.

$\eta$ -pairing-IV state.— This is the $\eta$ -pairing with a simple stripe alignment (Cyan rhombus in Fig. 7(b)). This $\eta$ -pairing is accompanied by CDW around $h=1.9$ shown in Fig. 7(c). As shown below, this stripe phase show an anisotropic behavior in linear response coefficients, while the other $\eta$ -pairing states are isotropic.

IV.2.2 Meissner response

Now we discuss the Meissner response. Figure 9(a,b,c) shows the Meissner kernels $K_{xx},~{}K_{yy}$ for the $\eta$ -pairing I, II and IV. The yellow-highlighted parts indicate the region where each $\eta$ -pairing becomes the ground state as identified from Fig. 7(b). The result for the $\eta$ -pairing III is not shown because it does not become a ground state at $U=-1.83$ . We confirm that the Meissner response is basically diamagnetic if the $\eta$ -pairing becomes the ground state as shown in Figs. 9(a,b,c). Thus the energetic stability and diamagnetic response are reasonably correlated. In the following, we discuss the properties of the Meissner kernel for each state.

The Meissner kernels for both $\eta$ -pairing I and $\eta$ -pairing II shown in Figs. 9(a) and (b) satisfy the relation $K_{xx}=K_{yy}$ , which means an isotropic linear response. For the $\eta$ -pairing I, the Meissner kernel becomes positive in the regions $h<1.2,~{}1.95<h<2.12$ , while the kernel becomes negative in the ground state region (Fig. 9(a)). Although the local current density is finite for the $\eta$ -pairing I state, it does not affect the expression of the Meissner kernel in Eq. (10) since the total current $\bm{j}(\bm{q}=0)$ is zero.

Next we disucuss the $\eta$ -pairing IV state. The Meissner kernel jumps at $h=1.8$ due to the emergence of the CDW order parameter as shown in Fig. 9(c1,c2). It is notable that the $\eta$ -pairing IV with the stripe pattern shows a difference between $x$ and $y$ directions as shown in Figs. 9(c1,c2), respectively. This characteristic behavior can be intuitively understood from Fig. 7(a), where the current along the $x$ -axis flows with experiencing a staggered pair potential, whereas the current in the $y$ -direction feels an uniform pair potential. In the Meissner response, $K_{xx}$ shows a characteristic behavior of the $\eta$ -pairing, while $K_{yy}$ is qualitatively the same as the kernel of BCS. Thus, as shown in Fig. 9(c1), the diamagnetic response in the $x$ -axis direction is related to to the odd-frequency pair, whereas the diamagnetic response in the $y$ -axis direction, shown in Fig. 9(c2), is related to even-frequency pair.

V Summary and Outlook

We have studied the square and the triangular lattice of the attractive Hubbard model by using the mean-field theory. Several types of $\eta$ -pairing have been found in the triangular lattice where a simple bipartite pattern is not allowed. Using the formulation of the Meissner kernel for a general tight-binding lattice, we have investigated the electromagnetic stability of $\eta$ -pairings. We have confirmed that the electromagnetic stability of the $\eta$ -pairing correlates with the internal energy. In a narrow parameter range, we also find that the $\eta$ -pairing state can show an unphysical paramagnetic response if we assume the two or three sublattice structure in the mean-field calculation. In this case, another solution with longer periodicity needs to be sought.

When the current path experiences the staggered phase of the superconducting order parameter, the odd-frequency component of the pair amplitude contributes to the diamagnetic response. This is in contrast to the conventional BCS case in which the even-frequency component of the pair amplitude contributes to the diamagnetism. We have further clarified that one of the $\eta$ -pairing states on the triangular lattice has a stripe pattern and shows an anisotropic Meissner response. In this case, the odd-frequency pair contributes diamagnetically or paramagnetically depending on the direction of current.

We comment on some issues which are not explored in this paper. We expect that the $\eta$ -pairing without a simple staggered phase will appear on pyrochlore, kagome and quasicrystalline lattice, whose phase-alignment could be qualitatively different from the triangular lattice. In addition, there is another model that shows $\eta$ -pairing in equilibrium. A two-channel Kondo lattice (TCKL) is an example of a model in which $\eta$ -pairing appears even in the absence of a Zeeman field [24]. Our preliminary calculation for the TCKL shows a number of ordered states which have similar energies. These additional studies provide more insight into the exotic superconductivity characteristic for the $\eta$ -pairing.

Acknowledgement

This work was supported by KAKENHI Grants No. 18H01176, No. 19H01842, and No. 21K03459.

Appendix A Self-consistent equations in mean-field theory

We derive self-consistent equations for the general interacting Hamiltonian. Let us begin with the Hamiltonian

\displaystyle\mathscr{H}

\displaystyle=\sum_{12}\varepsilon_{12}c_{1}^{\dagger}c_{2}+\sum_{1234}U_{1234}c_{1}^{\dagger}c_{2}^{\dagger}c_{4}c_{3}

(31)

where site-spin indices are written as $1=(i_{1},\sigma_{1})$ . The mean-field Hamiltonian is introduced as

\displaystyle\mathscr{H}_{\rm MF}

\displaystyle=\sum_{12}\quantity(E_{12}c_{1}^{\dagger}c_{2}+\Delta_{12}c_{1}^{\dagger}c_{2}^{\dagger}+\Delta_{12}^{*}c_{2}c_{1}).

(32)

We assume $\langle\mathscr{H}\rangle=\langle\mathscr{H}_{\rm MF}\rangle$ where the statistical average is taken with $\mathscr{H}_{\rm MF}$ . Then the self-consistent equation is obtained as

$\displaystyle E_{12}$	$\displaystyle=\frac{\partial\langle\mathscr{H}\rangle}{\partial\langle c_{1}^{\dagger}c_{2}\rangle}$
	$\displaystyle=\varepsilon_{12}+\sum_{34}\quantity(U_{1324}+U_{3142}-U_{1342}-U_{3124})\langle c_{3}^{\dagger}c_{4}\rangle$	(33)
$\displaystyle\Delta_{12}$	$\displaystyle=\frac{\partial\langle\mathscr{H}\rangle}{\partial\langle c_{1}^{\dagger}c_{2}^{\dagger}\rangle}=\sum_{34}U_{1234}\langle c_{4}c_{3}\rangle$	(34)

where the Wick’s theorem is used for the derivation. Although the variational principle for the free energy also gives the same equation, the above formalism gives a simple procedure to derive the self-consistent equations.

Appendix B Attractive Hubbard model on Cubic lattice

We analyze the $\eta$ -pairing on the cubic lattice, whose DOS does not have a van Hove singularity near zero energy. Here we choose the parameter $U=-1.375$ and the electron concentration is half-filled. As a result, the DOS for the $\eta$ -pairing around zero energy for each magnetic filed on the cubic lattice is smaller than the DOS of the normal state as shown in Fig. 10(a). For reference, we also show in Figs. 10(b) and (c) the pair amplitude similar to Fig. 4(b) in the main text. In addition, the odd-frequency pair amplitude increases when DOS near zero energy is enhanced as seen from Figs. 10(a) and (b).

References

[1] P. Fulde and R. A. Ferrell, Phys. Rev. $\bm{135}$ , A550 (1964).
[2] A. I. Larkin and Y. N. Ovchinnikov, Zh. Eksp. Teor. Fiz. $\bm{47}$ , 1136 (1964).
[3] C. N. Yang, Phys. Rev. Lett. $\bm{63}$ , 2144 (1989).
[4] For a review, see, D.F. Agterberg, J.C.S. Davis, S.D. Edkins, E. Fradkin, D.J. Van Harlingen, S.A. Kivelson, P.A. Lee, L. Radzihovsky, J.M. Tranquada, and Y. Wang, Ann. Rev. Condens. Matter Phys., 11, 231 (2020).
[5] R. R. P. Singh and R. T. Scalettar, Phys. Rev. Lett. $\bm{66}$ , 3203 (1991).
[6] P. Coleman, E. Miranda, and A. Tsvelik, Phys. Rev. Lett. 70, 2960 (1993).
[7] S. Hoshino and Y. Kuramoto, Phys. Rev. Lett. $\bm{112}$ , 167204 (2014).
[8] W. R. Czart, K. J. Kapcia, R. Micnas, S. Robaszkiewicz, Physica A: Statistical Mechanics and its Applications. $\bm{585}$ , 126403 (2022)
[9] T. Kaneko, T. Shirakawa, S. Sorella, and S. Yunoki, Phys. Rev. Lett. $\bm{122}$ , 077002 (2019).
[10] P. Werner, J. Li, D. Golez and M. Eckstein, Phys. Rev. B 100, 155130 (2019).
[11] J. Li, D. Golez, P. Werner, and M. Eckstein, Phys. Rev. B $\bm{102}$ , 165136 (2020).
[12] M. Nakagawa, N. Tsuji, N. Kawakami, M. Ueda, arXiv:2103.13624 (2021).
[13] X. M. Yang and Z. Song, Phys. Rev. B $\bm{105}$ , 195132 (2022).
[14] L. Jiajun, M. Müller, A. J. Kim, A. M. Läuchli, and P. Werner, arXiv:2202. $\bm{10176}$ (2022).
[15] Y. Motome, K. Nakamikawa, Y. Yamaji, and M. Udagawa, Phys. Rev. Lett. $\bm{105}$ , 036403 (2010).
[16] R. R. dos Santos Phys. Rev. B $\bm{48}$ , 3976 (1993).
[17] I. Affleck and J. B. Marston, Phys. Rev. B $\bm{37}$ , 3774 (1988).
[18] S. Chakravarty, R. B. Laughlin, D. K. Morr, and C. Nayak, Phys. Rev. B $\bm{63}$ , 094503 (2001).
[19] D. K. Morr, Phys. Rev. Lett. $\bm{89}$ , 106401 (2002).
[20] S. Zhou and Z. Wang, Phys. Rev. B $\bm{70}$ , 020501(R) (2004).
[21] H. Yokoyama, S. Tamura, and M. Ogata, J. Phys. Soc. Jpn. $\bm{85}$ , 124707 (2016).
[22] K. Kobayashi, H. Yokoyama, and Y. Toga, J. Phys.: Conf. Ser. $\bm{1054}$ , 012016 (2018).
[23] S. Fukuda, K. Yamazaki, H. Tsuchiura, and M. Ogata, J. Phys.: Conf. Ser. $\bm{1293}$ , 012026 (2019).
[24] S. Hoshino, Phys. Rev. B $\bm{90}$ , 115154 (2014).
[25] V.L. Berezinskii, ZhETF Pis Red. 20, 628 (1974), [JETP Lett. 20, 287 (1974)].
[26] T. R. Kirkpatrick and D. Belitz, Phys. Rev. Lett. 66, 1533 (1991).
[27] A. Balatsky and E. Abrahams, Phys. Rev. B 45, 13125 (1992).
[28] Y. Tanaka, M. Sato, and N. Nagaosa, J. Phys. Soc. Jpn. 81, 011013 (2012).
[29] J. Linder and A.V. Balatsky, Rev. Mod. Phys. 91, 045005 (2019).
[30] J. Cayao, C. Triola and M. Black-Schaffer, Eur. Phys. J.: Spec. Top. $\bm{229}$ , 545 (2020)
[31] H. Walter, W. Prusseit, R. Semerad, H. Kinder, W. Assmann, H. Huber, H. Burkhardt, D. Rainer, and J. A. Sauls, Phys. Rev. Lett. 80, 3598 (1998).
[32] Y. Tanaka, Y. Asano, A. A. Golubov, and S. Kashiwaya, Phys. Rev. B 72, 140503 (2005).
[33] F. S. Bergeret, A. F. Volkov, and K. B. Efetov, Rev. Mod. Phys. $\bm{77}$ , 1321 (2005).
[34] Ya. V. Fominov, Y. Tanaka, Y. Asano, and M. Eschrig, Phys. Rev. B $\bm{91}$ , 144514 (2015)
[35] Y. Tanaka and A. A. Golubov, Phys. Rev. Lett. $\bm{98}$ , 037003 (2007)
[36] Y. Tanaka, Y. Tanuma, and A. A. Golubov, Phys. Rev. $\bm{B}$ 76, 054522 (2007)
[37] S. Higashitani, J. Phys. Soc. Jpn. 66, 2556 (1997).
[38] T. Yokoyama, Y. Tanaka, and N. Nagaosa, Phys. Rev. Lett. 106, 246601 (2011).
[39] S.-I. Suzuki and Y. Asano, Phys. Rev. B 91, 214510 (2015).
[40] A. Di Bernardo, Z. Salman, X. L. Wang, M. Amado, M. Egilmez, M. G. Flokstra, A. Suter, S. L. Lee, J. H. Zhao, T. Prokscha, E. Morenzoni, M. G. Blamire, J. Linder, and J. W. A. Robinson, Phys. Rev. X 5, 041021 (2015).
[41] J. A. Krieger, A. Pertsova, S. R. Giblin, M. Dbeli, T. Prokscha, C. W. Schneider, A. Suter, T. Hesjedal, A. V. Balatsky, and Z. Salman, Phys. Rev. Lett. 125, 026802 (2020).
[42] R. Micnas, J. Ranninger, and S. Robaszkiewicz, Rev. Mod. Phys. $\bm{62}$ , 113 (1990).
[43] H. Shiba, Prog. Theor. Phys. $\bm{48}$ , 2171 (1972).
[44] Y. Claveau, B. Arnaud, S. D. Matteo, Eur. J. Phys. $\bm{35}$ , 035023 (2014).
[45] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
[46] D.J. Scalapino, S. R. White, and S. C. Zhang, Phys. Rev. Lett. 68, 2830 (1992).
[47] D. J. Scalapino, S. R. White, and S. C. Zhang, Phys. Rev. B 47, 7995 (1993).
[48] T. Kostyrko, R. Micnas, and K. A. Chao, Phys. Rev. B 49, 6158 (1994).
[49] S. Hoshino, K. Yada, and Y. Tanaka, Phys. Rev. B 93, 224511 (2016).
[50] A. Tsuruta, S. Hyodo, and K. Miyake, J. Phys. Soc. Jpn. $\bm{83}$ , 063706 (2014).
[51] D. Chakraborty and A. M. Black-Schaffer, Phys. Rev. B $\bm{106}$ , 024511 (2022).
[52] J. Pearl, App. Phys. Lett. 5, 65 (1964).

η\eta-pairing on bipartite and non-bipartite lattices