#  Toward fixed point and pulsation quantum search on graphs driven by quantum
walks with in- and out-flows: a trial to the complete graph

Yusuke Higuchi Mohamed Sabri Etsuo Segawa Department of Mathematics, Gakushuin
University, Tokyo 171-8588, Japan.Graduate School of Information Sciences,
Tohoku University, Aoba, Sendai 980-8579, Japan Graduate School of Environment
and Information Sciences, Yokohama National University, Hodogaya, Yokohama
240-8501, Japan

Abstract. We treat a quantum walk model with in- and out- flows at every time
step from the outside. We show that this quantum walk can find the marked
vertex of the complete graph with a high probability in the stationary state.
In exchange of the stability, the convergence time is estimated by $O(N\log
N)$, where $N$ is the number of vertices. However until the time step $O(N)$,
we show that there is a pulsation with the periodicity $O(\sqrt{N})$. We find
the marked vertex with a high relative probability in this pulsation phase.
This means that we have two chances to find the marked vertex with a high
relative probability; the first chance visits in the pulsation phase at short
time step $O(\sqrt{N})$ while the second chance visits in the stable phase
after long time step $O(N\log N)$. The proofs are based on Kato’s perturbation
theory.

Keywords: Quantum walk; Perturbation theory; Quantum search.

##  1 Introduction

Quantum search driven by quantum walks on graphs is one of the interesting
topics of studies on quantum walks with fruitful results (see for example [1,
2, 3] and its references). The time evolution operator of quantum walk is
driven by a unitary operator. Since the absolute values of the eigenvalues of
the unitary matrix are unit, each eigenvalue is rotated on the unit circle in
the unitary time iteration. Then the eigenspaces of the time evolution
operator which have assymptotically large overlap to the initial state with
respect to the system size are dominant. This is the reason for the
oscillation of the time course of the finding probability at the target vertex
in the quantum search algorithm. Thus we need to perform the sharpshooting to
the appropriate observation time, otherwise we would find the target vertex
with a very small probability. This oscillation is known as the soufflé
problem [4]. The soufflé problem has been solved by constructing the nested
time evolution which produces the convergence to a fixed point [5, 6].

In this paper, we also consider the dynamical system which converges to a
fixed point driven by a quantum walk but we set a time independent unitary
evolution operator and extend the domain of the initial state from the square
summable space to the uniformly bounded space. The initial state, which is no
longer square-summable, naturally provides the situation that the inflow
penetrates the internal graph at every time step while the outflow also goes
out from the internal graph. This quantum walk model is firstly introduced by
[7, 8] as a discrete-analogue of the stationary Schrödinger equation [9]. In
[10, 11], detailed aspect of the stationary state for the Szegedy walk with a
constant inflow are characterized by the electric circuit.

To see whether this quantum walk model with in- and out- flows finds the
marked vertex with a high relative probability, in this paper, we try to set
the complete graph with the vertex number $N$ as the internal graph and give a
perturbation to the marked vertex as a first trial. We show that (i) a high
relative probability at the marked vertex can be seen in the stable state;
(ii) the convergence time is estimated by $O(N\log N)$; (iii) there is a
pulsation with the periodicity $\Theta(\sqrt{N})$ until the time step $O(N)$.
This means that at time step $O(\sqrt{N})$, we can find the marked vertex with
a high probability. In addition, there is an insurance for missing this timing
of the observation; that is, the chance to find the marked vertex with a high
probability stationarily visits again after waiting for this dynamical system
$O(N\log N)$.

This paper is organized as follows. In Section 2, we set the definition of our
quantum walk model. In Section 3, we give the main results. In Section 4, the
proofs of the main results are shown. The proof is based on Kato’s
perturbation theory [12].

##  2 Setting

In this paper, we fix the internal graph as the complete graph with $N$
vertices; $K_{N}$. We set the semi-infinite length path and join its end
vertex to each vertex of the complete graph. Such resulting infinite graph
(like a hedgehog) is denoted by $\tilde{K}_{N}$. More precisely, let
$\left\\{\mathbb{P}_{j}:j=1,...,N\right\\}$ be the set of semi-infinite length
paths called the tails whose end vertex is denoted by $o(\mathbb{P}_{j})$ and
the vertex set of $K_{N}$ be $\\{u_{1},\dots,u_{N}\\}$. Then the tail
$\mathbb{P}_{j}$ is connected to the complete graph $K_{N}$ by identifying
$o(\mathbb{P}_{j})$ with $u_{j}\in V(K_{N})$. Next, let us introduce the
notion of the symmetric arc, which plays an important role to describe the
dynamics of the quantum walk. The symmetric arc set induced by $E(K_{N})$ and
$E(\tilde{K}_{N})$ by $A$ and $\tilde{A}$ respectively. For any
$a\in\tilde{A}$, there uniquely exists the inverse arc of $a$ denoted by
$\bar{a}$ in $\tilde{A}$. The origin and terminal vertices of $a\in\tilde{A}$
are denoted by $o(a)$ and $t(a)$, respectively. Then $o(a)=t(\bar{a})$ and
$t(a)=o(\bar{a})$ hold. The degree of vertex $u\in V(\tilde{K}_{N})$ is
denoted by $\deg(u):=\\#\\{a\in\tilde{A}\;|\;t(a)=u\\}$. Furthermore, we
choose a vertex $u_{*}$ as the marked vertex from $V(K_{N})$.

For any countable set $\Omega$, we describe $\mathbb{C}^{\Omega}$ by the
vector space whose standard basis are labeled by each element of $\Omega$. The
total state space of our quantum walk is described by
$\mathbb{C}^{\tilde{A}}$. The time evolution operator $U$ on
$\mathbb{C}^{\tilde{A}}$ is defined as follows.

$(U\Psi)(a)=(-1)^{\bm{1}_{\\{o(a)=u_{*}\\}}(a)}\left(-\Psi(\bar{a})+\frac{2}{\deg(u)}\sum_{b\in\tilde{A}:\;t(b)=o(a)}\Psi(b)\right)$
(1)

for any $a\in\tilde{A}$ and $\Psi\in\mathbb{C}^{\tilde{A}}$. Let $\Psi_{t}$ be
the $t$-th iteration of the time evolution operator $U$; that is,
$\Psi_{t+1}=U\psi_{t}$. This dynamics can be also explained as follows. Let
$\\{a_{1},\dots,a_{r}\\}\subset\tilde{A}$ be the set of all arcs whose
terminal vertices are arbitrary chosen $u\in V(\tilde{K}_{N})$. Note that the
value $r$ takes two possibilities; if $u\in V(K_{N})$, then $r=N$, while
$u\notin V(K_{N})$, then $r=2$. Then $\Psi_{t+1}$ and $\Psi_{t}$ locally
satisfies the following recursion.

$\begin{bmatrix}\Psi_{t+1}(\bar{a}_{1})\\\ \vdots\\\
\Psi_{t+1}(\bar{a}_{r})\end{bmatrix}=\begin{cases}-\mathrm{Gr}(r)\begin{bmatrix}\Psi_{t}(a_{1})\\\
\vdots\\\ \Psi_{t}(a_{r})\end{bmatrix}&\text{: $u=u_{*}$,}\\\ \\\
+\mathrm{Gr}(r)\begin{bmatrix}\Psi_{t}(a_{1})\\\ \vdots\\\
\Psi_{t}(a_{r})\end{bmatrix}&\text{: $u\neq u_{*}$.}\end{cases}$ (2)

Here $\mathrm{Gr}(r)$ is the $r$-dimensional Grover matrix; that is,
$\mathrm{Gr}(r)=(2/r)J_{r}-I_{r}$, where $J_{r}$ and $I_{r}$ are the all one
matrix and the identity matrix, respectively. This represents the dynamics
that a quantum walker on the arc $b\in\tilde{A}$ is locally scattered at each
vertex by the Grover matrix with the phase varied operator which is active
only on the marked vertex $u_{*}$. Note that if $r=2$; then the Grover matrix
is reduced to the Pauli matrix $\sigma_{X}$. This means that the walk on the
tails is the free. The initial state $\Psi_{0}\in\mathbb{C}^{\tilde{A}}$ is
set so that a quantum walker penetrates $K_{N}$ at every time step.

$\Psi_{0}(a)=\begin{cases}1&\text{: $a\in\tilde{A}\setminus A$,
$\mathrm{dist}(t(a),K_{N})<\mathrm{dist}(o(a),K_{N})$}\\\ 0&\text{:
otherwise}\end{cases}\;\;(a\in\tilde{A}),$ (3)

where $\mathrm{dist}(u,K_{N})$ is the shortest distance between the vertex of
the tail $u$ and $K_{N}$. Since the walk is free on the tails, the internal
graph $K_{N}$ receives the inflow from the tails at every time step $t$, on
the other hand, once a quantum walker goes out a tail, then it never go back,
which can be regarded as the outflow. It is shown in [10] that this dynamical
system converges to a fixed point; that is, there exists
$\Psi_{\infty}\in\mathbb{C}^{\infty}$ such that

$\lim_{t\to\infty}\Psi_{t}(a)=\Psi_{\infty}(a)$

for any $a\in\tilde{A}$. The relative probability of finding a vertex $u\in
V(K_{N})$ at time $t$ is given by

$\tilde{\nu}_{t}(u)=\sum\limits_{a\in A_{0}:t(a)=u}|\psi_{t}(a)|^{2}.$

Taking the normalization to $\tilde{\nu_{t}}$, we set

