# The decidability of the Riemann hypothesis

Kevin Broughan
Department of Mathematics and Statistics
University of Waikato
Hamilton, New Zealand
<EMAIL_ADDRESS>

###### Abstract

Using a result of recursive function theory and results of the complex
analysis of Takeuti which is based on a type theory and the work of Kreisel,
and which gives a conservative extension of first order Peano arithmetic (PA),
we show that RH is decidable in PA

Keywords: total partial recursive function, elementary complex analysis,
decidabilty, Peano arithmetic PA, Riemann Hypothesis, conservative extension

MSC2020 Subject Classification: 11M26

##  1 Introduction

It is apparent that many mathematical propositions are undecidable. See for
example [3, Chapter 10]. After over 160 years of strenuous endeavours by some
of the best mathematical minds of many generations but without a resolution of
the Riemann hypothesis (RH), one might suspect that it is undecidable also.
However, most who have thought about the issue believe it can be resolved, and
there has been recent significant progress in that direction. See for example
the work of Nicolas, Rogers and Tao, Dobner, Polymath15 and Ono et al,
recorded in [3].

In [3, Chapter 11] there are some proof ideas on showing RH is decidable.
However, the assumption that all critical line zeros are simple was made. It
is the purpose of this note to remove that assumption.

To do this two powerful theories are used: recursive function theory and
Takeuti’s extension of first order Peano arithmetic which is based on the work
of Kreisel and on a form of type theory. This extension not only provides a
great deal of real mathematical analysis and complex analysis, but also has
the property of being “conservative”. That means anything that can be proved
using its theory can also be proved using Peano arithmetic.

The result from recursive function theory is much simpler to state and
understand. It is in stark contrast to the theorem of Rice, which is often
used to demonstrate the ubiquity of undecidable propositions.

In Section 2 we give the main definitions, in Section 3 the supporting lemmas,
including the statement of the result from recursive function theory, give
some of the results we using from Takeuti’s system for real and complex
analysis, and derive in PAT an increasing enumeration of the distinct critical
zeros of $\zeta(s)$ which have positive imaginary part. Finally, in Section 4
we give the proof of the main result.

It goes without saying that a proof of the decidability of RH implies all of
the equivalent forms of the hypothesis are decidable, including those reported
in [2, 3]. In addition, the simplicity of the proof indicates that extensions
to families of L-functions such as the extended Selberg class should not be
too difficult.

This work is based principally on the methods developed by Gaisi Takeuti,
which follows that of Gehard Gentzen and George Kreisel. There is a summary of
some of this work, especially that of Gentzen and Takeuti, in [3]. See
Sections M.11.6, O.3 and O.4. To clarify confusion, Takeuti does not base his
development of real and complex analysis on the form of arithmetic shown by
Gentzen to be consistent, but simply first order Peano arithmetic PA.

##  2 Definitions

In this section we set out the definitions which are used in the note. These
are from recursive function theory, Takeuti’s elementary complex analysis from
mathematical logic, and from analytic number theory. There are no derivations
and the reader might need to consult the references.

###  2.1 Recursive functions

This is a study of a class of functions
$f:\mathbb{N}_{0}\rightarrow\mathbb{N}_{0}$ where
$\mathbb{N}_{0}=\\{0,1,2,3,\ldots\\}$ defined procedurally, and where the
domain is not necessarily the whole of $\mathbb{N}_{0}$ but is determined by
the structure of $f$. See [3, Sections N.2-N.4] and the references given
there, especially [4].

The class of partial recursive functions, denoted PRF, is the least class of
functions $f\colon\mathbb{N}_{0}^{k}\rightarrow\mathbb{N}_{0}$ for some $k\geq
1$, where the domain $A$ of $f$ need not be all of $\mathbb{N}_{0}^{k}$, but
is determined by the structure of the function’s definition, such that the
class is closed under the operations (a), (b) and (c), where

(a) PRF contains the zero function $Z(x)=0$ for all $x\in\mathbb{N}_{0}$, the
successor function $S(x)=x^{\prime}=x+1$ for all $x\in\mathbb{N}_{0}$, and for
each $j$ with $1\leq j\leq k$ the projection function
$\pi_{j}^{k}(x_{1},\ldots,x_{k})=x_{j}$. These three functions are everywhere
defined so are total.

(b) If $f(\mathbf{x})$ and $g_{1}(\mathbf{y}),\ldots,g_{k}(\mathbf{y})$ are in
PRF, so is the composite function
$h(\mathbf{y})=f(g_{1}(\mathbf{y}),\ldots,g_{k}(\mathbf{y}))$. If any of the
$g_{j}$ are undefined at a particular $\mathbf{y}$ then so is $f(\mathbf{y})$.

(c) PRF is closed under primitive recursion, namely if $f(\mathbf{x})$ and
$g(\mathbf{x},y,z)$ are functions in PRF then so is $h(\mathbf{x},y)$ which is
defined by

$\displaystyle h(\mathbf{x},0)$ $\displaystyle=f(\mathbf{x}),$ $\displaystyle
h(\mathbf{x},y+1)$
$\displaystyle=\begin{cases}&g(\mathbf{x},y,h(\mathbf{x},y)),\\\
&\uparrow\text{ if $h(\mathbf{x},y)$ is undefined}.\\\ \end{cases}$

The symbol $\uparrow$ means undefined and $\downarrow$ defined.

(d) The set of functions PRF is closed under the operation minimalization,
namely, if $g(\mathbf{x},y)$ is in PRF then so is
$h(\mathbf{x})=\mu_{y}(g(\mathbf{x},y)=0)$ where

$\displaystyle h(\mathbf{x})\text{ is the least $y$ such that
$g(\mathbf{x},y)=0$},$ and, for $0\leq j<y,~{}g(\mathbf{x},j)$, is defined and
non-zero $\displaystyle h(\mathbf{x})\text{ is undefined if no $y$ exists with
$g(\mathbf{x},y)=0$},$ $\displaystyle~{}~{}~{}~{}~{}\text{ or, if such a $y$
exists, there is a $j$ with $1\leq j<y$ with $g(\mathbf{x},j)$ undefined}.$

The functions in PRF which are everywhere defined are called total functions
and the corresponding subset denoted TPR.

We use the Church–Turing thesis to simplify arguments. A statement of this is
that a partial function $f:\mathbb{N}_{0}^{k}\rightarrow\mathbb{N}_{0}$ is in
PRF if and only if there is an algorithm which computes $f$. This is in the
sense that the domain of $f$ is all points where $f$ returns a value and the
complement of the domain is where the algorithm which computes the function
does not halt.

Let $A\subset\mathbb{N}_{0}^{k}$ be a subset. We say $R$ is decidable if its
characteristic function $C_{A}$ is a total partial recursive function TPR.

A subset $A\subset\mathbb{N}_{0}$ (or $\mathbb{N}_{0}^{k}$) is recursively
enumerable (denoted r.e.) if it is the domain of a partial recursive function
PRF. This is somewhat unintuitive, but the meaning is clearer in the case
$k=1$ where for a subset $A\subset\mathbb{N}_{0}$, the following are
equivalent: (a) $A$ is recursively enumerable, and (b) $A$ is the empty set or
the range of a total partial recursive function TPR on $\mathbb{N}_{0}$ [3,
Theorem N.4].

###  2.2 Mathematical logic

Very little is needed from this large field of ideas and techniques. For a
concise summary of the basic concepts see for example [3, Appendix M],
especially Section M.11.1 for first order Peano arithmetic PA, and Section
M.11.6 for Takeuti’s conservative extension of PA for real and complex
analysis which we call PAT (but Takeuti calls “elementary complex analysis”).
Consulting Takeuti’s well written [6] is worthwhile, especially Chapters 3 and
4. Motivation might be enhanced by commencing reading at the final section
Section 4.8, where we see after a description of the sort of results which
cannot be proved in PAT such as those relying on the axiom of choice,

> …classical analytic number theory, for example, the theory in [2] (which is
> Ingham’s 1971 text, The distribution of prime numbers) can be carried out in
> elementary complex analysis.

An extension $T_{2}$ of a theory $T_{1}$ is called conservative if any
proposition which can be proved in the theory $T_{2}$ can be proved in
$T_{1}$.

###  2.3 Analytic number theory

Very little is needed from Analytic number theory, only the two real valued
Riemann-Siegel functions $Z$ and $\vartheta$ which are defined as follows. See
[5, Section 6.5] or [7, Section 4.17].

Take $s=\textstyle{\frac{1}{2}}+it$ and $x=y=\sqrt{t/(2\pi)}$ in the
approximate functional equation for $\zeta(s)$, (see for example [3, Theorem
G.2]), to get

$\zeta(\textstyle{\frac{1}{2}}+it)=\sum_{n\leq
x}\frac{1}{n^{\textstyle{\frac{1}{2}}+it}}+\chi(\sigma+it)\sum_{n\leq
x}\frac{1}{n^{-\textstyle{\frac{1}{2}}+it}}+O\left(\frac{1}{t^{1/4}}\right),$

where $|\chi(\textstyle{\frac{1}{2}}+it)|=1$. Define

$\vartheta(t)=-\frac{\arg(\chi(\textstyle{\frac{1}{2}}+it))}{2}\text{ and
}Z(t)=e^{i\vartheta(t)}\zeta(\textstyle{\frac{1}{2}}+it).$

We get

$Z(t)=2\sum_{n\leq x}\frac{\cos(\vartheta(t)-t\log
n)}{\sqrt{n}}+O\left(\frac{1}{t^{1/4}}\right),$

and can write

$\vartheta(t)=\frac{1}{2}\log\left(\frac{t}{2\pi}\right)-\frac{t}{2}-\frac{\pi}{8}+O\left(\frac{1}{t}\right).$

##  3 Lemmas

The first lemma is from Takeuti’s real analysis FA [6, Chapter 3]:.

######  Lemma 1.

(Takeuti)[6, Proposition 3.5.1] The supremum of a real sequence which is
bounded above exists as a real number in PAT. The corresponding property of
the infimum can be stated in a similar manner.

######  Lemma 2.

[3, Theorem N.10][4, Theorem 7-2.14] A subset of the natural numbers
$\mathbb{N}_{0}$ is total partial recursive, and hence decidable, if and only
if it is the range of a strictly increasing total partial recursive function,
that is to say if it can be enumerated TPR in a strictly increasing order.

This next lemma is a variation on [3, Section 11.10 (3)]. All of the arguments
can be carried out in either Takeuti’s conservative extensions of PA called FA
for mathematical analysis or his “elementary complex analysis” for complex
analysis, i.e. in what we call PAT.

######  Lemma 3.

(Enumeration of the critical zeta zeros) The critical line zeros of $\zeta(s)$
with positive imaginary part are recursively enumerable TPR in increasing
order.

###### Proof.

(1) By Lemma 1 every real sequence which is bounded above has a real least
upper bound. This implies the Archimedian axiom for Takeuti’s system FA holds,
and thus for each pair of reals $x,y$ with $x<y$ there is a rational number
$r$ say with $x<r<y$. By [6, Theorem 4.5.2], the zeros of an entire function
are isolated. Thus the zeros of $\zeta(s)$ on the critical line are isolated
in PAT, and therefore, since

$Z(t)=exp(i\vartheta(t))\zeta(\textstyle{\frac{1}{2}}+it),~{}t\in\mathbb{R}$

so are those of $Z(t)$ in $\mathbb{R}$.

If $\gamma$ is such a zero and $x_{\gamma}<\gamma<y_{\gamma}$ is an interval
in $\mathbb{R}$ which contains the zeta zero with imaginary part $\gamma$ and
no other zero, let $r_{\gamma}$ be a rational with
$x_{\gamma}<r_{\gamma}<y_{\gamma}$. This shows that the zero’s of $Z(t)$ are
an enumerable set.

By the results of Titchmarsh [8], based on the Approximate Functional equation
[3, Theorem G.2], the positive zeros of $Z(t)$ have an infinite subset $h_{n}$
with $h_{n}\rightarrow\infty$.

Now let $\gamma_{1}$ be the real infimum of the set of zeros of $Z(t),~{}t>0$.
Since the zeros are isolated $Z(\gamma_{1})=0$. If $\gamma_{n}$ has been
defined, let $\delta>0$ be such that $Z(t)^{2}>0$ on
$(\gamma_{n},\gamma_{n}+\delta]$ and let $(a_{n})$ be an enumeration of the
zeros of $Z(t),~{}t>0$ which satisfy $a_{n}>\delta$. Let
$\gamma_{n+1}=\inf_{n}a_{n}$. Then $Z(\gamma_{n+1})=0$ and $Z(t)^{2}>0$ on
$(\gamma_{n},\gamma_{n+1})$. All this is true in FA and thus in PAT.

(2) Let $S:=\\{\gamma_{n}:n\in\mathbb{N}\\}$ and $T=\\{\gamma:Z(\gamma)=0\\}$
so $S\subset T$. If $S\neq T$ let $\gamma\in T\setminus S$. Note that for all
$n$, $\gamma_{n}<\gamma<\gamma_{n+1}$ is false, since by Step (1) we have
$Z(t)^{2}>0$ on $(\gamma_{n},\gamma_{n+1})$. Hence $\gamma_{n}<\gamma$ for all
$n$. By Lemma 1, the supremum $\gamma_{0}:=\sup\gamma_{n}$ exists as a real in
FA. By continuity we have $Z(\gamma_{0})=0$. But zeros are isolated so we must
have $\gamma_{0}=\gamma_{n}$ for some $n$, which as we have seen is false.
Hence $S=T$ and we can express the set of positive zeros of $Z(t)$ as

$\\{\gamma_{1}<\gamma_{2}<\gamma_{3}<\cdots\\}\text{ with
}\gamma_{n}\rightarrow\infty.$

This completes the proof.  $\sqcap$$\sqcup$

######  Lemma 4.

(Takeuti)[6, Theorem 4.5.6] Let $f$ be holomorphic in a circular disk $\Delta$
and be not identically zero in $\Delta$. Let $z_{j}$ be the zeros of $f(z)$,
each zero being counted as many times as its multiplicity. For every closed
curve $\gamma$ in $\Delta$ which does not pass through a zero we have

