# On the behaviour of a periodically forced and thermostatted harmonic chain

Pedro L. Garrido  Pedro L. Garrido
Universidad de Granada
Granada Spain<EMAIL_ADDRESS>,  Tomasz Komorowski  Tomasz Komorowski,
Institute of Mathematics, Polish Academy Of Sciences, Warsaw, Poland
and Institute of Mathematics, Maria Curie-Skłodowska University, Lublin,
Poland<EMAIL_ADDRESS>,  Joel L. Lebowitz  Joel L. Lebowitz,
Departments of Mathematics and Physics, Rutgers University
<EMAIL_ADDRESS>and  Stefano Olla  Stefano Olla, CEREMADE,
Université Paris-Dauphine, PSL Research University
and Institut Universitaire de France
and GSSI, L’Aquila<EMAIL_ADDRESS>
(Date: September 3, 2024File: no-flip-submitted-jsp.tex.)

###### Abstract.

We consider a chain consisting of $n+1$ pinned harmonic oscillators subjected
on the right to a time dependent periodic force ${\mathcal{F}}(t)$ while
Langevin thermostats are attached at both endpoints of the chain. We show that
for long times the system is described by a Gaussian measure whose covariance
function is independent of the force, while the means are periodic. We compute
explicitly the work and energy due to the periodic force for all $n$ including
$n\to\infty$.

###### Key words and phrases:

harmonic chain, periodic force, work into heat, resonance response

###### 2000 Mathematics Subject Classification:

80A19,80M22,82C22,82C70,70J35

We thank David Huse for helpful discussions, and the IAS for the hospitality
during part of this work. P.G. acknowledges the support of the Project I+D+i
Ref.No.PID2020-113681GB-I00, financed by MICIN/AEI/10.13039/501100011033 and
FEDER “A waytomakeEurope.” T.K. acknowledges the support of the NCN grant
2020/37/B/ST1/00426.

##  1\. Introduction

In this work we consider the conversion of work into heat in a simple model
system: a pinned harmonic chain of $n+1$ particles on which work is performed
by an external periodic force acting at one of the endpoints. The system is
also in contact with thermal reservoirs, placed at both of its endpoints,
which absorb the energy generated by the work. In the absence of the
reservoirs the response of the system to the external forcing depends entirely
on whether the frequency $\omega$ of the external force coincides with the
normal frequencies of the chain $\\{\omega_{j},j=0,\ldots,n\\}$. When
$\omega\neq\omega_{j}$ the system adjusts itself to be out of phase with the
force so that there is no work done on the average. If on the other hand the
system is in resonance with the force, i.e. $\omega=\omega_{j}$ for some $j$,
then the amplitude of the oscillation tends to infinity as time $t\to\infty$.

The situation is different in the presence of the thermostats. They cause the
oscillations at resonance to be damped and as a result the work done by the
force is strictly finite for all values of $\omega$.

There is still a strong dependence on $\omega$, as far as the magnitude of the
work is concerned, when $n$ gets large. This difference becomes qualitative
when $n\to\infty$ and the spectrum of the harmonic chain becomes dense in an
interval $\cal I$. The work done and the internal energy of the chain depend
strongly on whether $\omega$ lies in the interior of $\cal I$, or not.

Due to the linearity of the system there is a clear division, in the long time
properties of the system, between those due to temperatures of the thermal
reservoirs and those due to the external force. The energy flowing through the
system as a result of the presence of the thermal reservoirs we call _thermal
energy_. It is not influenced by the external force and its behavior is the
same as in [10] and [8]. The energy flow due to the work of the external force
we call _mechanical energy_. It is independent of the temperatures of the
reservoirs, and it is influenced only by the corresponding damping. For finite
$n$ and pure damping equal on both sites this was computed in section 4 of [9]
in terms of the Green function of the corresponding damped harmonic chain. The
main objective of the present note is the exact calculation of the asymptotic
behavior as $n\to\infty$ of the work and the mechanical energy. Calculation of
these quantities, turns out to be quite complicated, but leads to explicit
expressions for their asymptotics. In particular we show that, for forcing
frequency outside $\cal I$, the work, the mechanical energy and its flow
become negligible as $n\to\infty$. Inside $\cal I$ these quantities oscillate
fast and their asymptotic behavior can be described in terms of Young
measures.

The results of the present work remain also valid in the case of unpinned
harmonic chain. It suffices to set the pinning constant $\omega_{0}=0$ in our
formulas describing the work and energy functionals. Obviously we now always
have $\omega>\omega_{0}=0$ and consider the motion relative to the center of
mass positioned at zero.

For anharmonic interaction the situation is qualitatively very different. The
non-linearity produces many new effects described in [9]. We have also studied
the case of a harmonic chain with a random velocity flip [7].

##  2\. Description of the System

The configurations of our system, consisting of $n+1$ pinned harmonic
oscillators, are described by

$(\mathbf{q},\mathbf{p})=(q_{0},\dots,q_{n},p_{0},\dots,p_{n})\in\Omega_{n}:={\mathbb{R}}^{n+1}\times{\mathbb{R}}^{n+1}.$
(2.1)

We should think of the positions $q_{x}$ as the relative displacement of an
atom from a point $x$ belonging to the integer lattice interval
${{\mathbb{I}}_{n}=\\{0,\ldots,n\\}}$ and $p_{x}$ as its respective momentum.

The total energy of the chain is given by the Hamiltonian:

$\mathcal{H}_{n}(\mathbf{q},\mathbf{p}):=\sum_{x=0}^{n}e_{x}(\mathbf{q},\mathbf{p}),$
(2.2)

where the microscopic energy density at $x$ is given by

$e_{x}(\mathbf{q},\mathbf{p}):=\frac{p_{x}^{2}}{2}+\frac{1}{2}(q_{x}-q_{x-1})^{2}+\frac{\omega_{0}^{2}q_{x}^{2}}{2},\quad
x\in{\mathbb{I}}_{n}.$ (2.3)

Here we let $q_{-1}:=q_{0}$.

The microscopic dynamics of the process describing the total chain is given by

$\displaystyle\dot{q}_{x}(t)$ $\displaystyle=p_{x}(t),\qquad\qquad
x\in{\mathbb{I}}_{n},$ (2.4) $\displaystyle\dot{p}_{x}(t)$
$\displaystyle=\Delta_{x}q_{x}-\omega_{0}^{2}q_{x},\quad
x\in{\mathbb{I}}_{n}^{o}=\\{1,\dots,n-1\\}$

and at the boundaries by