$\nu_{t}(u):=\frac{\tilde{\nu}_{t}(u)}{\sum_{v\in
V(K_{N})}\tilde{\nu}_{t}(v)}.$

##  3 Main results

As a first interest, whether our quantum walker find a target vertex $u_{*}$
with a high probability in the long time limit. To see that, we define the
following distribution

######  Definition 1.

(The limit distribution) The normalized finding measure on $V(K_{N})$ is
defined by

$\mu_{N}(u)=\lim_{t\to\infty}\nu_{t}(u)$

for any $u\in V(K_{N})$.

######  Theorem 1.

(The finding probability of the marked vertex at the fixed point)
Let $\mu$ be the limit distribution of the finding probability defined as the
above. Then we have

$\lim_{N\to\infty}\mu_{N}(u_{*})=1/2.$

This theorem implies that we can find the marked vertex in the stable state
with probability $1/2$ while the celebrated quantum search algorithm has the
oscillation. Since the celebrated quantum search algorithm takes $O(\sqrt{N})$
if we measure the system at an appropriate time, now the next interest may be
that how does it take time for the convergence ? To answer the question, let
us introduce the finite time which we regard as convergence.

######  Definition 2.

($\ell^{2}$-mixing time of quantum walk) Let $\Psi_{t}$ be the $t$-th
iteration of the quantum walk and $\Psi_{\infty}$ be its stationary state. For
$\theta>0$, we set $t(\theta)$ by

$t(\theta):=\min\\{s>0\;:\;{}^{\forall}t>s,\;||\Psi_{\infty}-\Psi_{t}||_{K_{N}}<e^{-\theta}\\},$

where for $f\in\mathbb{C}^{\tilde{A}}$, the norm $||f||_{K_{N}}$ is the
$\ell^{2}$-norm with respect to the internal graph; that is,

$||f||_{K_{N}}^{2}:=\sum_{a\in A:\;t(a)\in V(K_{N})}|f(a)|^{2}.$

######  Theorem 2.

(The $\ell^{2}$-mixing time) Let $t(\theta)$ be the $\ell^{2}$-valued mixing
time defined as the above. Then we have

$t(\theta)\in\Theta(N\log N)$

for any fixed $\theta>0$.

Then unfortunately, in exchange for the stability, we need to wait for this
quantum walk until $\Theta(N\log N)$, which means a speed down, on the
contrary, this quantum walk takes more time than a classical algorithm.
However we obtain a good news as follows.

######  Theorem 3.

(Pulsation) Until the time step $O(N)$, the finding probability of the marked
vertex is estimated by

$\nu_{t}(u_{*})\sim\frac{1}{2}\frac{(\;1-c_{t}\cos[t\sqrt{2/N}]\;)^{2}}{1+c_{t}^{2}},$

where $c_{t}=e^{-(5/2)\cdot(t/N)}$ for large size $N$. In particular, at the
time step $t\sim\pi\sqrt{N/2}$, the finding probability of the marked vertex
is

$\nu_{t}(u_{*})\sim\frac{1}{2}\frac{(\;1+e^{-\frac{5\pi}{2\sqrt{2N}}}\;)^{2}}{1+e^{-\frac{5\pi}{\sqrt{2N}}}}>\frac{1}{2}.$

This theorem implies that until the time step $O(N)$, this walk has the
oscillation with the periodicity $\pi\sqrt{2N}$. Thus the locally maximal
value of the time sequence of $\\{\mu_{t}(u_{*})\\}_{t>0}$ can be observed
around $O(\sqrt{N})$ and its probability is higher than $1/2$.

In Fig. 1, we give a numerical simulation of the time courses of the finding
probabilities of each vertex for $N=100$.

Figure 1: The time courses of the finding probability at each positions for
$N=100$: The blue curve is the time course for the marked vertex. The others
are those of the non-marked vertices. The convergence to $1/2$ of the blue
curve derives from Theorem 1. The convergence time to $1/2$ is estimated by
Theorem 2; the final time in the numerical simulation is $[N\log N]=460$. The
first oscillation comes from Theorem 3; in the first half of the time course
of the target vertex, the peaks can be seen with the periodicity about
$\pi\sqrt{2N}\asymp 44$.

##  4 Proofs of theorems

###### Proof.

Step 1 (Reduction to an invariant subspace): Let us divide the arc set $A$
into the disjoint sets $A_{+}$, $A_{-}$, $A_{0}$, where

$\displaystyle A_{+}$ $\displaystyle=\\{a\in A\;:\;t(a)=u_{*}\\},$
$\displaystyle A_{-}$ $\displaystyle=\\{a\in A\;:\;o(a)=u_{*}\\},$
$\displaystyle A_{0}$ $\displaystyle=\\{a\in A\;:\;t(a),\;o(a)\neq u_{*}\\}.$

From the symmetricity of the time evolution, it is easy to notice that the
value $\psi_{t}(a)$ is invariant of the subsets $A_{+}$, $A_{-}$ and $A_{0}$.
Thus setting $a_{t},b_{t},c_{t}$ as the representative of $\psi_{t}(a)\;(a\in
A_{+})$, $\psi_{t}(b)\;(b\in A_{-})$ and $\psi_{t}(c)\;(c\in A_{0})$,
respectively, we obtain

$\begin{bmatrix}a_{t+1}\\\ b_{t+1}\\\
c_{t+1}\end{bmatrix}=\begin{bmatrix}0&-1+2/N&2-4/N\\\ -1+2/N&0&0\\\
0&2/N&1-4/N\end{bmatrix}\begin{bmatrix}a_{t}\\\ b_{t}\\\
c_{t}\end{bmatrix}+\frac{2}{N}\begin{bmatrix}1\\\ -1\\\ 1\end{bmatrix}$

from the definition of our quantum walk (2). Taking the normailization,

$\bm{\alpha}_{t}:=\begin{bmatrix}\sqrt{N-1}&0&0\\\ 0&\sqrt{N-1}&0\\\
0&0&\sqrt{(N-1)(N-2)}\end{bmatrix}\begin{bmatrix}a_{t}\\\ b_{t}\\\
c_{t}\end{bmatrix},$

we obtain the following master equation:

$\bm{\alpha}_{t+1}=T(\epsilon)\bm{\alpha}_{t}+\bm{b}_{\epsilon},\;\bm{\alpha}_{0}=0.$
(4)

Here

$T(\epsilon)=\begin{bmatrix}0&-1+\epsilon^{2}&\sqrt{2\epsilon^{2}(1-\epsilon^{2})}\\\
-1+\epsilon^{2}&0&0\\\
0&\sqrt{2\epsilon^{2}(1-\epsilon^{2})}&1-2\epsilon^{2}\end{bmatrix}$

and

$\bm{b}_{\epsilon}=\left[\epsilon\sqrt{2-\epsilon^{2}},\;-\epsilon\sqrt{2-\epsilon^{2}},\;\sqrt{4-6\epsilon^{2}+2\epsilon^{4}}\right]^{\top},$

where $\epsilon=\sqrt{2/N}$. Remark that the relative probability of the
marked vertex can be computed by $\nu_{t}(u_{*})=||\bm{\alpha}_{t}(1)||^{2}$.
Now let us consider the asymptotics for large $N$, which is equivalent to
considering (4) for small perturbation $\epsilon$.

Step 2 (Application of Kato’s perturbation theory): The matrix $T(\epsilon)$
can be expanded as follows:

$T(\epsilon)=T+\epsilon T^{(1)}+\epsilon^{2}T^{(2)}+\cdots,$

where the non-perturbed matrix is

$\displaystyle T$ $\displaystyle=\begin{bmatrix}0&-1&0\\\ -1&0&0\\\
0&0&2\end{bmatrix},$ (5)

which is symmetric and $T^{(1)}$ and $T^{(2)}$ are

$T^{(1)}=\begin{bmatrix}0&0&\sqrt{2}\\\ 0&0&0\\\
0&\sqrt{2}&0\end{bmatrix}\text{ and }T^{(2)}=\begin{bmatrix}0&1&0\\\ 1&0&0\\\
0&0&-2\end{bmatrix}.$

The set of eigenvalues of $T$ is

$\mathrm{Spec}(T)=\\{-1,+1,+1\\}.$ (6)

The eigenvalue $1$ has the multipilicity $2$, but it is semi-simple since $T$
is a symmetric matrix. Thus by a small perturbation $\epsilon$, the eigenvalue
$1$ is splitted into two parts $\lambda_{1,+}(\epsilon)$ and
$\lambda_{1,-}(\epsilon)$ while the simple eigenvalue $-1$ is just perturbed
which is denoted by $\lambda_{-1}(\epsilon)$. The eigenvalues of $T(\epsilon)$
can be expanded by

$\displaystyle\lambda_{-1}(\epsilon)$
$\displaystyle=-1+\epsilon\lambda_{-1}^{(1)}+\epsilon^{2}\lambda_{-1}^{(2)}+\cdots$
(7) $\displaystyle\lambda_{1,\pm}(\epsilon)$
$\displaystyle=1+\epsilon\lambda_{1,\pm}^{(1)}+o(\epsilon)+\cdots$ (8)

because $(-1)$-eigenvalue is simple while $(+1)$-eigenvalue is semi-simple.
Now let us compute the coefficients.

  1. 1.

Expansion of $\lambda_{-1}(\epsilon)$ : Since $-1$ is simple eigenvalue of the
non-perturbed matrix $T$, the weighted mean of the perturbed eigenvalues of
$T(\epsilon)$ around $-1$ is the perturbed eigenvalue $\lambda_{-1}(\epsilon)$
itself. Then following the formulas on the weighted average of the perturbed
eigenvalues [12], let us compute the first and second terms. According to
[(2.33) in Ch. II, Sect. 2.2 [12]], the first term and the second terms are
expressed by