$\sum_{j}n(\gamma,z_{j})=\frac{1}{2\pi
i}\oint_{\gamma}\frac{f^{\prime}(z)}{f(z)}~{}dz,$

where the sum has only a finite number of non-zero terms.

##  4 Result

Again, We use Takeuti’s conservative extension of Peano Arithmetic PAT, (see
Appendix M Section M.11.6), and claim, RH is decidable in PAT and therefore in
PA.

######  Theorem 5.

Assuming first order Peano arithmetic PA, the Riemann hypothesis is decidable.

###### Proof.

(1) Let $(\gamma_{n}:n\in\mathbb{N}_{0})$ be the increasing enumeration of the
zeros given by Lemma 3. For $n,j\in\mathbb{N}_{0}$, let $\mathbb{C}_{n,j}$ be
a circular contour center $\gamma_{n}$ and radius $1/j$. Then by Lemma 4 the
integral

$\frac{1}{2\pi i}\oint_{C_{n,j}}\frac{\xi^{\prime}(z)}{\xi(z)}~{}dz$

exists and is monotonically decreasing in $j$ with values in $\mathbb{N}$. By
Lemma 1, the limit

$\lim_{j\rightarrow\infty}\frac{1}{2\pi
i}\oint_{C_{n,j}}\frac{\xi^{\prime}(z)}{\xi(z)}~{}dz$

exists in PAT. It is the multiplicity of the zero $\gamma_{n}$. Let its value
be denoted by the positive integer $l_{n}$.

(2) Let $f(0)=\gamma_{0}=1$, and using Lemma 3, for all $n\geq 1$ let
$f(n)=\gamma_{n}$ be the strictly increasing enumeration of the imaginary
parts of the zeros of $\zeta(s)$ on the critical line. Let $J_{0}:=[i,1+i]$
and for each $n\geq 1$ set

$\delta_{n}:=\frac{f(n)+f(n-1)}{2}\text{ and
}J_{n}:=[i\delta_{n},1+i\delta_{n}].$

Then there is no zero of $\zeta(s)$ in the intervals $J_{0}$ or $J_{1}$. For
each zero of $\zeta(s)$ in the interval $J_{n}$, indent the contour with a
semi-circular path downwards, so the new zeros are associated with the $n$th
region, which is the interior of the indented $\Delta_{n}$ defined next.

(3) Now set $g(0)=0$, and for $n\in\mathbb{N}$ define the rectangular contour

$\Delta_{n}:=[i\delta_{n},1+i\delta_{n}]\cup[1+i\delta_{n},1+i\delta_{n+1}]\cup[1+i\delta_{n+1},i\delta_{n+1}]\cup[i\delta_{n+1},i\delta_{n}],$

indented as may be needed, and set

$M(n):=\frac{1}{2\pi i}\oint_{\Delta_{n}}\frac{\xi^{\prime}(z)}{\xi(z)}~{}dz.$

Recall $l_{n}$ is defined in Step (1). Let

$g(n):=\begin{cases}g(n-1)&\text{ if }M(n)=l_{n},\\\ n&\text{ otherwize}.\\\
\end{cases}$

Then $g$ is an non-decreasing function which, by PAT and the Church-Turing
thesis, is TPR. By Lemma 2, its range is a decidable set. But the range of $g$
indexes all of the (possibly indented) rectangles containing a zero of
$\zeta(s)$ off the critical line. Thus, RH is decidable.

$\sqcap$$\sqcup$

## Acknowledgements

The author acknowledges helpful communications with Tim Stokes, Yuri
Matiasevich and Alex Weiermann.

## References

  * [1]
  * [2] K. A. Broughan, Equivalents of the Riemann hypothesis, volumes one and two, Cambridge University Press, 2017.
  * [3] K. A. Broughan, Equivalents of the Riemann hypothesis, volume three, Cambridge University Press, 2023.
  * [4] N. Cutland, Computability: an introdution to recursive function theory
  * [5] H. M. Edwards, Riemann’s Zeta Function, Academic Press, 1974\. Reprinted by Dover, 2001.
  * [6] G. Takeuti, Two applications of logic to mathematics, Iwanami Shoten and Princeton University Press, 1978.
  * [7] E.C. Titchmarsh and D.R. Heath-Brown, The Theory of the Riemann Zeta-Function, 2nd edn, Oxford University Press, 1986.
  * [8] E. C. Titchmarsh, On van der Corput’s method and the zeta-function of Riemann (IV), Quart. J. Math. Oxford ser. 5 (1934), 98–105.


aainstitutetext: Max-Planck-Institut für Physik (Werner-Heisenberg-Institut),
Förhringer Ring 6, 80805 München, Germany. bbinstitutetext: State Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese
Academy of Science, Beijing 100190, People’s Republic of China

# Magnetic field induced lattice ground states from holography

Yan-Yan Bu a Johanna Erdmenger a Jonathan P. Shock a Migael Strydom
<EMAIL_ADDRESS>
###### Abstract

We study the holographic field theory dual of a probe $SU(2)$ Yang-Mills field
in a background $(4+1)$-dimensional asymptotically Anti-de Sitter space. We
find a new ground state when a magnetic component of the gauge field is larger
than a critical value. The ground state forms a triangular Abrikosov lattice
in the spatial directions perpendicular to the magnetic field. The lattice is
composed of superconducting vortices induced by the condensation of a charged
vector operator. We perform this calculation both at finite temperature and at
zero temperature with a hard wall cutoff dual to a confining gauge theory. The
study of this state may be of relevance to both holographic condensed matter
models as well as to heavy ion physics. The results shown here provide support
for the proposal that such a ground state may be found in the QCD vacuum when
a large magnetic field is present.

###### Keywords:

Gauge-gravity correspondence, Holography and condensed matter physics(AdS/CMT)

††preprint: MPP-2012-144

##  1 Introduction

The study of black hole instabilities is an important research topic that has
led to very interesting results. In particular, within gauge/gravity duality,
the study of Anti-de Sitter black hole solutions and their stability
properties is important for understanding thermal states on the gauge theory
side. Last year in Ammon:2011je , some of the authors of the present paper
studied an $SU(2)$ Einstein-Yang-Mills model at finite temperature in
asymptotically AdS space. They found that when a magnetic component of the
gauge field reaches a critical value in units of the temperature, the system
becomes unstable. Though the critical value of the magnetic field for onset of
the instability was calculated, the new ground state of the system was not
known. In the current work we calculate a ground state solution, using a
perturbative analysis similar to the one performed by Abrikosov in
Abrikosov:1956sx  for type II superconductors. In agreement with the work of
Abrikosov we find that the ground state is a triangular lattice. There have
been many attempts recently to model lattices holographically with the goal of
providing more realistic models for condensed matter systems Horowitz:2012ky ;
Flauger:2010tv , and this novel procedure for generating a lattice dynamically
adds to these developments. Moreover, our holographic model provides support
for recent QCD studies of $\rho$ meson condensation from a strong magnetic
field Chernodub:2010qx ; Chernodub:2010zw ; Chernodub:2011gs . The effect
described here is similar to the Nielsen-Olsen solution for gluon condensation
Nielsen:1978rm  and to magnetically catalysed W boson condensation
Ambjorn:1988tm ; Ambjorn:1989bd ; Ambjorn:1989sz .

Although this is the first holographic calculation to explicitly uncover an
Abrikosov lattice in 3+1 dimensions, it is not the first to examine spatially
inhomogeneous phases of strongly coupled field theories. In Horowitz:2012ky ,
the authors studied the holographic construction of an Einstein-Maxwell-scalar
theory at finite temperature and density. They looked at the gauge theory
optical conductivity, which is the conductivity in the direction of an applied
electric field. They broke the translational invariance explicitly by imposing
scalar field boundary conditions in the form of a lattice modulated in one of
the Minkowski spatial directions. A fully backreacted solution was found which
thus induces a spatially inhomogeneous black hole solution. This leads to an
extremely rich behaviour of the frequency dependent optical conductivity. At
low frequencies there appears a Drude peak. A Drude peak is a broadening of
the zero frequency delta peak in the conductivity. In real materials this is
due to impurities and finite temperature effects. The Drude peak is not
present when translational invariance is unbroken. The solution also exhibits
a power law behaviour at frequencies intermediate with respect to the
temperature, and a constant value in the high frequency regime. The power law
behaviour is the same as that found experimentally in cuprates, while the
constant value at high frequencies is expected from conformal invariance. The
setup in this context was a (2+1)-dimensional model where the lattice was
periodic in only one of the spatial dimensions. A more realistic lattice
structure would be highly desirable.

While the lattice in the approach of Horowitz:2012ky  was implemented in the
boundary conditions, there are a number of other mechanisms known that lead
dynamically to ground states without translational symmetry. One approach was
pioneered in Domokos:2007kt  by studying a Yang-Mills-Chern-Simons theory in
the gauge/gravity context. It was shown that the Chern-Simons term can induce
an instability which leads to a ground state with both translational and
rotational symmetry breaking. Such work was continued in Nakamura:2009tf ;
Chuang:2010ku ; Bergman:2011rf ; Bayona:2011ab ; Takeuchi:2011uk  where the
spatially homogeneous phase was found to be unstable in a variety of
gravitational contexts in the presence of Chern-Simons couplings. The
perturbative analysis of quasinormal modes that become tachyonic at finite
momentum gives a relatively simple computational tool for finding
instabilities to ground states without translational symmetry. These solutions
are found to induce a helical current Ooguri:2010kt ; Donos:2012gg ;
Donos:2012wi . Interestingly, a Chern-Simons term is not always enough to
induce such an instability. It was shown by Ammon:2011hz  that this type of
instability does not exist in the D3/D7 system.

Translational invariance can also be broken with a magnetic field or with
magnetic monopoles. The former was first studied by Gauntlett et al. in
Donos:2011qt  and Donos:2011bh  where the instability of the magnetically
charged black hole in a top-down framework was studied in detail. In the
latter work, an infinite family of solutions coming from $D=11$ supergravity
was shown to exhibit a magnetically catalysed instability. Such work is
important as it proves that these instabilities can also come from real string
theory constructions. The subject of magnetic monopoles in $(3+1)$-dimensional
AdS space was studied in Bolognesi:2010nb  and Sutcliffe:2011sr . These
magnetic monopoles are solutions to the scalar field in a Yang-Mills-Higgs
theory with gauge group $SU(2)$. In a certain limit where the monopole
magnetic charge becomes large and a “monopole wall” is formed, it was shown in
Bolognesi:2010nb  that there is a W boson instability. In Sutcliffe:2011sr  a
hexagonal lattice ground state of these monopole walls was found numerically.
In Allahbakhshi:2011nh  the holographic dual of a self-gravitating Julia-Zee
Dyon was constructed, and it was shown to contain a vortex condensate.

There are some holographic models exhibiting a superconducting phase
transition that results in a vortex lattice ground state. The first we mention
involves an $s$-wave superconductor. In Maeda:2009vf  a type II superconductor
was modelled using a $(3+1)$-dimensional gravitational setup. A type II
superconductor is one for which the external applied magnetic field has two
critical values. When the magnitude of the magnetic field increases beyond the
lower of the two critical values, the field starts to penetrate the
superconducting condensate. Some of the condensate remains until the magnitude
of the field is increased beyond the upper critical value, at which point
superconductivity is completely destroyed. Just before the upper critical
value is reached from below, the ground state of the system is a triangular
Abrikosov lattice Abrikosov:1956sx . The authors of Maeda:2009vf  constructed
a holographic superconductor modelling the behaviour of a type II
superconductor near the upper critical value of the magnetic field and found
the Abrikosov lattice ground state explicitly.111 The model of Maeda:2009vf
does not display the transition at the lower critical magnetic field value
because the gauge field is not dynamical. See Domenech:2010nf  for adding
dynamical gauge fields to holographic superconductors. The transition at the
upper critical value is present however because there the condensate is small
so the backreaction is negligible.  In Murray:2011gr  it was shown how to
construct a similar vortex lattice solution in a model describing a $p$-wave
superconductor. There the authors used a holographic model with an $SU(2)$
gauge field similar to the one described in the current paper. Both of theses
examples are different from our model, however, because here we find a
superconducting Abrikosov lattice ground state that is induced by an $SU(2)$
magnetic field, rather than being destroyed by it. Moreover, in contrast to
these models, we do not need a finite density. Our model is a cousin of
holographic p-wave superconductors where the condensation is induced by a
finite isospin density, holographically realised by a non-trivial temporal
component of the $SU(2)$ gauge field (see Gubser:2008wv  and Ammon:2008fc ;
Ammon:2009fe  as well as the recent Chunlen:2012zy ). Here, in contrast, a
spatial component of this gauge field has a non-trivial profile. Whereas in
Ammon:2008fc ; Ammon:2009fe , a Meissner effect is shown to occur by which a
magnetic field reduces the transition temperature, here it is again the
magnetic field which induces condensation at zero density.

In addition to being interesting in the broader context of holographic
lattices, the model we discuss serves as supporting evidence for a phenomenon
first described by Chernodub et al. in Chernodub:2010qx ; Chernodub:2010zw .
There it was proposed that the QCD$\times$QED vacuum may itself be susceptible
to a superconducting transition when a magnetic field of the order of the QCD
scale is present. Such extreme conditions are rare but they may be present for
a few femtoseconds during highly off-centre heavy ion collisions. The
discovery of this phase came about through the study of an effective field
theory description (the DSGS model proposed by Djukanovic, Schindler, Gegelia
and Scherer in Djukanovic:2005ag ) of $\rho$ mesons interacting with a
magnetic field. A destabilisation of the vacuum was shown that would clearly
lead to the condensation of charged and neutral $\rho$ mesons. This breaks the
$U(1)$ gauge symmetry and leads to a superconductor with the quark-antiquark
pairs in the mesons acting as Cooper pairs. The instability was also found
using an extended Nambu–Jona-Lasinio model with $SU(3)$ colour and $SU(2)$
flavour in Chernodub:2011mc . Lattice gauge theory studies were then performed
looking at QCD in strong magnetic fields and these indicate the same
instability. Moreover, using the DSGS model and guided by the Ginzburg-Landau
model of type II superconductors, a solution was found in which the $\rho$
meson condensate forms an Abrikosov lattice made up of superconducting
vortices Chernodub:2011gs . This may be relevant experimentally. Evidence has
mounted at both RHIC and the ALICE experiment at CERN that strong magnetic
fields may contribute to the physics of the strongly coupled quark gluon
plasma as charges are quickly accelerated during the interaction period
Skokov:2009qp ; Bzdak:2011yy . The importance of these effects remains a
contested topic because the time scales involved are small. However, given
that strong magnetic fields may be present, it is interesting to ask if traces
of this $\rho$ meson condensate could be detected.