$\displaystyle\mathrm{d}p_{0}(t)$
$\displaystyle=\;\Big{(}q_{1}(t)-q_{0}(t)-\omega_{0}^{2}q_{0}\Big{)}\mathrm{d}t-2\gamma_{-}p_{0}(t)\mathrm{d}t+\sqrt{4\gamma_{-}T_{-}}\mathrm{d}\widetilde{w}_{-}(t),\vphantom{\Big{(}}$
(2.5) $\displaystyle\mathrm{d}p_{n}(t)$
$\displaystyle=\;\Big{(}q_{n-1}(t)-q_{n}(t)-\omega_{0}^{2}q_{n}(t)\Big{)}\mathrm{d}t+\;{\mathcal{F}}(t/\theta)\mathrm{d}t-2\gamma_{+}p_{n}(t)\mathrm{d}t+\sqrt{4\gamma_{+}T_{+}}\mathrm{d}\widetilde{w}_{+}(t).\vphantom{\Big{(}}$

Here $\Delta q_{x}=q_{x+1}+q_{x-1}-2q_{x}$, $x\in{\mathbb{Z}}$ is the
laplacian on the integer lattice ${\mathbb{Z}}$, $\omega_{0}>0$ is a pinning
constant, $\widetilde{w}_{-}(t)$ and $\widetilde{w}_{+}(t)$ are two
independent standard one dimensional Wiener processes and $\gamma_{\pm}$ are
non-negative constants that describe the respective strengths of the Langevin
thermostats.

We assume the force ${\mathcal{F}}(t)$ to be a smooth periodic function of
period $1$ and parameter $\theta$ rescales the period. We will suppose,
without losing generality, that

$\int_{0}^{1}{\mathcal{F}}(t)\mathrm{d}t=0,\qquad\int_{0}^{1}{\mathcal{F}}(t)^{2}\mathrm{d}t>0.$
(2.6)

The generator of the dynamics is given by

$\mathcal{G}_{t}=\mathcal{A}_{t}+2\sum_{\iota\in\\{-,+\\}}\gamma_{\iota}S_{\iota},$
(2.7)

where

$\mathcal{A}_{t}=\sum_{x=0}^{n}p_{x}\partial_{q_{x}}+\sum_{x=0}^{n}(\Delta
q_{x}-\omega^{2}_{0}q_{x})\partial_{p_{x}}+{\mathcal{F}}(t/\theta)\partial_{p_{n}}.$
(2.8)

By convention we let $q_{n+1}:=q_{n}$ and $q_{-1}=q_{0}$. Furthermore

$S_{-}=T_{-}\partial_{p_{0}}^{2}-p_{0}\partial_{p_{0}},\qquad
S_{+}=T_{+}\partial_{p_{n}}^{2}-p_{n}\partial_{p_{n}}$ (2.9)

The energy currents are

$\displaystyle\mathcal{G}_{t}e_{x}=j_{x-1,x}-j_{x,x+1},$ (2.10) $\displaystyle
j_{x,x+1}:=-p_{x}(q_{x+1}-q_{x}),\qquad\text{if }x\in\\{0,...,n-1\\},$

and at the boundaries

$j_{-1,0}:=2{\gamma}_{-}\left(T_{-}-p_{0}^{2}\right)\qquad
j_{n,n+1}:=-2{\gamma}_{+}\left(T_{+}-p_{n}^{2}\right)-{\mathcal{F}}(t/\theta)p_{n}.$
(2.11)

We are interested in the long time behavior of the system. In the absence of
the external forcing, ${\mathcal{F}}(t)\equiv 0$, this is just the model
considered in [10], with $\omega_{0}=0$, and in [8] for $\omega_{0}>0$. In the
case when ${\mathcal{F}}(t)\equiv 0$, starting with any initial configuration
$(\mathbf{q}(0),\mathbf{p}(0))$ (or any initial probability distribution
$\mu_{0}(\mathrm{d}\mathbf{q},\mathrm{d}\mathbf{p})$) the system approaches a
stationary Gaussian distribution $\mu_{\rm
stat}(\mathrm{d}\mathbf{q},\mathrm{d}\mathbf{p})$, in which the expectation
values of $q_{x}$ and $p_{x}$ vanish, i.e. $\overline{q}_{x}(t)=0$ and
$\overline{p}_{x}(t)=0$, while the covariances between components of
$(\mathbf{q},\mathbf{p})$ are given explicitly.

In particular the expectation of the energy current $\overline{j}_{x,x+1}$
between sites $x$ and $x+1$, that is independent of $x$ and $t$, is given by

$\overline{j}_{x,x+1}=(c+o(1))(T_{-}-T_{+}),\quad\mbox{as }n\gg 1,$ (2.12)

with

$c=\frac{\gamma}{1+4\gamma^{2}+2\gamma\omega_{0}(\gamma\omega_{0}+\sqrt{1+4\gamma^{2}+(\gamma\omega_{0})^{2}})},$

when $\gamma_{-}=\gamma_{+}=\gamma$, see [8, formula (37), p. 240]. In the
case $\omega_{0}=0$ the term $o(1)$ in the formula (2.12) can be omitted (no
dependence on $n$) and we have $c=\frac{\gamma}{1+4\gamma^{2}}$, see [8,
formula (40), p. 241].

Eq. (2.12) implies that the thermal conductivity is proportional to $n$ \- the
size of the system - and becomes infinite in the limit $n\to+\infty$, see also
[10]. In fact the ”temperature” $T_{x}$, defined as the variance of
$p_{x}^{2}$, is independent of $x$, except near the boundary points $x=0,n$.
Adding now the periodic force of period $\theta$ leads, as $t\to+\infty$, to a
Gaussian, periodic stationary state $\\{\mu_{t}^{P},t\in[0,+\infty)\\}$, whose
covariances are the same as in the case when no force is applied. For any
functions $F=F(\bf q,\bf p)$ and $G=G(t)$ define

$\overline{F}(t)=\int_{\Omega_{n}}F\mathrm{d}\mu_{t}^{P}\quad\mbox{and}\quad\langle\langle
G\rangle\rangle=\frac{1}{\theta}\int_{0}^{\theta}G(t)\mathrm{d}t.$ (2.13)

The periodic stationary state has the property that
$\langle\langle\overline{\mathcal{G}F}\rangle\rangle=0$ for any $F$ in the
domain of $\mathcal{G}_{t}$.

The expectation values of the position and momentum $\overline{q}_{x}(t)$ and
$\overline{p}_{x}(t)$ are now $\theta$-periodic and independent of the
temperature of the reservoirs. They are given by

$\left(\begin{array}[]{c}\overline{\bf q}(t)\\\ \overline{\bf
p}(t)\end{array}\right)=\int_{-\infty}^{t}e^{-A(t-s)}\;{\mathcal{F}}(s/\theta){\rm
e}_{p,n+1}\mathrm{d}s.$

Here $A$ is a $2\times 2$ block matrix made of $(n+1)\times(n+1)$ matrices of
the form

$A=\left(\begin{array}[]{cc}0&-{\rm Id}_{n+1}\\\ -\Delta_{\rm
N}+\omega_{0}^{2}&\Gamma\end{array}\right),$

where ${\rm Id}_{n+1}$ is the $(n+1)\times(n+1)$ identity matrix, $\Delta_{\rm
N}$ is the Neumann laplacian on ${\mathbb{I}}_{n}$:

$\Delta_{\rm N}f_{x}:=\Delta f_{x},\quad
x\in{\mathbb{I}}_{n}^{o}\quad\mbox{and }\quad\Delta_{\rm
N}f_{0}=f_{1}-f_{0},\quad\Delta_{\rm N}f_{n}=f_{n-1}-f_{n}.$ (2.14)

Furthermore $\Gamma$ is the diagonal matrix

$\Gamma=2\left(\begin{array}[]{ccccc}\gamma_{-}&0&\ldots&0&0\\\
0&0&\ldots&0&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\ 0&0&\ldots&0&0\\\
0&0&\ldots&0&\gamma_{+}\end{array}\right).$

The column vector ${\rm e}_{p,n+1}$ is given by ${\rm
e}_{p,n+1}^{T}=[\underbrace{0,\ldots,0}_{2n+1-\rm{times}},1].$ Notice that the
first of the conditions (2.4) implies that
$\langle\langle\overline{p}_{x}\rangle\rangle=0$, while the second gives
$\langle\langle\overline{q}_{x}\rangle\rangle=0$.

The expected value of energy, averaged over a period, breaks up into the
mechanical part, coming from the averaged position $\overline{\bf q}(t)$ and
momentum $\overline{\bf p}(t)$, which is independent of the temperature of the
reservoirs, and the thermal part, which is independent of the external force.
More precisely

$\langle\langle e_{x}\rangle\rangle=\langle\langle e_{x}^{\rm
mech}\rangle\rangle+\langle\langle e_{x}^{\rm th}\rangle\rangle,$ (2.15)

where the mechanical component of the energy is given by

$e_{x}^{\rm
mech}(t):=\frac{1}{2}\Big{[}\overline{p}_{x}^{2}(t)+\omega_{0}^{2}\overline{q}_{x}^{2}(t)+\big{(}\overline{q}_{x}(t)-\overline{q}_{x-1}(t)\big{)}^{2}\Big{]},\quad
x\in{\mathbb{I}}_{n}.$ (2.16)

and the thermal part is

$e_{x}^{\rm
th}(t):=\frac{1}{2}{\mathbb{E}}\Big{[}p_{x}^{\prime}(t)^{2}+\omega_{0}^{2}q_{x}^{\prime}(t)^{2}+\big{(}q_{x}^{\prime}(t)-q_{x-1}^{\prime}(t)\big{)}^{2}\Big{]},\quad
x\in{\mathbb{I}}_{n}.$ (2.17)

where $q_{x}^{\prime}(t)=q_{x}(t)-\overline{q}_{x}(t)$ and
$p_{x}^{\prime}(t)=p_{x}(t)-\overline{p}_{x}(t)$ and ${\mathbb{E}}$ denotes
the average with respect to the initial data and the realizations of the
Wiener processes in (2.5). As before, we adopt the convention
$\overline{q}_{-1}(t):=\overline{q}_{0}(t)$ and likewise
$q_{-1}^{\prime}(t):=q_{0}^{\prime}(t)$.

As already mentioned in the Introduction one of the goals of the present paper
is to describe the work done by the force on the system. It is given by

$W(n)=\frac{1}{\theta}\int_{0}^{\theta}{\mathcal{F}}(t/\theta)\overline{p}_{n}(t)\mathrm{d}t.$
(2.18)

$W(n)$ is always positive, generates energy fluxes into the two heat
reservoirs. Furthermore, we describe the time average of the mechanical energy
functional given by eq. (2.16). Its thermal counterpart does not depend on
time and has been described in [8, 10]. We mention here also that the case
$n=0$, i.e. a single oscillator in contact with a heat bath and driven by an
external unbiased time-periodic force, has been fully characterized in [11].

##  3\. Results

In what follows we will use the dispersion relation of the infinite chain
given by

$\omega(r)=\sqrt{\omega_{0}^{2}+4\sin^{2}\left(\frac{\pi r}{2}\right)},\qquad
r\in[0,1]$ (3.1)

and its inverse defined for $\omega\in\cal
I:=[\omega_{0},\sqrt{\omega_{0}^{2}+4}]$ by the formula

$r(\omega)=\frac{2}{\pi}\arcsin\left(\frac{1}{2}\sqrt{\omega^{2}-\omega_{0}^{2}}\right),$
(3.2)

###  3.1. Work done by the force on the system

The work $W(n)$ performed by the force on the system, see (2.18), depends on
the period $\theta$. Considering for simplicity the simple mode case when

${\mathcal{F}}(t/\theta)=F\cos(\omega t),\quad\omega:=\frac{2\pi}{\theta}$
(3.3)

the work done is given by (see Appendix):

$W(\omega,n)=\big{(}\omega F\big{)}^{2}\frac{N(\omega,n)}{D(\omega,n)}.$ (3.4)

Here

$\displaystyle N(\omega,n)=\ $
$\displaystyle\gamma_{-}G^{1}(\omega,n)^{2}+\gamma_{+}G^{0}(\omega,n)^{2}+4{\gamma_{-}^{2}\gamma_{+}}\omega^{2}\left(G^{0}(\omega,n)^{2}-G^{1}(\omega,n)^{2}\right)^{2}$
$\displaystyle D(\omega,n)=\ $ $\displaystyle
1+8\gamma_{-}\gamma_{+}{\omega^{2}}G^{1}(\omega,n)^{2}+4\omega^{2}G^{0}(\omega,n)^{2}(\gamma_{-}^{2}+\gamma_{+}^{2})$
$\displaystyle+16\gamma_{-}^{2}\gamma_{+}^{2}\omega^{4}\left(G^{0}(\omega,n)^{2}-G^{1}(\omega,n)^{2}\right)^{2},$
(3.5)

where

$G^{s}(\omega,n)=G_{0,sn}(\omega,n),\quad s=0,1,$ (3.6)

and

$G_{x,y}(\omega,n)=\frac{1}{n+1}\cdot\frac{1}{\omega_{0}^{2}-\omega^{2}}+\frac{2}{n+1}\sum_{j=1}^{n}\frac{\cos(\frac{\pi
j(2x+1)}{2(n+1)})\cos(\frac{\pi
j(2y+1)}{2(n+1)})}{\omega_{j}^{2}-\omega^{2}},\quad x,y\in{\mathbb{I}}_{n}$
(3.7)

is the Green’s functions of $-\Delta_{N}+\omega_{0}^{2}-\omega^{2}$, and
$\pm\omega_{j}$, $j=0,\ldots,n$ are the eigenvalues of
$-\Delta_{N}+\omega_{0}^{2}$ defined by
$\omega_{j}=\omega\left(\frac{j}{n+1}\right)$ where $\omega(r)$ is given by
(3.1).

It is easy to see from (3.1) that $4\omega^{2}\gamma_{-}N\leqslant
D+\gamma_{-}^{2}G^{1}(\omega,n)^{2}$. Therefore, the following bound can be
found

$W(\omega,n)\leqslant\frac{\big{(}\omega
F\big{)}^{2}}{4}\Big{(}\frac{1}{\gamma_{-}}+\frac{1}{\gamma_{+}}\Big{)},\quad
n=1,2,\ldots.$ (3.8)

The functions $G^{s}(\omega,n)$ can be computed explicitly:

$G^{s}(\omega,n)=\frac{1}{n+1}\cdot\frac{1}{\omega_{0}^{2}-\omega^{2}}+\frac{2}{n+1}\sum_{j=1}^{n}\frac{(-1)^{js}\cos^{2}(\frac{\pi
j}{2(n+1)})}{\omega_{j}^{2}-\omega^{2}}\quad,\quad s=0,1,$ (3.9)

There are very different behaviors of $W(\omega,n)$ depending on whether
$\omega$ is in the spectrum of the harmonic chain, or not, see Figure 1.

Figure 1. Behavior of the work for $\omega_{0}=1$, $n=50$ with
$\gamma_{\pm}=1$ (left figure) and $\gamma_{-}=1$, $\gamma_{+}=1/10$ (right
figure). The red points are the values of work computed at the points
$\omega_{j}$ of the harmonic spectra using equation (3.11). Note the larger
magnitude of the vertical scale on the right diagram.

In particular, the formula (3.4) cannot be applied directly when
$\omega=\omega_{j}$ as then both $G^{s}(\omega,n)$, $s=0,1$ are divergent.
However, we can still use the formula to find $W(\omega_{j},n)$, because both
$N(\omega,n)$ and $D(\omega,n)$ have the same order of magnitude in the
neighborhood of $\omega_{j}$ and, due to the cancellation, the work remains
finite.

More precisely, assume that given $j$ we have
$\omega^{2}=\omega^{2}_{j}+\epsilon$ for some $\epsilon\ll 1/n$. The Green’s
functions can be then written in the form

$G^{s}(\omega,n)=\frac{2(-1)^{js}}{n+1}\cos^{2}\left(\frac{\pi
j}{2(n+1)}\right)\epsilon^{-1}+\overline{G}^{s}(\omega,n)$ (3.10)

where $\overline{G}^{s}(\omega,n)$ is of order $O(1)$ for $\epsilon\ll 1/n$.
We obtain then

$W(\omega_{j},n)=\frac{F^{2}}{4}\cdot\frac{\gamma_{+}+\gamma_{-}+16\gamma_{+}\gamma_{-}^{2}\omega_{j}^{2}S(\omega_{j},n)}{(\gamma_{+}+\gamma_{-})^{2}+16\gamma_{+}^{2}\gamma_{-}^{2}\omega_{j}^{2}S(\omega_{j},n)}$
(3.11)

where

$S(\omega_{j},n)=\overline{G}^{0}(\omega_{j},n)-(-1)^{j}\overline{G}^{1}(\omega_{j},n).$
(3.12)

In particular, (3.11) implies that

$W(\omega_{j},n)\sim\frac{F^{2}}{4\gamma_{-},},\quad\mbox{ as
}\gamma_{+}\rightarrow 0\quad\mbox{ and }\quad
W(\omega_{j},n)\sim\frac{F^{2}}{4\gamma_{+},},\quad\mbox{ as
}\gamma_{-}\rightarrow 0.$

If $\gamma_{+}=0$ the formula (3.4) for the work simplifies to

$W(\omega,n)=\left(\omega
F\right)^{2}\frac{\gamma_{-}G^{1}(\omega,n)^{2}}{1+4\gamma_{-}^{2}\omega^{2}G^{0}(\omega,n)^{2}}$
(3.13)

that gives $W(\omega,n)\to 0$, as $\gamma_{-}\to 0$, if
$\omega\neq\omega_{j}$. This means that outside the resonance frequences, no
work is done on the system if dissipation is absent. Recall also that when
$\gamma_{+}=\gamma_{-}=0$ and $\omega-\omega_{j}$ the stationary periodic
state does not exist as the energy keeps accumulating inside the system.

####  3.1.1. Work in the case $n\to+\infty$ when $\omega$ lies outside the
harmonic chain spectrum

Consider now the case $n\gg 1$. The spectrum becomes then the interval ${\cal
I}:=[\omega_{0},\sqrt{\omega_{0}^{2}+4}].$ For $\omega$ outside ${\cal I}$ the
right hand side of the formula for the Green’s function, see (3.9), does not
contain any singular term and $G^{s}(\omega,n)$ can be approximated by:

$\begin{split}G^{0}(\omega,n)&={2\int_{0}^{1}\frac{\cos^{2}\left(\frac{\pi
r}{2}\right)}{\omega_{0}^{2}-\omega^{2}+4\sin^{2}\left(\frac{\pi
r}{2}\right)}\;dr+O\Big{(}\frac{1}{n}\Big{)}}\\\
&=-\frac{1}{2}+\frac{1}{2\pi}\left(\omega_{0}^{2}+4-\omega^{2}\right)\int_{0}^{\pi}\,\frac{\mathrm{d}q}{\omega_{0}^{2}+2-\omega^{2}-2\cos
q}+O\Big{(}\frac{1}{n}\Big{)}.\end{split}$ (3.14)

Using [6, formula 2.553.3] we get 111Note that formula (3.15) makes also sense
in case $\omega_{0}=0$, as then any $\omega$ outside ${\cal I}$ satisfies
$\omega^{2}>4$.

$\overline{G}^{0}(\omega)=\lim_{n\rightarrow\infty}G^{0}(\omega,n)=-\frac{1}{2}+\frac{|\omega_{0}^{2}+4-\omega^{2}|}{2\sqrt{(\omega_{0}^{2}-\omega^{2})(\omega_{0}^{2}+4-\omega^{2})}}\quad\omega\notin{\cal
I}.$ (3.15)

Likewise, we can show

$\overline{G}^{1}(\omega)=\lim_{n\rightarrow\infty}G^{1}(\omega,n)=0,\quad\omega\notin{\cal
I}.$ (3.16)

Combining the above the work corresponding to $\omega$ outside the harmonic
spectra is given by

$\overline{W}(\omega)=\lim_{n\rightarrow\infty}W(\omega,n)=\frac{\gamma_{+}\Big{(}FH\Big{)}^{2}\Big{[}1+4\Big{(}\gamma_{-}H\Big{)}^{2}\Big{]}}{1+4(\gamma_{+}^{2}+\gamma_{-}^{2})H^{2}+16\Big{(}\gamma_{+}\gamma_{-}H^{2}\Big{)}^{2}},$
(3.17)

where

$\displaystyle H(\omega):=\omega\overline{G}^{0}(\omega).$ (3.18)

Observe that $\overline{W}(\omega)$ tends to $0$, when $\gamma_{+}\to 0$.
Likewise $\overline{W}(\omega)\to 0$, when either $\omega\to\infty$ or
$\omega\to 0$. Notice that there is still a strictly positive work done even
if $\omega\notin\cal I$, as long as there is dissipation on the point where
work is applied ($\gamma_{+}>0$) and $\omega$ is finite. We will see in
section 4 that this work flows directly into the right reservoir while the
current of mechanical energy through the system vanishes as $n\to\infty$. In
particular, it follows from (3.17) that

$\begin{split}\lim_{\omega\uparrow\omega_{0}}\overline{W}(\omega)&=\overline{W}(\omega_{0}):=\frac{F^{2}}{4{\gamma_{+}}},\\\
\lim_{\omega\downarrow\sqrt{\omega_{0}^{2}+4}}\overline{W}(\omega)&=\overline{W}(\sqrt{\omega_{0}^{2}+4})\\\
&=\frac{\gamma_{+}F^{2}(\omega_{0}^{2}+4)}{4}\cdot\frac{1+\gamma_{-}^{2}(\omega_{0}^{2}+4)}{1+(\gamma_{-}^{2}+\gamma_{+}^{2})(\omega_{0}^{2}+4)+\gamma_{-}^{2}\gamma_{+}^{2}(\omega_{0}^{2}+4)^{2}}.\end{split}$
(3.19)

This helps to understand the different scales on vertical lines in Figure 1
depending on the value of $\gamma_{-}$.

####  3.1.2. The case $n\to+\infty$ and $\omega$ is inside of the harmonic
chain spectrum

The computation of the $n\rightarrow\infty$ limit for the Green’s functions
when $\omega$ is inside the harmonic spectral interval ${\cal I}$ is more
complicated because there are singularities at the harmonic frequencies
$\omega_{j}$ and the distance between singularities is of order $1/n$.

Fix $\omega$ inside of ${\cal I}$. To describe the behavior of $W(\omega,n)$
near the selected frequency $\omega$ we introduce a function
$\overline{W}(r,u)$, see formula (B.14). This function is $1$-periodic in both
variables and satisfies
$W(\omega,n)=\overline{W}\big{(}r(\omega),(n+1)r(\omega))+o(1)$, as
$n\to+\infty$. The description of $W(\omega,n)$ in terms of the associated
family of Young measures is given in (B.15) below. The work $W(\omega,n)$ in
the limit, when $n$ is large, is plotted in Figure 2.

Figure 2. Behavior of the work functional. First row:
$(\gamma_{-},\gamma_{+})=(1,1)$. Second row:
$(\gamma_{-},\gamma_{+})=(1,1/10)$. Left column: work inside the harmonic
spectrum computed using limiting expression (3.4) for $n\rightarrow\infty$.
Black dotted curve represents $W(\omega,n)$ with $n=50$. Red dashed lines
stand for the limit of the harmonic spectrum. Blue, cyan and orange lines
indicate the harmonic frequencies $\omega=1.0478$, $1.41421$ and $2.101$,
respectively. Right column: diagrams of $\overline{W}\big{(}r,u\big{)}$,
$u\in[-1,1]$ around the harmonic frequencies $\omega=1.0478$ ($r=0.1$, blue),
$1.41421$ ($r=0.66$, cyan) and $2.101$ ($r=0.75$, orange). Note the larger
magnitude of the vertical scale in the second row.

####  3.1.3. The case of a general periodic force

Finally, we remark that in the general case of a $\theta$-periodic force of
the form

${\mathcal{F}}(t/\theta)=\sum_{\ell=1}^{+\infty}F_{\ell}\cos(\omega(\ell)t),\quad\mbox{where
}\omega(\ell):=\frac{2\pi\ell}{\theta}$ (3.20)

whose real valued Fourier coefficients satisfy $\sum_{\ell=1}^{+\infty}(\ell
F_{\ell})^{2}<+\infty$, the work performed by the force can be determined from
the formula:

$W(n)=\sum_{\ell=1}^{+\infty}\big{(}\omega(\ell)F_{\ell}\big{)}^{2}\frac{N(\omega(\ell),n)}{D(\omega(\ell),n)}.$
(3.21)

Therefore its behavior, as $n$ gets large, can be determined from the term by
term analysis of the series appearing on the right hand side of (3.21).

###  3.2. Energy

As in Section 3.1 we assume that the periodic force ${\cal F}(t)$ is given by
(3.3). The time average of the expectation of the total energy energy of the
chain $E(\omega,n)$ breaks up into the sum of thermal component $E_{\rm
th}(\omega,n)=\sum_{x\in{\mathbb{I}}_{n}}\langle\langle e_{x}^{\rm
th}\rangle\rangle$ and the mechanical one $E_{\rm
mech}(\omega,n)=\sum_{x\in{\mathbb{I}}_{n}}\langle\langle e_{x}^{\rm
mech}\rangle\rangle$, with $e_{x}^{\rm th}$ and $e_{x}^{\rm mech}$ defined in
(2.17) and (2.16), respectively.

Considering the behavior of the thermal energy functional, defined in (2.15),
it has been shown in [10], that in the case $\omega_{0}=0$ and
$\gamma_{-}=\gamma_{+}$ we have $\langle\langle e_{x}^{\rm
th}\rangle\rangle=\frac{1}{2}(T_{-}+T_{+})$ for all $x=1,\ldots,n-1$. If
$\omega_{0}>0$ and $\gamma_{-}=\gamma_{+}$, then [8, formulas (38) and (42)]
give

$\langle\langle e_{x}^{\rm
th}\rangle\rangle=\frac{1}{2}(T_{-}+T_{+})(1+o_{x}),\quad\mbox{where
}|o_{x}|\leqslant\frac{C}{g^{x\wedge(n+1-x)}}$

for some constants $C>0$, $g>1$ independent of $n$. As a result we have
$E^{\rm th}(\omega,n)\sim n$, as $n\to+\infty$.

####  3.2.1. Formula for the total mechanical energy functional for a single
mode oscillating force

In what follows we consider the behavior of the mechanical component of the
energy. Again, assume that the force is given by (3.3). It turns out, see
Section C of the Appendix, that the time average over the period of the
microscopic mechanical energy density equals

$\langle\langle e_{x}^{\rm
mech}\rangle\rangle=\frac{F^{2}}{2}\cdot\frac{M_{x}(\omega,n)}{D(\omega,n)},$
(3.22)

where $D(\omega,n)$ is given by (3.1) and

$M_{x}(\omega,n)=G^{1}_{x}(\omega,n)^{2}(\omega^{2}+\omega_{0}^{2})+(\nabla^{\star}G^{1}_{x})(\omega,n)^{2}+(2\omega\gamma_{-})^{2}\Big{[}{\cal
G}_{x}(\omega,n)^{2}+(\nabla^{\star}{\cal G}_{x})(\omega,n)^{2}\Big{]},$

with (see (3.7))

$G^{0}_{x}(\omega,n)=G_{0,x}(\omega,n)\quad\mbox{ and }\quad
G^{1}_{x}(\omega,n)=G_{x,n}(\omega,n).$ (3.23)

Using (3.7) we get

$\displaystyle
G^{s}_{x}(\omega,n)=\frac{1}{n+1}\cdot\frac{1}{\omega_{0}^{2}-\omega^{2}}+\frac{2}{n+1}\sum_{j=1}^{n}\frac{(-1)^{js}\cos(\frac{\pi
j}{2(n+1)})\cos(\frac{\pi
j(2x+1)}{2(n+1)})}{\omega_{j}^{2}-\omega^{2}}\quad,\quad s=0,1$

and (recall $G^{s}=G^{s}_{0}$, $s=0,1$)

${\cal
G}_{x}(\omega,n)=G^{0}(\omega,n)G^{1}_{x}(\omega,n)-G^{1}(\omega,n)G^{0}_{x}(\omega,n).$

The explicit formula for the total mechanical energy functional, obtained by
summing over all $x$ expression (3.22), is presented in (C.1) below.

####  3.2.2. Energy in the case $\omega$ lies outside harmonic chain spectrum

Analogously as in the case of the work functional the behavior $E_{\rm
mech}(\omega,n)$ depends on whether the force frequency belongs to the inside
or outside of the spectrum of the harmonic chain. If $\omega\not\in{\cal I}$
the asymptotics of $E_{\rm mech}(\omega,n)$, as $n\to+\infty$, can be obtained
by a Riemann sum approximation. Then,

$\lim_{n\to+\infty}E(\omega,n)=\frac{F^{2}\left\\{1+4\big{(}\gamma_{-}H\big{)}^{2}\right\\}}{4[1+4(\gamma_{-}^{2}+\gamma_{+}^{2})H^{2}+16(\gamma_{+}\gamma_{-}H^{2})^{2}]}\Big{[}K_{0}\left(\omega^{2}+\omega_{0}^{2}\right)+K_{1}\Big{]}.$
(3.24)

Here $H$ is given by (3.18) and

$\displaystyle K_{0}=\frac{\mathrm{d}H}{\mathrm{d}\omega^{2}},\qquad
K_{1}=\frac{\mathrm{d}}{\mathrm{d}\omega^{2}}\left(\Gamma_{0}(\omega)-\Gamma_{2}(\omega)\right),$

where $\Gamma_{x}(\omega)$ is the Green’s function of the lattice
${\mathbb{Z}}$ laplacian. It is given by

$\displaystyle\Gamma_{x}(\omega)=\left\\{\Big{[}\omega_{0}^{2}-\omega^{2}\Big{]}\Big{[}4+\omega_{0}^{2}-\omega^{2}\Big{]}\right\\}^{-1/2}$
$\displaystyle\times\left\\{1+\frac{1}{2}\Big{[}\omega_{0}^{2}-\omega^{2}\Big{]}+\frac{1}{2}\left\\{\Big{[}\omega_{0}^{2}-\omega^{2}\Big{]}\Big{[}4+\omega_{0}^{2}-\omega^{2}\Big{]}\right\\}^{1/2}\right\\}^{-|x|},\quad\mbox{when
$\omega_{0}^{2}>\omega^{2}$ }$

and

$\displaystyle\Gamma_{x}(\omega)=-\left\\{\Big{[}\omega^{2}-\omega_{0}^{2}\Big{]}\Big{[}\omega^{2}-\omega_{0}^{2}-4\Big{]}\right\\}^{-1/2}$
$\displaystyle\times\left\\{1-\frac{1}{2}\Big{[}\omega^{2}-\omega_{0}^{2}\Big{]}-\frac{1}{2}\left\\{\Big{[}\omega^{2}-\omega_{0}^{2}\Big{]}\Big{[}\omega^{2}-\omega_{0}^{2}-4\Big{]}\right\\}^{1/2}\right\\}^{-|x|},\quad\mbox{when
$\omega_{0}^{2}+4<\omega^{2}$.}$

Note that when $\gamma_{+}\to 0$, the formula (3.24) simplifies and we have

$\lim_{n\to+\infty}E(\omega,n)=\frac{F^{2}}{4}\Big{[}K_{0}\left(\omega^{2}+\omega_{0}^{2}\right)+K_{1}\Big{]}.$
(3.25)

####  3.2.3. The case when $\omega$ is inside of the harmonic chain spectrum

If, $\omega$ is inside of ${\cal I}$, the time average of $E_{\rm
mech}(\omega,n)$ is proportional to the size of the system. After
normalization we obtain, see Section C of the Appendix,

$\displaystyle\frac{1}{n}E(\omega,n)=\overline{E}\big{(}r,(n+1)r\big{)}+o(1)$

as $n\to+\infty$, where $\overline{E}\big{(}r,u\big{)}$ is $1$-periodic in the
first and $2$-periodic in the second variable. It is described by formulas
(C.5) and (C.6). Here $r$ is determined from $\omega$ by formula (3.1).

Behavior of the energy functional is illustrated in Fig. 3.

Figure 3. Behavior of the energy for $n=50$ with $\gamma_{\pm}=1$ (left) and
$\gamma_{+}=1$, $\gamma_{-}=1/10$ (right).

##  4\. Current of mechanical energy

The currents of the mechanical energy are given by

$\begin{split}\overline{j_{x,x+1}^{\text{mech}}}(t)&=-\overline{p}_{x}(t)\left(\overline{q}_{x+1}(t)-\overline{q}_{x}(t)\right),\qquad
x=0,\dots,n-1\\\
\overline{j_{-1,0}^{\text{mech}}}(t)&=-\gamma_{-}\overline{p}_{0}(t)^{2},\qquad\overline{j_{n,n+1}^{\text{mech}}}(t)=\gamma_{+}\overline{p}_{n}(t)^{2}-\mathcal{F}(t/\theta)\overline{p}_{n}(t).\end{split}$
(4.1)

They have all the same time average over the period:

$J^{\text{mech}}(n):=\langle\langle\overline{j_{x,x+1}^{\text{mech}}}\rangle\rangle.$
(4.2)

Note that $W^{-}(n):=-J^{\text{mech}}(n)$ is the amount of work that goes into
the left reservoir. Of course when $\gamma_{+}=0$ we have $W^{-}(n)=W(n)$. If
however $\gamma_{+}>0$, then some of the work, denoted by
$W^{+}(n)=W(n)-W^{-}(n)$, goes into the right reservoir.

We compute first $W^{-}(n)$, using $\overline{j_{-1,0}^{\text{mech}}}(t)$, as
it involves simpler formulas. From (A.5) we have

$\begin{split}&\overline{p}_{0}(t)=\text{Re}\left(i\omega{e^{i\omega
t}}\widetilde{q}_{0}(\omega)\right)\\\
&=\frac{FG^{1}(\omega,n)\left[\text{Re}(\widetilde{D}(\omega,n))\omega\cos(\omega
t)-\text{Im}(\widetilde{D}(\omega,n))\omega\sin(\omega
t)\right]}{\left|\widetilde{D}(\omega,n)\right|^{2}}\end{split}$ (4.3)

and, recalling that $\omega=\frac{2\pi}{\theta}$,

$\begin{split}J^{\text{mech}}(n)=-\frac{\gamma_{-}}{\theta}\int_{0}^{\theta}\overline{p}_{0}(t)^{2}\mathrm{d}t=-\gamma_{-}\Big{(}F\omega
G^{1}(\omega,n)\Big{)}^{2}\end{split}$ (4.4)

As a result, combining with (3.4), we get

$\displaystyle W^{-}(n)=\gamma_{-}\Big{(}F\omega
G^{1}(\omega,n)\Big{)}^{2},\qquad W^{+}(n)=W(n)-W^{-}(n).$ (4.5)

Notice that if $\omega\notin\cal I$, since
$G^{1}(\omega,n)\mathop{\longrightarrow}_{n\to\infty}0$, we have
$J^{\text{mech}}(n)\mathop{\longrightarrow}_{n\to\infty}0$. Comparing with
(3.17) we deduce that if $\omega\notin\cal I$, all the work goes to the right
thermostat as $n\to\infty$.

If $\omega\in(\omega_{0},\sqrt{\omega_{0}^{2}+4})$ , then
$W^{-}(n)=\overline{W}\vphantom{1}^{-}\Big{(}r(\omega),(n+1)r(\omega)\Big{)}+o(1),$
where the formula for $\overline{W}\vphantom{1}^{-}(r,u)$ can be obtained from
(4.5) by replacing $G^{1}(\omega,n)$ by the function $\overline{G}^{1}(r,u)$
defined in (B.13). We also have
$W^{+}(n)=\overline{W}\vphantom{1}^{+}\Big{(}r(\omega),(n+1)r(\omega)\Big{)}+o(1),$
where
$\overline{W}\vphantom{1}^{+}(r,u)=\overline{W}(r,u)-\overline{W}\vphantom{1}^{-}(r,u)$,
where $\overline{W}(r,u)$ is given by (B.14).

##  Appendix A Time harmonics of the position and momenta averages

Recall that ${\cal F}(t/\theta)={\rm Re}\Big{(}Fe^{i\omega t}\Big{)}$.
Consider the Fourier coefficients of the means of the positions and momenta

$\widetilde{p}_{x}=\frac{1}{\theta}\int_{0}^{\theta}e^{-i\omega
t}\overline{p}_{x}(t)\mathrm{d}t,\quad\widetilde{q}_{x}=\frac{1}{\theta}\int_{0}^{\theta}e^{-i\omega
t}\overline{q}_{x}(t)\mathrm{d}t,\quad x\in{\mathbb{I}}_{n}.$ (A.1)

We have $\overline{p}_{x}(t)={\rm Re}\Big{(}\widetilde{p}_{x}e^{i\omega
t}\Big{)}$ and $\overline{q}_{x}(t)={\rm Re}\Big{(}\widetilde{q}_{x}e^{i\omega
t}\Big{)}$.

From (2.4) and (2.5) we obtain $\widetilde{p}_{x}=i\omega\widetilde{q}_{x}$
and

$i\omega\widetilde{p}_{x}=\Big{(}\Delta_{\rm
N}-\omega_{0}^{2}-2i\gamma_{x}\omega\Big{)}\widetilde{q}_{x}+F\delta_{x,n},\;\quad
x\in{\mathbb{I}}_{n}.$ (A.2)

Here $\gamma_{x}=\gamma_{-}\delta_{0,x}+\gamma_{+}\delta_{n,x}$. Substituting
into (A.2) for $\widetilde{p}_{x}$ we get the equation

$0=\Big{(}\Delta_{x}+\omega^{2}-\omega_{0}^{2}-2i\omega\gamma_{x}\Big{)}\widetilde{q}_{x}+F\delta_{x,n},\;\quad
x\in{\mathbb{I}}_{n}.$ (A.3)

Hence, using the notation of (3.23), we can write

$\Big{(}F-2i\gamma_{+}\omega\widetilde{q}_{n}\Big{)}G^{1}_{x}(\omega,n)-2i\gamma_{+}\omega\widetilde{q}_{0}G^{0}_{x}(\omega,n)=\widetilde{q}_{x},\quad
x\in{\mathbb{I}}_{n}.$ (A.4)

For $x=0,n$ we get a closed system of $2$ equations for $\widetilde{q}_{0}$
and $\widetilde{q}_{n}$ that can be solved explicitly and we obtain

$\widetilde{q}_{0}=\frac{FG^{1}(\omega,n)}{\widetilde{D}(\omega,n)},\quad\widetilde{q}_{n}=\frac{F\widetilde{N}(\pm\omega,n)}{\widetilde{D}(\omega,n)},$
(A.5)

where, using the notation of (3.6), we have

$\displaystyle\widetilde{N}(\omega,n)=G^{0}(\omega,n)+2i\omega\gamma_{-}\left(G^{0}(\omega,n)^{2}-G^{1}(\omega,n)^{2}\right)$
$\displaystyle\widetilde{D}(\omega,n)=1-4\gamma_{+}\gamma_{-}\omega^{2}\left(G^{0}(\omega,n)^{2}-G^{1}(\omega,n)^{2}\right)+2i\omega(\gamma_{+}+\gamma_{-})G^{0}(\omega,n).$
(A.6)

Substituting back into (A.4) we conclude that

$\begin{split}&\widetilde{q}_{x}=F\Big{(}aG^{0}_{x}(\omega,n)+bG^{1}_{x}(\omega,n)\Big{)},\quad\mbox{where},\\\
&a=1-2i\omega\gamma_{+}\frac{\widetilde{N}(\omega,n)}{\widetilde{D}(\omega,n)},\quad
b=-2i\omega\gamma_{-}\frac{G^{1}(\omega,n)}{\widetilde{D}(\omega,n)}.\end{split}$
(A.7)

Using (2.18) and the fact that $\overline{p}_{n}(t)=-\omega{\rm
Im}\Big{(}\widetilde{q}_{n}e^{i\omega t}\Big{)}$ we obtain (3.4).

##  Appendix B Time average of work functional when $\omega$ is inside ${\cal
I}$ and $n\to+\infty$

We consider now $\omega\in(\omega_{0},\sqrt{\omega_{0}^{2}+4})$. We will
parametrize the spectrum using $r(\omega)\in(0,1)$, defined by (3.2), and we
study here the asymptotic behaviour of $\widetilde{W}(r,n)=W(\omega(r),n)$.
Similarly we define $\widetilde{G}^{s}(r,n),s=0,1$.

Denote $j(r)=[(n+1)r]$ (where $[x]$ denotes the integer part of $x$) and

$u(r)=(n+1)r-[(n+1)r]\in(0,1).$ (B.1)

Since we are choosing $\omega\neq\omega_{j}$, we have that $u(r)\in(0,1)$. .

To compute $\widetilde{G}^{0}(r,n)$ we start with extracting the singular term
at $\omega_{j}$. From (3.9) we get

${\widetilde{G}^{0}(r,n)}=\frac{1}{n+1}\cdot\frac{1}{\omega_{0}^{2}-\omega(r)^{2}}+I_{-}(r;0,j-1)+I_{+}(r;1,n-j),$
(B.2)

where

$I_{\pm}(r;m,k)=\frac{2}{n+1}\sum_{\ell=m}^{k}\frac{\cos^{2}(\frac{\pi(j(r)\pm\ell)}{2(n+1)})}{\omega\left(\frac{j(r)\pm\ell}{n+1}\right)^{2}-\omega(r)^{2}},$
(B.3)

For any $1\leqslant k_{0}\leqslant k$ we break $I_{\pm}(\omega;1,k)$ in two
terms: one with the first $k_{0}$ terms and the other with the remaining
$k-k_{0}$ ones. The idea is to assume that $k$ is of order $n$ and $k_{0}$ is
of order $n^{a}$, with $a\in(0,1)$, when $n\rightarrow\infty$. The first term
can be summed up explicitly and for the second we can use the Riemann sum
approximation, since we are far away from the singularity that occurs at
$\omega_{j}$. More precisely we can write
$I_{+}(r;0,k)=I_{+}(r;0,k_{0})+I_{+}(r;k_{0}+1,k)$. Using the formula

$\cot(\pi x)=\frac{1}{\pi
x}-\frac{2x}{\pi}\sum_{j=1}^{+\infty}\frac{1}{j^{2}-x^{2}},$ (B.4)

for $k_{0}\sim n^{a}$, $a<1$, and large $n$ we have

$\begin{split}I_{+}(r;0,k_{0})=\frac{1}{n+1}\sum_{\ell=0}^{k_{0}}\frac{\cos^{2}(\frac{\pi(j(r)+\ell)}{2(n+1)})}{\cos(\pi
r)-\cos\left(\pi\frac{j(r)+\ell}{n+1}\right)}\\\
=\frac{1}{n+1}\sum_{\ell=0}^{k_{0}}\frac{\cos^{2}(\frac{\pi(j(r)+\ell)}{2(n+1)})}{-\sin\left(\frac{j(r)+\ell}{n+1}\right)\frac{\pi}{n+1}\left((n+1)r-j(r)-\ell\right)}+o(1)\\\
=-\frac{1}{2\pi}\sum_{\ell=0}^{k_{0}}\frac{\cot\left(\frac{\pi(j(r)+\ell)}{2(n+1)}\right)}{\left(u(r)-\ell\right)}+o(1)=\frac{\cot(\pi
r/2)}{2\pi}\sum_{\ell=0}^{k_{0}}\frac{1}{\ell-u(r)}+o(1).\end{split}$ (B.5)

The sum in the last expression diverges, when $k_{0}\rightarrow\infty$.
However, in the expression (B.2) for $\widetilde{G}^{0}$ we have also

$I_{-}(r;1,k_{0})=-\frac{\cot(\pi
r/2)}{2\pi}\sum_{\ell=1}^{k_{0}}\frac{1}{\ell+u(r)}+o(1)$ (B.6)

and, as a result of the cancelation, the sum of them has a finite limit as
$k_{0}\rightarrow\infty$. It can be computed and the result is:

$\lim_{k_{0}\rightarrow\infty}\big{(}I_{+}(r;0,k_{0})+I_{-}(r;1,k_{0})\big{)}=-\frac{1}{2}\cot(\pi
r/2)\cot(\pi u(r)).$ (B.7)

Now we compute the remaining expressions $I_{\pm}(r;k_{0}+1,k)$ by using the
Riemann sum approximation:

$\begin{split}&I_{\pm}(r;k_{0}+1,k)=\frac{2}{n+1}\sum_{\ell=k_{0}+1}^{k}\frac{\cos^{2}(\frac{\pi(j(r)\pm\ell)}{2(n+1)})}{\omega\left(\frac{j(r)\pm\ell}{n+1}\right)^{2}-\omega^{2}(r)}\\\
&=2\int_{v_{0}}^{v}\frac{\cos^{2}\left(\frac{\pi}{2}(r\pm\overline{v})\right)}{\omega(r\pm\overline{v})^{2}-\omega(r)^{2}}d\overline{v}+o(1)=\frac{1}{2}\int_{v_{0}}^{v}\,\frac{\big{[}1+\cos\big{(}\pi(r\pm\overline{v})\big{)}\big{]}}{\cos(\pi
r)-\cos(\pi(r\pm\overline{v}))}dv+o(1),\end{split}$ (B.8)

where $v_{0}=k_{0}/(n+1)$, $v=k/(n+1)$. The last integral has a logarithmic
singularity when $v_{0}\rightarrow 0$ (i.e. $k_{0}\ll n$). Nevertheless, when
putting together the two terms, we obtain the principal value of the integral
at the singular point and, as a result,

$\displaystyle
I(r)=\lim_{n\to+\infty}\Big{(}I_{+}(r;k_{0},j)+I_{-}(r;k_{0},n-j)\Big{)}=\frac{1}{2}{\rm
p.v.}\int_{-r}^{1-r}\,\frac{\big{[}1+\cos\big{(}\pi(r+\overline{v})\big{)}\big{]}\mathrm{d}\overline{v}}{\cos(\pi
r)-\cos(\pi(r+\overline{v}))}$ $\displaystyle=-\frac{1}{2}+\cos^{2}\big{(}\pi
r/2\big{)}{\rm p.v.}\int_{0}^{1}\,\frac{\mathrm{d}\overline{v}}{\cos(\pi
r)-\cos(\pi\overline{v})}$ (B.9)

(we recall that $k_{0}\simeq n^{a}$ with $a\in(0,1)$, and $j=[r(n+1)]$). Using
[6, formula 2.551.3, p. 171]

$\int\,\frac{\mathrm{d}\overline{v}}{\cos(\pi
r)-\cos(\pi\overline{v})}=\frac{1}{\pi\sin(\pi r)}\log\Big{|}\frac{(1+\sin(\pi
r))/\cos(\pi
r)-\tan\left(\frac{\pi}{2}\big{(}\overline{v}+\frac{1}{2}\big{)}\right)}{(1-\sin(\pi
r))/\cos(\pi
r)-\tan\left(\frac{\pi}{2}\big{(}\overline{v}+\frac{1}{2}\big{)}\right)}\Big{|}.$

we conclude that that the principal value of the integral on the utmost right
hand side of (B) equals null. Hence $I(r)=-1/2$.

Finally, putting together (B.7) and (B) we find:

$\displaystyle{\widetilde{G}^{0}(r,n)=\overline{G}^{0}\big{(}r,(n+1)r\big{)}}+o(1),\quad\mbox{where}$
$\displaystyle\overline{G}^{0}(r,u)=-\frac{1}{2}\Big{(}\cot(\pi r/2)\cot(\pi
u)+1\Big{)}.$ (B.10)

We will consider $u\in{\mathbb{R}}$ and extend periodically the function
$\overline{G}^{0}(r,u)$.

We compute $\overline{G}^{1}(r)=\lim_{n\to+\infty}\widetilde{G}^{1}(r,n)$ by
using formula:

$\frac{1}{2}\Big{(}\widetilde{G}^{0}(r,n)+\widetilde{G}^{1}(r,n)\Big{)}=\frac{1}{n+1}\cdot\frac{1}{\omega_{0}^{2}-\omega(r)^{2}}+\frac{2}{n+1}\sum_{\ell=1}^{n/2}\frac{\cos^{2}(\frac{\pi\ell}{n+1})}{4\sin^{2}(\frac{\pi\ell}{n+1})+\omega_{0}^{2}-\omega(r)^{2}}.$
(B.11)

That is a very similar expression to the original one for
$\widetilde{G}^{0}(r,n)$, see (3.9) with $s=0$, but with factors $2$ not
present in the denominators of fractions appearing in the infinite sum.
Following analogous arguments to the ones used before we find, in the limit
$n\rightarrow+\infty$

$\begin{split}&\frac{1}{2}\Big{(}\widetilde{G}^{0}(r,n)+\widetilde{G}^{1}(r,n)\Big{)}=\overline{H}\Big{(}r,(n+1)r\Big{)}+o(1),\quad\mbox{where}\\\
&\overline{H}(r,u)=-\frac{1}{4}\Big{(}\cot(\pi r/2)\cot(\pi
u/2)+1\Big{)}.\end{split}$ (B.12)

Therefore we get
$\widetilde{G}^{1}(r,n)=\overline{G}^{1}\big{(}r,(n+1)r\big{)}+o(1)$, as
$n\to+\infty$, where

$\overline{G}^{1}(r,u)=(-1)^{[u]+1}\frac{\cot(\pi r/2)}{2\sin(\pi u)}.$ (B.13)

We have shown therefore that

$\begin{split}&\widetilde{W}(r,n)=\overline{W}(r,(n+1)r)+o(1),\quad\mbox{where}\\\
&{\overline{W}(r,u)=F^{2}\omega(r)^{2}\frac{\overline{N}(r,u)}{\overline{D}(r,u)}}.\end{split}$
(B.14)

The functions $\overline{N}(r,u)$ and $\overline{D}(r,u)$ are given by
analogues of (3.1), with $G^{s}(\omega,n)$ replaced by
$\overline{G}^{s}(r,u)$, respectively for $s=0,1$. As in (3.8) we get

$0\leqslant\overline{W}(r,u)\leqslant\frac{F^{2}}{4}\left(\frac{1}{\gamma_{-}}+\frac{1}{\gamma_{+}}\right).$

Equality (B.14) can be used to find the family of Young measures that is
associated with $\widetilde{W}(r,n)$. This is defined, see e.g. [5, Section
1.E.3, p. 16], as a family of measures $\mu(r,\mathrm{d}v)$,
$v\in{\mathbb{R}}$ such that for any test function $\varphi\in L^{1}[0,1]$ and
a bounded continuous function $\Phi\in C_{b}({\mathbb{R}})$

$\displaystyle\lim_{n\to+\infty}\int_{0}^{1}\Phi\Big{(}\widetilde{W}(r,n)\Big{)}\varphi(r)\mathrm{d}r=\int_{0}^{1}\overline{\Phi}(r)\varphi(r)\mathrm{d}r,$
(B.15)

with $\overline{\Phi}(r):=\int_{\mathbb{R}}\Phi(v)\mu(r,\mathrm{d}v).$ Thanks
to (B.14) we conclude that the probability measures $\mu(r,\mathrm{d}v)$
obtained by transporting the Lebesgue measure ${\rm m}(\mathrm{d}u)$ on
$(0,1)$ by the mapping $u\mapsto\overline{W}\big{(}r,u\big{)}$ constitute the
family of Young measures associated with the sequence $\widetilde{W}(r,n)$. We
have $\mu(r,A)={\rm m}\Big{[}u:\overline{W}\big{(}r,u\big{)}\in A\Big{]}$ for
any Borel measurable subset $A$ of ${\mathbb{R}}$. Since
$\overline{W}\big{(}r,u\big{)}$ is bounded, piecewise $C^{1}$-smooth and
$\\{u:\overline{W}\big{(}r,u\big{)}=v\\}$ is finite for each $r$, the Young
measures $\mu(r,\mathrm{d}v)$, have compactly supported densities. Using the
frequency domain in the description of the Young measures, we conclude from
(B.15) that

$\displaystyle\lim_{n\to+\infty}\int_{\cal
I}\Phi\Big{(}W(\omega,n)\Big{)}\varphi(\omega)\mathrm{d}\omega=\int_{\cal
I}\overline{\Phi}(\omega)\varphi(\omega)\mathrm{d}\omega,$ (B.16)

for any function $\varphi\in C_{b}({\cal I})$. Here
$\overline{\Phi}(\omega):=\int_{\mathbb{R}}\Phi(v)\mu(r(\omega),\mathrm{d}v).$
We have

$\displaystyle\lim_{\omega\downarrow\omega_{0}}\mu(r(\omega),\mathrm{d}v)=\delta_{\overline{W}(\omega_{0})}(\mathrm{d}v),$
$\displaystyle\lim_{\omega\uparrow\sqrt{\omega_{0}^{2}+4}}\mu(r(\omega),\mathrm{d}v)=\delta_{\overline{W}(\sqrt{\omega_{0}^{2}+4})}(\mathrm{d}v),$

where, as we recall $\overline{W}(\omega_{0})$ and
$\overline{W}(\sqrt{\omega_{0}^{2}+4})$ are given in (3.19). The limit holds
in the sense of the weak convergence of measures.

##  Appendix C Time average of energy in case $\mathbf{\omega}$ is inside of
${\cal I}$

Formula (3.22) is a direct consequence of (2.16) and formula
(LABEL:021205-21e). Summing over all $x$ we conclude that

$\displaystyle E_{\rm mech}(\omega,n)$
$\displaystyle=F^{2}\Biggl{\\{}(|a|^{2}+|b|^{2})\biggl{[}(\omega^{2}+\omega_{0}^{2})I_{0}(\omega,n)+J_{0}(\omega,n)\biggl{]}$
(C.1) $\displaystyle+{\rm
Re}\big{(}a^{*}b)\biggl{[}(\omega^{2}+\omega_{0}^{2})I_{1}(\omega,n)+J_{1}(\omega,n)\biggl{]}\Biggl{\\}},$
(C.2)

where

$\displaystyle
I_{s}(\omega,n)=\frac{\mathrm{d}}{\mathrm{d}\omega^{2}}G_{s}(\omega,n),\quad
J_{s}(\omega,n)=\frac{\mathrm{d}}{\mathrm{d}\omega^{2}}R_{s}(\omega,n)$ (C.3)

and

$R_{s}=\frac{2}{n+1}\sum_{j=1}^{n}\frac{(-1)^{js}\sin^{2}(\frac{\pi
j}{(n+1)})}{4\sin^{2}(\frac{\pi j}{2(n+1)})+\omega_{0}^{2}-\omega^{2}},\quad
s=0,1.$ (C.4)

The remaining terms have been defined in Section A.

The four functions: $I_{s}(\omega,n)$ and $J_{s}(\omega,n)$, $s=0,1$,
appearing in (C.3) diverge, as $n\rightarrow\infty$, for $\omega$ inside of
${\cal I}$. Computations involving these functions use the same technique as
in the case of the asymptotics of the work functional considered in Section B
of the Appendix. We obtain

$\displaystyle e(\omega,n)=\frac{E_{\rm
mech}(\omega,n)}{n}=\overline{e}\big{(}r,(n+1)r\big{)}+o(1),\quad\mbox{where}$
(C.5)
$\displaystyle\overline{e}\big{(}r,u\big{)},=F^{2}\left[(\omega^{2}+\omega_{0}^{2})\left[(|\overline{a}|^{2}+|\overline{b}|^{2})I_{0}\big{(}r,u\big{)}+(\overline{a}^{*}\overline{b}+\overline{a}\overline{b}^{*})I_{1}\big{(}r,u\big{)}\right]\right.$
$\displaystyle+\left.(|\overline{a}|^{2}+|\overline{b}|^{2})J_{0}\big{(}u\big{)}+(\overline{a}^{*}\overline{b}+\overline{a}\overline{b}^{*})J_{1}\big{(}u\big{)}\right]+o(1),$

with the formulas for terms $\overline{a}$ and $\overline{b}$ given by
analogues of (LABEL:021205-21e), where the Green’s functions $G^{s}(\omega,n)$
are replaced by $\overline{G}^{s}(r,u)$, $s=0,1$, defined in (B.10) and
(B.13). Here $r$ is determined from $\omega$ by eqt. (3.1) and

$\displaystyle I_{s}(r,u)=\frac{(-1)^{s[u]}}{8}\cdot\frac{\cos(s\pi
u)}{\sin(\pi r/2)\sin^{2}(\pi u)},$ $\displaystyle
J_{s}(u)=\frac{(-1)^{s[u]}}{2}\cdot\frac{\cos(s\pi u)}{\sin^{2}(\pi
u)},\quad\quad s=0,1.$ (C.6)

Figure 4 illustrates the behaviour of $e(\omega)$.

Figure 4. Behavior of the energy per-oscillator. First row:
$(\gamma_{-},\gamma_{+})=(1,1)$. Second row:
$(\gamma_{-},\gamma_{+})=(1,1/10)$. Left column: energy computed with the
limiting expressions for $n\rightarrow\infty$. The oscillating part is
obtained directly using the Green’s function expressions with $n=40$ and black
dotted curves inside the harmonic spectra zone are computed used the energy
expression for $n\rightarrow\infty$. Red dashed lines define the limits of the
harmonic spectra. Blue, cyan and orange lines indicate the harmonic
frequencies $\omega=1.0478$, $1.41421$ and $2.101$. Right column: Scaled
energy $\overline{e}(r,u)$ around the harmonic frequencies $\omega=1.0478$
(blue), $1.41421$ (cyan) and $2.101$ (orange).

## Declarations

### Data Availability

Data sharing is not applicable to this article as no datasets were generated
or analysed during the current study.

### Conflict of interest

In addition, the authors have no conflicts of interest to declare that are
relevant to the content of this article.

## References

  * [1] C. Bernardin, S. Olla, Transport Properties of a Chain of Anharmonic Oscillators with Random Flip of Velocities, J. Stat Phys (2011) 145:1224-1255 DOI 10.1007/s10955-011-0385-6
  * [2] F. Bonetto, J. L. Lebowitz and J. Lukkarinen Fourier’s Law for a Harmonic Crystal with Self-Consistent Stochastic Reservoirs. J. of Stat. Physics, Vol. 116, 2004
  * [3] J.L. Lebowitz, P.G. Bergmann, _Irreversible Gibbsian Ensembles_ , Annals of Physics, Vol. 1, N.1, 1-23, 1957\. https://doi.org/10.1016/0003-4916(57)90002-7
  * [4] P. Carmona, _Existence and uniqueness of an invariant measure for a chain of oscillators in contact with two heat baths_ , Stochastic Processes and their Applications 117, (2007), no. 8, 1076–1092.
  * [5] Evans, L. C., Weak convergence methods for nonlinear partial differential equations. CBMS Regional Conf. Ser. in Math., 74, by the American Mathematical Society, Providence, RI, 1990.
  * [6] I.S. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series, and Products, Seventh Edition, Academic Press-Elsevier 2007.
  * [7] T. Komorowski, J.L. Lebowitz, S. Olla, _Heat flow in a periodically forced, thermostatted chain_ , Comm.Math.Phys., 400, 2181–2225 (2023), https://doi.org/10.1007/s00220- 023-04654-4
  * [8] H. Nakazawa On the Lattice Thermal Conduction, Supplement of the Progress of Theoretical Physics, No. 45, (1970), pp 231-262.
  * [9] Abhinav Prem, Vir B. Bulchandani, and S. L. Sondhi, _Dynamics and transport in the boundary-driven dissipative Klein-Gordon chain_ , Phys. Rev. B 107, 104304, 2023
  * [10] Rieder, Z., Lebowitz, J.L., Lieb, E.: _Properties of harmonic crystal in a stationary non-equilibrium state_. J. Math. Phys. 8, 1073–1078 (1967).
  * [11] Yaghoubi M., Foulaadvand M. E., Bérut, A., Łuczka, J., Energetics of a driven Brownian harmonic oscillator, J. Stat. Mech. (2017) 113206


††institutetext: Department of Physics, Hanyang University, Seoul 04763, Korea

# Classification of the Mott gap

Debabrata Ghorai  Taewon Yuk  Young-Kwon Han  Sang-Jin Sin
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we demonstrate the classification of the gap in a holographic
setup by studying the density of states. A gap can be classified into order
gap and Mott gap depending on the presence of the order due to the symmetry
breaking or not. A Mott insulating gap appears in the fermion spectrum due to
the strong Coulomb interaction between the electrons. We then classify all
Mott gaps as well as order gaps in one-flavor and two-flavor fermions. We also
identified possible non-minimal interactions that may produce a flatband.

##  1 Introduction

The understanding of strongly correlated systems still remains mysterious
because of the lack of tools to calculate strange behaviors of materials
RevModPhys.78.17  although many useful methods were developed including the
dynamical mean-field theory (DMFT) RevModPhys.68.13  to get the modification
from the many-body effect. One of the celebrated features of strongly
interacting systems is the Mott gap, which can not be described by the gap
induced by the order and symmetry breaking. It is induced by the electron-
electron on site interaction. Mott physics may be essential in understanding
the physics of high-temperature superconductivity since the latter is
considered as a doped Mott insulator PhilipNature . Its physics can be best
represented by the Hubbard Hamiltonian RevModPhys.40.677  that involves the
competition between hopping and on-site repulsion. Although the Hubbard model
is not solved completely in two and higher dimensions, the Mott gap is known
to be well explained by DMFT RevModPhys.68.13 . The gap appears when onsite
repulsion is dominating the hopping.

Gauge/gravity duality Maldacena:1997re ; Gubser:1998bc ; Witten:1998qj
provides an alternative tool to study a strongly coupled system with the help
of weakly interacting theory of gravity in one higher spatial dimension HKMS ;
Seo:2016vks ; Oh:2021xbe ; Song:2019asj ; Oh:2020cym ; Oh:2018wfn ;
Seo:2017oyh ; Seo:2017yux . The application of this duality in condensed
matter physics has been widely used in the past decade. Since gravity models
map strongly coupled boundary field theory, one obvious inquiry is to find a
holographic setup. Since DMFT is used to explain the Mott gap, another inquiry
is to find an alternative mean field theory for a strongly coupled system.
Recently, an alternative mean field theory for a strongly coupled system has
been proposed using a holographic approach in Sukrakarn:2023ncp . The main
task is to find a gap-like feature without any symmetry breaking in the
holographic setup. This is only possible when we consider non-minimal coupling
between the gauge and fermionic fields. Because such holographic theory aims
at the general theory for strongly interacting system, it would be interesting
to exlain the Mott physics in holographic mean field theory and compare with
the DMFT. Along this line of thinking, the Pauli interaction term has been
utilized in the holographic literature PhysRevLett.106.091602 ;
PhysRevD.83.046012 ; PhysRevD.90.044022 ; philip3 ; Seo_2018 . However, the
spectral function of the holographic fermion with Pauli interaction has higly
asymmetric gap in the sense that uppper side of the Fermi sea has gap but
lower side of FS is touched by a spectral peak line, which is clear in
spectral function as well as the density of states (DoS). From the DoS in the
presence of the dipole coupling, the Mott gap in this theory is ‘soft’ as well
as asymmetric. On the other hand, the typical DMFT calculation shows that the
Mott gap is symmetric. This motivates us to consider other non-minimal
interactions between the gauge field and fermion. Finding a symmetric Mott gap
in the DoS analysis of the holographic fermion is the main motivation of this
paper. In doing this we can classify Mott gaps as well as the ordered gaps.

In this paper, we first reproduced all the spectral functions for the dipole
interaction. We then analyzed the DoS corresponding to the spectral functions,
which has not been done in the literature. We then propose our holographic
setup with different interactions. From the holographic mean field theory
Sukrakarn:2023ncp , we know that only (pseudo) scalar type Yukawa interaction
can give a proper gap, while we also know that there should not be any order
parameter field involved in the description of the Mott gap. Then the only
field we can use is the gauge field describing the density effect and the only
way to form a gauge invariant scalar out of the gauge field is
$F_{\mu\nu}F^{\mu\nu}$ so that we should try a few version of the
$F^{2}{\bar{\psi}}\Gamma^{scalar}\psi$ term. This is the idea of the paper and
as a consequence of adding such interaction, we get symmetric Mott gap from
some of them. Then, the dependence on temperature, chemical potential,
coupling constant, and the effect of fermion mass have been discussed in
detail. From the boundary point of view, the possible gamma scalar are
$\Gamma^{scalar}=\mathbb{I}_{4},\Gamma^{z},i\Gamma^{5},i\Gamma^{5z}$. This
analysis is extended to the two-flavor fermion case. For completeness, we
revisit all ordered gaps in holographic setups. Finally, we classify all
interactions in the holographic set from the gap point of view.

This paper is organized as follows. In section 2, we have revisited the
holographic setup with the dipole interaction and proposed our setup with
different interactions. The DoS analysis for Mott gaps is presented in section
3. In section 4, we present the classifications of all interactions in terms
of Mott gap, ordered gap and flatband. We summarize our findings in section 5.

##  2 Basic setup

###  2.1 Pauli interaction term for Mott gap

Before discussing our proposal for holographic Mott gap, we would like to
revisit the previous model for Mott gap PhysRevLett.106.091602 . The
holographic Mott gap model is based on the Pauli or dipole interaction term.
The Lagrangian is given by

$\displaystyle\mathcal{L}_{f}=i\bar{\psi}\left(\Gamma^{\mu}D_{\mu}-m_{f}-i\frac{p}{2}F_{\mu\nu}\Gamma^{\mu\nu}\right)\psi$
(1)

The above interaction term can not be mapped with Hubbard interaction term for
Mott gap 27. The background geometry was considered in PhysRevLett.106.091602
as follows:

$\displaystyle
ds^{2}=\frac{1}{z^{2}}\left[-f(z)dt^{2}+\frac{dz^{2}}{f(z)}+dx^{2}+dy^{2}\right]~{}~{}\text{with}~{}~{}f(z)=1-Mz^{3}+Q^{2}z^{4}~{}~{}.$
(2)

The above metric represents Reissner-Nordström (RN) AdS black hole geometry.
The zero temperature limit of the boundary theory implies to the black hole
extremal limit. In the extremal limit, the mass and the charge of the black
hole are fixed with values $M=4$ and $Q=\sqrt{3}$. For the fixed values of the
other parameters ($p=\pm 4.5,m_{f}=0,q_{f}=1,r_{h}=1,\mu=\sqrt{3}$), we have
reproduced all previous spectral function. In the previous investigation
PhysRevLett.106.091602 , the energy density curve
$``A(\omega)~{}\text{vs}~{}\omega"$ for fixed $k$ has been used to descibe the
gap feature in the fermionic spectral function. For $p=-4.5$, the spectral
function (figure 1(a)) exhibits a gap feature while Fermi level $(\omega=0)$
appears to be touching the valence band. On the other hand, for $p=+4.5$, the
spectral function flips, suggesting a gap in the negative energy region,
although the Fermi level touches the conduction band Oh:2020cym . Using the
definition of density of state $(DoS=\frac{1}{2\pi}\int A(k,\omega)kdk)$, we
have shown the DoS corresponding to their spectral function plot. The
identification of the Mott gap in the DoS is not clear, although the spectral
function for $p=-4.5$ shows a gap-like feature in figure 1(a). This Mott gap
is soft and asymmetrical in the DoS. This motivates us to consider other
possible non-minimal interactions between the gauge field and fermion.

###  2.2 Our proposal

We propose the following fermionic action

$\displaystyle S_{\psi}=\int
d^{4}x\sqrt{-g}\left[i\bar{\psi}(\Gamma^{\mu}D_{\mu}-m_{f}-\eta
F^{2}\Gamma^{scalar})\psi\right]~{}~{}\text{and}~{}~{}S_{bdy}=i\int
d^{3}x\sqrt{-h}\bar{\psi}\psi~{}.$ (3)

Here, $\eta$ represents the coupling constant and
$F^{2}=F_{\mu\nu}F^{\mu\nu}$. The spinor’s covariant derivative is denoted by
$D_{\mu}=\partial_{\mu}+\frac{1}{4}\omega_{\mu\bar{\alpha\beta}}\Gamma^{\bar{\alpha\beta}}-iq_{f}A_{\mu}$.
The possible gamma matrices are
$\Gamma^{scalar}=\mathbb{I}_{4},i\Gamma^{5},\Gamma^{z},i\Gamma^{5z}$.
Comparing with the interaction term in Hubbard Hamiltonian (see appendix A and
eq.(23)), we can identify $\Gamma^{scalar}=\mathbb{I}_{4}$, which will show a
symmetric Mott gap in the DoS. The bulk gamma matrices for this study are the
following:

$\displaystyle\Gamma^{\underline{t}}=\sigma_{1}\otimes
i\sigma_{2},\quad\Gamma^{\underline{x}}=\sigma_{1}\otimes\sigma_{1},\quad\Gamma^{\underline{y}}=\sigma_{1}\otimes\sigma_{3},\quad\Gamma^{\underline{z}}=\sigma_{3}\otimes\sigma_{0},$
(4)

where underline indices represent tangent space indices. We obtain the Dirac
equation

$\displaystyle(\Gamma^{\mu}D_{\mu}-m_{f}-\eta F^{2}\Gamma^{scalar})\psi=0~{}.$
(5)

To simplify the analysis, we express the fermionic field as follows:

$\displaystyle\psi(t,x,y,z)=\frac{1}{(-gg^{zz})^{1/4}}e^{-i\omega
t+ik_{x}x+ik_{y}y}\Psi(z)~{}~{}.$ (6)

This form allows us to eliminate the spin connection term in the spinor
equation of motion. By substituting the above expression into the Dirac
equations, we derive the following expressions:

$\displaystyle\
\left[\Gamma^{\underline{z}}\partial_{z}-i\sqrt{\frac{g^{tt}}{g^{zz}}}(\omega+q_{f}A_{t})\Gamma^{\underline{t}}-\frac{m_{f}}{\sqrt{g^{zz}}}\right]\Psi$
$\displaystyle\
+i\left(\sqrt{\frac{g^{xx}}{g^{zz}}}k_{x}\Gamma^{\underline{x}}+\sqrt{\frac{g^{yy}}{g^{zz}}}k_{y}\Gamma^{\underline{y}}\right)\Psi-\frac{\eta
F^{2}}{\sqrt{g^{zz}}}\Gamma^{scalar}\Psi=0.$ (7)

##  3 Mott gaps in DoS

In this section, we will calculate the density of states by solving the Dirac
equation. To solve the Dirac equation, we express the four-component spinor as
$\Psi(z)=\left(\Psi_{+}(z),\Psi_{-}(z)\right)^{T}$ where
$\Psi_{\pm}=\left(\Psi_{\pm 1},\Psi_{\pm 2}\right)$. First, we focus on the
gamma matrix $\Gamma^{scalar}=\mathbb{I}_{4}$, which can be mapped to the
interaction term in Hubbard Hamiltonian. The Dirac equation becomes

$\displaystyle\left[\partial_{z}\mp\frac{m_{f}}{\sqrt{g^{zz}}}\right]\Psi_{\pm}(z)=\pm\left[iK_{j}\gamma^{j}\Psi_{\mp}(z)\right]\pm\frac{\eta
F^{2}}{\sqrt{g^{zz}}}\Psi_{\pm}$ (8)

where
$K_{j}=\left(\sqrt{\frac{g^{tt}}{g^{zz}}}(\omega+q_{f}A_{t}),-\sqrt{\frac{g^{xx}}{g^{zz}}}k_{x},-\sqrt{\frac{g^{yy}}{g^{zz}}}k_{y}\right)$,
$\gamma^{j}=\left(i\sigma_{2},\sigma_{1},\sigma_{3}\right)$. In the asymptotic
limit as $z\rightarrow 0$, we consider $g^{\mu\nu}\rightarrow
z^{2}\eta^{\mu\nu}$, where $\eta^{\mu\nu}$ is the Minkowski metric. In the
asymptotic limit, the source and condensation are given for
$|m_{f}|<\frac{1}{2}$ by Ghorai:2023wpu

$\displaystyle\Psi_{+}(z)\overset{z\rightarrow
0}{=}\mathbf{A}z^{m_{f}},\quad\Psi_{-}(z)\overset{z\rightarrow
0}{=}\mathbf{D}z^{-{m_{f}}}~{}~{}~{}.$ (9)

Following the same procedure in Ghorai:2023vuo , we can write down the
boundary action in the following form

$\displaystyle S_{bdy}=\int
d^{3}x\left[\Psi_{-}^{\dagger}(z)\tilde{\Gamma}\Psi_{+}(z)+\Psi_{+}^{\dagger}(z)\tilde{\Gamma}\Psi_{-}(z)\right]$
(10)

where the boundary gamma matrix $\tilde{\Gamma}=-\sigma_{2}$. Recasting the
Dirac equation, the flow equation for bulk Green’s function $\mathbb{G}(z)$
has been derived in Appendix B. Using the horizon value of the $\mathbb{G}(z)$
and solving the flow equation, we can numerically calculate the bulk Green’s
function $\mathbb{G}(z)$. From this bulk Green’s function, the retarded
Green’s function $\mathbb{G}_{R}$ is obtained using the following relation

$\displaystyle\mathbb{G}_{R}=\lim_{z\to 0}U(z)\mathbb{G}(z)U(z)$ (11)

where $U(z)=\text{diag}(z^{m_{f}},z^{m_{f}})$. The spectral function is
defined as

$\displaystyle A(k_{x},k_{y},\omega)=\text{Tr}(\text{Im}(\mathbb{G}_{R}))~{}.$
(12)

From this, the density of state (DoS) is defined in the following way:

$\displaystyle\text{DoS}=\frac{1}{(2\pi)^{2}}\int_{k_{cut}}dk_{x}dk_{y}A(k_{x},k_{y},\omega)=\frac{1}{2\pi}\int_{k_{cut}}A(k_{x},k_{y},\omega)kdk~{}$
(13)

where $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$ and ${k_{cut}}$ is the momentum cutoff
region in which we are counting the degree of freedom of the system. We have
considered the gauge field ansatz $A=A_{t}(z)dt$ to compute the spectral
function, which behaves as $A_{t}=\mu(1-z/z_{h})$, where $z_{h}$ is the
horizon and $\mu$ is the chemical potential of the boundary theory. Using this
gauge field solution, we obtain the RN-AdS black hole as background spacetime
which has form

$\displaystyle
ds^{2}=\frac{1}{z^{2}}\left[-f(z)dt^{2}+\frac{dz^{2}}{f(z)}+dx^{2}+dy^{2}\right]~{}~{}\text{with}~{}~{}f(z)=1-\frac{z^{3}}{z^{3}_{h}}+\frac{\mu^{2}}{2}\left(\frac{z^{4}}{z^{2}_{h}}-\frac{z^{3}}{z_{h}}\right)~{}.$
(14)

The spectral function with DoS at finite temperature $T=0.025\mu$ for
$q_{f}=0,1$ is shown in figure 1. A clear Mott gap feature is observed in the
spectral function as well as in the density of states.

(a) $F_{\mu\nu}\Gamma^{\mu\nu}$ at $T=0,q_{f}=1$

(b) $F^{2}$ at $T=0.025\mu,q_{f}=0$

(c) $F^{2}$ at $T=0.025\mu,q_{f}=1$

Figure 1: Spectral function and density of state: (a) For dipole interaction
at zero temperature, Mott gap is soft and asymmetric. (b) For density
interaction ($\eta F^{2}\bar{\psi}\psi$) with $q_{f}=0$, Mott gap is hard and
symmetric. (c) With minimal interaction included also, the Mott gap is hard
and asymmetric due to the shifted Fermi surface.

We now explore all possible interactions that can generate a gap without any
symmetry breaking. As we mentioned earlier, we have also examined other three
possible cases $\Gamma^{scalar}=i\Gamma^{5},\Gamma^{z},i\Gamma^{5z}$. From the
spectral function and DoS, we have found that $F^{2}\bar{\psi}i\Gamma^{5}\psi$
generates Mott gap while the other two interactions do not show any gap
feature in the spectral function (table 1). One noticing point is that the gap
size for the $F^{2}\bar{\psi}\psi$ interaction is larger than the gap size for
the $F^{2}\bar{\psi}i\Gamma^{5}\psi$ interaction. The comparison of the DoS
with the same parameter values for different interactions at $T=0.025\mu$ is
shown in figure 2. This is clear evidence that $F^{2}\bar{\psi}\psi$
interaction is more suitable for describing the Mott gap in holography. The
holographic flatband can also be realized using non-minimal gauge coupling
with fermion ($F_{\mu\nu}\Gamma^{\mu\nu}\Gamma^{5}$). The fermion’s finite
charge bends and shifts the flatband.

Gauge | $-i\eta\bar{\psi}L_{int}\psi$
---|---
Field | Gapless | Gap | Flatband | Effect of $q_{f}$
$A_{t}(z)dt$ |  $\Gamma^{z}F^{2}$, $i\Gamma^{5z}F^{2}$ |  $F^{2}(\eta>0),iF^{2}\Gamma^{5}$, $-iF_{\mu\nu}\Gamma^{\mu\nu}$ | $F_{\mu\nu}\Gamma^{\mu\nu}\Gamma^{5}$ | Shifting & Bending
Table 1: Mott gaps from different gauge field interactions Figure 2: DoS for
Mott gaps with different interactions at the same paratemeters.

We would like to investigate the effect of temperature, coupling strength,
fermion mass in the spectral function as well as in the density of state. We
know that the effect of the charge is to shift the position of the Fermi
surface. In the presence of the charge $q_{f}=1$, we have also found the same
feature here (figure 1(c)). The gap size $\Delta_{M}$ is measured from the
DoS, where DoS$\leq 0.001$ region is considered as the gap region. It depends
on the chemical potential and the temperature. Since the Mott gap is generated
due to non-minimal gauge field interaction, the gap size is proportional to
the chemical potential. The effect of the coupling strength on the gap size is
shown in figure 3(a). The dependency of gap size on the temperature (figure
3(b)) shows a phase transition from Mott insulator to metal transition. The
gap size decreases as the bulk fermion mass increases from zero to
$\frac{1}{2}$. The nature of the spectral function changes to pole type when
$m_{f}=\frac{1}{2}$ which is consistent with the previous investigation
zaanen2015holographic . Since gap generation is due to interacting term, and
$m_{f}=\frac{1}{2}$ in holography leads to non-interacting theory, the gap
tends to vanish as $m_{f}=\frac{1}{2}$. However, we observe that there is a
finite gap size at $m_{f}=\frac{1}{2}$ when coupling strength is turned on. We
have shown the effect of coupling strength, fermion mass and temperature on
the gap size in figure 3.

(a) Gap vs Coupling

(b) Gap vs Temperature

(c) Gap vs Fermion mass

Figure 3: Mott gap size $\Delta_{M}$ from DoS where DoS$\leq 0.001$ value is
considered as cutoff for gap. (a) At $T=0.025\mu,m_{f}=0$ (b) For
$\mu=2,\eta=1,m_{f}=0$, (c) For $T=0.025\mu,\eta=1$.

##  4 Classification of gaps in two flavor fermions

In this section, we promote our flavor analysis to two flavor fermions related
to sublattice symmetry in materials. We will examine whether the dipole
interaction can create a symmetric Mott gap. Besides this, we will also
investigate $F^{2}\bar{\psi}_{1}\Gamma^{scalar}\psi_{2}$ interaction. In the
two-flavor fermions setup, the Lagrangian density becomes

$\displaystyle\mathcal{L}=\sum_{j=1}^{2}i\bar{\psi}_{j}\left(\Gamma^{\mu}D_{\mu}-m_{f}\right)\psi_{j}-i\bar{\psi_{1}}L_{int}\psi_{2}-i\bar{\psi_{2}}L_{int}\psi_{1}~{}~{}.$
(15)

The corresponding coupled equation of motion of two flavor fermions read

$\displaystyle\left(\Gamma^{\mu}D_{\mu}-m_{f}\right)\psi_{1}-L_{int}\psi_{2}=0,~{}~{}~{}\text{and}~{}~{}~{}\left(\Gamma^{\mu}D_{\mu}-m_{f}\right)\psi_{2}-L_{int}\psi_{1}=0~{}.$
(16)

The procedure of Green’s function derivation for two flavor fermion is almost
the same with one flavor case, which is given in detail in Oh:2020cym ;
Sukrakarn:2023ncp . In the two flavor scenario, we mainly focus on the
standard-standard (SS) and standard-alternative (SA) quantization. One
noticeable point is that $G^{\text{SS}}=\Gamma^{5}G^{\text{SA}}$ relates SS
and SA quantizations in terms of the output of the spectral function. The
spectral function for two flavor case for $L_{int}=i\eta F^{2}\Gamma^{5}$ with
SS quantization is shown in figure 4. The effect of the charge and the fermion
mass on the spectral function is shown in figure 4 and 5 respectively. The
role of bulk fermion mass is to determine the singularity structure of the
Green’s function zaanen2015holographic . It changes branch-type singularity to
pole type singularity when $m_{f}\rightarrow\frac{1}{2}$. Unlike one flavor
case, the Mott gap in two flavors almost vanishes when $m_{f}=\frac{1}{2}$,
which is shown in the spectral function plot (see figure 5). Although the
interaction term is turned on, the system tends to be non-interacting when
$m_{f}=\frac{1}{2}$. In the two-flavor fermion with SS-quantization, the
$L_{int}=i\eta F^{2}\Gamma^{5}$ interaction only shows the Mott gap feature in
the spectral function. The dipole interaction term in two flavor fermion shows
no gap feature. We have summarized all non-minimal gauge field interactions in
Table 2.

(a) For charge $q_{f}=0$

(b) For charge $q_{f}=1$

Figure 4: Spectral function along with DoS for $L_{int}=i\eta F^{2}\Gamma^{5}$
at $T=0.025\mu$.

(a) For $m_{f}=0.1$

(b) For $m_{f}=0.4$

(c) For $m_{f}=0.5$

Figure 5: Effect of the fermion mass on SF at $T=0.025\mu$. Gauge |  $-i\bar{\psi}_{1}L_{int}\psi_{2}$ with SS-quantization, $SS=\Gamma^{5}SA$
---|---
$A=A_{t}(z)dt$ | Gapless | Gap | Flatband | Effect of $q$
$F^{2}$-term |  $\eta F^{2}$, $i\eta F^{2}\Gamma^{z}$ | $i\eta F^{2}\Gamma^{5}$ |  | Shifting FS
Dipole-term |  $i\eta F_{\mu\nu}\Gamma^{\mu\nu}$, $\eta F_{\mu\nu}\Gamma^{\mu\nu}\Gamma^{z}$ |  | $\eta F_{\mu\nu}\Gamma^{\mu\nu}\Gamma^{5}$ | Shifting & Bending
Table 2: Mott gaps and flatband from non-minimal interactions for two flavor
fermions.

For completeness, we now focus on the ordered gap generated by the symmetric
breaking. The superconducting gap is classified as an ordered gap. There are
three ordered gap types: $s$-, $p$\- and $d$-wave ordered gap. In holographic
models, these three ordered gaps have been realized from charged scalar
Hartnoll:2008vx ; Hartnoll:2008kx ; Horo2009 ; Horowitz:2009ij ;
ghorai2016higher ; Ghorai:2016tvk ; Ghorai:2021uby , vector Gubser:2008wv ;
Vegh:2010fc ; Ghorai:2022gzx , and symmetric tensor Benini_2010 ; Chen:2010mk
; Zeng:2010fx ; Gao:2011aa ; Krikun:2013iha ; Nishida:2014lta ; Krikun:2015tga
fields. The action for holographic superconductors (bosonic sector) are given
as follows:

$\displaystyle S_{b}=\int
d^{4}x\frac{\sqrt{-g}}{2\kappa^{2}}\left[R+\frac{6}{L^{2}}+2\kappa^{2}\left(-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\mathcal{L}_{matter}\right)\right]$
(17)

where $\mathcal{L}_{matter}$ for scalar, vector and tensor fields are given in
Yuk:2022lof , Ghorai:2023wpu  and Ghorai:2023vuo  respectively. In these
references, the spectral function analysis for $s,p,d$-wave superconductors
has also been done in detail, which is summarized in table 3. We set
$\kappa=1,L=1,q=2$. Given the value of $\frac{T}{\mu}$, we can solve
numerically all coupled bosonic field equations. In the presence of these
bosonic fields, we observe different ordered gaps in the fermionic spectral
function for the scalar field, vector field, and symmetric tensor field. To
incorporate the particle-hole symmetry, we have to consider Nambu-Gorkov (NG)
representation Yuk:2022lof , where conjugate $\psi_{c}=\psi^{*}$ is considered
as independent degree of freedom.

Ordered |  $\eta\bar{\psi}L_{int}\psi_{c}$ with Nambu-Gorkov spinor Ghorai:2023wpu
---|---
Gap | Sc. Gap | Flatband
$s$-wave | $\phi i\Gamma^{z}$ |
$p$-wave | $V_{\mu}\Gamma^{\mu},V_{\mu}i\Gamma^{5}\Gamma^{\mu}(R),V_{\mu\nu}\Gamma^{\mu\nu}$ |  $V_{\mu}\Gamma^{\mu}\Gamma^{z}$ (1-dim.)
$d$-wave | $B_{\mu\nu}\Gamma^{\mu}D^{\nu}$ |
Table 3: The ordered gap: $\phi$ is scalar field, $V_{\mu}$ is vector field
and $B_{\mu\nu}$ is symmetric tensor field. $V_{\mu\nu}$ is the covariant
derivative of the vector field and $R$ presents the rotation of Fermi arc. The
spatial component of vector field and symmetric tensor field show $p$-wave and
$d$-wave superconductivity respectively.

Here, we compare the order gap and the Mott gap from the density of states.
With the same parameters at $T=0.025\mu$, the Mott gap size is larger than
superconducting gap, as shown in figure 6. Given the values of $T$ and $\mu$,
the Mott gap is generated directly through the gauge coupling, whereas the
superconducting gap in fermion spectral function is triggered by the
condensation value, which is generated by the spontaneous symmetry breaking.
The dependency of parameters for Mott gaps and superconducting gaps show the
same behavior. However, the mechanisms of these two gap generations are
completely different. The superconducting gap in a fermionic spectral function
is associated with $U(1)$-symmetry breaking and minimal coupling between the
fermion and the order parameter, whereas the Mott gaps are associated with
non-minimal coupling of the gauge field and fermion. In the NG representation,
the one dimensional flatband is also realized from a spatial vector field
interaction with fermion ($V_{\mu}\Gamma^{\mu}\Gamma^{z}$), which shown in
table 3.

Figure 6: Compare of DoS for Mott gap and order gap at $T=0.025\mu$ with same
$\eta$.

##  5 Summary

In this paper, we have addressed the analysis of the density of states for
Mott gaps and superconducting gaps in holographic setups. Without manifest
symmetry breaking, the gap generation is known as the Mott gap, whereas the
superconducting gap is generated because of symmetry breaking. Using
gauge/gravity duality, these two gaps are explained in the literature, where
Mott physics is elucidated using dipole interaction PhysRevLett.106.091602 .
It is a non-minimal gauge coupling with fermions that generates a gap in the
fermionic spectral function without a symmetry-breaking mechanism. In the
previous investigation PhysRevLett.106.091602 , the Mott gap in the spectral
function is asymmetric, and one of the bands seems to touch the Fermi level at
$\omega=0$. To gain a clearer view of the Mott gap size, we have calculated
the density of states (DoS), which shows a soft and asymmetric Mott gap,
whereas DMFT results show differently. This motivates us to consider other
non-minimal gauge couplings with fermions. First, we have considered an
density-density coupling ($F^{2}\bar{\psi}\psi$) term, which can be mapped to
the interaction part of the Hubbard Hamiltonian. This mapping also justifies
our proposal. Then, we have calculated the spectral function along with the
DoS, which clearly shows a strong and symmetric Mott gap feature. Other
possible non-minimal gauge field interactions (26) have been examined, and the
Mott gap size in the DoS (figure 2) has been compared. The DoS figure also
supports our claim about the proper Mott gap in the fermionic spectral. We
have classified Mott gaps, namely, soft gap and hard gap. For the dipole
interaction case, we can treat the Mott gap as a soft and asymmetric gap,
whereas the $F^{2}$ interactions produce a hard and symmetric Mott gap.
The effects of the fermion charge, coupling strength, temperature, and fermion
mass are important for further investigation. To observe the effects of these
parameters on the Mott gaps, we measured the gap size $\Delta_{M}$ from the
DoS, where DoS$\leq 0.001$ region is a gapped region. The gap size is plotted
in figure 3. We observed that the Mott gap is vanishing at $T=0.13\mu$.
Therefore, this can be treated as the critical temperature for the Mott
insulator-metal transition. Finite fermion mass decreases the Mott gap size.
In the limit of the fermion mass $m_{f}=\frac{1}{2}$, we found that the
structure of singularity type changes from a branch-cut type to a simple pole
structure, which matches with previous findings zaanen2015holographic . We
also have found that the dipole interaction term with $\Gamma^{5}$ shows a
flatband in the spectral function. This is the non-minimal gauge interaction
which creates a flatband in a holographic setup.
Next, we have classified interactions in terms of Mott and superconducting
gaps. To achieve this, we first examined all non-minimal gauge field
interactions in two-flavor fermions. We found that in two-flavor fermions,
only the density interaction ($F^{2}i\Gamma^{5}$) creates a Mott gap, while
the dipole term doesn’t produce any gap feature in the spectral function. For
the fermion mass $m_{f}=\frac{1}{2}$, the system nearly becomes non-
interacting, despite the interaction strength being high $(\eta=1)$, which is
not seen in the one-flavor fermion case. We have replicated three types of
order gaps, which are summarized in table 3. When comparing the gap size in
DoS for the same parameters, we observe that the Mott gap size is larger than
the superconducting gap which is consistent with literature. Our future
direction is to study the Fermi-Dirac distribution function from the
holographic fermion.

###### Acknowledgements.

This work is supported by Mid-career Researcher Program through the National
Research Foundation of Korea grant No. NRF-2021R1A2B5B02002603,
RS-2023-00218998 and NRF-2022H1D3A3A01077468. We thank the APCTP for the
hospitality during the focus program, where part of this work was discussed.

##  Appendix A Mott Physics

Hubbard Hamiltonian (HH) is given by

$\displaystyle
H=-\sum_{<ij>,\sigma}t_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}+H_{I}~{}.$ (18)