$\displaystyle\lambda_{-1}^{(1)}=\tr(T^{(1)}P_{-1}),$ (9)
$\displaystyle\lambda_{-1}^{(2)}=\tr(T^{(2)}P_{-1}-T^{(1)}S_{-1}T^{(1)}P_{-1}),$
(10)

where $P_{\lambda}$ is the eigenprojection of $\ker(\lambda-T)$ and
$S_{\lambda}$ is the reduced resolvent at $\lambda$
($\lambda\in\mathrm{Spec}(T)$) which can be computed by

$P_{-1}=\begin{bmatrix}1/2&1/2&0\\\ 1/2&1/3&0\\\ 0&0&0\end{bmatrix}$

and

$\displaystyle S_{-1}$
$\displaystyle=-\lim_{\xi\to-1}\sum_{\lambda\in\mathrm{Spec}(T)\setminus\\{-1\\}}(\xi-\lambda)^{-1}P_{\lambda}=\frac{1}{2}(I-P_{-1})$
$\displaystyle=\frac{1}{2}\begin{bmatrix}1/2&-1/2&0\\\ -1/2&1/2&0\\\
0&0&1\end{bmatrix}.$

Then we have

$\lambda_{-1}^{(1)}=0,\text{ and }\lambda_{-1}^{(2)}=1/2.$

As a consequence, we obtain

$\lambda_{-1}(\epsilon)=-1+\frac{1}{2}\epsilon^{2}+O(\epsilon^{3}).$ (11)

  2. 2.

Expansion of $\lambda_{1,\pm}(\epsilon)$: Since $(+1)$-eigenvalue is semi-
simple, but has the multiplicity, we choose a different method from the case
for $\lambda_{-1}(\epsilon)$. Such a case, we consider

$\tilde{T}(\epsilon):=\frac{1}{\epsilon}(T(\epsilon)-1)P_{1}(\epsilon)$

instead of $T(\epsilon)$, which is so called the reduction process [12]. Here
$P_{1}(\epsilon)$ is the total projection for the $+1$-group. The new matrix
$\tilde{T}(\epsilon)$ can be expanded by

$\tilde{T}(\epsilon)={\tilde{T}}^{(1)}+\epsilon{\tilde{T}}^{(2)}+\cdots$

because $+1$ is semi-simple eigenvalue of $T$. Here $\tilde{T}^{(1)}$ and
$\tilde{T}^{(2)}$ are given by [(2.20) in Ch. II, Sect. 2.2 [12]]:

$\displaystyle\tilde{T}^{(1)}$
$\displaystyle=P_{1}T^{(1)}P_{1}=\frac{1}{\sqrt{2}}\begin{bmatrix}0&0&-1\\\
0&0&1\\\ 1&-1&0\end{bmatrix},$ $\displaystyle\tilde{T}^{(2)}$
$\displaystyle=P_{1}T^{(1)}P_{1}-P_{1}T^{(1)}P_{1}T^{(1)}S_{1}-P_{1}T^{(1)}S_{1}T^{(1)}P_{1}-S_{1}T^{(1)}P_{1}T^{(1)}P_{1}.$

The coefficient $\lambda_{1,\pm}^{(1)}$ can be obtained by [(2.40) in Ch. II,
Sect. 2-3]:

$\lambda_{1,\pm}^{(1)}\in\mathrm{Spec}\left(\tilde{T}^{(1)}|_{\mathrm{Ran}{(}P_{1})}\right)=\\{\pm
i\\}.$

Since $\tilde{T}^{(1)}$ is a skew-Hermitian matrix, the eigenprojections of
$\pm i$ can be directly computed by its normalized eigenvectors $\bm{v}_{1,\pm
i}$ by

$\bm{v}_{1,\pm i}=(1/2)\;[-i,\;+i,\;\sqrt{2}]^{\top}.$

Indeed we can check that

$P_{1}=\bm{v}_{1,+i}\bm{v}_{1,+i}^{*}+\bm{v}_{1,-i}\bm{v}_{1,-i}^{*}.$

We set $P_{1,\pm}:=v_{1,\pm}v_{1,\pm i}^{*}$ as a standard decomposition of
$\ker(1-T)$. Moreover, since the eigenvalues $\pm i$ are simple, the term
$o(\epsilon)$ of the expansion of $\lambda_{1,\pm}(\epsilon)$ can be expressed
by

$o(\epsilon)=\epsilon^{2}\lambda_{1,\pm}^{(2)}+O(\epsilon^{3})$

and the coefficient $\lambda_{1,\pm}^{(2)}$ can be obtained by using the
formula of the weighed mean of the perturbed eigenvalues again; that is,

$\lambda_{1,\pm}^{(2)}=\tr(\tilde{T}^{(1)}P_{1,\pm})=-5/4.$

After all, we conclude that

$\lambda_{1,\pm}(\epsilon)=1\pm
i\epsilon-\frac{5}{4}\epsilon^{2}+O(\epsilon^{3}).$ (12)

Step 3 (Feedback of the expansions to the expression of $\bm{\alpha}_{t}$):
By (4), we have

$\bm{\alpha}_{t}=\left(I+T(\epsilon)+T^{2}(\epsilon)+\cdots+T^{t-1}(\epsilon)\right)\bm{b}_{\epsilon}.$

Since $T$ is semi-simple, the perturbed matrix $T(\epsilon)$ can be decompose
as

$T(\epsilon)=\lambda_{-1}(\epsilon)P_{1}(\epsilon)+\lambda_{1,+}(\epsilon)P_{1,+}(\epsilon)+\lambda_{1,-}(\epsilon)P_{1,-}(\epsilon),$

where $P_{j}(\epsilon)$ is the eigenprojection onto the each perturbed
eigenvalues and can be expanded by

$P_{j}(\epsilon)=P_{j}+O(\epsilon).$

This implies

$\bm{\alpha}_{t}=\left(\frac{1-\lambda_{-1}^{t}(\epsilon)}{1-\lambda_{-1}(\epsilon)}P_{-1}(\epsilon)+\frac{1-\lambda_{1,+}^{t}(\epsilon)}{1-\lambda_{1,+}(\epsilon)}P_{1,+}(\epsilon)+\frac{1-\lambda_{1,-}^{t}(\epsilon)}{1-\lambda_{1,-}(\epsilon)}P_{1,-}(\epsilon)\right)\bm{b}_{\epsilon}.$
(13)

Inserting the expansions of $\lambda_{-1}(\epsilon)$,
$\lambda_{1,\pm}(\epsilon)$ in (11), (12) into (13), we have another
expression for $\bm{\alpha}_{t}$ as follows:

$\displaystyle\bm{\alpha}_{t}$
$\displaystyle=\frac{1}{2}\left\\{1-(-1)^{t}e^{-\frac{\epsilon^{2}t}{2}\;(1+O(\epsilon))}\right\\}e^{\frac{\epsilon^{2}}{4}(1+O(\epsilon))}\bm{b}_{\epsilon}^{(-1)}$
$\displaystyle+\frac{i}{\epsilon}\left\\{1-e^{i\epsilon
t}e^{-\frac{5}{4}\epsilon^{2}t(1+O(\epsilon))}\right\\}e^{\frac{5}{4}i\epsilon(1+O(\epsilon))}\bm{b}_{\epsilon}^{(1,+)}$
$\displaystyle-\frac{i}{\epsilon}\left\\{1-e^{-i\epsilon
t}e^{-\frac{5}{4}\epsilon^{2}t(1+O(\epsilon))}\right\\}e^{-\frac{5}{4}i\epsilon(1+O(\epsilon))}\bm{b}_{\epsilon}^{(1,-)}.$
(14)

Here $\bm{b}_{\epsilon}^{(j)}:=P_{j}(\epsilon)\bm{b}_{\epsilon}$. We notice
that $P_{-1}\bm{b}_{0}=\bm{0}$, which implies that there exists
$\bm{b}^{\prime}$ such that
$\bm{b}_{\epsilon}^{(-1)}=\epsilon(\bm{b}^{\prime}+O(\epsilon))$.

Step 4 (Finishing the proofs):
Proof of Theorem 1. Taking $t\to\infty$ in (13) and inserting the expansions
of $\lambda_{-1}(\epsilon)$ and $\lambda_{1,\pm}(\epsilon)$ in (11) and (12),
we have

$\displaystyle\bm{\alpha}_{\infty}$
$\displaystyle:=\lim_{t\to\infty}\bm{\alpha}_{t}$
$\displaystyle=\frac{1}{2}e^{\epsilon^{2}(1+O(\epsilon))/4}P_{-1}(\epsilon)\bm{b}_{\epsilon}+\frac{i}{\epsilon}e^{5i\epsilon(1+O(\epsilon))/4}P_{1,+}(\epsilon)\bm{b}_{\epsilon}-\frac{i}{\epsilon}e^{-5i\epsilon(1+O(\epsilon))/4}P_{1,-}(\epsilon)\bm{b}_{\epsilon}.$

Then

$\displaystyle\lim_{\epsilon\to 0}\epsilon\bm{\alpha}_{\infty}$
$\displaystyle=\lim_{\epsilon\to
0}(ie^{5i\epsilon/4}P_{1,+}(\epsilon)\bm{b}_{\epsilon}-ie^{-5i\epsilon/4}P_{1,-}(\epsilon)\bm{b}_{\epsilon})$
$\displaystyle=\lim_{\epsilon\to
0}2\mathrm{Re}\left(ie^{5i\epsilon/4}P_{1,+}\bm{b}_{\epsilon}\right)$
$\displaystyle=-2\mathrm{Im}(P_{1,+}\bm{b}_{0})$
$\displaystyle=\frac{-1}{\sqrt{2}}\begin{bmatrix}-i\\\ i\\\ 0\end{bmatrix}.$