In the current paper we find a possible ground state of the system in
Ammon:2011je . As mentioned above, this system was shown to be unstable under
the imposition of a large $SU(2)$ magnetic field.222It was shown in
Callebaut:2011ab  that the same sort of instability occurs in the Sakai-
Sugimoto model, but there the ground state has also not been found. We show
that it has very similar properties to the ground state of a type II
superconductor near the upper critical magnetic field as well as to the ground
state in the model of Chernodub et al. In other words, the ground state is a
triangular Abrikosov lattice. We take here a very simple model of a strongly
coupled finite temperature quantum field theory in $(3+1)$-dimensions with a
global $SU(2)$ symmetry. The dual gravity theory is an $SU(2)$ Einstein-Yang-
Mills theory in $(4+1)$-dimensions with a magnetic component of the $SU(2)$
switched on. We work entirely within the probe approximation, which means that
the Yang-Mills term is small compared to the Einstein-Hilbert term in the
action. We also fix the gauge in such a way that the gauge theory condensate
is transformed under a $U(1)$ subgroup of the global $SU(2)$ symmetry. There
appears to be a certain universality to the triangular lattice ground state.
Here we show that it forms in both the AdS Schwarzschild background (dual to a
finite temperature field theory) as well as the hard wall cutoff model (dual
to a confining field theory). It would be interesting to uncover exactly how
universal these results are.

The two holographic models that we study have several important differences
from QCD. In the finite temperature model there is no confinement or chiral
symmetry breaking and so there are no goldstone bosons (pions) present which
are the normal decay modes of the $\rho$ meson in QCD. The hard wall model has
its conformal symmetry broken only by an IR boundary condition which sets a
confinement scale. However, the phenomenology of these two models appears to
be close enough to that of QCD to compare qualitatively with the models of
Chernodub et al.

In section 2 we provide the details of the holographic setup. There we also
explain the strategy behind the perturbative expansion of the $SU(2)$ gauge
field near the critical magnetic field. Since we follow the philosophy of
Abrikosov’s calculation of the ground state in type II superconductors, which
was done in the Ginzburg-Landau model, in section 3 we give a brief outline of
this approach and then follow it to solve perturbatively up to third order. In
section 4 we discuss the numerical results and analyse the free energy of the
different lattice solutions, showing that the triangular lattice has the
lowest free energy of all the Abrikosov solutions studied. It is important to
note that we are not able to show conclusively that we have found _the_ ground
state but we are able to find a state with lower free energy than the
translationally invariant state and that has lowest energy within a large
class of lattice solutions. In section 5 we conclude and give an outline of
important future work.

##  2 Holographic setup

###  2.1 The finite temperature and hard wall backgrounds

The system we study is an Einstein-Yang-Mills theory on the (Poincaré patch
of) an asymptotically AdS5 geometry with an $SU(2)$ gauge field. The action is

$S=\int d^{5}x\sqrt{-g}~{}\left\\{\frac{1}{16\pi
G_{N}}\left(R+\frac{12}{L^{2}}\right)-\frac{1}{4\hat{g}^{2}}\mathrm{tr}\left(F_{\mu\nu}F^{\mu\nu}\right)\right\\}~{},$
(1)

where $\hat{g}$ is the Yang-Mills coupling, $G_{N}$ is the 5D gravitational
constant and $L$ is the AdS5 radius. $R$ and $F$ are the Ricci scalar and
Yang-Mills field strength respectively.

We consider the probe approximation, where the Yang-Mills term is small
compared to the Einstein-Hilbert term, so that the backreaction of the gauge
fields on the geometry can be neglected. We thus choose a fixed 5-dimensional
background metric, given by

$ds^{2}=\frac{L^{2}}{u^{2}}\left(-f(u)dt^{2}+dx^{2}+dy^{2}+dz^{2}+\frac{du^{2}}{f(u)}\right)~{},$
(2)

where the asymptotically AdS region is at $u\rightarrow 0$. We study two
different models. The first is a finite temperature model where the background
is AdS Schwarzschild, first proposed in Witten:1998qj . In this case,
$f(u)=1-\frac{u^{4}}{u_{H}^{4}}$, where $u_{H}$ is the location of the planar
black hole horizon. The Hawking temperature of the black hole is $T=1/\pi
u_{H}$. The second model is the hard wall cutoff model, proposed in
Erlich:2005qh ; DaRold:2005zs , where $f(u)=1$ and the geometry terminates at
a radial distance $u_{C}$. This model corresponds to a zero temperature theory
($u_{H}=\infty$), but it still has a scale $u_{C}$ which corresponds to a
confinement scale in the gauge theory. The intrinsic scales in these theories
allow us to form a dimensionless magnitude for the magnetic field. This will
be the parameter that we tune in order to find the instability of the
spatially invariant ground state. Without loss of generality we can choose
units where $u_{H}=1$ in the finite temperature theory and $u_{C}=1$ in the
confining theory. Factors of $u_{H}$ and $u_{C}$ can then be restored through
dimensional analysis. In the following the exact form of the metric is not
important until we come to solving the numerical equations in the radial
direction of AdS.

###  2.2 The Yang-Mills action

The relevant part of the action simplifies to

$S=-\frac{1}{4\hat{g}^{2}}\int
d^{5}x\sqrt{-g}~{}\mathrm{tr}\left(F_{\mu\nu}F^{\mu\nu}\right)~{},$ (3)

with the equations of motion

$\nabla^{\mu}F^{a}_{\mu\nu}+\epsilon^{abc}\mathcal{A}^{b\mu}F^{c}_{\mu\nu}=0~{}.$
(4)

The $SU(2)$ gauge field is
$\mathcal{A}=\mathcal{A}_{\mu}^{a}\tau^{a}dx^{\mu}$, for $a=1\dots 3$. We use
the convention where the Lie algebra basis is given by
$\tau^{a}=\frac{\sigma^{a}}{2i}$, with $\sigma^{a}$ the Pauli matrices, and
the structure constants $f^{abc}$ are defined by
$[\tau^{a},\tau^{b}]=\epsilon^{abc}\tau^{c}$ so that $f^{abc}=\epsilon^{abc}$.
With these definitions, the components of the field-strength tensor
$F=d\mathcal{A}+\mathcal{A}\wedge\mathcal{A}$ become

$F^{a}_{\mu\nu}=\partial_{\mu}\mathcal{A}^{a}_{\nu}-\partial_{\nu}\mathcal{A}^{a}_{\mu}+\epsilon^{abc}\mathcal{A}^{b}_{\mu}\mathcal{A}^{c}_{\nu}~{}.$
(5)

It will be important to understand how gauge transformations affect the
system. Under a gauge transformation $e^{i\Lambda(x^{\mu})}$, $\mathcal{A}$
transforms as

$\mathcal{A}_{\mu}\rightarrow\mathcal{A}_{\mu}+\delta\mathcal{A}_{\mu}=e^{i\Lambda}\mathcal{A}_{\mu}e^{-i\Lambda}-i\partial_{\mu}e^{i\Lambda}e^{-i\Lambda}~{}.$
(6)

When $\Lambda(x^{\mu})$ is an infinitesimal transformation, this becomes

$\delta\mathcal{A}^{a}_{\mu}=\mathcal{D}_{\mu}\Lambda^{a}=\partial_{\mu}\Lambda^{a}+\epsilon^{abc}\mathcal{A}^{b}_{\mu}\Lambda^{c}~{}.$
(7)

The gauge transformations give us the freedom to fix the gauge
$\mathcal{A}^{a}_{u}=0$. We work in this gauge from now on.

In this paper we look at the effect of a strong (flavour-)magnetic field given
by $F^{3}_{xy}=B$, with all other components of $F^{a}_{\mu\nu}$ vanishing. As
we will see, when $B$ becomes large333Since we have chosen the units where
$u_{H}=1$ or $u_{C}=1$, $B$ is a dimensionless quantity. Restoring the units,
the statement is that $Bu_{H}^{2}=B/(\pi T)^{2}$ or $Bu_{C}^{2}\sim
B/\Lambda_{QCD}^{2}$ is large, or that $B$ is large compared to the radial
scale of the background., other components of $F$ become non-zero dynamically.
To get a consistent set of equations we therefore consider a gauge field
$\mathcal{A}$ of the form

$\mathcal{A}=\sum_{a=1,2,3,\mu=x,y}\mathcal{A}^{a}_{\mu}(x,y,u)\tau^{a}dx^{\mu}~{}.$
(8)

It turns out that we can turn off the $t$ and $z$ dependence of the gauge
field and still have consistent equations. This simplifies the equations.
Turning off the $t$ dependence guarantees a static solution. Turning off the
$z$ dependence, where the $z$ direction is parallel to the magnetic field,
yields a lattice in the $x,y$-plane.

The action 3 has an $SU(2)$ gauge freedom. Choosing the solution
$F^{3}_{xy}=B$, with all other components vanishing, breaks this symmetry.
Only $U(1)$ transformations of the form $\Lambda=\Lambda^{3}\tau^{3}$ leave it
invariant. For $B$ large enough, all the components in 8 become nonzero due to
the dynamics. We thus claim to have a superconductor, because the $U(1)$
symmetry is broken dynamically. Note however that it is technically a
superfluid because the $U(1)$ gauge symmetry in the bulk theory gets mapped to
a global symmetry in the field theory. Taking the linear combinations
$\mathcal{E}^{\pm}_{\mu}=\mathcal{A}^{1}_{\mu}\pm i\mathcal{A}^{2}_{\mu}$
gives fields that transform in the fundamental of the remaining gauge
symmetry. It can be checked from 7 that $\mathcal{E}^{\pm}_{\mu}\rightarrow\mp
i\Lambda^{3}\mathcal{E}^{\pm}_{\mu}$ whenever $\Lambda=\Lambda^{3}\tau^{3}$.
Later on we work only with the fields $\mathcal{E}^{+}$, which we rename to
$\mathcal{E}$.

###  2.3 Perturbative expansion of the gauge fields

Substituting the ansatz 8 into equation 4 yields nine coupled partial
differential equations in the variables $x$, $y$ and $u$. Of these nine
equations of motion, six are dynamical equations for each field
$\mathcal{A}^{1,2,3}_{x,y}$, and three equations are constraints. The
constraint equations arise from the equations of motion for the components
$\mathcal{A}^{1,2,3}_{u}$, which were chosen to be zero using gauge symmetry.

In solving the PDE’s, we follow the strategy of Abrikosov:1956sx ;
AbrikosovBook , which works as follows. When the magnetic field $B$ is smaller
than some critical value $B_{c}$, the field configuration
$\mathcal{A}^{3}_{y}=xB$, $\mathcal{A}^{3}_{x}=0$ and
$\mathcal{A}^{1,2}_{x,y}=0$ solves the equations of motion. This is the normal
phase of the superconductor. As shown in Ammon:2011je , the system enters a
new phase when the magnetic field is increased beyond some critical value
$B_{c}$. In this phase, the superconducting phase, the ground state has a non-
trivial profile for all fields in the ansatz equation 8. We look for this
configuration at some value of $B$ infinitesimally above $B_{c}$, where the
condensate is still small. This lets us do a perturbative expansion in a small
parameter $\varepsilon\sim\frac{B-B_{c}}{B_{c}}$. For notational convenience
we leave this parameter $\varepsilon$ explicit when studying the expansion.
However, it will be absorbed into the definition of the perturbative
corrections to the fields when we come to minimising the energy. We thus write
an ansatz for the expansion in the form

$\displaystyle\mathcal{A}^{3}_{y}$ $\displaystyle=xB_{c}+\varepsilon
A^{3}_{y}+\varepsilon^{2}a^{3}_{y}+\dots,$ (9)
$\displaystyle\mathcal{A}^{a}_{\mu}$ $\displaystyle=\varepsilon
A^{a}_{\mu}+\varepsilon^{2}a^{a}_{\mu}+\dots~{}~{}~{}~{}~{}\mathrm{for~{}}(a,\mu)\neq(3,y)~{},$
(10)

and solve the equations order by order in $\varepsilon$, as detailed in
section 3.

###  2.4 Gauge field boundary conditions

The holographic dictionary relates field theory operators to gravity theory
fields through the relation

$\displaystyle e^{-W_{\mathrm{CFT}}[\mathcal{A}^{(0)}]}=\langle
e^{\int_{\partial
AdS}\mathcal{A}^{(0)}_{\mu}J^{\mu}}\rangle=e^{-S_{\mathrm{on-shell}}}~{}.$
(11)

The minus sign on the right-hand side is because we are in Euclidean space for
simplicity. Here $\mathcal{A}^{(0)}$ is the value of the gauge field
$\mathcal{A}$ at the AdS boundary. It acts as a source in the boundary field
theory. In our setup, the only source we want in the field theory comes from
the component $\mathcal{A}^{3}_{y}=xB$, producing the magnetic field. For the
other components in 8, there should be no explicit source because we want to
model spontaneous symmetry breaking. The spontaneous symmetry breaking results
in a vev444We also need to take holographic renormalisation into account to
yield a finite on-shell action.,