where interaction Hamiltonian $H_{I}=U\sum_{i}n_{i\uparrow}n_{i\downarrow}$
with the fermionic number operator
$n_{i\uparrow}=c^{\dagger}_{i\uparrow}c_{i\uparrow}$. In the Mean field
approximation, we can write

$\displaystyle n_{i\uparrow}n_{i\downarrow}$ $\displaystyle\approx$
$\displaystyle n_{i\uparrow}\langle
n_{i\downarrow}\rangle+n_{i\downarrow}\langle n_{i\uparrow}\rangle-\langle
n_{i\uparrow}\rangle\langle n_{i\downarrow}\rangle$ (19)

which gives the interaction term in the HH as follows:

$\displaystyle H_{I}=U\sum_{i}\left[n_{i\uparrow}\langle
n_{i\downarrow}\rangle+n_{i\downarrow}\langle n_{i\uparrow}\rangle-\langle
n_{i\uparrow}\rangle\langle n_{i\downarrow}\rangle\right]~{}.$ (20)

The interpretation of this expression is that the up-spin fermions interact
with the average density of the down-spin fermions, and similarly the down-
spin fermions interact with the average density of the up-spin fermions. This
HH governs Mott metal-insulator transition. We need to find a similar kind of
interaction term in the holographic setup. The number density can be expressed
in terms of charge density ($\rho$) which maps to the gauge field in the bulk
theory. First, we need to assume that $\langle n_{i\uparrow}\rangle=\langle
n_{i\downarrow}\rangle\propto\rho$. Therefore, the interaction part now
becomes