Thus the finding probability at the target vertex $u_{*}$ is $1/2$ because the
relative probability is given by $|\bm{\alpha}_{\infty}(1)|^{2}$.

Proof of Theorem 2. By (4), we have

$\displaystyle||\bm{\alpha}_{\infty}-\bm{\alpha}_{t}||$
$\displaystyle=\frac{e^{-\frac{\epsilon^{2}t}{2}(1+o(1))}}{\epsilon}\;\bigg{|}\bigg{|}\frac{\epsilon^{2}}{2}(-1)^{t}(\bm{b}^{\prime}+o(1))+ie^{-\frac{3}{4}\epsilon^{2}t}(\bm{b}_{0}^{(1,+)}+o(1))-ie^{-\frac{3}{4}\epsilon^{2}t}(\bm{b}_{0}^{(1,-)}+o(1))\bigg{|}\bigg{|}.$
(15)

There is a conflict between $\epsilon^{2}$ and $e^{-3\epsilon^{2}t/4}$ in the
RHS. Note that $\epsilon^{2}\ll e^{-3\epsilon^{2}t/4}$ if and only if
$t<8|\log\epsilon|/(3\epsilon^{2})$. Let us see the lower bound of $t$ such
that $||\bm{\alpha}_{\infty}-\bm{\alpha}_{t}||<e^{-\theta}$ if
$t<8|\log\epsilon|/(3\epsilon^{2})$. By (15), there exists $m>0$ such that

$\displaystyle||\bm{\alpha}_{\infty}-\bm{\alpha}_{t}||$
$\displaystyle=(m+o(1))\frac{e^{-\frac{5\epsilon^{2}t}{4}(1+o(1))}}{\epsilon}.$
(16)

Then solving

$(m+o(1))\frac{e^{-\frac{5\epsilon^{2}t}{4}(1+o(1))}}{\epsilon}<e^{-\theta},$

we obtain the lower bound of $t$ by

$t>(\theta+\log(m+o(1))+|\log\epsilon|)\frac{4}{5\epsilon^{2}(1+o(1))}.$

Thus the lower bound of such a $t$ must belong to at least
$\Theta(|\log\epsilon|/\epsilon^{2})$. On the other hand, if $t\geq
8|\log\epsilon|/(3\epsilon)^{2}$, then (15) is rewritten by

$\displaystyle||\bm{\alpha}_{\infty}-\bm{\alpha}_{t}||$
$\displaystyle=(m^{\prime}+o(1))\epsilon
e^{-\frac{\epsilon^{2}t}{2}(1+o(1))},$ (17)

which implies

$t>(\theta+\log(m^{\prime}+o(1))-|\log\epsilon|)\frac{2}{\epsilon^{2}(1+o(1))}.$

Since $t>8|\log\epsilon|/(3\epsilon)^{2}$, this inequality is satisfied for
any fixed $\theta$ if $\epsilon$ is sufficiently small. Therefore we can
conclude that $t(\theta)\in\Theta(|\log\epsilon|/\epsilon^{2})=\Theta(N\log
N)$.

Proof of Theorem 3. Since $\bm{b}_{\epsilon}^{(-1)}\in O(\epsilon)$, if $t\ll
1/\epsilon^{2}$, the second and third terms in (13) are main terms; that is,

$\bm{\alpha}_{t}\sim\frac{i}{\epsilon}(1-e^{i\epsilon
t}e^{-\frac{5}{4}\epsilon
t^{2}})\bm{b}_{\epsilon}^{(1,+)}-\frac{i}{\epsilon}(1-e^{-i\epsilon
t}e^{-\frac{5}{4}\epsilon t^{2}})\bm{b}_{\epsilon}^{(1,-)}$

for sufficiently small $\epsilon$. Note that

$\displaystyle\bm{b}_{\epsilon}^{(1,+)}$
$\displaystyle=P_{1,+}\bm{b}_{\epsilon}=\langle\bm{v}_{1,i},\bm{b}_{0}\rangle\bm{v}_{1,i}+O(\epsilon)$
$\displaystyle=[-i/\sqrt{2},i/\sqrt{2},1]^{\top}+O(\epsilon).$
$\displaystyle\bm{b}_{\epsilon}^{(1,-)}$
$\displaystyle=P_{1,-}\bm{b}_{\epsilon}=\langle\bm{v}_{1,-i},\bm{b}_{0}\rangle\bm{v}_{1,-i}+O(\epsilon)$
$\displaystyle=[i/\sqrt{2},-i/\sqrt{2},1]^{\top}+O(\epsilon).$

Therefore until $t\in O(1/\epsilon)$, we have

$\displaystyle\bm{\alpha}_{t}\sim\frac{i}{\epsilon}\left[\;-\sqrt{2}(1-c_{t}\cos\epsilon
t),\;\sqrt{2}(1+c_{t}\cos\epsilon t),\;2\sin\epsilon t\;\right]^{\top},$

where $c_{t}=e^{-(5/4)\epsilon^{2}t}$. Then the normalized constant of the
relative probability can be computed by $4(1+c_{t}^{2})$. Then we have

$\nu_{t}(u_{*})\sim\frac{2(1-c_{t}\cos\epsilon t)^{2}}{4(1+c_{t}^{2})},$

for $t\in O(1/\epsilon)$. In particular, if $t=\pi/\epsilon$, then the finding
probability at the marked vertex takes the local maximal value

$\nu_{t}(u_{*})=\frac{1}{2}\frac{(1+e^{-5\pi/(2\sqrt{2N})})^{2}}{1+e^{-5\pi/\sqrt{2N}}}>1/2.$

∎

##  5 Summary

In this paper, we considered a quantum search driven by quantum walk with in-
and out-flows. We mark one-vertex in the complete graph. The internal graph
receives the inflow and releases the out flow at every time step. This quantum
walk model converges to a fixed point. We find that there are two phases of
the time evolution of this quantum walk: the first phase is the pulsation
phase until $O(N)$ and the send phase is the stationary phase after $O(N\log
N)$. We show that (i) the high probability of the marked vertex in the
stationary phase (Theorem 1); (ii) the mixing time is estimated by $O(N\log
N)$ (Theorem 2); (iii) until $O(N)$, there is a pulsation with the periodicity
$O(\sqrt{N})$ and we find the marked vertex with a high probability in this
pulsation phase (Theorem 3). By extending this model to general graphs, to
classify the graphs where such an aspect of a kind of a phase transition
appear is one of the interesting future’s problems.

Acknowledgments: M.S. acknowledges financial supports from JST SPRING (Grant
No. JPMJSP2114). Yu.H. acknowledges financial supports from the Grant-in-Aid
of Scientific Research (C) Japan Society for the Promotion of Science (Grant
No. 18K03401). E.S. acknowledges financial supports from the Grant-in-Aid of
Scientific Research (C) Japan Society for the Promotion of Science (Grant No.
19K03616) and Research Origin for Dressed Photon. The authors would like to
thank prof. Hajime Tanaka for the invaluable discussions and the comments
which were very useful to carry out this study.

## References

  * [1] A. Ambainis, Quantum walks and their algorithmic applications, International Journal of Quantum Information 1 (2003), pp.507-518.
  * [2] A. M. Childs, Universal computation by quantum walk, Physical Review Letter 102 (2009), 180501.
  * [3] R. Portugal, Quantum Walk and Search Algorithms, 2nd Ed., Springer Nature, Switzerland (2018).
  * [4] G. Brassard, Searching a quantum phone book, Science 275 (1997), 627.
  * [5] L. K. Grover, Fixed-point quantum search, Physical Review Letter 95 (2005), 15050.
  * [6] T. J. Yoder, G. H. Low and I. L. Chuang, Fixed-point quantum search with an optimal number of queries. Physical Review Letter 21 (2014), 210501.
  * [7] E. Feldman and M. Hillery, Quantum walks on graphs and quantum scattering theory, Coding Theory and Quantum Computing, edited by D. Evans, J. Holt, C. Jones, K. Klintworth, B. Parshall, O. Pfister, and H. Ward, Contemporary Mathematics 381 (2005), pp.71-96.
  * [8] E. Feldman and M. Hillery, Modifying quantum walks: A scattering theory approach, Journal of Physics A: Mathematical and Theoretical 40 (2007), 11319.
  * [9] S. Albeverio, F. Gesztesy, R. Høegh-Krohn, H. Holden and P. Exner, “Solvable Model in Quantum Mechanics”, AMS Chelsea publishing (2004).
  * [10] Yu. Higuchi and E. Segawa, Dynamical system induced by quantum walks, Journal of Physics A: Mathematical and Theoretical 52 (2009), 395202.
  * [11] S. Mohamed, Yu. Higuchi and E. Segawa, Electric circuit induced by quantum walks, Journal of Statistical Physics 181 (2020), pp.603–617.
  * [12] T. Kato, A Short Introduction to Perturbation Theory for Linear Operator, Springer-Verlag, New York (1982).


# High Temperature Superconductivity: A Simple Model Exploiting Hydrogen Bonds

Daniel Kaplan,1∗ Yoseph Imry,1

1Dept. of Condensed Matter Physics, Faculty of Physics, Weizmann Institute of
Science
Rehovot, 76100, Israel