$\displaystyle\langle J^{\mu}\rangle=\left.\frac{\delta
W_{\mathrm{CFT}}}{\delta\mathcal{A}^{(0)}_{\mu}}\right|_{\mathcal{A}^{(0)}_{\mu}=0}=\left.\frac{\delta
S_{\mathrm{on-
shell}}}{\delta\mathcal{A}^{(0)}_{\mu}}\right|_{\mathcal{A}^{(0)}_{\mu}=0}=\left.-\int
d^{4}x\frac{\partial\mathcal{L}}{\partial\left(\partial_{u}\mathcal{A}_{\mu}\right)}\right|_{u=0}$
(12)

The second equality is a generalisation to the radial coordinate of one of the
steps in deriving the Hamilton-Jacobi equation. It relates the variation of
final value of a generalised coordinate with respect to the on-shell action
and the conjugate momentum at the final time.

It is interesting to note that the on-shell action can be written as

$\displaystyle S_{\mathrm{on-shell}}=$
$\displaystyle-\frac{1}{2\hat{g}^{2}}\int_{\partial
AdS}d^{d}x\sqrt{-\gamma}n_{\mu}A^{a}_{\nu}F^{a\mu\nu}+\frac{1}{4\hat{g}^{2}}\int_{AdS}d^{d+1}x\sqrt{-g}\epsilon^{abc}A^{a}_{\mu}A^{b}_{\nu}F^{c\mu\nu}~{},$
(13)

where we integrated by parts and substituted in the equations of motion. The
second term on the right-hand side, the bulk term, is not present in non-
interacting theories. In our case, however, it is present and nonzero even
after using ansatz 8. This bulk term should seemingly influence the
calculation of the condensate when varying with respect to the boundary value.
It turns out that, due to the formula at the right of equality 12, it makes no
contribution.

Equations 11 and 12 imply that in an expansion of the gauge fields near the
AdS boundary, the leading term is the source and the subleading term is
proportional to the vev. The field $\mathcal{A}_{x,y}^{3}$ has a boundary
expansion given by

$\left.\mathcal{A}_{x,y}^{3}\right|_{u\rightarrow
0}=s_{x,y}^{(3)}+v_{x,y}^{(3)}u^{2}+\dots~{},$ (14)

where $s_{x,y}^{(3)}$ is the value of the source, which in this case is the
externally applied magnetic field potential. $v_{x,y}^{(3)}$ is proportional
to the vev corresponding to the magnetisation. We set the boundary conditions
so that the applied magnetic field is not corrected by the higher order
perturbations in $\varepsilon$, whereas the magnetisation will obtain a non-
zero value.

Similarly, the fields $A_{x,y}^{1,2}$ have a boundary expansion given by

$\left.{\mathcal{A}_{x,y}^{1,2}}\right|_{u\rightarrow
0}=s_{x,y}^{(1,2)}+v_{x,y}^{(1,2)}u^{2}+\dots~{},$ (15)

where $s^{(1,2)}_{x,y}$ corresponds to the source of the operator that will
condense to break the $U(1)$ symmetry. We adjust the boundary conditions in
such a way that this vanishes. This means that the symmetry breaking is
spontaneous. $v_{x,y}^{(1,2)}$ is proportional to the vacuum expectation value
of this operator, which we read off to find the resulting supercurrent in the
superconducting phase.

Boundary conditions are also imposed on the fields in the IR. In the case of
the black hole background, we impose regularity at the horizon and in the case
of the hard wall model we impose Neumann boundary conditions.

###  2.5 The gauge theory ground state energy

In finding the ground state, it is important to be able to calculate the
energy of the field theory solution from the action. We would like to compare
the solutions in the normal phase to those in the superconducting phase. The
energy $\mathcal{F}$ of the gauge theory solution is found by using the
holographic dictionary. In the case of the finite temperature solution, we are
in the canonical ensemble and we calculate the free energy, which is
$\mathcal{F}/T=-\ln\mathcal{Z}=-S_{cl}$ with our conventions. Here
$S_{cl}=-\frac{1}{4\hat{g}^{2}}\int d^{5}x\sqrt{-g}F^{a}_{\mu\nu}F^{a\mu\nu}$
is the classical action. In the hard wall case, we are simply calculating the
energy of the field configuration, which is defined in terms of the classical
action in the same way. Since we are only interested in whether the energy of
a particular superconducting solution is lower than that of the normal phase
solution, we can simply calculate the difference
$\Delta\mathcal{F}=\mathcal{F}_{s}-\mathcal{F}_{n}$ and thus do not need to
implement holographic renormalisation. Here $\mathcal{F}_{s}$ is the energy of
the superconducting phase, while $\mathcal{F}_{n}$ is the normal phase energy
with $\mathcal{A}^{3}_{y}=xB$ and all other components zero. We also need to
take care of the fact that $S_{cl}$ diverges when we perform the integral over
the Minkowski directions. This is easy to fix by considering the energy
density555We divide the free energy by $T$ in the finite temperature model to
get a dimensionless $\Omega$. This means that in both models, our total
dimensionless energy is simply -$S_{cl}$. $\Omega$, which is obtained by
integrating $S_{cl}$ only over the world volume of one lattice cell and
dividing by its volume. Having explained how to calculate the energy of a
field configuration, in the next section we turn to the problem of solving the
equations of motion to find the ground state.

##  3 Solving the equations

###  3.1 The comparison with Ginzburg-Landau theory

Before turning to the equations of motion, it helps to first look at the
Ginzburg-Landau equations for an analogy. In some suitable units defined in
Abrikosov:1956sx ; AbrikosovBook , they are

$\displaystyle\left(-i\nabla-\vec{A}\right)^{2}\psi-\psi+|\psi|^{2}\psi=0~{},$
(16)
$\displaystyle~{}\nabla\times\nabla\times\vec{A}=-i\left(\bar{\psi}\nabla\psi-\psi\nabla\bar{\psi}\right)-|\psi|^{2}\vec{A}~{}.$
(17)

Only the structure of these equations is important, so we have ignored
constant factors. Here $\psi$ is the wave function of Cooper pairs, and
$\vec{A}$ is the electromagnetic vector potential. The nine equations of
motion in our system can be split into two groups that roughly correspond to
the two equations above.

The first of the two groups, hereafter called the condensate equations,
contains the six equations for the fields $\mathcal{A}^{1,2}_{x,y,u}$. The
superconducting condensate of the dual field theory, which is like $\psi$
above, is found by differentiating the on-shell action with respect to the
boundary values of $\mathcal{A}^{1,2}_{x,y}$, as in equation 12. Of the six
equations in this group, the dynamical equations are for
$\mathcal{A}^{1,2}_{x,y}$ and the constraint666Recall that we have set
$\mathcal{A}^{a}_{u}=0$. However, its equations of motion still impose
constraints on the other fields.  equations are for $\mathcal{A}^{1,2}_{u}$.
So this first group is analogous to equation 16. The analogy can be made more
clear. As mentioned above, we can make the field definitions
$\mathcal{E}_{x,y}=\mathcal{A}^{1}_{x,y}+i\mathcal{A}^{2}_{x,y}$. Doing so
allows us to combine the six real equations into three complex equations, two
dynamical and one constraint. The constraint equation relates
$\mathcal{E}_{x}$ and $\mathcal{E}_{y}$ such that there is only one complex
degree of freedom left, which is analogous to the state $\psi$. All this is
hard to see at the non-perturbative level, but it illustrates the strategy we
follow for solving the equations at each order: we use the constraint equation
to reduce the two dynamical equations into one, and then solve it.

The second group of equations, which we call the magnetic field equations, is
for the fields $\mathcal{A}^{3}_{x,y,u}$, corresponding to $\vec{A}$ in
equation 17 above. There are three such equations, one of which is a
constraint. At each order we will be able to use the constraint to separate
the equations into one for $\mathcal{A}^{3}_{x}$ and one for
$\mathcal{A}^{3}_{y}$.

###  3.2 The gauge field perturbative expansion in more detail

Having defined the ansatz for our gauge potential in equation 9 we can learn
more about the perturbative expansion by studying the non-linear structure of
the equations of motion. The equation for $\mathcal{A}^{3}_{u}$ is

$\displaystyle-\mathcal{A}^{2}_{x}\partial_{u}\mathcal{A}^{1}_{x}-\mathcal{A}^{2}_{y}\partial_{u}\mathcal{A}^{1}_{y}+\mathcal{A}^{1}_{x}\partial_{u}\mathcal{A}^{2}_{x}+\mathcal{A}^{1}_{y}\partial_{u}\mathcal{A}^{2}_{y}+\partial_{y}\partial_{u}\mathcal{A}^{3}_{y}+\partial_{x}\partial_{u}\mathcal{A}^{3}_{x}=0~{}.$
(18)

We see that the magnetic field components appear in the linear terms, while
the condensate components appear in quadratic terms. This suggests that a
contribution to the condensate components that is first order in the
perturbative expansion influences a second order contribution in the magnetic
field components. More generally, an odd order contribution to the condensate
components influences an even order contribution to the magnetic field
components.

This structure is common throughout all the equations of motion. It turns out
that terms in the perturbative expansion of the magnetic field components that
have an odd order vanish. The even order terms in the condensate components
can then also be set to zero. We can thus constrain the expansion ansatz of
equation 9 to

$\displaystyle\mathcal{E}_{x,y}$ $\displaystyle=\varepsilon
E_{x,y}+\varepsilon^{3}e_{x,y}+\mathcal{O}(\varepsilon^{5})~{},$
$\displaystyle\mathcal{A}^{3}_{y}$
$\displaystyle=xB_{c}+\varepsilon^{2}a^{3}_{y}+\mathcal{O}(\varepsilon^{4})~{},$
(19) $\displaystyle\mathcal{A}^{3}_{x}$
$\displaystyle=\varepsilon^{2}a^{3}_{x}+\mathcal{O}(\varepsilon^{4})~{}.$

Here the calligraphic letters denote the non-perturbative fields. $E_{x,y}$
and $e_{x,y}$ are first and third order contributions to the condensate
components, respectively, while $a^{3}_{x,y}$ are second order corrections to
$\mathcal{A}^{3}_{x,y}$.

Because of this convenient expansion of the fields, the condensate components
and the magnetic components decouple at each order. That means that at each
order, we only need to work with fields we have already solved at previous
orders. Our strategy is thus to solve for the fields in the following
sequence:

$\textstyle{{\mathcal{E}_{x,y}}=~{}~{}~{}~{}~{}~{}~{}}$$\textstyle{\varepsilon
E_{x,y}}$$\textstyle{+}$$\textstyle{\varepsilon^{3}e_{x,y}}$$\textstyle{+~{}~{}\mathcal{O}(\varepsilon^{5})~{},}$$\textstyle{{\begin{array}[]{c}\mathcal{A}^{3}_{y}\\\
{}\hfil\\\ \mathcal{A}^{3}_{x}\end{array}}{\begin{array}[]{c}=\\\ {}\hfil\\\
=\end{array}}{\begin{array}[]{c}xB_{c}\\\ {}\hfil\\\
{}\hfil\end{array}}}$$\textstyle{\begin{array}[]{c}+\\\ {}\hfil\\\
{}\hfil\end{array}}$$\textstyle{\begin{array}[]{c}\varepsilon^{2}a^{3}_{y}\\\
{}\hfil\\\ \varepsilon^{2}a^{3}_{x}\end{array}}$
$\textstyle{{\begin{array}[]{c}+~{}~{}\mathcal{O}(\varepsilon^{4})~{},\\\
{}\hfil\\\ +~{}~{}\mathcal{O}(\varepsilon^{4})~{}.\end{array}}}$

In the next section we start with the linear order solution, which will shed
more light on the procedure that must be implemented at higher orders.

###  3.3 Solving the equations to linear order

Using the expansion 19 and keeping terms to linear order, we find that there
are six remaining equations given (in complex form) by

$\displaystyle 0$
$\displaystyle=-iB_{c}x\partial_{u}E_{y}-\partial_{y}\partial_{u}E_{y}-\partial_{x}\partial_{u}E_{x}~{},$
(20) $\displaystyle 0$
$\displaystyle=B_{c}^{2}x^{2}E_{x}-iB_{c}E_{y}+\left(\frac{f}{u}-f^{\prime}\right)\partial_{u}E_{x}-f\partial_{u}^{2}E_{x}-2iB_{c}x\partial_{y}E_{x}$
$\displaystyle~{}-\partial_{y}^{2}E_{x}+iB_{c}x\partial_{x}E_{y}+\partial_{x}\partial_{y}E_{y}~{},$
(21) $\displaystyle 0$
$\displaystyle=2iB_{c}E_{x}+\left(\frac{f}{u}-f^{\prime}\right)\partial_{u}E_{y}-f\partial_{u}^{2}E_{y}+iB_{c}x\partial_{x}E_{x}+\partial_{x}\partial_{y}E_{x}-\partial_{x}^{2}E_{y}~{}.$
(22)

Here, as above, $f(u)=1-u^{4}$ for the AdS Schwarzschild model and $f(u)=1$
for the hard wall model.

We can solve these equations by following Abrikosov AbrikosovBook . The
solution is given by

$\displaystyle E_{y}$ $\displaystyle=-iE_{x}~{},$ (23) $\displaystyle E_{x}$
$\displaystyle=\sum_{n=-\infty}^{\infty}C_{n}e^{-inky-\frac{1}{2}B_{c}\left(x-\frac{nk}{B_{c}}\right)^{2}}U(u)~{}.$
(24)

$U(u)$ is determined by solving

$U^{\prime\prime}+\left(\frac{f^{\prime}(u)}{f(u)}-\frac{1}{u}\right)U^{\prime}+\frac{B_{c}}{f(u)}U=0~{},$
(25)