$\displaystyle H_{I}\propto\rho
U\sum_{i}\left[c^{\dagger}_{i\uparrow}c_{i\uparrow}+c^{\dagger}_{i\downarrow}c_{i\downarrow}\right]-U\rho^{2}$
(21)

which translates in momentum space as

$\displaystyle H_{I}\propto
U\rho\sum_{k}\left[c^{\dagger}_{k\uparrow}c_{k\uparrow}+c^{\dagger}_{k\downarrow}c_{k\downarrow}\right]-U\rho^{2}$
$\displaystyle\approx$ $\displaystyle
U\rho\left[c^{\dagger}_{+\uparrow}c_{+\uparrow}+c^{\dagger}_{+\downarrow}c_{+\downarrow}+c^{\dagger}_{-\uparrow}c_{-\uparrow}+c^{\dagger}_{-\downarrow}c_{-\downarrow}\right]-U\rho^{2}~{}~{}~{}~{}~{}~{}~{}~{}$
(22)

The last term in the right hand side of the above equation gives the shifting
of the energy. We need to find a interaction terms in the holographic setup,
which can map to
$\left[c^{\dagger}_{+\uparrow}c_{+\uparrow}+c^{\dagger}_{+\downarrow}c_{+\downarrow}+c^{\dagger}_{-\uparrow}c_{-\uparrow}+c^{\dagger}_{-\downarrow}c_{-\downarrow}\right]$.
We can identify the bulk fermion field
$\psi=(c_{+\uparrow},c_{+\downarrow},c_{-\uparrow},c_{-\downarrow})^{T}$ where
we are embedding $H_{I}$ into a higher dimension. Now the fermionic
annihilation $(c_{\pm\uparrow\downarrow})$ and creation
$(c^{\dagger}_{\pm\uparrow\downarrow})$ operators are function of radial
coordinate and momentum Yuk:2022lof . Therefore, the suitable non-minimal
coupling term for gauge field is $F^{2}\psi^{\dagger}\psi$ which is the non-
relativistic(NR) interaction. We can promot this NR interaction term to the
relativistic bulk interaction term in the following way:

$\displaystyle\mathcal{L}_{int}=-i\bar{\psi}\eta F^{2}\psi$ (23)

where $\eta$ is the coupling strength and $F^{2}=F_{\mu\nu}F^{\mu\nu}$. Note
that $i\Gamma^{\underline{t}}$ is the dressing factor for the promoting
relativistic form from the NR interaction since the kinetic term of the
Lagrangian also contains this dressing factor. Since potential term in a
Lagrangian is opposite to the Hamiltonian, we need to consider the negative
sign in the interaction Lagrangian. Therefore, the most suitable In general,
we can consider the following interaction Lagrangian

$\displaystyle\mathcal{L}_{int}=\left\\{\begin{array}[]{ c l
}-i\eta\bar{\psi}F^{2}\Gamma^{scalar}\psi&\quad\textrm{Density type
interaction}\\\ -i\eta
F_{\mu\nu}\bar{\psi}\Gamma^{\mu\nu}\Gamma^{scalar}\psi&\quad\textrm{Pauli
interaction}\end{array}\right.~{}~{}.$ (26)

where the possible gamma matrices are
$\Gamma^{scalar}=\mathbb{I}_{4},\Gamma^{z},i\Gamma^{5},i\Gamma^{5z}$ although
other interaction terms with different $\Gamma^{scalar}$ are not proportional
to
$\left[c^{\dagger}_{+\uparrow}c_{+\uparrow}+c^{\dagger}_{+\downarrow}c_{+\downarrow}+c^{\dagger}_{-\uparrow}c_{-\uparrow}+c^{\dagger}_{-\downarrow}c_{-\downarrow}\right]$.
For the Pauli interaction $(\Gamma^{scalar}=\mathbb{I}_{4})$ in the non-
relativistic limit, the interaction term
$F_{\mu\nu}\bar{\psi}\Gamma^{\mu\nu}\psi$ becomes
$\psi^{\dagger}F_{\mu\nu}\Gamma^{\mu\nu}\psi$. For the gauge field ansatz
$A=A_{t}(z)dt$, the dipole interaction term can be expressed as

$\displaystyle\mathcal{L}_{int}^{dipole}\propto\psi^{\dagger}F_{zt}\Gamma^{zt}\psi=z^{2}(\partial_{z}A_{t})\left[c^{\dagger}_{+\uparrow}c_{-\downarrow}-c^{\dagger}_{+\downarrow}c_{-\uparrow}-c^{\dagger}_{-\uparrow}c_{+\downarrow}+c^{\dagger}_{-\downarrow}c_{+\uparrow}\right]~{}.$
(27)

The above expression can not be mapped to
$\left[c^{\dagger}_{+\uparrow}c_{+\uparrow}+c^{\dagger}_{+\downarrow}c_{+\downarrow}+c^{\dagger}_{-\uparrow}c_{-\uparrow}+c^{\dagger}_{-\downarrow}c_{-\downarrow}\right]$
term. Therefore, we can argue that the reasonable interaction for Mott gap is
the density type interaction (23). Since some interactions in eq.(26) produce
gap feature in the spectral function, we will examine all possible interactios
to find Mott gap feature in the DoS.

##  Appendix B Derivation of Green’s function

Rearranging all components of equations, we can recast the Dirac equations in
the following structure

$\displaystyle\partial_{z}\Psi_{+}+\mathbb{M}_{1}\Psi_{+}+\mathbb{M}_{2}\Psi_{-}$
$\displaystyle=$ $\displaystyle 0~{},$ (28)
$\displaystyle\partial_{z}\Psi_{-}+\mathbb{M}_{3}\Psi_{-}+\mathbb{M}_{4}\Psi_{+}$
$\displaystyle=$ $\displaystyle 0$ (29)

where $2\times 2$-matrix $\mathbb{M}_{i},~{}i=1,2,3,4$ are given by

$\displaystyle\mathbb{M}_{1}=$ $\displaystyle\ -\frac{m_{f}+\eta
F^{2}}{z\sqrt{f}}\mathbb{I}_{2},$ $\displaystyle\mathbb{M}_{2}=$
$\displaystyle\
\frac{i}{\sqrt{f}}\begin{pmatrix}k_{y}&-\frac{(\omega+q_{f}A_{t})}{\sqrt{f}}+k_{x}\\\
\frac{(\omega+q_{f}A_{t})}{\sqrt{f}}+k_{x}&-k_{y}\end{pmatrix}~{},$
$\displaystyle\mathbb{M}_{3}=$ $\displaystyle\
-\mathbb{M}_{1},\quad\mathbb{M}_{4}=-\mathbb{M}_{2}~{}.~{}~{}~{}$ (30)

There are two independent solutions since $\Psi_{+}$ and $\Psi_{-}$ are two
components spinor. The general solution can be written in a linear combination
of two solutions as

$\displaystyle\Psi_{+}=\sum_{i=1}^{2}c_{i}\Psi_{+}^{(i)}=\mathbb{S}(z)\mathbf{c},$
$\displaystyle\Psi_{-}=\sum_{i=1}^{2}c_{i}\Psi_{-}^{(i)}=\mathbb{C}(z)\mathbf{c}$
(31)

The $2\times 2$-matrix $\mathbb{S}(z),\mathbb{C}(z)$ are constructed from the
solution, where the constant column vector $\mathbf{c}$ is constructed from
the two coefficients of the linear combination. Substituting above eq.(31) in
eq.(s)(28, 29), we obtain

$\displaystyle\partial_{z}\mathbb{S}(z)+\mathbb{M}_{1}\mathbb{S}(z)+\mathbb{M}_{2}\mathbb{C}(z)$
$\displaystyle=$ $\displaystyle 0~{},$ (32)
$\displaystyle\partial_{z}\mathbb{C}(z)+\mathbb{M}_{3}\mathbb{C}(z)+\mathbb{M}_{4}\mathbb{S}(z)$
$\displaystyle=$ $\displaystyle 0~{}.$ (33)

The boundary solution tells us that we need to define
$U(z)=diag(z^{m_{f}},z^{m_{f}})$ to get normalized boundary Green’s function.
Then we can write the boundary solution from eq.(31)

$\displaystyle\Psi_{+}(z)\overset{z\rightarrow
0}{\approx}U(z)\mathbb{S}_{0}\mathbf{c},~{}~{}~{}~{}\Psi_{-}\overset{z\rightarrow
0}{\approx}U(z)^{-1}\mathbb{C}_{0}\mathbf{c}~{},$ (34)

where $\mathbb{S}_{0},\mathbb{C}_{0}$ are the $z$-independent boundary
$2\times 2$-matrix. We can define

$\displaystyle\mathcal{J}=\mathbb{S}_{0}\mathbf{c},~{}~{}~{}~{}~{}\mathcal{C}=\mathbb{C}_{0}\mathbf{c}$
(35)

which translate the boundary solution (eq.(34)) as

$\displaystyle\Psi_{+}\overset{z\rightarrow
0}{\approx}U(z)\mathcal{J},~{}~{}~{}~{}\Psi_{-}\overset{z\rightarrow
0}{\approx}U(z)^{-1}\mathcal{C}~{}.$ (36)

Comparing eq.(36) with eq.(9), we find

$\displaystyle\mathcal{J}=\mathbf{A},~{}~{}~{}~{}\mathcal{C}=\mathbf{D}$ (37)

We can also get the relation between $\mathcal{C}$ and $\mathcal{J}$ from
eq.(35)

$\displaystyle\mathcal{C}=\mathbb{C}_{0}\mathbb{S}^{-1}_{0}\mathcal{J}~{}.$
(38)

From the boundary action (eq.(10)), we can write

$\displaystyle S_{bdy}=\int
d^{3}x\mathcal{J}^{\dagger}\tilde{\Gamma}\mathcal{C}+h.c.$ (39)

Using eq.(38), the boundary action now becomes

$\displaystyle S_{bdy}=\int
d^{3}x\mathcal{J}^{\dagger}\tilde{\Gamma}\mathbb{C}_{0}\mathbb{S}^{-1}_{0}\mathcal{J}+h.c.=\int
d^{3}x\mathcal{J}^{\dagger}\mathbb{G}_{R}\mathcal{J}+h.c.~{}~{}~{}~{}~{}$ (40)

where the retarded Green’s function
$\mathbb{G}_{R}=\tilde{\Gamma}\mathbb{C}_{0}\mathbb{S}^{-1}_{0}$. We can
promote this boundary Green’s function into bulk Green’s function by
considering the $z$-dependent Green’s function as follows:

$\displaystyle\mathbb{G}=\tilde{\Gamma}\mathbb{C}(z)\mathbb{S}^{-1}(z)$ (41)

where $\mathbb{C}(z),\mathbb{S}(z)$ is defined in eq.(31). Taking the
derivative of the above equation, we get

$\displaystyle\partial_{z}\mathbb{G}(z)=\tilde{\Gamma}\left[\partial_{z}\mathbb{C}(z)\mathbb{S}^{-1}(z)-\mathbb{C}(z)\mathbb{S}^{-1}(z)(\partial_{z}\mathbb{S}(z))\mathbb{S}^{-1}(z)\right]~{}~{}~{}~{}~{}$
(42)

Using eq.(s)(32,33), we have found

$\displaystyle\partial_{z}\mathbb{G}(z)+\tilde{\Gamma}\mathbb{M}_{3}\tilde{\Gamma}\mathbb{G}(z)-\mathbb{G}(z)\mathbb{M}_{1}-\mathbb{G}(z)\mathbb{M}_{2}\tilde{\Gamma}\mathbb{G}(z)+\tilde{\Gamma}\mathbb{M}_{4}=0~{}~{}~{}~{}~{}$
(43)

This is the desired flow equation to know the bulk Green’s function
$\mathbb{G}(z)$. From eq.(34), we can express

$\displaystyle\mathbb{S}(z)\overset{z\rightarrow
0}{\approx}U(z)\mathbb{S}_{0}~{}~{}\text{and}~{}~{}\mathbb{C}(z)\overset{z\rightarrow
0}{\approx}U(z)^{-1}\mathbb{C}_{0}~{}.$ (44)

By substituting the above relations, we can now map the boundary Green’s
function with bulk Green’s function near the boundary in the following way

$\displaystyle\mathbb{G}(z)\overset{z\rightarrow
0}{\approx}U(z)^{-1}\mathbb{G}_{R}U(z)^{-1}$ (45)

where we have used the fact $\tilde{\Gamma}U(z)^{-1}\tilde{\Gamma}=U(z)^{-1}$.
To solve the flow equation, we need to know the horizon value of the Green’s
function which is $\mathbb{G}(z_{h})=i\mathbb{I}_{2}$.

## References

  * (1) P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a mott insulator: Physics of high-temperature superconductivity, Rev. Mod. Phys. 78 (Jan, 2006) 17–85.
  * (2) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68 (Jan, 1996) 13–125.
  * (3) P. W. Phillips, L. Yeo, and E. W. Huang, Exact theory for superconductivity in a doped mott insulator, Nature Physics 16 (2020) 1175–1180.
  * (4) N. F. MOTT, Metal-insulator transition, Rev. Mod. Phys. 40 (Oct, 1968) 677–683.
  * (5) J. M. Maldacena, The Large N limit of superconformal field theories and supergravity, Int.J.Theor.Phys. 38 (1999) 1113–1133, [hep-th/9711200].
  * (6) S. S. Gubser, I. R. Klebanov, and A. M. Polyakov, Gauge theory correlators from noncritical string theory, Phys. Lett. B428 (1998) 105–114, [hep-th/9802109].
  * (7) E. Witten, Anti-de Sitter space and holography, Adv. Theor. Math. Phys. 2 (1998) 253–291, [hep-th/9802150].
  * (8) S. A. Hartnoll, P. K. Kovtun, and M. S. Müller, Theory of the Nernst effect near quantum phase transitions in condensed matter and in dyonic black holes, Phys. Rev. B 76 (Oct., 2007) 144502, [arXiv:0706.3215].
  * (9) Y. Seo, G. Song, P. Kim, S. Sachdev, and S.-J. Sin, Holography of the Dirac Fluid in Graphene with two currents, Phys. Rev. Lett. 118 (2017) 036601, [arXiv:1609.03582].
  * (10) E. Oh, T. Yuk, and S.-J. Sin, The emergence of strange metal and topological liquid near quantum critical point in a solvable model, JHEP 11 (2021) 207, [arXiv:2103.08166].
  * (11) G. Song, J. Rong, and S.-J. Sin, Stability of topology in interacting Weyl semi-metal and topological dipole in holography, JHEP 10 (2019) 109, [arXiv:1904.09349].
  * (12) E. Oh, Y. Seo, T. Yuk, and S.-J. Sin, Ginzberg-Landau-Wilson theory for Flat band, Fermi-arc and surface states of strongly correlated systems, arXiv:2007.12188.
  * (13) E. Oh and S.-J. Sin, Entanglement String and Spin Liquid with Holographic Duality, arXiv:1811.07299.
  * (14) Y. Seo, G. Song, and S.-J. Sin, Strong Correlation Effects on Surfaces of Topological Insulators via Holography, Phys. Rev. B96 (2017), no. 4 041104, [arXiv:1703.07361].
  * (15) Y. Seo, G. Song, C. Park, and S.-J. Sin, Small Fermi Surfaces and Strong Correlation Effects in Dirac Materials with Holography, JHEP 10 (2017) 204, [arXiv:1708.02257].
  * (16) S. Sukrakarn, T. Yuk, and S.-J. Sin, Mean field theory for strongly coupled systems: Holographic approach, arXiv:2311.01897.
  * (17) M. Edalati, R. G. Leigh, and P. W. Phillips, Dynamically generated mott gap from holography, Phys. Rev. Lett. 106 (Mar, 2011) 091602.
  * (18) M. Edalati, R. G. Leigh, K. W. Lo, and P. W. Phillips, Dynamical gap and cupratelike physics from holography, Phys. Rev. D 83 (Feb, 2011) 046012.
  * (19) G. Vanacore and P. W. Phillips, Minding the gap in holographic models of interacting fermions, Phys. Rev. D 90 (Aug, 2014) 044022.
  * (20) G. Vanacore, S. Ramamurthy, and P. Phillips, Evolution of holographic fermi arcs from a mott insulator, Journal of High Energy Physics 2018 (2018), no. 9.
  * (21) Y. Seo, G. Song, Y.-H. Qi, and S.-J. Sin, Mott transition with holographic spectral function, Journal of High Energy Physics 2018 (Aug., 2018).
  * (22) D. Ghorai, T. Yuk, and S.-J. Sin, Fermi arc in p-wave holographic superconductors, JHEP 10 (2023) 003, [arXiv:2304.14650].
  * (23) D. Ghorai, T. Yuk, and S.-J. Sin, Order parameter and spectral function in d-wave holographic superconductors, Phys. Rev. D 109 (2024), no. 6 066004, [arXiv:2309.01634].
  * (24) J. Zaanen, Y. Liu, Y.-W. Sun, and K. Schalm, Holographic duality in condensed matter physics.  Cambridge University Press, 2015.
  * (25) S. A. Hartnoll, C. P. Herzog, and G. T. Horowitz, Building a Holographic Superconductor, Phys.Rev.Lett. 101 (2008) 031601, [arXiv:0803.3295].
  * (26) S. A. Hartnoll, C. P. Herzog, and G. T. Horowitz, Holographic Superconductors, JHEP 12 (2008) 015, [arXiv:0810.1563].
  * (27) G. T. Horowitz and M. M. Roberts, Holographic superconductors with various condensates, Physical Review D 78 (2008), no. 12 126008.
  * (28) G. T. Horowitz and M. M. Roberts, Zero temperature limit of holographic superconductors, JHEP 11 (2009) 015.
  * (29) D. Ghorai and S. Gangopadhyay, Higher dimensional holographic superconductors in born–infeld electrodynamics with back-reaction, The European Physical Journal C 76 (2016), no. 3 1–12.
  * (30) D. Ghorai and S. Gangopadhyay, Holographic free energy and thermodynamic geometry, Eur. Phys. J. C 76 (2016), no. 12 702, [arXiv:1607.05187].
  * (31) D. Ghorai and S. Gangopadhyay, Analytical study of holographic superconductor with backreaction in 4d Gauss-Bonnet gravity, Phys. Lett. B 822 (2021) 136699, [arXiv:2105.09423].
  * (32) S. S. Gubser and S. S. Pufu, The Gravity dual of a p-wave superconductor, JHEP 11 (2008) 033, [arXiv:0805.2960].
  * (33) D. Vegh, Fermi arcs from holography, arXiv:1007.0246.
  * (34) D. Ghorai, Y.-S. Choun, and S.-J. Sin, Momentum dependent gap in holographic superconductors revisited, JHEP 09 (2022) 098, [arXiv:2205.02514].
  * (35) F. Benini, C. P. Herzog, R. Rahman, and A. Yarom, Gauge gravity duality for d-wave superconductors: prospects and challenges, Journal of High Energy Physics 2010 (nov, 2010).
  * (36) J.-W. Chen, Y.-J. Kao, D. Maity, W.-Y. Wen, and C.-P. Yeh, Towards A Holographic Model of D-Wave Superconductors, Phys. Rev. D 81 (2010) 106008, [arXiv:1003.2991].
  * (37) H.-B. Zeng, Z.-Y. Fan, and H.-S. Zong, Characteristic length of a Holographic Superconductor with $d$-wave gap, Phys. Rev. D 82 (2010) 126014, [arXiv:1006.5483].
  * (38) D. Gao, Vortex and droplet in holographic D-wave superconductors, Phys. Lett. A 376 (2012) 1705–1709, [arXiv:1112.2422].
  * (39) A. Krikun, Charge density wave instability in holographic d-wave superconductor, JHEP 04 (2014) 135, [arXiv:1312.1588].
  * (40) M. Nishida, Phase Diagram of a Holographic Superconductor Model with s-wave and d-wave, JHEP 09 (2014) 154, [arXiv:1403.6070].
  * (41) A. Krikun, Phases of holographic d-wave superconductor, JHEP 10 (2015) 123, [arXiv:1506.05379].
  * (42) T. Yuk and S.-J. Sin, Flow equation and fermion gap in the holographic superconductors, JHEP 02 (2023) 121, [arXiv:2208.03132].