∗To whom correspondence should be addressed; E-mail:
<EMAIL_ADDRESS>
> Lately, there has been much interest in high temperature superconductors,
> and more recently hydrogen-based superconductors. This work offers a simple
> model which explains the behavior of the superconducting gap based on BCS
> theory, and reproduces most effects seen in experiments, including the
> isotope effect and $T_{c}$ enhancement as a function of pressure. We show
> that this is due to a combination of the factors appearing in the gap
> equation: the matrix element between the proton states, and the level
> splitting of the proton.

## Introduction

Recently, there has been a surge of interest concerning the discovery of high
temperature superconductivity in hydrogen compounds (?, ?, ?). These showed
that metallization and superconductivity in such materials can be obtained by
changing the pressure of the system. Furthermore, the discovery of the inverse
isotope effect in hydrides (notably, PdH) (?) helped fuel interest in these
materials. Historically, arguments have been put forward that the maximal
value of $T_{c}$ in BCS theory is roughly $T_{c}\sim 30K$ (?); but since then
higher values of $T_{c}$ were demonstrated, e.g. in the discovery of the
cuperates, pnictides and $MgB_{2}$ (?). Suggested first by N. Ashcroft (?),
metallic hydrogen (and hydrogen compounds) can be a platform for high
temperature superconductivity. This is accomplished by maximizing the $N_{0}U$
parameter in the BCS gap equation, while increasing the energy transfer to the
phonon $\Delta E$. In this paper, we present a model which accounts for these
properties in hydrogen sulfide, and gives a qualitative explanation for the
enhancement of $T_{c}$, the isotope effect in these compounds, and the reason
for pressure dependence in enhancing (or weakening) $T_{c}$.

## The Harmonic oscillator case

We first consider the emergence of a superconducting gap by assuming that
electrons interact with a harmonic oscillator,

$\displaystyle H=\frac{p^{2}}{2M}+\frac{M\omega_{0}^{2}x^{2}}{2}$ (1)

With well-known solutions for wavefunctions and eigen-energies. We consider
the usual virtual transitions between the ground-state and a phonon-excited
state, which have an energy difference $\Delta E=\hbar\omega_{0}$, and another
eletron de-excites that state. Following (?) for the calculation of the full
attractive interaction matrix element.

$\displaystyle\left\langle
1|H_{int}|2\right\rangle=2\frac{|W|^{2}_{q}}{\hbar}\left(\frac{\omega_{q}}{\omega^{2}-\omega_{q}^{2}}\right)$
(2)

Here the states $\left.|1\right\rangle$, $\left.|2\right\rangle$ correspond to
two Cooper pairs. One with incoming momentum pairs $(k,-k)$, and the other
$(k^{\prime},-k^{\prime})$. $\vec{q}=\vec{k}-\vec{k}^{\prime}$. $W_{q}$ is the
relevant electron-phonon matrix element, for the emission or absorption of a
phonon with wave-vector $\vec{q}$. $\omega_{q}$ is the energy of the phonon,
and $\omega=|E(1)-E(2)|$ is the energy difference between states (1 and 2,
respectively). We will consider the particle occupying the H.O. as “the
proton”. In the lattice, each energy level becomes a narrow band, with states
characterized by a wave-vector $\vec{q}$.

### Electron-proton coupling

We explore two forms of this term. For simplicity, we show this in 1D but the
approach is easily generalized to 3D. Firstly, we approximate the electron-
proton interaction by a delta function, with the form $V_{int}=-g\delta(x-y)$,
such that $x$ is the coordinate of the oscillator, and $y$ is the coordinate
of the electron. Suppose that an electron excites this H.O. from its ground
state to the first excited state. The initial state of the
electron,$\left.|i\right\rangle=\frac{1}{\sqrt{2\pi}}\exp(-iky)$ and the final
state (after interacting with the
oscillator),$\left.|f\right\rangle=\frac{1}{\sqrt{2\pi}}\exp(-ik^{\prime}y)$.
All in all the total state of the system is expressed by the notation
$\left.|{0,i}\right\rangle$ (where the $0$ refers to the ground state of the
oscillator) and $\left.|{1,f}\right\rangle$, where the $1$ now indicates the
first excited state. Then the matrix element is,

$\displaystyle\begin{matrix}\displaystyle W_{q}=\left\langle
0,i|V_{int}|1,f\right\rangle=-\frac{g}{2\pi}\iint
dxdy\psi_{0}^{*}(x)e^{\left(-i(k-k^{\prime})y\right)}\delta(x-y)\psi_{1}(x)=\\\
\displaystyle-\frac{g}{2\pi}\int
dx\psi_{0}^{*}(0)\psi_{1}(x)e^{-iqx}\end{matrix}$ (3)

Where we denoted $k-k^{\prime}=q$, i.e., the momentum transfer to the phonon.
$\psi_{0},\psi_{1}$ are the ground-state, first excited state wavefunctions,
respectively. For a harmonic oscillator this result is known exactly. We have,

$\displaystyle W_{q}=\frac{igqa}{2\pi}e^{-\frac{q^{2}a^{2}}{4}}$ (4)

Where $a^{2}=\frac{\hbar}{M\omega_{0}}$ is the length scale of the oscillator.
In our case, $qa\ll 1$, and a reasonable value for $g$ above is $g\sim
2eV\AA$. We also consider a screened-Coulomb interaction. For this form,
$V_{int}=-\frac{Ze^{2}}{4\pi}\frac{e^{-k_{s}r}}{r}$, where $k_{s}$ is the
screening length, the matrix element of Eq. 3 becomes,

$\displaystyle\left\langle
0,i|V_{int}|1,f\right\rangle=-\frac{Ze^{2}}{8\pi^{2}}\int
dx\psi_{0}^{*}(0)\psi_{1}(x)\frac{e^{-iqx}}{q^{2}+k_{s}^{2}}$ (5)

For H2S, in the metallic phase, $k_{s}\sim 1\frac{1}{\AA}$ (?).

### The BCS superconducting gap

We focus on a particular longitudinal phonon mode, which, by work in (?), has
been shown to be particularly dominant in these materials. Consider the matrix
element of Eq. 2. Substituting $W_{q}$ from Eq. 4 one finds,

$\displaystyle U=\left\langle
1|H_{int}|2\right\rangle=\frac{g^{2}q^{2}e^{-\frac{q^{2}\hbar}{2m\omega_{0}}}}{2m\omega_{0}}\times\frac{\omega_{q}}{\omega^{2}-\omega_{q}^{2}}$
(6)

We make the simplifying assumption that the processes we’re interested in
occur at $\omega=0$. We base this on the fact the probability of inelastic
processes where $\omega\neq 0$ is relatively small (which is implied by the
$qa\ll 1$ limit) (?). Finally, we are left with,

$\displaystyle
U=-2\frac{g^{2}q^{2}a^{2}e^{-\frac{q^{2}a^{2}}{2}}}{\hbar\omega_{0}}$ (7)

This form of the coupling is also found in (?). We then replaced
$\omega_{q}=\omega_{0}$ since the phonon in this case is represented by the
H.O, and the relevant branch has a very small dispersion in $\vec{q}$. The
plot below shows $U$ as a function of $\omega_{0}$ for a given $q$. Typical
$q$ is $q\sim 0.1\frac{1}{A}$. $g\sim 2eV$, M is the proton mass.. The
relevant dimensionless parameter is the BCS interaction parameter, which we
take to be $N_{0}U$, where $N_{0}$ is the electronic density at the Fermi
energy. For $H_{2}S$, we have $N_{0}\sim 0.3eV^{-1}$ (?).

Figure 1: $N_{0}|U|$ as a function of $\omega_{0}$

The BCS gap has the form (?) (?),

$\displaystyle\Delta=\hbar\omega_{0}e^{-\frac{1}{N_{0}|U|}}$ (8)

Note that the prefactor to the exponential is an increasing function of
$\omega_{0}$, whereas the matrix element, per Fig. 1 is a decreasing one.
Clearly, a maximum might exist, This is shown below.

Figure 2: $\Delta$ as a function of $\omega_{0}$, for $q\sim 0.1\frac{1}{A}$,
according to Eq. 8

Using results of the gap equation at finite temperatures (see below), the
maximum in Fig. 2 appears at $\omega_{0}=22meV$, and corresponds to
$T_{c}\approx 92K$.

### The Isotope effect

Having found $\Delta$ it is possible to investigate the dependence of the gap
on the Mass of the proton. We calculate,

$\displaystyle\frac{1}{\Delta}\frac{d\Delta}{dM}\sim-\frac{\omega_{0}^{2}}{q^{2}}\left(1-\frac{q^{2}a^{2}}{2}\right)$
(9)

Recall that $a^{2}\sim\frac{1}{M\omega_{0}}$, this shows, as suggested by Fig.
2, one may obtain either a positive, or negative isotope effect, depending on
the sign of $1-q^{2}a^{2}/2$. In the range $qa\ll 1$ the gap changes
negatively (i.e. decreases) with mass. This is confirmed by (?, ?), and
experimentally observed.

### Gap equation at finite temperature

For completeness, we include the gap equation at finite temperature which
permits an accurate determination of $T_{c}$. BCS theory puts
$T_{c}\sim\Delta$ (up to a numerical factor of order 1) (?). The BCS result is
(using our simplified approach),

$\displaystyle 1=N_{0}U\int\limits_{0}^{\Delta
E}\frac{dE}{E^{2}+\Delta^{2}}\left(1-2f(\sqrt{E^{2}+\Delta^{2}})\right)$ (10)