subject to the constraints $U(0)=0$ and $U^{\prime}(1)=0$. For the AdS
Schwarzschild model ($f(u)=1-u^{4}$), the latter constraint comes from
imposing regularity at the horizon. It is possible to calculate $B_{c}$ by
numerically finding the value at which $U(u)$ satisfies these constraints.
There is an infinite tower of solutions to $B_{c}$, but we are only interested
in the lowest one, which is where the phase transition occurs. For further
details on solving this equation in the AdS Schwarzschild model, see
Ammon:2011je . For the hard wall model, $f(u)=1$ so the equation simplifies to
the extent that it can be solved analytically, the solution being a Bessel
function. Qualitatively the solutions for $U(u)$ in both models look very
similar and we are only interested in their numerical form. For the AdS
Schwarzschild model, we get $B_{c}\approx 5.1$, while we get $B_{c}\approx
5.8$777This is the zero of the Bessel function of the first kind
$J_{0}(\sqrt{B})$. for the hard wall model.

It should be noted that the solution 24 for $E_{x}$ agrees precisely (except
for the factor of $U(u)$) with the linear order solution for the order
parameter close to the upper critical magnetic field $H_{c2}$ in the theory of
type II superconductors, as seen in AbrikosovBook . It is also the result
found by Chernodub et al. in Chernodub:2011gs . Depending on the values of the
parameters $C_{n}$ and $k$ (to be determined by the higher order equations in
the perturbative expansion), $E_{x}$ corresponds to different inhomogeneous
functions in the $x,y$-plane. We are particularly interested in finding those
with lattice symmetries that represent evenly spaced vortices running in the
$z$ direction in the gauge theory.

###  3.4 The Abrikosov lattice solution

Before going beyond linear order, we discuss the possible solutions we can
expect. The number of coefficients specifying a configuration can make the
problem of finding the lowest energy solution unmanageable without making use
of some symmetries. We can argue that, since nothing in the setup is
explicitly breaking translational invariance in the $x,y$-directions, the
solution should be a highly symmetric lattice. A nice review of how lattices
can be formed from the Abrikosov solution (24) is given in RosensteinLi2010 .
There the authors explain that in order for $|E_{x}|$ to be a lattice
solution, the coefficients $C_{n}$ must have the same magnitude $|C_{n}|$ and
moreover be periodic in some integer $P$, that is, $C_{n}=C_{n+P}$.

In Abrikosov:1956sx , Abrikosov first studied the simplest solution, a square
lattice. In this case, $P=1$, implying that $C_{n}=C$ for all $n$, and
$k=\sqrt{2\pi B_{c}}$. Later Kleiner et al. in Kleiner:1964  generalised the
analysis by looking at $P=2$, with $C_{1}=\pm iC_{0}=\pm iC$. This choice of
coefficients specifies a general rhombic lattice, with the shape of the
rhombus controlled by varying $k$. In particular, a square lattice can be
obtained by choosing $k=\sqrt{\pi B_{c}}$. This square lattice is the same as
Abrikosov’s solution with $P=1$, but it is rotated by $\pi/4$ and translated.
A triangular lattice is obtained by choosing $k=3^{\frac{1}{4}}\sqrt{\pi
B_{c}}$.

To show how this works, we first substitute $P=2$ and $C_{1}=iC_{0}=iC$ into
the solution for $E_{x}$, which simplifies to

$\displaystyle E_{x}$
$\displaystyle=C\sum_{n=-\infty}^{\infty}e^{i\frac{\pi}{2}n^{2}-inky-\frac{1}{2}B_{c}\left(x-\frac{nk}{B_{c}}\right)^{2}}U(u)~{}.$
(26)

It is then easy to see the symmetries
$|E_{x}(x+[m+\frac{1}{2}q]L_{x},y+[n+\frac{1}{2}q]L_{y})|=|E_{x}(x,y)|$ for
integers $m$, $n$ and $q$. $L_{x}$ and $L_{y}$ are the lengths of the lattice
cell in the $x$ and $y$ directions, and are given by $L_{x}=2k/B_{c}$ and
$L_{y}=2\pi/k$. See figure 1.

Figure 1: A lattice cell, illustrating the meanings of $L_{x}$ and $L_{y}$ for
a fixed area cell.

We follow the approach of Kleiner et al, which is to compute the energy
density of the lattice for a range of values of the ratio
$L_{x}/L_{y}=k^{2}/\pi B_{c}$. This essentially means that we vary $k$. The
energy is computed numerically from the analytic expressions we obtain at each
order in the following sections. What we find agrees with their result that
the triangular lattice has the lowest energy of the $P\leq 2$ solutions. When
doing this, magnetic flux conservation is an important constraint. The total
applied magnetic field per unit area is constant, and each lattice cell
corresponds to a vortex with a single quantum of magnetic flux. This means
that when comparing the energy of different lattices, we should make sure that
they have the same magnetic flux per unit area, which in turn means that their
lattice cells have the same area. Fortunately with this ansatz that is always
the case since the area $L_{x}L_{y}=4\pi/B_{c}$ is independent of $k$.

In the following sections we calculate analytic expressions for the higher
order corrections to the gauge field. We keep $P$ and the coefficients $C_{n}$
general, except for imposing the periodicity condition $C_{n}=C_{n+P}$.

###  3.5 Higher order contributions to the energy

In order to find the ground state solution we must calculate the energy of the
superconducting solutions and compare them to the normal phase solution. We
can study the form of the energy as defined in section 2.5 to see how far we
must go in the perturbative expansion of the gauge fields. The energy has
terms that are quadratic and quartic in the gauge potential. The quartic term
ensures that the energy is bounded below, because it has a positive
coefficient. The quartic terms have lowest perturbative contributions of order
$\varepsilon^{4}$. One might expect contributions of order $\varepsilon^{3}$
coming from the zeroth order magnetic field contribution multiplied by three
first order corrections. However, from equation 19 it can be shown that such
terms do not arise. Thus we should expect to expand to third order in
$\mathcal{A}_{x,y}^{1,2}$ and fourth order in $\mathcal{A}_{x,y}^{3}$ .
However, it turns out that going to fourth order is not necessary because
inserting ansatz 19 into the action of equation 3, we find that the only
fourth order terms from $\mathcal{A}_{x,y}^{3}$ that appear at the fourth
order of the action are proportional to
$\sim\partial_{y}a^{(4)3}_{x}-\partial_{x}a^{(4)3}_{y}$. Here $a^{(4)3}_{x}$
and $a^{(4)3}_{y}$ are the fourth order corrections to $\mathcal{A}^{3}_{x}$
and $\mathcal{A}^{3}_{y}$, respectively. This term respects the lattice
symmetries, thus on performing the integration over the lattice cell to get
the free energy density, it vanishes by Stokes’ theorem.

We saw above that the parameters $k$ and $C_{n}$ in the solution 24 are not
fixed by the equations of motion to linear order. This is due to the fact that
to linear order, the different vortices do not interact. We can therefore not
expect to fix any of the coefficients $C_{n}$ or the spacing parameter $k$ at
this order. In fact, trying to calculate $\Delta\Omega$ to this order, which
has no quartic terms in $\mathcal{A}$, one finds that the free energy density
is not bounded below; increasing the overall magnitude of the condensate
always decreases $\Delta\Omega$. To see which configuration, that is, which
set of values for $C_{n}$ and $k$, is energetically favourable, we clearly
have to go beyond linear order.

###  3.6 Solving the equations to higher orders

In this section we solve the equations of motion up to third order in the
perturbation parameter.

The second order corrections to the gauge fields contribute to the potentials
$\mathcal{A}^{3}_{x}$ and $\mathcal{A}^{3}_{y}$, that is, $a^{3}_{x}$ and
$a^{3}_{y}$ in 19. These fields source the external magnetic field and the
magnetisation. We impose that these corrections must vanish at the AdS
boundary, so that the dual field theory has a constant applied magnetic field.
We find however that they do not vanish throughout the bulk. In particular
they develop non-vanishing subleading terms in the boundary expansion,
representing a magnetisation in the field theory.

In appendix A we explain how the equations for the Fourier modes of the fields
$a^{3}_{x}$ and $a^{3}_{y}$ can be decoupled. This yields the following
equations

$\displaystyle
u\partial_{u}\left(\frac{f}{u}\partial_{u}\hat{a}^{3}_{x,y}(m,n,u)\right)-\left(k^{2}n^{2}+\frac{4B_{c}^{2}m^{2}\pi^{2}}{k^{2}P^{2}}\right)\hat{a}^{3}_{x,y}(m,n,u)$
$\displaystyle+T_{x,y}e^{-\frac{k^{2}n^{2}}{4B_{c}}+\frac{inm\pi}{P}-\frac{B_{c}m^{2}\pi^{2}}{k^{2}P^{2}}}\left(\sum_{l=0}^{P-1}e^{\frac{2ilm\pi}{P}}\bar{C}_{l}C_{l+n}\right)U^{2}=0~{},$
(27)

where

$T_{x}=-i\frac{\sqrt{B_{c}\pi}}{P}n~{},~{}~{}~{}T_{y}=2i\frac{\pi^{3/2}B_{c}^{3/2}}{k^{2}P^{2}}m~{},$
(28)

and

$a^{3}_{x,y}(x,y,u)=\sum_{m}\sum_{n}e^{-i\frac{2\pi
mB_{c}}{Pk}x-inky}~{}\hat{a}^{3}_{x,y}(m,n,u)~{}.$ (29)

As before, $P$ defines the periodicity in the $C_{n}$. The parameters $m$ and
$n$ correspond to the Fourier space levels of these fields. In order to
calculate the solution $a_{x,y}^{3}(x,y,u)$ we will in theory need to solve
these equations for all values of $m$ and $n$. However, it will turn out to be
sufficient to only study the first few Fourier modes. The numerical procedure
for solving these will be explained in section 3.7

At third order we are studying the perturbative corrections to the condensate.
Here we calculate the corrections $e_{x}$ and $e_{y}$. It is reasonable to
assume that the answer is of the form

$\displaystyle\varepsilon E_{x}+\varepsilon^{3}e_{x}$
$\displaystyle=\varepsilon\sum_{n=-\infty}^{\infty}\left(C_{n}U(u)+\varepsilon^{2}c_{x,n}(u)\right)e^{-inky-\frac{1}{2}B_{c}\left(x-\frac{nk}{B_{c}}\right)^{2}}~{},$
(30) $\displaystyle\varepsilon E_{y}+\varepsilon^{3}e_{y}$
$\displaystyle=\varepsilon\sum_{n=-\infty}^{\infty}\left(-iC_{n}U(u)+\varepsilon^{2}c_{y,n}(u)\right)e^{-inky-\frac{1}{2}B_{c}\left(x-\frac{nk}{B_{c}}\right)^{2}}~{},$
(31)

where we have made use of equation 23 to relate the first order terms
$C_{n}U(u)$ in $\mathcal{E}_{x}$ and $\mathcal{E}_{y}$. $c_{(x,y),n}(u)$ is
the first perturbative correction to the condensate where the $u$ dependence
is a function of $n$ in contrast to the first order term.

We can write $e_{x}$ and $e_{y}$ in Fourier space, then use the three
condensate equations discussed in section 3 to calculate these corrections.
The one constraint equation can be used to decouple the other two equations.
We then have one equation for $c_{x,n}(u)$ and one for $c_{y,n}(u)$. Further
details are provided in appendix B.

###  3.7 Numerical solutions

Having separated the equations into ordinary differential equations in $u$ by
the method outlined in the appendices, we can now solve them numerically. Both
the second and third order equations take the same general form, given by

$\displaystyle
u\partial_{u}\left(\frac{f}{u}\partial_{u}\phi\right)+G(m,n)\phi+H(m,n,u)=0~{}.$
(32)

This equation can be solved numerically by picking some parameters for $C_{n}$
and $k$ that give a particular lattice and then using a shooting method to
integrate from $u=1$ (the horizon/hard wall cutoff) to $u=0$ (the AdS
boundary). It is an inhomogeneous second order differential equation, so there
are two integration constants. The first is fixed by imposing regularity at
the horizon or Neumann boundary conditions at the hard wall cutoff. This fixes
the value of $\partial_{u}\phi(1)$. The second constant is obtained by
demanding that $\phi(0)=0$, so that the fields vanish at the AdS boundary.
This vanishing corresponds to both the magnetic field strength corrections and
the source for the condensate being set to zero. We fulfil this boundary
condition by adjusting $\phi(1)$. Unlike in the case of the first order
equations, the equations here are not homogeneous and thus the source sets a
scale with which the value $\phi(1)$ can be compared. Changing $\phi(1)$ in
this case thus acts as more than just a scaling for the solution and so is
used as the tuning parameter to satisfy the UV constraint.

For all of the equations, we can implement this procedure for arbitrary
integers $m$ and $n$, corresponding to the different Fourier modes of the
gauge fields. This will then give a Fourier coefficient
$\hat{a}^{3}_{x,y}(m,n,u)$ that can be used to determine $a^{3}_{x,y}(x,y,u)$.
Fortunately we do not have to do the calculation for many different values of
$m$ and $n$, because as the values get large, the source term gets suppressed
exponentially. This can be seen in equation 27 for the second order terms and
is true also for the third order equation. For a vanishing source, the
equations for $\hat{a}^{3}_{x}$ or $\hat{a}^{3}_{y}$ have only the trivial
solution. This means that $\hat{a}^{3}_{x,y}(m,n,u)$ is negligibly small for
large $m$ or $n$, and we can therefore truncate the Fourier series for
$a^{3}_{x,y}$ beyond $m,n\approx 3$.

##  4 Results

###  4.1 Finding the minimum energy state

As explained above, we wish to find the values of the parameters $k$, $P$ and
$C_{n}=C_{n+P}$ that give the minimum energy state. These parameters define
the shape of the lattice. Our analysis is only valid for $B$ slightly above
$B_{c}$, where $B_{c}$ was determined in section 3.3. The first step is thus
to pick a value for $B$ in this vicinity. We then choose a set of lattice
parameters that give us the lattice solution we wish to consider. As mentioned
in RosensteinLi2010 , for lattice solutions all the $C_{n}$ must have the same
magnitude $C$. We can therefore fix $C_{n}$ up to the normalisation $C$, along
with a value of $k$, according to the discussion in section 3.4. We then
substitute these values into the energy density that was defined in section
2.5. It takes the form $\Delta\Omega=a_{1}\varepsilon
C+a_{2}\varepsilon^{2}C^{2}+\dots~{}$. At this point we see that we can
redefine $C$ by absorbing a factor of $\varepsilon$, which we call
$C_{\varepsilon}$. $C_{\varepsilon}$ is the only parameter left unfixed up to
this point in the analysis. Here the $a_{i}$ are values that are calculated
numerically from substituting the solutions to the equations of motion into
the expression for the energy derived in appendix C. $\Delta\Omega$ forms a
Mexican hat potential, which is easy to minimise numerically. An illustration
of this procedure is shown in figure 2.

Figure 2: The change in energy density in units of temperature as a function
of $C_{\varepsilon}$, the overall condensate scale. The leftmost curve
corresponds to $B=B_{c}$, which is never negative for nonzero condensate.
Curves for $B<B_{c}$ are similar. Increasing $B$ beyond $B_{c}$ yields the
curves to the right, and we see the formation of a clear minimum of the energy
that is lower than the energy of the normal phase. The dashed line traces out
the minimum of each of these curves, which corresponds to the energetically
preferred size of the condensate as a function of $B$. This plot was generated
in the AdS Schwarzschild model for $P=2$ and $k=3^{\frac{1}{4}}\sqrt{\pi B}$,
corresponding to a triangular lattice. $B$ takes the values $B\approx
B_{c},1.04B_{c},1.07B_{c},1.1B_{c}$ from left to right. Changing $P$ and $k$
to correspond to different lattices or using the hard wall model yields
qualitatively similar results.

The plot in figure 3 shows the energy-minimising value of $C_{\varepsilon}$ as
a function of magnetic field near the phase transition at $B_{c}$. It shows
that $C_{\varepsilon}\sim(B-B_{c})^{\frac{1}{2}}$, so the condensate888Note
that only the combination $\varepsilon C$ is physically relevant, not $C$ or
$\varepsilon$ independently. has a critical exponent of $1/2$. A fit to the
numerical data for the triangular lattice gives that
$C_{\varepsilon}=0.58(B-B_{c})^{\frac{1}{2}}$ in the AdS Schwarzschild model
and $C_{\varepsilon}=0.53(B-B_{c})^{\frac{1}{2}}$ in the hard wall model.

Figure 3: $C_{\varepsilon}\sim$ the overall condensate size for the AdS
Schwarzschild solution in units of the temperature, as a function of the
external magnetic field $B$. For $B<B_{c}$, the condensate is zero, and for
$B$ slightly above $B_{c}$, we see a $(B-B_{c})^{\frac{1}{2}}$ scaling
behaviour. This plot was generated for $P=2$ and $k=3^{\frac{1}{4}}\sqrt{\pi
B}$, corresponding to a triangular lattice. The plot for different lattices in
both the AdS Schwarzschild and hard wall models is the same, up to a scaling
of the $B$ and $C_{\varepsilon}$ axes. For the triangular lattice, the AdS
Schwarzschild model has scaling behaviour
$C_{\varepsilon}=0.58(B-5.1)^{\frac{1}{2}}$ and the hard wall model has
$C_{\varepsilon}=0.53(B-5.8)^{\frac{1}{2}}$.

Having minimised with respect to $C_{\varepsilon}$ for a given value of $B$
and a given lattice configuration, we can plot the difference in the energy
between the normal and superconducting states. Figure 4 shows $\Delta\Omega$,
the difference between the energy density in the superconducting and normal
phases, as a function of external magnetic field for two different lattices.
The first lattice is square, and the second is triangular. Both are described
in section 3.4.

Figure 4: The change in energy density (compared to the normal phase) for the
triangular and square lattices as the external applied magnetic field is
varied. The phase transition happens at $B_{c}\approx 5.1$, which is where the
coordinate axes are centred.
$\Delta\Omega_{\text{square}}-\Delta\Omega_{\text{triangle}}$ is so small that
the two plots are almost on top of each other. This is for the AdS
Schwarzschild model, but the plots for the hard wall model are identical
except for the scale on the axes. In the hard wall model, $B_{c}\approx 5.8$.

The curves in figure 4 are the result of calculations in the AdS Schwarzschild
model, but we get the same results up to a rescaling of the axes for the hard
wall model. In the AdS Schwarzschild model, the critical magnetic field
$B_{c}\approx 5.1$, while in the hard wall model $B_{c}\approx 5.8$. Each
curve shows that the free energy density is proportional to
$\left(B-B_{c}\right)^{2}$. This shows that the phase transition is second
order, as expected if one looks at the analogous case in Ginzburg-Landau
theory. There one can show (TinkhamBook ) that the free energy is proportional
to $\left(T-T_{c}\right)^{2}$, where $T_{c}$ is the phase transition critical
temperature.

###  4.2 An analysis of $P=2$ solutions

We now specialise to the case where the periodicity of the $C_{n}$ is $P=2$.
This describes a general rhombic lattice solution which includes both the
triangular and square lattices. The $P=1$ square lattice can be found within
the $P=2$ solutions up to translation and rotation. We here perform the
analysis done in Kleiner:1964  as described at the end of section 3.4.

The energy difference as a function of $R$ is plotted in figure 5. By looking
at the form of equation 26, it is possible to see that the triangular lattice
occurs for $R=L_{x}/L_{y}=\sqrt{3}$ and $R=1/\sqrt{3}$. In general, $R$ and
$1/R$ give the same lattice but with the $x$ and $y$ directions flipped. This
is why figure 5 displays the symmetry $\Delta\Omega(R)=\Delta\Omega(1/R)$. The
triangular lattice corresponds to a global minimum of the energy as a function
of $R$, as seen from the figure. There is a local maximum for the square
lattice, which is when $R=1$. As $R\rightarrow\infty$ (or $R\rightarrow 0$),
the free energy increases. Intuitively one can understand this by making use
of the properties of Abrikosov vortices that we understand from type II
superconductors. These vortices repel. Since $R\rightarrow\infty$ and
$R\rightarrow 0$ correspond to elongating the rhombic lattice cell (while
keeping the area constant) neighbouring vortices are squeezed together, and
since they repel, this is energetically unfavourable.

Figure 5: The change in free energy density as a function of $R=L_{x}/L_{y}$,
the ratio of side lengths of a constant area lattice cell. This plot is for
the AdS Schwarzschild model, but the plot for the hard wall model is identical
up to a rescaling of the axes. When $R=1$, the lattice is square and the free
energy achieves a local maximum. When $R=\sqrt{3}$ and $1/{\sqrt{3}}$, the
lattice is triangular and the free energy is at a global minimum. Note that
the plot has the symmetry $\Delta\Omega(R)=\Delta\Omega(1/R)$, which simply
corresponds to swapping the $x,y$-axes.

We can calculate the condensate in the minimum energy state using equation 12.
The result, to linear order in $\varepsilon$, is

$\langle J^{+}_{x}\rangle\equiv\frac{\delta S_{\mathrm{on-shell}}}{\delta
E^{(0)}_{x}}=\frac{L}{2\hat{g}^{2}}U_{\mathrm{sub}}C_{\varepsilon}\sum_{n=-\infty}^{\infty}e^{-i\frac{\pi}{2}n^{2}+inky-\frac{1}{2}B_{c}\left(x-\frac{nk}{B_{c}}\right)^{2}}$
(33)

The AdS radius can be related to field theory quantities through the relation
$L^{4}=2\lambda{\alpha^{\prime}}^{2}$, where $\lambda$ is the ’t Hooft
coupling and $\alpha^{\prime}$ the string tension. The factor
$U_{\mathrm{sub}}$ is equal to the subleading term in the boundary expansion
of $U(u)$. Using equation 25 it is possible to show that

$U_{\mathrm{sub}}=B_{c}\int_{0}^{u_{H}}\frac{U(u)}{u}du~{},$ (34)

so it can be determined numerically. In figure 6 we present the contour plot
of
$3^{\frac{1}{4}}\sqrt{8}\frac{\hat{g}^{4}}{L^{2}U_{\mathrm{sub}}^{2}C_{\varepsilon}^{2}}\left|\langle
J^{+}_{x}\rangle\right|^{2}$, the modulus squared of the condensate in the
$x,y$-plane for the minimum energy solution corresponding to the triangular
lattice. The factors are chosen so that the maximum value is 1. Substituting
in the numerical values, we find that the maximum value the condensate takes
is $|\langle
J^{+}_{x}\rangle|=1.0\frac{L}{\hat{g}^{2}}\left(B-B_{c}\right)^{\frac{1}{2}}$
for the AdS Schwarzschild model, where $B_{c}\approx 5.1$, and $|\langle
J^{+}_{x}\rangle|=1.3\frac{L}{\hat{g}^{2}}\left(B-B_{c}\right)^{\frac{1}{2}}$
for the hard wall model, where $B_{c}\approx 5.8$.

Figure 6: A contour plot of
$3^{\frac{1}{4}}\sqrt{8}\frac{\hat{g}^{4}}{L^{2}U_{\mathrm{sub}}^{2}C_{\varepsilon}^{2}}\left|\langle
J^{+}_{x}\rangle\right|^{2}$, the modulus squared of the field theory
condensate dual to $E_{x}$ in the ground state triangular lattice. At the
center of the dark vortices, the condensate vanishes.

We could also plot the magnetisation of the ground state, which is found from
the normalisable term in the boundary value expansion of
$\partial_{x}a^{3}_{y}-\partial_{y}a^{3}_{x}$. However, it takes the same form
as the condensate and the numerics indicate that it differs only up to a
scale.

##  5 Conclusion

In this work we have found a likely ground state for the black hole Yang-Mills
instability analysed in Ammon:2011je . The solution, being of a lattice form,
clearly has much potential for analysis in condensed matter models, where the
breaking of translational invariance has already been shown to be very
important in getting realistic phenomenology. As we have explained, it also
has possible implications for heavy ion collider physics. There are a number
of interesting areas where we could apply similar techniques and perform
further calculations to elucidate the phenomenology of the ground state found
here.

It is certainly important to understand exactly how universal this result is.
We have seen that the difference between the normal phase and superconducting
phase energy density of this solution is, up to a scale, independent of the
background geometry in the two models that we have studied here. It would be
very useful to understand precisely where this universality stems from and
find how much we can deform the gravity solutions until the Abrikosov lattice
is no longer the ground state.

In the present work we have analysed lattices with $P=1$ and $P=2$,
corresponding to square and rhombic forms, respectively. Going to $P=3$
requires a large increase in computational power. While this would be an
interesting further calculation, the analogous cases of type II
superconductivity and the model of Chernodub et al point to the triangular
lattice being the true ground state. We thus expect higher $P$ lattices to be
energetically disfavoured.

Having found this solution, there are some extensions that can be made. It
will now be possible to study time dependent fluctuations about the ground
state. In order to do this we would have to introduce a second perturbative
parameter in addition to the parameter $\varepsilon$ used in the current work.
This would be analogous to the parameter $\alpha^{\prime}$ in a D-brane
construction. This would allow us to study the transport properties of the
lattice ground state, by looking at current-current correlation functions. If
we wish to study the effect of the lattice on the shear viscosity to entropy
density ratio, we would have to introduce gravitational back reaction in our
model. Clearly this will be a much more involved calculation, with many non-
linear couplings, but if we want to study the theory in a more realistic
scenario, where the stress energy tensor also has a lattice structure, such a
calculation would clearly be important.

It is expected that if the QCD vacuum is unstable to $\rho$ meson condensation
in extremely off-centre heavy ion collisions, then the timescale of the
instability would not be enough to form a well-defined lattice. Abrikosov
vortices may form, but the magnetic field would likely drop below the critical
value before they had time to arrange themselves into a lattice. It would be
very interesting to perform a real-time calculation in order study the
formation of the vortices and their movements as the magnetic field increased
and decreased through the lifetime of a single off-centre collision.

## Acknowledgements

Y.Y.B. was supported by MPG-CAS Doctoral Promotion Programme. The work of J.S.
was supported by the European Union through a Marie Curie Fellowship. The work
of J.E. and M.S. was supported in part by the DFG cluster of excellence
‘Origin and Structure of the Universe’ (www.universe-cluster.de). We thank
Maxim Chernodub for useful discussions.

##  Appendix A Deriving the equations for $a^{3}_{x,y}$

We substitute the ansatz 19 into the full equations of motion and neglect
terms beyond quadratic order in $\varepsilon$. To get rid of all appearances
of $E_{y}$, we use the relation that $E_{y}=-iE_{x}$ from 23. Then we find
that there are only three equations in which the fluctuations $a^{3}_{x,y}$
appear. We focus on those three equations.

The simplest equation of the three is the constraint equation, which came from
the equation of motion for $A^{3}_{u}$. To quadratic order, this equation is
simply

$\partial_{u}\partial_{x}a^{3}_{x}+\partial_{u}\partial_{y}a^{3}_{y}=0~{}.$
(35)

The first thing to do is integrate by $u$. This gives an integration constant,
but by the fact that both $a^{3}_{x}$ and $a^{3}_{y}$ must vanish at $u=0$,
this integration constant vanishes. So the even simpler constraint is

$\partial_{x}a^{3}_{x}+\partial_{y}a^{3}_{y}=0~{}.$ (36)

This is all we need to decouple the other two equations in $a^{3}_{x}$ and
$a^{3}_{y}$. These equations now become

$\displaystyle 0=$
$\displaystyle\frac{3}{2}E_{x}\partial_{y}\bar{E}_{x}+\frac{3}{2}\bar{E}_{x}\partial_{y}E_{x}-\frac{1}{2}iE_{x}\partial_{x}\bar{E}_{x}+\frac{1}{2}i\bar{E}_{x}\partial_{x}E_{x}$
$\displaystyle+u\partial_{u}\left(\frac{f}{u}\partial_{u}a^{3}_{x}\right)+\partial^{2}_{y}a^{3}_{x}+\partial^{2}_{x}a^{3}_{x}~{},$
(37) $\displaystyle 0=$ $\displaystyle-
B_{c}x\bar{E}_{x}E_{x}-\frac{1}{2}iE_{x}\partial_{y}\bar{E}_{x}+\frac{1}{2}i\bar{E}_{x}\partial_{y}E_{x}$
$\displaystyle-\frac{3}{2}E_{x}\partial_{x}\bar{E}_{x}-\frac{3}{2}\bar{E}_{x}\partial_{x}E_{x}+u\partial_{u}\left(\frac{f}{u}\partial_{u}a^{3}_{y}\right)+\partial^{2}_{y}a^{3}_{y}+\partial^{2}_{x}a^{3}_{y}~{},$
(38)