Where $f$ is the Fermi-Dirac distribution. At the transition temperature
$T=T_{c}$, the gap $\Delta=0$, yielding,

$\displaystyle T_{c}\sim\frac{\Delta}{1.75}$ (11)

## A simple new model of the superconducting transition in $H_{2}S$

Numerical as well as experimental work on the structure of $H_{2}S$ have
confirmed that the material experiences a transition of crystalline-structure
upon becoming superconducting (?). In the $Im\bar{3}m$ phase of hydrogen
sulfide, shown below, the hydrogen bond between $S$ atoms is symmetrized such
that the body-centered phase allows for superconductivity.

Figure 3: The body-centered structure of $H_{3}S$, taken from (?)

We model the potential seen by the proton using a sum of two Lennard Jones
(L-J) potentials
$V=4\epsilon\left(\left(\frac{\sigma}{r}\right)^{6}-\left(\frac{\sigma}{r}\right)^{12}\right)$.
Consider placing two $S$ atoms at co-linear points separated by $x_{0}$ in
space (for simplicity, the dynamics are on the $x$ axis). This was observed,
experimentally (in hydrogen-bonded compounds) (?). Then, the potential an $H$
atom experiences is,

$\displaystyle
V=4\epsilon\left(\left(\frac{\sigma}{(x-x_{0})}\right)^{6}-\left(\frac{\sigma}{(x-x_{0})}\right)^{12}\right)+4\epsilon\left(\left(\frac{\sigma}{(x+x_{0})}\right)^{6}-\left(\frac{\sigma}{(x+x_{0})}\right)^{12}\right)$
(12)

The parameter governing the form of the potential is then $x_{0}$ – such that
$2x_{0}$ is the separation of $S$ atoms. Below are several forms of the
potential for different values of $x_{0}$.

Figure 4: Various forms of the potential in Eq. 12. Note the transition from
single-well to double-well, occurring at $x_{0}\approx 1.2\sigma$

Typical values for $H_{2}S$ are $\epsilon\sim 300K$,$\sigma\sim 2.0A$ (?).
After diagonalizing the Hamiltonian to find the energy levels and
eigenfunctions, required by the procedure in Eq. 3. The figure below shows the
first two levels’ splitting, as a function of $x_{0}$,

Figure 5: $\Delta E_{12}$ (energy difference between ground state and first
excited state). $\sigma\sim 2.0A$, $\epsilon\sim 300K$

Using Eqs. 2,8, we can now obtain the superconducting gap. We replace above,
in Eq. 8, $\omega_{0}$ with the energy difference $\Delta E\equiv\Delta
E(x_{0})$. The usual solution is then,

$\displaystyle\Delta=(\Delta E)e^{-\frac{1}{N_{0}U}}$ (13)

Where $N_{0}\sim 0.3\frac{1}{eV}$. Finally, the plot below shows all the
relevant factors that enter Eq. 13. The red line depicts the energy difference
$\Delta E$. The brown line is the matrix element squared as it appears in Eq.
3. The combination of these, detailed above gives the gap, which is plotted in
blue.

Figure 6: Numerical results for the L-J potential. Red is the energy
difference between the first two levels, brown is the electron-proton
interaction matrix element squared, blue is the resultant gap

The figure above clearly illustrates that $\Delta$ possesses a maximum. This
is obtained at $x_{0}\approx 1.82A$, with a value of $\Delta\approx 22meV$,
giving a $T_{c}\approx 150K$. This is a considerable enhancement over the
harmonic case (shown before), and suggests that (relatively) high $T_{c}$ can
be obtained in this model. This follows from the large matrix element,
resulting from the separation of the two minima. The change in the separation
may be caused by pressure on the system which modifies the positions of the
$S$ atoms.

## Conclusions

This work demonstrates qualitatively that $T_{c}$ may be changed by altering
the hydrogen bonds, in hydrogen compounds, using a simple model of $H_{2}S$.
We have shown that changes in $x_{0}$ (half $S$ atoms separation in the
model), due to changes in pressure affect the value of the superconducting
gap. This gap is optimized at a specific value of $x_{0}$, which is determined
by $P$ (pressure of the system), since $\Delta$ decreases for too high values
of $x_{0}$. This is both theoretically and experimentally confirmed (?, ?).
Our model gives a qualitative explanation of the isotope effect (both normal,
and inverse), depending on the value of $x_{0}$, or $\omega_{0}$ (in the
H.O.). More accurate determination of the parameters of the $H-S$ potential,
can, we believe, allow to fully account for $\Delta$ and its dependence on
physical parameters. At present, very high pressures (of the order of 100 GPa)
are needed to bring about the superconductivity in these compounds. However,
judicious Chemistry with larger and smaller ions might create internal
pressures in the material to simulate such conditions.

## Acknowledgements

Y.I. acknowledges the late Prof. I Pelah for conversations on hydrogen bonds,
and the late Prof. W Kohn on superconductivity and DFT. D.K. thanks Or Ben Zvi
for his immeasurable support and attention. Financial support from the
Weizmann Institute of Science is gratefully acknowledged.

## References

  * 1. I. Errea, M. Calandra, C. J. Pickard, J. R. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Nature 532, 81 (2016).
  * 2. A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, and S. I. Shylin, Nature 525, 73 (2015), arXiv:1412.0460 .
  * 3. J. A. Flores-Livas, A. Sanna, and E. K. Gross, The European Physical Journal B  89, 63 (2016), arXiv:1501.06336 .
  * 4. T. Skoskiewicz, A. W. Szafranski, W. Bujnowski, and B. Baranowski, Journal of Physics C: Solid State Physics 7, 2670 (1974).
  * 5. J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Physical Review 108, 1175 (1957).
  * 6. P. W. Anderson, “bibfield journal “bibinfo journal Journal of Physics: Conference Series“ “textbf “bibinfo volume 449,“ “bibinfo pages 012001 (“bibinfo year 2013).
  * 7. N. W. Ashcroft, “bibfield journal “bibinfo journal Physical Review Letters“ “textbf “bibinfo volume 21,“ “bibinfo pages 1748 (“bibinfo year 1968).
  * 8. P.-G. De Gennes,  _Superconductivity of Metals and Alloys_ , Vol. 1 (1966).
  * 9. S. Havriliak, R. W. Swenson, and R. H. Cole, “bibfield journal “bibinfo journal The Journal of Chemical Physics“ “textbf “bibinfo volume 23,“ “bibinfo pages 134 (“bibinfo year 1955).
  * 10. Y. Imry, in “emph “bibinfo booktitle Tunneling Phenomena in Solids (Springer US, Boston, MA, 1969) pp. 563–576.
  * 11. N. A. Kudryashov, A. A. Kutukov, and E. A. Mazur, “bibfield journal “bibinfo journal Physics of Metals and Metallography“ “textbf “bibinfo volume 118,“ “bibinfo pages 113 (“bibinfo year 2017).
  * 12. J. Bardeen and D. Pines, Physical Review 99, 1140 (1955).
  * 13. Y. Ge, F. Zhang, and Y. Yao, “bibfield journal “bibinfo journal Physical Review B“ “textbf “bibinfo volume 93,“ “bibinfo pages 224513 (“bibinfo year 2016), arXiv:1507.08525 .
  * 14. Y. Imry, I. Pelah, and E. Wiener, “bibfield journal “bibinfo journal The Journal of Chemical Physics“ “textbf “bibinfo volume 43,“ “bibinfo pages 2332 (“bibinfo year 1965).
  * 15. G. Galliero and C. Boned, “bibfield journal “bibinfo journal Fluid Phase Equilibria“ “textbf “bibinfo volume 269,“ “bibinfo pages 19 (“bibinfo year 2008).
  * 16. A. F. Goncharov, S. S. Lobanov, V. B. Prakapenka, and E. Greenberg, “bibfield journal “bibinfo journal Physical Review B“ “textbf “bibinfo volume 95,“ “bibinfo pages 140101 (“bibinfo year 2017).
  * 17. M. I. Eremets and A. P. Drozdov, “bibfield journal “bibinfo journal Physics-Uspekhi“ “textbf “bibinfo volume 59,“ “bibinfo pages 1154 (“bibinfo year 2016).


[a]E. Garnacho Velasco

# Anomalous transport phenomena on the lattice

B. B. Brandt  F. Cuteri  G. Endrődi  G. Markó

###### Abstract

The interrelation between quantum anomalies and electromagnetic fields leads
to a series of non-dissipative transport effects in QCD. In this work we study
anomalous transport phenomena with lattice QCD simulations using improved
staggered quarks with physical masses in the presence of a background magnetic
field. In particular, we focus on the chiral separation effect and calculate
the corresponding conductivity both for free quarks and in the interacting
case, analyzing the dependence on several parameters, such as the temperature
and the quark mass.

##  1 Introduction

The topological sector of Quantum Chromodynamics (QCD) has been the subject of
intense study for more than 40 years, yielding a series of very interesting
phenomena associated with it. Some of these effects are called anomalous
transport phenomena, as they arise from the interplay between quantum
anomalies and electromagnetic fields or vorticities. Topology and quantum
anomalies are related via the index theorem, so the manifestation of these
non-dissipative transport effects provides an exciting opportunity to probe
the non-trivial topological structure of QCD.

One of the most celebrated among these phenomena is the Chiral Magnetic Effect
(CME) [1]. The CME has been detected in condensed matter systems [2] and, what
is more relevant for our discipline, is actively sought for in heavy ion
collision experiments. The most recent result from the STAR collaboration [3]
could not find a signal of CME in a dedicated run in RHIC with isobar
collisions. Although the results are still under discussion [4], understanding
this apparent suppression of such effects is one of the motivations behind
this contribution.