which are partial differential equations with sources that come from the
linear order solutions. These two equations only differ by their source terms,
so we will focus on $a^{3}_{x}$. $a^{3}_{y}$ should be similar. Using the
expression 24, we can see that the source term is periodic with $y\sim
y+\frac{2\pi}{k}$. $a^{3}_{x}$ must have the same periodicity, so we can write
it as a Fourier series,

$\displaystyle
a^{3}_{x}(x,y,u)=\sum_{n=-\infty}^{\infty}e^{-inky}\tilde{a}^{3}_{x}(x,n,u)~{}.$
(39)

The equation becomes

$\displaystyle\sum_{m}-ie^{-\frac{1}{2}B_{c}\left(-\frac{km}{B_{c}}+x\right)^{2}-\frac{1}{2}B_{c}\left(-\frac{k(n+m)}{B_{c}}+x\right)^{2}}kn\bar{C}_{m}C_{n+m}U^{2}$
$\displaystyle-k^{2}n^{2}\tilde{a}^{3}_{x}+u\partial_{u}\left(\frac{f}{u}\partial_{u}\tilde{a}^{3}_{x}\right)+\partial^{2}_{x}\tilde{a}^{3}_{x}$
$\displaystyle=0~{}.$ (40)

We notice that the source term in this equation is periodic in the
$x$-direction; $x\sim x+\frac{Pk}{B_{c}}$. This lets us expand
$\tilde{a}^{3}_{x}$ as a Fourier series in $x$ as well:

$\displaystyle\tilde{a}^{3}_{x}=\sum_{m}e^{-i\frac{2\pi
mB_{c}}{Pk}x}\hat{a}^{3}_{x}(m,n,u)~{}.$ (41)

Writing the source term as a series lets us then obtain the equation 27 for
the coefficients $\hat{a}^{3}_{x}(m,n,u)$.

Calling the source term $S(x)$, the naïve way of finding its Fourier
coefficients is to use the formula

$\displaystyle\tilde{S}_{n}=\frac{B_{c}}{Pk}\int_{0}^{\frac{Pk}{B_{c}}}e^{i\frac{2\pi
nB_{c}}{Pk}x}S(x)~{}.$ (42)

However, the source terms contains Gaussians, and those are much easier to
integrate when the domain of integration is the entire real line. So we do the
following trick. Doing a continuous Fourier transform on a periodic function
gives a sum of $\delta$-functions,

$\displaystyle\int dx~{}e^{ipx}S(x)$ $\displaystyle=\int
dx~{}e^{ipx}\sum_{m}e^{-i\frac{2\pi mB_{c}}{Pk}x}\tilde{S}_{n}$
$\displaystyle=2\pi\sum_{m}\tilde{S}_{n}\delta\left(p-\frac{2\pi
mB_{c}}{Pk}\right)~{}.$ (43)

The coefficients in front of the $\delta$-functions are what we are looking
for. We get

$\displaystyle\int dx~{}e^{ipx}S(x)$
$\displaystyle=-\sqrt{\frac{\pi}{B_{c}}}\sum_{m,n}ie^{-\frac{k^{2}n^{2}}{4B_{c}}+\frac{ikmp}{B_{c}}+\frac{iknp}{2B_{c}}-\frac{p^{2}}{4B_{c}}}kn\bar{C}_{m}C_{m+n}U^{2}~{}.$
(44)

Using

$\displaystyle\sum_{m=-\infty}^{\infty}f(m)=\sum_{m=-\infty}^{\infty}\sum_{l=0}^{P-1}f(Pm+l)$
(45)

and then using the symmetry $C_{i+P}=C_{i}$, the only $m$-dependence remaining
in the sum comes from $e^{\frac{ikPmp}{B_{c}}}$. Making use of the identity

$\sum_{m=-\infty}^{\infty}e^{imq}=2\pi\sum_{m=-\infty}^{\infty}\delta(q-2\pi
m)$ (46)

and $\delta(\alpha x)=\frac{\delta(x)}{|\alpha|}$ gives us the sum over
$\delta$-functions from 43. Then we can simply read off the coefficients
$\tilde{S}_{n}$. This gives us the equation 27.

##  Appendix B Deriving the equations for $c_{x,n}$, $c_{y,n}$

The third order equations of motion are

$\displaystyle 0$
$\displaystyle=ia^{3}_{x}\partial_{u}E_{x}+a^{3}_{y}\partial_{u}E_{x}-iE_{x}\partial_{u}a^{3}_{x}-E_{x}\partial_{u}a^{3}_{y}+iB_{c}x\partial_{u}e_{y}+\partial_{y}\partial_{u}e_{y}+\partial_{x}\partial_{u}e_{x}~{},$
(47) $\displaystyle 0$
$\displaystyle=-iB_{c}xa^{3}_{x}E_{x}-2B_{c}xa^{3}_{y}E_{x}-\bar{E}_{x}E_{x}^{2}-a^{3}_{x}\partial_{y}E_{x}+2ia^{3}_{y}\partial_{y}E_{x}-a^{3}_{y}\partial_{x}E_{x}$
$\displaystyle~{}-2E_{x}\partial_{y}a^{3}_{x}+iE_{x}\partial_{y}a^{3}_{y}+E_{x}\partial_{x}a^{3}_{y}+iB_{c}e_{y}-iB_{c}x\partial_{x}e_{y}-\partial_{x}\partial_{y}e_{y}$
$\displaystyle~{}-B_{c}^{2}x^{2}e_{x}+2iB_{c}x\partial_{y}e_{x}+\partial^{2}_{y}e_{x}+u\partial_{u}\left(\frac{f}{u}\partial_{u}e_{x}\right)$
(48) $\displaystyle 0$
$\displaystyle=B_{c}xa^{3}_{x}E_{x}+i\bar{E}_{x}E_{x}^{2}-ia^{3}_{x}\partial_{y}E_{x}+2a^{3}_{x}\partial_{x}E_{x}-ia^{3}_{y}\partial_{x}E_{x}$
$\displaystyle~{}+iE_{x}\partial_{y}a^{3}_{x}+E_{x}\partial_{x}a^{3}_{x}-2iE_{x}\partial_{x}a^{3}_{y}-2iB_{c}e_{x}-iB_{c}x\partial_{x}e_{x}-\partial_{x}\partial_{y}e_{x}+\partial^{2}_{x}e_{y}+u\partial_{u}\left(\frac{f}{u}\partial_{u}e_{y}\right)~{}.$
(49)

The first of these is the constraint equation. We use it to relate $e_{x}$ and
$e_{y}$. In order to do this, we first simplify it by noticing that, since

$\displaystyle e_{y}$
$\displaystyle=\sum_{n=-\infty}^{\infty}c_{y,n}(u)e^{-inky-\frac{1}{2}B_{c}\left(x-\frac{nk}{B_{c}}\right)^{2}}~{},$
(50)

we have that
$iB_{c}x\partial_{u}e_{y}+\partial_{y}\partial_{u}e_{y}=-i\partial_{x}\partial_{u}e_{y}$.
We can then integrate the entire equation by $u$, imposing vanishing boundary
conditions at the AdS boundary. The constraint equation then simplifies to

$\displaystyle 0$
$\displaystyle=-2i\frac{E_{x}}{U}J_{x}-2\frac{E_{x}}{U}J_{y}+ia^{3}_{x}E_{x}+a^{3}_{y}E_{x}+\partial_{x}e_{x}-i\partial_{x}e_{y}~{},$
(51)

where

$\displaystyle
J_{x,y}(x,y,u)=\int_{0}^{u}U(\tilde{u})\partial_{\tilde{u}}a^{3}_{x,y}(x,y,\tilde{u})d\tilde{u}~{}.$
(52)

This allows us to eliminate $e_{x}$ in equation 49 (after differentiating it
by $x$). We write each function as a Fourier series in $y$ and find an
equation for the coefficients $c_{y,n}$. At this point the equation still has
an $x$ dependence, which can be eliminated by multiplying the equation by
$(nk-B_{c}x)$ to make it an even function in $x$ and then integrating
$\int_{-\infty}^{\infty}dx$. In doing so we use the solution for $E_{x}$ and
the form for $e_{y}$ given by 50, as well as the Fourier series representation
of the other functions. Once this is done, we are left with an equation for
$e_{y}$ in the form 32.

The resulting equation for $c_{y,n}$ is

$\displaystyle 0$
$\displaystyle=\sum_{q,r=-\infty}^{\infty}\left\\{e^{-\frac{2\pi
q\left(ik^{2}P(n-r)+B_{c}\pi
q\right)}{k^{2}P^{2}}}\left[\frac{C_{n-r}\left(-2\left(k^{2}Pr+2iB_{c}\pi
q\right)\hat{J}_{x,qr}\right)}{kP}\right.\right.$
$\displaystyle\left.+\frac{C_{n-r}\left(\left(2ik^{2}Pr-4B_{c}\pi
q\right)\hat{J}_{y,qr}+\left(2iB_{c}\pi
q\hat{a}^{3}_{x}+\left(-ik^{2}Pr+4B_{c}\pi
q\right)\hat{a}^{3}_{y}\right)U\right)}{kP}\right]$
$\displaystyle\left.-\frac{ie^{-\frac{k^{2}\left(3r^{2}-3rq+q^{2}\right)}{3B_{c}}}\left(3B_{c}+2k^{2}q(-2r+q)\right)\bar{C}_{n+q}C_{n+r}C_{n-r+q}U^{3}}{3\sqrt{3}B_{c}}\right\\}$
$\displaystyle-
B_{c}c_{y,n}+u\partial_{u}\left(\frac{f}{u}\partial_{u}c_{y,n}\right)~{},$
(53)

where

$\displaystyle\hat{J}_{i,qr}(u)$
$\displaystyle=\int_{0}^{u}U(\tilde{u})\partial_{\tilde{u}}\hat{a}^{3}_{i}(q,r,\tilde{u})d\tilde{u}~{},$
(54)

for $i=x,y$.

A similar procedure gives the constraint equation in terms of the
coefficients,

$\displaystyle 0$ $\displaystyle=c_{x,n}(u)-ic_{y,n}(u)$
$\displaystyle~{}~{}~{}+\frac{1}{PkB_{c}}\sum_{q,r=-\infty}^{\infty}\Biggl{\\{}e^{-\frac{2\pi
q\left(ik^{2}P(n-r)+B_{c}\pi q\right)}{k^{2}P^{2}}}\left(-ik^{2}Pr+2B_{c}\pi
q\right)C_{n-r}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\left(2\hat{J}_{x,qr}-2i\hat{J}_{y,qr}-(\hat{a}^{3}_{x,qr}-i\hat{a}^{3}_{y,qr})U(u)\right)\Biggr{\\}}~{}.$
(55)

Once the coefficients $c_{y,n}$ are found, we use this to calculate $c_{x,n}$.

##  Appendix C Calculating the energy

The difference between the energy of the superconducting phase and that of the
normal phase is

$\displaystyle\Delta\mathcal{F}$ $\displaystyle=\frac{1}{4\hat{g}^{2}}\int
d^{5}x\sqrt{-g}\left(\left.F^{a}_{\mu\nu}F^{a\mu\nu}\right|_{superconducting}-\left.F^{a}_{\mu\nu}F^{a\mu\nu}\right|_{normal}\right)~{}.$
(56)

Note that for the AdS Schwarzschild model we implicitly divided by the
temperature to make the energy dimensionless. We calculate the energy density
by averaging over the domain $0\leq y<\frac{2\pi}{k}$, $0\leq
x<\frac{Pk}{B_{c}}$, $0\leq u\leq 1$ and $t,z\in\mathbb{R}$. Since the
integrand is independent of $t$ and $z$, the averaging amounts to simply
dropping the integration over those variables. In the following expression we
use

$\displaystyle\mathcal{E}_{x,y}=\mathcal{A}^{1}_{x,y}+i\mathcal{A}^{2}_{x,y}=\sum_{n}\mathcal{C}_{(x,y),n}(u)e^{-ikny-\frac{1}{2}B_{c}\left(x-\frac{nk}{B_{c}}\right)^{2}}~{},$
(57)

we write $\mathcal{A}^{3}_{x}=a^{3}_{x}$ and
$\mathcal{A}^{3}_{y}=xB+a^{3}_{y}$, and call the averaged energy
$\Delta\Omega$. The result is

$\displaystyle 4\hat{g}^{2}\Delta\Omega$ $\displaystyle=\int
du\left\\{\Omega_{1}(u)+\sum_{m,n=-\infty}^{\infty}\left[\Omega_{2}(m,n,u)+\Omega_{3}(m,n,u)+\Omega_{4}(m,n,u)\right]\right.$
$\displaystyle~{}\left.\sum_{m,n,p,q=-\infty}^{\infty}\Omega_{5}(m,n,q,r,u)\right\\}~{},$
(58)

where

$\displaystyle\Omega_{1}=$ $\displaystyle\frac{\sqrt{\pi
B}}{kPu}\sum_{l=0}^{P-1}\frac{B}{2}\left(\sum_{j=x,y}\left(f\partial_{u}\bar{\mathcal{C}}_{j,l}\partial_{u}\mathcal{C}_{j,l}+\bar{\mathcal{C}}_{j,l}\mathcal{C}_{j,l}\right)+3(i\bar{\mathcal{C}}_{y,l}\mathcal{C}_{x,l}-i\bar{\mathcal{C}}_{x,l}\mathcal{C}_{y,l})\right)~{},$
(59) $\displaystyle\Omega_{2}=$
$\displaystyle\frac{1}{u}\left\|kn\hat{a}^{3}_{x}(m,n,u)-\frac{2Bm\pi}{kP}\hat{a}^{3}_{y}(m,n,u)\right\|^{2}+\frac{f}{u}\sum_{j=x,y}\left\|\partial_{u}\hat{a}^{3}_{j}(m,n,u)\right\|^{2}~{},$
(60) $\displaystyle\Omega_{3}=$ $\displaystyle\frac{\sqrt{\pi
B}}{2k^{2}P^{2}u}\sum_{l=0}^{P-1}e^{-\frac{k^{2}m^{2}}{4B}-\frac{i(2l+m)n\pi}{P}-\frac{Bn^{2}\pi^{2}}{k^{2}P^{2}}}\Bigl{(}\left(3k^{2}mP+2iBn\pi\right)\hat{a}^{3}_{x}(n,m,u)\bar{\mathcal{C}}_{x,l+m}\mathcal{C}_{y,l}$
$\displaystyle+\hat{a}^{3}_{x}(n,-m,u)\bar{\mathcal{C}}_{y,l}\left(\left(3k^{2}mP+2iBn\pi\right)\mathcal{C}_{x,l+m}-2ik^{2}mP\mathcal{C}_{y,l+m}\right)$
$\displaystyle+\hat{a}^{3}_{y}(n,-m,u)\mathcal{C}_{x,l+m}\left(-4iBn\pi\bar{\mathcal{C}}_{x,l}+\left(ik^{2}mP+6Bn\pi\right)\bar{\mathcal{C}}_{y,l}\right)$
$\displaystyle+\hat{a}^{3}_{y}(n,m,u)\bar{\mathcal{C}}_{x,l+m}\mathcal{C}_{y,l}\left(-ik^{2}mP-6Bn\pi\right)\Bigr{)}~{},$
(61) $\displaystyle\Omega_{4}=$ $\displaystyle-\frac{1}{4kPu}\sqrt{\frac{\pi
B}{2}}e^{-\frac{k^{2}\left(m^{2}+n^{2}\right)}{2B}}\times$
$\displaystyle\sum_{l=0}^{P-1}\left(\bar{\mathcal{C}}_{y,l+m}\bar{\mathcal{C}}_{y,l+n}\mathcal{C}_{x,l}\mathcal{C}_{x,l+m+n}-2\bar{\mathcal{C}}_{x,l+m}\bar{\mathcal{C}}_{y,l+n}\mathcal{C}_{x,l+m+n}\mathcal{C}_{y,l}+\bar{\mathcal{C}}_{x,l}\bar{\mathcal{C}}_{x,l+m+n}\mathcal{C}_{y,l+m}\mathcal{C}_{y,l+n}\right)~{},$
(62) $\displaystyle\Omega_{5}=$ $\displaystyle\frac{\sqrt{\pi
B}}{Pku}\sum_{l=0}^{P-1}e^{-\frac{k^{2}m^{2}}{4B}-\frac{i(2l+m)n\pi}{P}-\frac{Bn^{2}\pi^{2}}{k^{2}P^{2}}}\times$
$\displaystyle\left(\hat{a}^{3}_{y}(n-q,-(m+r),u)\hat{a}^{3}_{y}(q,r,u)\bar{\mathcal{C}}_{x,l}\mathcal{C}_{x,l+m}-\hat{a}^{3}_{x}(n-q,r,u)\hat{a}^{3}_{y}(q,m-r,u)\bar{\mathcal{C}}_{x,l+m}\mathcal{C}_{y,l}\right.$
$\displaystyle\left.-\hat{a}^{3}_{x}(n-q,r,u)\hat{a}^{3}_{y}(q,-(m+r),u)\bar{\mathcal{C}}_{y,l}\mathcal{C}_{x,l+m}+\hat{a}^{3}_{x}(n-q,-(m+r),u)\hat{a}^{3}_{x}(q,r,u)\bar{\mathcal{C}}_{y,l}\mathcal{C}_{y,l+m}\right)~{}.$
(63)

In these expressions, $\mathcal{C}_{x,n}$ and $\mathcal{C}_{y,n}$ are
functions of $u$. Their complex conjugates are given by
$\bar{\mathcal{C}}_{x,n}$ and $\bar{\mathcal{C}}_{y,n}$, respectively. All the
infinite sums in the energy C can be terminated at a small finite value
because of exponential suppression in the $\Omega_{1\dots 5}$ terms.

In deriving this, it helps to make use of the formulae

$\displaystyle\int_{0}^{L}dx$
$\displaystyle\sum_{m=-\infty}^{\infty}e^{-\frac{B_{c}}{2}\left(x-mL\right)^{2}}=\int_{-\infty}^{\infty}dxe^{-\frac{B_{c}}{2}x^{2}}~{},$
(64) $\displaystyle\int_{0}^{L}dx$
$\displaystyle\sum_{m,n=-\infty}^{\infty}e^{-\frac{B_{c}}{2}\left(x-\frac{mL}{P}\right)^{2}-\frac{B_{c}}{2}\left(x-\frac{nL}{P}\right)^{2}}h(x,m,n)$
(65)
$\displaystyle=\int_{-\infty}^{\infty}dx\sum_{l=0}^{P-1}\sum_{m=-\infty}^{\infty}e^{-\frac{B_{c}}{2}\left(x-\frac{mL}{P}\right)^{2}-\frac{B_{c}}{2}\left(x-\frac{lL}{P}\right)^{2}}h(x,m,l)~{},$
(66)

where the latter is valid whenever $h(x,m,n)=h(x+L,m+P,n+P)$.

## References

  * (1) M. Ammon, J. Erdmenger, P. Kerner, and M. Strydom, Black Hole Instability Induced by a Magnetic Field, Phys.Lett. B706 (2011) 94–99, [arXiv:1106.4551].
  * (2) A. Abrikosov, On the Magnetic properties of superconductors of the second group, Sov.Phys.JETP 5 (1957) 1174–1182.
  * (3) G. T. Horowitz, J. E. Santos, and D. Tong, Optical Conductivity with Holographic Lattices, JHEP 1207 (2012) 168, [arXiv:1204.0519].
  * (4) R. Flauger, E. Pajer, and S. Papanikolaou, A Striped Holographic Superconductor, Phys.Rev. D83 (2011) 064009, [arXiv:1010.1775].
  * (5) M. Chernodub, Superconductivity of QCD vacuum in strong magnetic field, Phys.Rev. D82 (2010) 085011, [arXiv:1008.1055].
  * (6) M. Chernodub, Electromagnetically superconducting phase of QCD vacuum induced by strong magnetic field, AIP Conf.Proc. 1343 (2011) 149–151, [arXiv:1011.2658].
  * (7) M. Chernodub, J. Van Doorsselaere, and H. Verschelde, Electromagnetically superconducting phase of vacuum in strong magnetic field: structure of superconductor and superfluid vortex lattices in the ground state, arXiv:1111.4401.
  * (8) N. Nielsen and P. Olesen, An Unstable Yang-Mills Field Mode, Nucl.Phys. B144 (1978) 376.
  * (9) J. Ambjorn and P. Olesen, On Electroweak Magnetism, Nucl.Phys. B315 (1989) 606.
  * (10) J. Ambjorn and P. Olesen, A Condensate Solution Of The Electroweak Theory Which Interpolates Between The Broken And The Symmetric Phase, Nucl.Phys. B330 (1990) 193.
  * (11) J. Ambjorn and P. Olesen, Electroweak Magnetism: Theory And Application, Int.J.Mod.Phys. A5 (1990) 4525–4558.
  * (12) S. K. Domokos and J. A. Harvey, Baryon number-induced Chern-Simons couplings of vector and axial-vector mesons in holographic QCD, Phys.Rev.Lett. 99 (2007) 141602, [arXiv:0704.1604].
  * (13) S. Nakamura, H. Ooguri, and C.-S. Park, Gravity Dual of Spatially Modulated Phase, Phys.Rev. D81 (2010) 044018, [arXiv:0911.0679].
  * (14) W.-y. Chuang, S.-H. Dai, S. Kawamoto, F.-L. Lin, and C.-P. Yeh, Dynamical Instability of Holographic QCD at Finite Density, Phys.Rev. D83 (2011) 106003, [arXiv:1004.0162].
  * (15) O. Bergman, N. Jokela, G. Lifschytz, and M. Lippert, Striped instability of a holographic Fermi-like liquid, JHEP 1110 (2011) 034, [arXiv:1106.3883].
  * (16) C. B. Bayona, K. Peeters, and M. Zamaklar, A Non-homogeneous ground state of the low-temperature Sakai-Sugimoto model, JHEP 1106 (2011) 092, [arXiv:1104.2291].
  * (17) S. Takeuchi, Modulated Instability in Five-Dimensional U(1) Charged AdS Black Hole with R**2-term, JHEP 1201 (2012) 160, [arXiv:1108.2064].
  * (18) H. Ooguri and C.-S. Park, Holographic End-Point of Spatially Modulated Phase Transition, Phys.Rev. D82 (2010) 126001, [arXiv:1007.3737].
  * (19) A. Donos and J. P. Gauntlett, Helical superconducting black holes, Phys.Rev.Lett. 108 (2012) 211601, [arXiv:1203.0533].
  * (20) A. Donos and J. P. Gauntlett, Black holes dual to helical current phases, Phys.Rev. D86 (2012) 064010, [arXiv:1204.1734].
  * (21) M. Ammon, J. Erdmenger, S. Lin, S. Muller, A. O’Bannon, et al., On Stability and Transport of Cold Holographic Matter, JHEP 1109 (2011) 030, [arXiv:1108.1798].
  * (22) A. Donos, J. P. Gauntlett, and C. Pantelidou, Spatially modulated instabilities of magnetic black branes, JHEP 1201 (2012) 061, [arXiv:1109.0471].
  * (23) A. Donos and J. P. Gauntlett, Holographic striped phases, JHEP 1108 (2011) 140, [arXiv:1106.2004].
  * (24) S. Bolognesi and D. Tong, Monopoles and Holography, JHEP 1101 (2011) 153, [arXiv:1010.4178].
  * (25) P. Sutcliffe, Monopoles in AdS, JHEP 1108 (2011) 032, [arXiv:1104.1888].
  * (26) D. Allahbakhshi, On Holography of Julia-Zee Dyon, JHEP 1109 (2011) 085, [arXiv:1105.3677].
  * (27) K. Maeda, M. Natsuume, and T. Okamura, Vortex lattice for a holographic superconductor, Phys.Rev. D81 (2010) 026002, [arXiv:0910.4475].
  * (28) O. Domenech, M. Montull, A. Pomarol, A. Salvio, and P. J. Silva, Emergent Gauge Fields in Holographic Superconductors, JHEP 1008 (2010) 033, [arXiv:1005.1776].
  * (29) J. M. Murray and Z. Tesanovic, Isolated Vortex and Vortex Lattice in a Holographic p-wave Superconductor, Phys.Rev. D83 (2011) 126011, [arXiv:1103.3232].
  * (30) S. S. Gubser and S. S. Pufu, The Gravity dual of a p-wave superconductor, JHEP 0811 (2008) 033, [arXiv:0805.2960].
  * (31) M. Ammon, J. Erdmenger, M. Kaminski, and P. Kerner, Superconductivity from gauge/gravity duality with flavor, Phys.Lett. B680 (2009) 516–520, [arXiv:0810.2316].
  * (32) M. Ammon, J. Erdmenger, M. Kaminski, and P. Kerner, Flavor Superconductivity from Gauge/Gravity Duality, JHEP 0910 (2009) 067, [arXiv:0903.1864].
  * (33) S. Chunlen, K. Peeters, P. Vanichchapongjaroen, and M. Zamaklar, Instability of N=2 gauge theory in compact space with an isospin chemical potential, JHEP 1301 (2013) 035, [arXiv:1210.6188].
  * (34) D. Djukanovic, M. R. Schindler, J. Gegelia, and S. Scherer, Quantum electrodynamics for vector mesons, Phys.Rev.Lett. 95 (2005) 012001, [hep-ph/0505180].
  * (35) M. Chernodub, Spontaneous electromagnetic superconductivity of vacuum in strong magnetic field: evidence from the Nambu–Jona-Lasinio model, Phys.Rev.Lett. 106 (2011) 142003, [arXiv:1101.0117].
  * (36) V. Skokov, A. Y. Illarionov, and V. Toneev, Estimate of the magnetic field strength in heavy-ion collisions, Int.J.Mod.Phys. A24 (2009) 5925–5932, [arXiv:0907.1396].
  * (37) A. Bzdak and V. Skokov, Event-by-event fluctuations of magnetic and electric fields in heavy ion collisions, Phys.Lett. B710 (2012) 171–174, [arXiv:1111.1949].
  * (38) N. Callebaut, D. Dudal, and H. Verschelde, Holographic rho mesons in an external magnetic field, arXiv:1105.2217.
  * (39) E. Witten, Anti-de Sitter space and holography, Adv.Theor.Math.Phys. 2 (1998) 253–291, [hep-th/9802150].
  * (40) J. Erlich, E. Katz, D. T. Son, and M. A. Stephanov, QCD and a holographic model of hadrons, Phys.Rev.Lett. 95 (2005) 261602, [hep-ph/0501128].
  * (41) L. Da Rold and A. Pomarol, Chiral symmetry breaking from five dimensional spaces, Nucl.Phys. B721 (2005) 79–97, [hep-ph/0501218].
  * (42) A. Abrikosov, Fundamentals of the Theory of Metals.  North-Holland, Amsterdam, 1988.
  * (43) B. Rosenstein and D. Li, Ginzburg-landau theory of type ii superconductors in magnetic field, Rev. Mod. Phys. 82 (Jan, 2010) 109–168.
  * (44) W. H. Kleiner, L. M. Roth, and S. H. Autler, Bulk solution of ginzburg-landau equations for type ii superconductors: Upper critical field region, Phys. Rev. 133 (Mar, 1964) A1226–A1227.
  * (45) M. Tinkham, Introduction to Superconductivity.  Robert E. Krieger Publishing Company, Malabar, Florida, 1980.