However, in this work we will focus on another anomalous transport phenomenon,
the Chiral Separation Effect (CSE) [5, 6]. The objective is to use lattice QCD
simulations to determine the CSE conductivity $C_{\text{CSE}}$ and analyze its
dependence on relevant parameters, like the temperature $T$ or the mass $m$ of
the quarks. This conductivity has been calculated analytically only for free
quarks, so the physical setup involving full QCD simulations will certainly
help understand how this effect appears in realistic situations. The present
effort will also be useful for studying further anomalous transport effects
like the CME, for which most lattice simulations so far [7, 8, 9, 10, 11, 12]
were either based on indirect approaches or are yet to be performed in full
QCD at the physical point.

This contribution is organized as follows: in section 2 we review the main
features of anomalous transport phenomena. In section 3 we present our setup
on the lattice and the techniques used to calculate the conductivity. In
section 4 we present our results for $C_{\text{CSE}}$ both in the free case
and in full QCD. Finally, we present our conclusions in section 5.

##  2 Anomalous transport phenomena

As it is widely known, QCD is observed to possess CP-symmetry, or
equivalently, the $\theta$ parameter is experimentally bound to be practically
zero. This fine-tuning issue is known as the strong CP-problem. One of the
most interesting characteristics of anomalous transport effects is that they
represent a local or event-by-event CP-violation in QCD.

In this theory, anomalous transport phenomena arise from the
$\mathrm{U}_{\text{A}}(1)$ anomaly. The non-trivial topological solutions of
QCD can transfer chirality to quarks through this anomaly, and it is this
origin that explains its CP-odd nature. Globally, the expectation value of
these anomalous currents has to vanish since the topology has to be trivial in
the system with zero topological charge
$\quantity(\expectationvalue{Q_{\text{top}}}=0)$. However, localized regions
with non-trivial topology may exist (for example instantons). In these
topological environments anomalous phenomena can arise.

Figure 1: Schematic representation of the CSE in a system with only $u$ and
$\bar{u}$ massless quarks.

As already mentioned above, the prime example of these effects is the CME: the
generation of a vector current in the presence of a chiral imbalance and a
magnetic field $B$. Another example is the CSE, which can be thought of as a
“dual” of the CME. The CSE is the emergence of an axial current in the
presence of finite density and a magnetic field. To first order, the current
is linear in the magnetic field (assumed to point in the $z$ direction,) and
in the baryon chemical potential $\mu$

$J_{3}^{5}=\sigma_{\text{CSE}}\,eB=C_{\text{CSE}}\,\mu\,eB+\mathcal{O}(\mu^{3})$
(1)

with $J_{3}^{5}$ the $z$-component of the axial current

$J_{3}^{5}=\int\differential^{4}x\,\bar{\psi}(x)\gamma_{3}\gamma_{5}\psi(x).$
(2)

This expression holds for a system consisting of a single fermion with charge
$e$ and no color.

Since this is the effect we will be focusing on, let us understand intuitively
how this phenomenon works. In Fig. 1 we can see a sketch of how the CSE
current can appear. Consider a situation with only massless $u,\bar{u}$ quarks
for simplicity (0). When a strong magnetic field is applied (1), all the
particles are in the lowest Landau level and their spins align with the
direction of $B$ (anti-align if the charge of the particle is negative). At
finite density (2), there is a net surplus of particles (supposing $\mu>0$)
and right-handed and left-handed particles would flow in different directions,
separating particles with different chiralities and thus creating an axial
current, which is precisely the CSE current (3).

For non-interacting quarks, an analytical treatment of this problem gives [5,
6]

$\sigma_{\text{CSE}}=\dfrac{1}{4\pi^{2}}n_{m}(T,\mu)$ (3)

where

$n_{m}(T,\mu)=\int^{\infty}_{-\infty}\differential
p_{3}\quantity[n\quantity(\sqrt{p_{3}^{2}+m^{2}})+n\quantity(-\sqrt{p_{3}^{2}+m^{2}})]$
(4)

and $n(E)$ is the usual Fermi-Dirac distribution

$n(E)=\dfrac{1}{\exp{(E-\mu)/T}+1}.$ (5)

In the case of massless quarks (or, equivalently, asymptotically high
temperatures), this expression simplifies to

$C_{\text{CSE}}=\dfrac{1}{2\pi^{2}}.$ (6)

One of the long-standing questions raised about the CSE (and in general about
anomalous transport effects) is how this coefficient is modified by
interactions. In the first years after the proposal of anomalous transport
phenomena, it was believed that the conductivity did not change in the full
theory, as it was topologically protected and hence was not modified by
gluonic interactions. However, later works are pointing to corrections of this
coefficient in the interacting theory. Since these corrections are non-
perturbative, lattice QCD provides an optimal tool to study these effects [7,
8, 9, 10, 11, 12].

Figure 2: Numerical value of $C_{\text{CSE}}$ as a function of $m/T$ using Eq.
(3). $C_{\text{CSE}}$ is the numerical coefficient of the first order
contribution in $B$ and $\mu$ to the current, so in terms of Eq. (3) it can be
defined as $\frac{1}{4\pi^{2}}\evaluated{\derivative{n_{m}}{\mu}}_{\mu=0}$.

##  3 Lattice setup

Simulations at finite real $\mu$ suffer from the infamous sign problem. That
is why former efforts to calculate $C_{\text{CSE}}$ in the full theory with
lattice QCD relied on special setups where this issue does not appear, for
example in the quenched theory [13], where no significant corrections were
found, or in two-color QCD [14], where at high $T$ the conductivity approached
the free case result and for low $T$ CSE was found to be suppressed.

In this work, we take a different approach. Rather than simulating at finite
$\mu$ and $B$, we can measure derivatives of the CSE current with respect to
the baryon chemical potential. Using Eq. (1), the derivative yields

$\evaluated{\derivative{\expectationvalue{J^{5}_{3}}}{\mu}}_{\mu=0}=C_{\text{CSE}}\,eB.$
(7)

Employing this, leading-order Taylor expansion of the current only requires
simulations at $\mu=0$, free of the sign problem. Then we can take a numerical
derivative (linear fit) of the result with respect to $eB$ to obtain
$C_{\text{CSE}}$.

As mentioned above, the last expression holds for a fermion of charge $e$ and
no color. For a system with several quark flavors with charge
$q_{f}=\tilde{q}_{f}e$ and $N_{c}$ colors, the formula is modified to

$\evaluated{\derivative{\expectationvalue{J^{5}_{3}}}{\mu}}_{\mu=0}=C_{\text{CSE}}\,N_{c}\sum_{f}\tilde{q}^{2}_{f}\,eB\equiv
C_{\text{CSE}}C_{\text{dof}}\,eB$ (8)

Figure 3: Derivative of the CSE (axial) current with respect to $\mu$ as a
function of the magnetic field for a $24^{3}\times 6$ lattice in the free case
(left) and in full QCD with $2+1$ flavors and physical quark masses at $T=305$
MeV (right). As expected, the behavior is linear in both cases and the slope
of the fit gives the value of $C_{\text{CSE}}$.

with

$J_{3}^{5}=\sum_{f}\tilde{q}_{f}\int\differential^{4}x\,\bar{\psi}_{f}(x)\gamma_{3}\gamma_{5}\psi_{f}(x).$
(9)

These are overall factors that can always be restored, so from this point we
rescale all our results by $C_{\text{dof}}$.

##  4 Results for CSE

The measurement of the current derivative in the (rooted) staggered
formulation involves a usual connected and a disconnected term, plus an
additional term coming from the derivative of the staggered Dirac matrices
$\Gamma_{\mu}$ with respect to $\mu$

$\begin{split}\evaluated{\derivative{\expectationvalue{J^{5}_{z}}}{\mu}}_{\mu=0}=\dfrac{T}{V}\Bigg{[}\dfrac{1}{4}&\expectationvalue{\text{Tr}\quantity(\Gamma_{4}M^{-1})\text{Tr}\quantity(\Gamma_{3}\Gamma_{5}M^{-1})}_{\mu=0}\\\
-\dfrac{1}{16}&\expectationvalue{\text{Tr}\quantity(\Gamma_{4}M^{-1}\Gamma_{3}\Gamma_{5}M^{-1})}_{\mu=0}\\\
+\dfrac{1}{4}&\expectationvalue{\text{Tr}\quantity(\dfrac{\partial(\Gamma_{3}\Gamma_{5})}{\partial\mu}M^{-1})}_{\mu=0}\Bigg{]}\end{split}$
(10)

where $M=\not{D}+m$, $V=(aN_{s})^{3}$ is the spatial volume, $T=(aN_{t})^{-1}$
is the temperature, with $N_{s}$ and $N_{t}$ the number of spatial and
temporal points of the lattice respectively.111Note that the last term in (10)
arises because the staggered discretization of the Dirac matrices (see, e.g.,
Ref. [15]) contain links that explicitly depend on the chemical potential. We
have measured the three terms at physical quark masses [16] in an already
existing ensemble of configurations for different magnetic fields [17, 18]. In
Fig. 3 we show an example of the obtained results. In this plot, we can see
that the expected linear behavior is confirmed. However, we are mostly
interested in the quantitative result of the slope, so that is what we will be
presenting next, both in the free case and in the interacting theory.

Figure 4: Results for $C_{\text{CSE}}$ with free staggered quarks at different
values of $m/T$ and $LT$. The red dashed line represents the numerical value
of $C_{\text{CSE}}$ at the given $m/T$ according to Fig. 2, while the black
dashed line represents $1/2\pi^{2}$, the value for the analytical calculation
of the same coefficient in the massless limit.

Let us start with the free case. In this particular scenario, the sign problem
is not present, but we will still use the same approach presented before since
we can check the consistency of our setup by comparing to the analytical
prediction from Eq. (3). In Fig. 4 we show the free case results for different
$m/T$ values, in each plot a continuum limit is taken for different aspect
ratios $LT$ by increasing $N_{s}$ and $N_{t}$ (with $LT$ and $m/T$ kept
constant). For small values of $m/T$, we can see a divergence when $LT$ goes
to zero. This behavior disappears when we go towards higher values of $m/T$.
This tells us finite size effects at $LT\rightarrow 0$ are sizeable if $m/T$
is not large enough. It is also worth noting the importance of the continuum
limit, since, as can be seen in the plots, non-continuum extrapolated results
would underestimate $C_{\text{CSE}}$ if the value of $N_{t}$ is not large
enough. Finally, the main result is that $C_{\text{CSE}}$ approaches the value
given by Eq. (3) when $LT\rightarrow\infty$ for every value of $m/T$. This
serves as a cross-check of our setup since we can reproduce the analytical
result for $C_{\text{CSE}}$ with our lattice simulations in the free case.

Figure 5: Results for $C_{\text{CSE}}$ with $2+1$ flavors of staggered quarks
at physical masses for a wide range of temperatures and four different lattice
sizes. The black dashed line represents the analytical prediction for the free
case with massless quarks.

Finally, we present our main result, the conductivity $C_{\text{CSE}}$ in full
QCD, in particular for $N_{f}=2+1$ flavors of staggered fermions at physical
quark masses. This is the first fully non-perturbative result for
$C_{\text{CSE}}$ at the physical point. Fig. 5 demonstrates the dependence of
the conductivity on the temperature for several finite-temperature lattice
ensembles $24^{3}\times 6$, $24^{3}\times 8$, $28^{3}\times 10$ as well as a
zero-temperature ensemble $24^{3}\times 32$. For temperatures above the QCD
transition temperature $T_{c}$, $C_{\text{CSE}}$ is found to approach the free
case prediction. This is in accordance with the expectation based on
asymptotic freedom, i.e. that at high $T$, QCD approaches a gas of quasi-free
quarks and gluons. At temperatures below $T_{c}$, the coefficient decreases
until reaching zero. This indicates a suppression of the CSE at low
temperatures, which is consistent with a previous study in two-color QCD [14]
and can be understood, at a qualitative level, using chiral effective theories
[19]. Although a proper continuum and thermodynamic limit is yet to be taken,
our results include different lattices spacings and different volumes and the
dependence is observed to be minor, so the (qualitative) behavior of
$C_{\text{CSE}}$ is not expected to change after taking these limits.

##  5 Summary and Outlook

In this proceedings article, we have presented a study of the CSE using
lattice QCD. In particular, we have calculated the conductivity
$C_{\text{CSE}}$ using dynamical staggered quarks both in the absence of
gluonic interactions and in full QCD at the physical point. In the free case,
we have recovered the expected value given in Eq. (3) once the continuum limit
is taken and finite size effects are under control. Having cross-checked our
setup, we moved on to full QCD, where we have determined the dependence of
$C_{\text{CSE}}$ on the temperature at physical quark masses. At temperatures
higher than $T_{c}$, the coefficient approaches the free case prediction as
expected, while at temperatures below $T_{c}$, the conductivity goes to zero,
showing that CSE might be very suppressed at low $T$. This is the first study
of this dependence at physical quark masses and it agrees with some of the
expectations based on investigations in two-color QCD. The next steps in this
project will be to analyze the mass dependence of $C_{\text{CSE}}$ in the
interacting case, as well as to obtain continuum extrapolated results for all
the cases. Finally, our technique can be generalized to study the CME, which
can contribute to a better theoretical understanding of anomalous transport
phenomena, as a counterpart to the experimental efforts being made to detect
this effect.

## Acknowledgments

This research was funded by the DFG (Collaborative Research Center CRC-TR 211
“Strong-interaction matter under extreme conditions” - project number
315477589 - TRR 211) and by the Helmholtz Graduate School for Hadron and Ion
Research (HGS-HIRe for FAIR). The authors are grateful for inspiring
discussions with Pavel Buividovich, Kenji Fukushima, Dirk Rischke, Sören
Schlichting, Igor Shovkovy and Lorenz von Smekal.

## References

  * [1] K. Fukushima, D. E. Kharzeev, and H. J. Warringa, “The Chiral Magnetic Effect,” Phys. Rev. D 78 (2008) 074033, arXiv:0808.3382 [hep-ph].
  * [2] Q. Li, D. E. Kharzeev, C. Zhang, Y. Huang, I. Pletikosic, A. V. Fedorov, R. D. Zhong, J. A. Schneeloch, G. D. Gu, and T. Valla, “Observation of the chiral magnetic effect in ZrTe5,” Nature Phys. 12 (2016) 550–554, arXiv:1412.6543 [cond-mat.str-el].
  * [3] STAR Collaboration, M. Abdallah et al., “Search for the chiral magnetic effect with isobar collisions at $\sqrt{s_{NN}}$=200 GeV by the STAR Collaboration at the BNL Relativistic Heavy Ion Collider,” Phys. Rev. C 105 no. 1, (2022) 014901, arXiv:2109.00131 [nucl-ex].
  * [4] D. E. Kharzeev, J. Liao, and S. Shi, “Implications of the isobar run results for chiral magnetic effect in heavy ion collisions,” arXiv:2205.00120 [nucl-th].
  * [5] D. T. Son and A. R. Zhitnitsky, “Quantum anomalies in dense matter,” Phys. Rev. D 70 (2004) 074018, arXiv:hep-ph/0405216.
  * [6] M. A. Metlitski and A. R. Zhitnitsky, “Anomalous axion interactions and topological currents in dense matter,” Phys. Rev. D 72 (2005) 045011, arXiv:hep-ph/0505072.
  * [7] P. V. Buividovich, M. N. Chernodub, E. V. Luschevskaya, and M. I. Polikarpov, “Numerical evidence of chiral magnetic effect in lattice gauge theory,” Phys. Rev. D 80 (2009) 054503, arXiv:0907.0494 [hep-lat].
  * [8] P. V. Buividovich, M. N. Chernodub, E. V. Luschevskaya, and M. I. Polikarpov, “Quark electric dipole moment induced by magnetic field,” Phys. Rev. D 81 (2010) 036007, arXiv:0909.2350 [hep-ph].
  * [9] A. Yamamoto, “Lattice study of the chiral magnetic effect in a chirally imbalanced matter,” Phys. Rev. D 84 (2011) 114504, arXiv:1111.4681 [hep-lat].
  * [10] P. V. Buividovich, “Anomalous transport with overlap fermions,” Nucl. Phys. A 925 (2014) 218–253, arXiv:1312.1843 [hep-lat].
  * [11] G. S. Bali, F. Bruckmann, G. Endrődi, Z. Fodor, S. D. Katz, and A. Schäfer, “Local CP-violation and electric charge separation by magnetic fields from lattice QCD,” JHEP 04 (2014) 129, arXiv:1401.4141 [hep-lat].
  * [12] N. Astrakhantsev, V. V. Braguta, M. D’Elia, A. Y. Kotov, A. A. Nikolaev, and F. Sanfilippo, “Lattice study of the electromagnetic conductivity of the quark-gluon plasma in an external magnetic field,” Phys. Rev. D 102 no. 5, (2020) 054516, arXiv:1910.08516 [hep-lat].
  * [13] M. Puhr and P. V. Buividovich, “Numerical Study of Nonperturbative Corrections to the Chiral Separation Effect in Quenched Finite-Density QCD,” Phys. Rev. Lett. 118 no. 19, (2017) 192003, arXiv:1611.07263 [hep-lat].
  * [14] P. V. Buividovich, D. Smith, and L. von Smekal, “Numerical study of the chiral separation effect in two-color QCD at finite density,” Phys. Rev. D 104 no. 1, (2021) 014511, arXiv:2012.05184 [hep-lat].
  * [15] S. Durr, “Taste-split staggered actions: eigenvalues, chiralities and Symanzik improvement,” Phys. Rev. D 87 no. 11, (2013) 114501, arXiv:1302.0773 [hep-lat].
  * [16] S. Borsányi, G. Endrődi, Z. Fodor, A. Jakovác, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabó, “The QCD equation of state with dynamical quarks,” JHEP 11 (2010) 077, arXiv:1007.2580 [hep-lat].
  * [17] G. S. Bali, F. Bruckmann, G. Endrődi, Z. Fodor, S. D. Katz, S. Krieg, A. Schäfer, and K. K. Szabó, “The QCD phase diagram for external magnetic fields,” JHEP 02 (2012) 044, arXiv:1111.4956 [hep-lat].
  * [18] G. S. Bali, F. Bruckmann, G. Endrődi, Z. Fodor, S. D. Katz, and A. Schäfer, “QCD quark condensate in external magnetic fields,” Phys. Rev. D 86 (2012) 071502, arXiv:1206.4205 [hep-lat].
  * [19] A. Avdoshkin, A. V. Sadofyev, and V. I. Zakharov, “IR properties of chiral effects in pionic matter,” Phys. Rev. D 97 no. 8, (2018) 085020, arXiv:1712.01256 [hep-ph].