# The NumericalCertification package in Macaulay2

Kisun Lee

###### Abstract

The package NumericalCertification implements methods for certifying numerical
approximations of solutions for a given system of polynomial equations. For
certifying regular solutions, the package implements Smale’s $\alpha$-theory
and Krawczyk method. For a singular solution, we implement soft verification
using the iterative deflation method. We demonstrate the functionalities of
the package focusing on interaction with current numerical solvers in
Macaulay2.

##  1 Introduction

Systems with polynomial equations arise in many fields in mathematics and
applied science. Specially, the interest on polynomial systems turns into a
problem in algebraic geometry that finding all isolated solutions for a given
polynomial system. Due to a recent development in numerical algebraic geometry
(e.g. see [SW05]), a family of numerical algorithms called the homotopy
continuation gains popularity as a way to find solutions for a polynomial
system. There are several known implementations Bertini [BHSW],
HomotopyContinuation.jl [BT18], Hom4PS-3 [CLL14], NumericalAlgebraicGeometry
[Ley11] and PHCpack [Ver99] which are all widely used.

One remark for the homotopy continuation algorithm is that its output is not
certified. It means that numerical approximations obtained by the algorithm
might not satisfy the users depending on their purposes. We say that a
numerical approximation is certified if a compact region that contains a
unique solution can be obtained from the given approximation by applying a
sort of algorithm. For a certified approximation, the unique solution
contained in the compact region is called an associated solution and the given
approximation is called an approximate solution. We call this series of
algorithms numerical certification.

As implementations for numerical certification for a polynomial system, we
point out alphaCertified [HS11] and a function certify in the software
HomotopyContinuation.jl [BRT20]. alphaCertified implements Smale’s
$\alpha$-theory [Sma86] as a way for numerical certification [HS12]. On the
other hand, HomotopyContinuation.jl exploits Krawczyk method [Kra69] using
interval arithmetic [Moo77] as a tool for certification.

The package NumericalCertification in Macaulay2[GS02] executes regular
solution certification using both $\alpha$-theory and Krawczyk method. A
preferred method can be chosen as an option by users. As an improved version
of the software presented in [Lee19], it includes soft verification for a
singular solution using the idea of the deflation method [LVZ06]. Finally, the
package provides an interface to the software alphaCertifed.

The rest of the paper consists of two sections. In the next section, we
discuss the required preliminaries for certification. The implementation
details are given in the last section.

##  2 Preliminaries

In this section, we review the concepts used for numerical certification.
Smale’s $\alpha$-theory and a combination of Krawczyk method and interval
arithmetic are used as methods for regular solution certification. We
introduce the deflation method for a notion providing an idea for singular
solution certification.

###  2.1 Smale’s $\alpha$-theory

Consider an $n\times n$-square system $F$, i.e. a system with $n$ polynomial
equations with $n$ variables. For a point $x\in\mathbb{C}^{n}$, recall the
Newton operator $N_{F}(x)$ defined like the following :

$N_{F}(x)=\left\\{\begin{array}[]{ll}x-F^{\prime}(x)^{-1}F(x)&\text{if
}F^{\prime}(x)\text{ is invertible},\\\ x&\text{otherwise}\end{array}\right.$

We say that a sequence $\\{N_{F}^{k}(x)\\}_{k=1}^{\infty}$ converges
quadratically to an associated solution $x^{\star}$ of $F$ if for every
$k\in\mathbb{Z}_{>}$,

$\left\|N_{F}^{k}(x)-x^{\star}\right\|\leq\left(\frac{1}{2}\right)^{2^{k}-1}\|x-x^{\star}\|.$

In this case, $x$ is an approximate solution for $F$. When $F^{\prime}(x)$ is
not invertible, we say $x$ is an approximate solution if and only if $F(x)=0$.
The $\alpha$-theory provides a certificate for the quadratic convergence of a
given point. The certificate is obtained from the three auxiliary parameters :

$\begin{array}[]{ccl}\alpha(F,x)&:=&\beta(F,x)\gamma(F,x)\\\
\beta(F,x)&:=&\|x-N_{F}(x)\|=\|F^{\prime}(x)^{-1}F(x)\|\\\
\gamma(F,x)&:=&\sup\limits_{k\geq
2}\left\|\frac{F^{\prime}(x)^{-1}F^{(k)}(x)}{k!}\right\|^{\frac{1}{k-1}}\end{array}$

where $F^{(k)}(x)$ in the definition of $\gamma(F,x)$ is a symmetric tensor
whose components are the $k$-th partial derivatives of $F$, see [Lan83,
Chapter 5]. The norm in $\beta(F,x)$ is the usual Euclidean norm and the norm
in $\gamma(F,x)$ is the operator norm on $S^{k}\mathbb{C}^{n}$ (for details,
see [HS11]). When $F^{\prime}$ is not invertible at $x$, we define
$\alpha(F,x)=\beta(F,x)=\gamma(F,x)=\infty$. The contents of $\alpha$-theory
are summarized below.

######  Theorem 2.1 (c.f. [BCSS12, HS12]).

Let $F:\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$ be a square polynomial system
with a point $x\in\mathbb{C}^{n}$. Then,

  1. 1.

if $\alpha(F,x)<\frac{13-3\sqrt{17}}{4}$, then $x$ is an approximate solution
for $F$, and $\|x-x^{\star}\|\geq 2\beta(F,x)$ where $x^{\star}$ is an
associated solution to $x$.

  2. 2.

if $\alpha(F,x)<0.03$ and $\|x-y\|<\frac{1}{20\gamma(F,x)}$ for a point $y$,
then $x$ and $y$ are both approximate solutions for $F$ to the same solution
$x^{\star}$, and $\|x-x^{\star}\|\leq\frac{1}{20\gamma(F,x)}$.

  3. 3.

if $\|x-\overline{x}\|>4\beta(F,x)$ for the conjugate $\overline{x}$ of $x$,
then $x^{\star}$ is not real.

For an implementation of $\alpha$-theory, the step for computing (or bounding)
$\gamma$ is required. For a degree $d$ polynomial $f=\sum_{|\nu|\leq
d}a_{\nu}x^{\nu}$, we recall that the Bombieri-Weyl norm is defined as

$\|f\|^{2}=\frac{1}{d!}\sum\limits_{|\nu|\leq d}\nu!(d-|\nu|)!|a_{\nu}|^{2}.$

For a system of polynomials $F=\\{f_{1},\dots,f_{n}\\}$, we define a norm for
the system

$\|F\|^{2}=\sum_{i=1}^{n}\|f_{i}\|^{2}.$

Let $d_{i}=\deg f_{i}$ for each $i=1,\dots,n$ and $d=\max d_{i}$. For a point
$x\in\mathbb{C}$, define $\|(1,x)\|^{2}=1+\sum_{i=1}^{n}|x_{i}|^{2}$, and we
let $\Delta_{F}(x)$ be the diagonal matrix with entries
$\Delta_{F}(x)_{ii}:=\sqrt{d_{i}}\|(1,x)\|^{d_{i}-1}$. Combining all these, a
bound for $\gamma(F,x)$ is given as follows:

######  Proposition 2.2.

[HS12, Proposition 5] Let $F$ be a square system of polynomials and
$x\in\mathbb{C}^{n}$ be a point. Suppose that $F^{\prime}(x)$ is nonsingular.
Define

$\mu(F,x):=\max\left\\{1,\|F\|\|F^{\prime}(x)^{-1}\Delta_{F}(x)\|\right\\}$
(1)

where the norm in $\|F^{\prime}(x)^{-1}\Delta_{F}(x)\|$ is the operator norm.
Then,

$\gamma(F,x)\leq\frac{\mu(F,x)d^{\frac{3}{2}}}{2\|(1,x)\|}.$ (2)

###  2.2 Interval arithmetic and Krawczyk method

Interval arithmetic is introducing arithmetic operators between intervals to
achieve conservative results on numerical computations. For an operator
$\odot$ with intervals $[a,b]$ and $[c,d]$ over $\mathbb{R}$, we define
$[a,b]\odot[c,d]=\\{x\odot y\mid x\in[a,b],y\in[c,d]\\}$. The real interval
arithmetic can be extended over $\mathbb{C}$. For describing an interval over
the complex, we use two intervals to construct an interval box
$[a_{1},b_{1}]+i[a_{2},b_{2}]$ containing numbers in $\mathbb{C}$. Then,
complex interval arithmetic can be done similarly using arithmetic over the
complex numbers. A set of complex intervals is denoted by $\mathbb{IC}$.
Likewise, a set of $n$-dimensional complex interval boxes is denoted by
$\mathbb{IC}^{n}$. For a function $F:\mathbb{C}^{n}\rightarrow\mathbb{C}^{m}$
and an interval box $I\in\mathbb{IC}^{n}$, we define an interval extension of
$F$ as a set containing the image of all $F$ on $I$, and it is denoted by
$\square F(I)$.

Krawczyk method is a combination of interval arithmetic and generalized
Newton’s method to get certificates for the existence and uniqueness of a
solution for a square system in a given interval. Suppose that a square
differentiable system $F:\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$ is given
with an interval $I\in\mathbb{IC}^{n}$. Let $Y$ be an $n\times n$-invertible
matrix and $x$ be a point in $I$. Then, we define Krawczyk operator centered
at $x$ like the following :

$K_{x,Y}(I):=x-YF(x)+(I_{n}-Y\square F^{\prime}(I))(I-x)$

where $I_{n}$ is the $n\times n$ identity matrix. Then, the following theorem
summarizes the propositions required for interval arithmetic-based
certification.

######  Theorem 2.3 (c.f. [Kra69, BLL19, BRT20]).

Suppose that $F:\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$ is a square
differentiable system with a given interval extension $\square F(I)$ on an
interval $I$. For an $n\times n$-invertible matrix $Y$ and a point $x$,

  1. 1.

if a root $x^{\star}$ of $F$ is in $I$, then $x^{\star}\in K_{x,Y}(I)$.

  2. 2.

if $K_{x,Y}\subset I$, then $I$ contains a root $x^{\star}$ of $F$.

  3. 3.

if $I$ contains a root of $F$ and $\sqrt{2}\|I_{n}-Y\square
F^{\prime}(I)\|<1$, then the root $x^{\star}$ in $I$ is unique.

  4. 4.

if $I$ contains a root of $F$, $\sqrt{2}\|I_{n}-Y\square F^{\prime}(I)\|<1$
and a set of conjugates $\\{\overline{y}\mid y\in K_{x,Y}(I)\\}$ for
$K_{x,Y}(I)$ is contained in $I$, then the root $x^{\star}$ in $I$ is unique
and real.

Here, $\|I_{n}-Y\square F^{\prime}(I)\|$ is the maximum operator norm of the
interval matrix $I_{n}-Y\square F^{\prime}(I)$ under the max-norm.

Note that the invertible matrix $Y$ is chosen for minimizing $\|I_{n}-Y\square
F^{\prime}(I)\|$. In an actual implementation, a natural choice for $Y$ can be
$F^{\prime}(m(I))^{-1}$ where $m(I)$ is the midpoint of the box $I$.

######  Remark 2.4.

In general, interval arithmetic certification allows working with less
precision than $\alpha$-theory. On the other hand, $\alpha$-theory shows a
better convergence rate to an actual solution. An example in [BLL19, Section
5.1] shows a comparison between two methods.

###  2.3 The deflation method

A deflation is a series of method to reinstate the quadratic convergence of
Newton iteration for an isolated singular solution of a system of equations.
The basic idea is introducing more equations to construct an augmented system
with reduced singularity (e.g. multiplicity). For an isolated singular
solution $x^{\star}$ of a square system $F=\\{f_{1},\dots,f_{n}\\}$, define
$\kappa=\dim\ker F^{\prime}(x^{\star})$. Then, for a randomly chosen vector
$B=\begin{bmatrix}b_{1},\dots,b_{n}\end{bmatrix}^{\top}\in\mathbb{C}^{n\times
1}$ from the kernel of $F^{\prime}(x^{\star})$, an augmented system

$G=\begin{bmatrix}F\\\ F^{\prime}B\end{bmatrix}$ (3)

has a solution $x^{\star}$ with a lower multiplicity than that of $F$. It is
known that a singular solution is regularized within finitely many iterations
by applying iterative deflation [LVZ06].

##  3 Implementation details

The package NumericalCertification is designed to interplay with other
numerical solvers in Macaulay2, for example, NumericalAlgebraicGeometry
[Ley11], Bertini [BGLR13] or PHCpack [GPV11]. Hence, the package supports
PolySystem and AbstractPoint types of input.

The most direct way to use the package is certifySolutions. It takes a
polynomial system and a list of numerical solutions as input.

      i1 : needsPackage "NumericalCertification"
      i2 : R = CC[x1,x2,y1,y2];
      i3 : F = polySystem {3*y1 + 2*y2 -1, 3*x1 + 2*x2 -3.5, x1^2 + y1^2 -1, x2^2 + y2^2 -1};
      i4 : sols = solveSystem f; -- a list of numerical solutions
      i5 : c = certifySolutions(F,sols);Ψ

It supports three strategies as options, alphaTheory, intervalArithmetic and
alphaCertified, and the default value is alphaTheory. The function returns
MutableHashTable and we can peek the output to see the certification results.
The option alphaTheory returns the values of $\alpha(F,x)$ for each numerical
solution (in an order of the input), the list of certified distinct, real,
regular and singular solutions. The list of non-certified solutions also
returned and it may be certified again after refinement.

      i6 : peek c
      o6 = MutableHashTable{alphaValues => {2.07811e-30, 1.97421e-40}
             certifiedDistinct => {{.652548, .771177, .757747, -.63662},
                  {.95437, .318445, -.298627, .947941}}
             certifiedReal => {{.652548, .771177, .757747, -.63662},
                  {.95437, .318445, -.298627, .947941}}
             certifiedRegular => {{.652548, .771177, .757747, -.63662},
                  {.95437, .318445, -.298627, .947941}}
             certifiedSingular => {}
             nonCertified => {}
             }

The option intervalArithmetic returns the list of Krawczyk operators for
certified real, regular and singular solutions. For the list of non-certified
solutions, it returns the input interval boxes used for certification.

      i7 : c = certifySolutions(F,sols, Strategy=>"intervalArithmetic");
      i8 : peek c
      o8 = MutableHashTable{certifiedReal => {|  [.95437,.95437] +
                  [-1.33962e-27,1.33962e-27]*ii ... }
             certifiedRegular => {|  [.95437,.95437] +
                  [-1.33962e-27,1.33962e-27]*ii ... }
             certifiedSingular => {}
             nonCertified => {}
             }

For the option alphaCertified, we need to set a path to the software
installed. It might require reloading the package. When we run
certifySolutions with the option, it runs the software and creates files in
the directory where alphaCertified is installed.

      i9 : loadPackage("NumericalCertification",
               Configuration=>{"ALPHACERTIFIEDexec"=>"path/to/alphaCertified/"},
               Reload=>true)Ψ
      i10 : certifySoltions(F,sols, Strategy=>"alphaCertified")

        alphaCertified v1.3.0 (October 16, 2013)
        Jonathan D. Hauenstein and Frank Sottile
               GMP v6.2.1 & MPFR v4.1.0

                          ...

What follows is the implementation details of each strategy used in the
package.

###  3.1 $\alpha$-theory certification

In this section, we look into functions executed in certification with
$\alpha$-theory. These functions can also be used separately.

The function computeConstants takes a polynomial system $F$ and a numerical
point $x$, and computes three parameter values $\alpha(F,x),\beta(F,x)$ and
$\gamma(F,x)$. For the value of $\gamma(F,x)$, the upper bound given in (2) is
used instead. We use the Frobenius norm to bound the matrix norm used in
$\mu(F,x)$ (1). The polynomial system and its Jacobian are evaluated by a
straight-line program (see [BCS13, Section 4.1]) implemented in the package
SLPexpressions [CDLS] for a faster evaluation. Note that the polynomial system
and the point must be in the same coefficient ring.

      i11 : x = point{{.652548, .771177, .757747, -.63662_CC}};
      i12 : computeConstants(F,x)
      o12 = (1.16708e-10, 5.22384e-13, 223.414)
      o12 : Sequence

The function certifyRegularSolution certifies the given solution by checking
the inequality $\alpha(F,x)<\frac{13-3\sqrt{17}}{4}$. It returns true if input
satisfies the inequality, false otherwise.

      i13 : certifyRegularSolution(F,x)
      o13 = true

The function certifyDistinctSolutions takes a polynomial system and two points
as input. It returns false if given two points converge to the same actual
solution of the system, otherwise true.

      i14 : y = point{{.95437, .318445, -.298627, .947941_CC}};
      i15 : certifyDistinctSolutions(F,x,y)
      o15 = true

A given solution may converge to a solution over the real numbers even though
it is a complex-valued solution. The function certifyRealSolution checks if a
given solution corresponds to a real solution or not.

      i16 : x = point{{.652548, .771177, .757747, -.63662+0.001*ii}};
      i17 : certifyRealSolution(F,x)
      o17 = true

The function alphaTheoryCertification takes a polynomial system and a list of
numerical solutions and runs all aforementioned functions at once according to
the algorithm established in [HS12, Section 2.2]. Unlike certifySolutions, it
does not execute singular solution certification.

      i18 : sols = {x,y};
      i19 : c = alphaTheoryCertification(F,sols);
      i20 : peek c
      o20 = MutableHashTable{alphaValues => {.000223414, 1.04693e-10}}
              certifiedDistinct => {x, y}
              certifiedReal => {x, y}
              certifiedRegular => {x, y}

Finally, the package supports the exact arithmetic over the rational numbers
or Gaussian rationals for $\alpha$-theory certification. For example,
constants $\alpha(F,x),\beta(F,x)$ and $\gamma(F,x)$ can be computed over the
Gaussian rationals as follows :

      i21 : CR = QQ[i]/ideal(i^2+1); -- a ring of Gaussian rationals
      i22 : R = CR[x1,x2,y1,y2];
      i23 : F = polySystem {3*y1 + 2*y2 -1, 3*x1 + 2*x2 -7/2, x1^2 + y1^2 -1, x2^2 + y2^2 -1};
      i24 : x = point(sub(matrix{{5/9,3/4,3/4,-1/2}},CR));
      i25 : computeConstants(F,x)
              73052652544805089   9731461   60054828392
      o25 =  (-----------------, ---------, -----------)
               8695980754208352  303595776   229146291
      o25 : Sequence

###  3.2 Interval arithmetic certification

For interval arithmetic certification, the package provides a type of
intervals over the complex numbers. The function intervalCCi returns a complex
interval from a pair of real intervals representing real and imaginary part
respectively.

      i26 : I1 = intervalCCi(interval(.8,.9),interval(-0.1,0.1))
      o26 = [.8,.9] + [-.1,.1]*ii
      o26 : CCi

When only one real interval is given as input, it returns a complex interval
with the zero interval for its imaginary part.

      i27 : I2 = intervalCCi(interval(.2,.3))
      o27 = [.2,.3] + [0,-0]*ii
      o27 : CCi

The package supports a basic interval arithmetic for the complex intervals and
matrices with complex interval entries.

      i28 : I1 + I2
      o28 = [1,1.2] + [-.1,.1]*ii
      o28 : CCi
      i29 : I1 * I2
      o29 = [.16,.27] + [-.03,.03]*ii
      o29 : CCi
      i30 : I1^3
      o30 = [.486,.756] + [-.244,.244]*ii
      o30 : CCi
      i31 : M = matrixCCi{{I1,I2},{I2,I1}}
      o31 = |  [.8,.9] + [-.1,.1]*ii [.2,.3] + [0,-0]*ii |
            |  [.2,.3] + [0,-0]*ii [.8,.9] + [-.1,.1]*ii |
      o31 : CCiMatrix
      i32 : M^2
      o32 = |  [.67,.91] + [-.18,.18]*ii [.32,.54] + [-.06,.06]*ii |
            |  [.32,.54] + [-.06,.06]*ii [.67,.91] + [-.18,.18]*ii |
      o32 : CCiMatrix

The function pointToInterval construct an interval box from a given point.
This function helps to make interval input without defining them separately.
There are two ways to use the function. The first is inputting a point and a
desired radius for an interval box. Then, it returns an interval box with the
given radius centered at the given point.

      i33 : x = point{{-1.6,-1.3*ii}};
      i34 : I = pointToInterval(x,1)
      o34 = |  [-2.6,-.6] + [-1,1]*ii [-1,1] + [-2.3,-.3]*ii |
      o34 : CCiMatrix

In many cases, a proper radius can be different depending on the polynomial
system or the accuracy of the approximation. Running the function
pointToInterval with a polynomial system and a point as input returns an
interval box with a radius estimating the distance between the point and the
convergence limit by using Newton-Kantorovich theorem.

      i35 : R = CC[x,y];
      i36 : F = polySystem {x^2 + y^2 -1, x - y^2};
      i37 : x = point{{-1.61803, -1.27202*ii}};
      i38 : I = pointToInterval(F,x)
      o38 = |  [-1.61803,-1.61803] + [-8.75505e-15,8.75291e-15]*ii
                  [-1.80849e-14,1.80894e-14] + [-1.27202,-1.27202]*ii |
      o38 : CCiMatrix

The function krawczykOperator computes Krawczyk operator from a given
polynomial system and an interval box or a point. When a point is given as
input, it computes Krawczyk operator from the interval obtained by
pointToInterval(F,x).

      i39 : krawczykOperator(F,I)
      o39 = |  [-1.61803,-1.61803] + [-2.2629e-25,2.2629e-25]*ii
                  [-1.65999e-25,1.65999e-25] + [-1.27202,-1.27202]*ii |
      o39 : CCiMatrix
      i40 : krawczykOperator(F,x)
      o40 = |  [-1.61803,-1.61803] + [-2.2629e-25,2.2629e-25]*ii
                  [-1.65999e-25,1.65999e-25] + [-1.27202,-1.27202]*ii |
      o40 : CCiMatrix

The function krawczykTest checks if Krawczyk operator satisfies 2 and 3 of
Theorem 2.3. It returns true if it is (hence the given interval is certified
to contain a solution uniquely), false otherwise.

      i41 : krawczykTest(F,I)
      o41 = true

Finally, the function krawczykRealnessTest certifies if a given interval
corresponds to a real solution to the system or not. It returns true if the
given interval contains a unique real solution to the system, false otherwise.

      i42 : krawczykRealnessTest(F,I)
      o42 = false
      i43 : y = point{{.618034, -.786151}}; -- a real solution
      i44 : krawczykRealnessTest(F,y)
      o44 = true

As shown in i44, both functions krawczykTest and krawczykRealnessTest also
take a point as input.

###  3.3 Singular solution certification

The method of iterated deflation is exploited for singular solution
certification. For a polynomial system and a numerical approximation, we
consider a subsystem

$\hat{F}=F+F^{\prime}B=\left\\{f_{i}+\sum_{j=1}^{n}\frac{\partial
f_{i}}{\partial x_{j}}b_{j}\right\\}_{i=1,\dots,n}$

of the overdetermined system $G$ given in (3) obtained by the deflation. The
$\alpha$-theory or interval arithmetic certification is applied on the square
subsystem and the given numerical solution. If the numerical solution is still
singular, then we construct an augmented system and take a square subsystem
repeatedly. Since applying the deflation on a singular solution must terminate
within finitely many iterations, the given numerical approximation becomes an
approximation of a regular solution of the square subsystem.

It is possible to produce a false positive result as $\hat{F}$ is obtained by
squaring-up the overdetermined system from a randomly chosen vector $B$.
However, it can recover the quadratic convergence of Newton iteration for a
singular solution with probability $1$, and so it can be used as soft
verification of a singular solution.

Singular solution certification is done by the function
certifySingularSolution to a given polynomial system and a numerical solution.

      i45 : F = polySystem {x^2 + y, x^3 - y^2};
      i46 : x = point{{1e-7,2e-7*ii}};
      i47 : certifySingularSolution(F,x)
      o47 = true

Both strategies alphaTheory and intervalArithmetic are available as options.
The function executes the iterated deflation until the given singular solution
is regularized. Therefore, the function might not terminate if a poor
approximation is given. To prevent this, if a user knows the number of
iterations required in advance, it can be provided as input.

      i48 : certifySingularSolution(F,x,1) -- an insufficient number of iterations
      o48 = false
      i49 : certifySingularSolution(F,x,2)
      o49 = true

## Acknowledgements

The author is grateful to Anton Leykin for encouraging the project and
suggesting various improvements. The author would also like to thank Michael
Burr and Thomas Yahl for help on the package, and the organizers of the
Macaulay2 workshop at Cleveland State University where the project was greatly
improved.

## References

  * [BCS13] Peter Bürgisser, Michael Clausen, and Mohammad A Shokrollahi.  Algebraic complexity theory, volume 315.  Springer Science & Business Media, 2013.
  * [BCSS12] Lenore Blum, Felipe Cucker, Michael Shub, and Steve Smale.  Complexity and real computation.  Springer Science & Business Media, 2012.
  * [BGLR13] Daniel J Bates, Elizabeth Gross, Anton Leykin, and Jose Israel Rodriguez.  Bertini for macaulay2.  arXiv preprint arXiv:1310.3297, 2013.
  * [BHSW] Daniel J. Bates, Jonathan D. Hauenstein, Andrew J. Sommese, and Charles W. Wampler.  Bertini: Software for numerical algebraic geometry.  Available at bertini.nd.edu with permanent doi: dx.doi.org/10.7274/R0H41PB5.
  * [BLL19] M. Burr, K. Lee, and A. Leykin.  Effective certification of approximate solutions to systems of equations involving analytic functions.  In Proceedings of the 2019 on International Symposium on Symbolic and Algebraic Computation, pages 267–274, 2019.
  * [BRT20] Paul Breiding, Kemal Rose, and Sascha Timme.  Certifying zeros of polynomial systems using interval arithmetic.  arXiv preprint arXiv:2011.05000, 2020.
  * [BT18] Paul Breiding and Sascha Timme.  HomotopyContinuation.jl: A package for homotopy continuation in Julia.  In International Congress on Mathematical Software, pages 458–465. Springer, 2018.
  * [CDLS] Justin Chen, Timothy Duff, Anton Leykin, and Mike Stillman.  SLPexpressions.m2: a Macaulay2 package for straight line programs and expressions for evaluation circuits.
  * [CLL14] Tianran Chen, Tsung-Lin Lee, and Tien-Yien Li.  Hom4PS-3: a parallel numerical solver for systems of polynomial equations based on polyhedral homotopy continuation methods.  In International Congress on Mathematical Software, pages 183–190. Springer, 2014.
  * [GPV11] Elizabeth Gross, Sonja Petrović, and Jan Verschelde.  PHCpack in Macaulay2.  arXiv preprint arXiv:1105.4881, 2011.
  * [GS02] D. R. Grayson and M. E. Stillman.  Macaulay2, a software system for research in algebraic geometry.  Available at http://www.math.uiuc.edu/Macaulay2/, 2002.
  * [HS11] Jonathan D Hauenstein and Frank Sottile.  alphaCertified: Software for certifying numerical solutions to polynomial equations.  Available at http://math.tamu.edu/~sottile/research/stories/alphaCertified, 2011.
  * [HS12] J. D. Hauenstein and F. Sottile.  Algorithm 921: alphaCertified: certifying solutions to polynomial systems.  ACM Transactions on Mathematical Software (TOMS), 38(4):28, 2012\.
  * [Kra69] Rudolf Krawczyk.  Newton-algorithmen zur bestimmung von nullstellen mit fehlerschranken.  Computing, 4(3):187–201, 1969.
  * [Lan83] Serge Lang.  Real analysis.  Addison-Wesley Publishing Company, Advanced Book Program, Reading, MA, second edition, 1983.
  * [Lee19] Kisun Lee.  Certifying approximate solutions to polynomial systems on Macaulay2.  ACM Communications in Computer Algebra, 53(2):45–48, 2019.
  * [Ley11] A. Leykin.  Numerical algebraic geometry.  Journal of Software for Algebra and Geometry, 3(1):5–10, 2011.
  * [LVZ06] A. Leykin, J. Verschelde, and A. Zhao.  Newton’s method with deflation for isolated singularities of polynomial systems.  Theoretical Computer Science, 359(1-3):111–122, 2006.
  * [Moo77] R. E. Moore.  A test for existence of solutions to nonlinear systems.  SIAM Journal on Numerical Analysis, 14(4):pp. 611–615, 1977.
  * [Sma86] S. Smale.  Newton’s method estimates from data at one point.  The Merging of Disciplines: New Directions in Pure, Applied, and Computational Mathematics, 1986.
  * [SW05] A. Sommese and C. Wampler.  The Numerical Solution of Systems of Polynomials Arising in Engineering and Science.  World Scientific, 2005.
  * [Ver99] Jan Verschelde.  Algorithm 795: PHCpack: A general-purpose solver for polynomial systems by homotopy continuation.  ACM Transactions on Mathematical Software (TOMS), 25(2):251–276, 1999.

Department of Mathematics, University of California San Diego,
---
9500 Gilman Drive, La Jolla, CA 92093
E-mail address<EMAIL_ADDRESS>
URL: https://klee669.github.io


# QCD Phase Diagram by Holography

Yi Yang<EMAIL_ADDRESS>Department of Electrophysics, National Chiao
Tung University, Hsinchu, ROC  Physics Division, National Center for
Theoretical Sciences, Hsinchu, ROC  Pei-Hung Yuan<EMAIL_ADDRESS>School of physics, University of Chinese Academy of Sciences, Beijing 100049,
China  Kavli Institute for Theoretical Sciences, University of Chinese Academy
of Sciences, Beijing 100049, China

###### Abstract

We explore QCD phase diagram by constructing a holographic QCD model using the
Einstein-Maxwell-Scalar system. The chiral transition is investigated by
adding a probe scalar and confinement transition is studied by adding a probe
string into the system. By interpreting the black hole phase transition in the
bulk spacetime as the quarkyonic transition in the dual QCD theory and
introducing the bypass mechanism for deconfinement transition, we give an
explanation why chiral symmetry breaking and deconfinement transition lines
coincide with each other despite their different physical origins.

Introduction It is widely believed that Quantum chromodynamics (QCD) is in the
confinement and chiral symmetry breaking ($\chi$SB) phase at low temperature
and small chemical potential region, while in the deconfinement and chiral
symmetry restored ($\chi$S) phase at high temperature and large chemical
potential territory. It is then conjectured that there is a phase transition
between these two phases Cabibbo . However, deconfinement and $\chi$SB have
quite different physical origins. In principle, the deconfinement transition
is defined in the quench limit $(m_{q}\rightarrow\infty)$ with the Polyakov
loop as its order parameter, while the $\chi$SB is defined in the chiral limit
$(m_{q}\rightarrow 0)$ with the quark condensation as its order parameter. In
practice, it is difficult to determine the QCD phase diagram due to the strong
interaction near the phase transition region where the conventional
perturbation method does not work.

For a few decades, lattice QCD simulation (LQCD) is the only reliable method
to attack this problem. Around zero chemical potential $\mu\sim 0$, LQCD
showed that the transition temperatures for the deconfinement and $\chi$SB are
almost the same 0609068 , and the phase transition reduces to a crossover for
the physical quark mass 1111.4953  that implies a critical end point (CEP) at
certain chemical potential $\mu$. However, LQCD suffers the sign problem for
finite chemical potentials.

(a) QCD phase diagram in Large $N_{c}$ limit

(b) Conjectured QCD phase diagram for finite $N_{c}$

Figure 1: The blue, purple and red lines label the deconfinement, $\chi$SB and
quarkyonic transition, respectively.

Moreover, large $N_{c}$ QCD provides an expressive picture for the QCD phase
diagram as shown in FIG.1(a). In the large $N_{c}$ limit, the transition
temperature of the deconfinement phase transition $T_{d}$ (blue line) is
independent of the chemical potential since quarks do not affect the gluons in
this limit. Above $T_{d}$, QCD is in the quark gluon plasma (QGP) phase with
deconfinement and $\chi$S. Below $T_{d}$, QCD is in the hadronic phase with
confinement and $\chi$SB at small $\mu$, and in the quarkyonic phase at large
$\mu$. There is a quarkyonic phase transition (purple line) between the two
phases below $T_{d}$ with the baryon density as its order parameter 0706.2191
. The hadronic, quarkyonic and QGP phases encounter at a tri-point
$(\mu_{t},T_{t})$. Furthermore, it has been argued that the $\chi$SB
transition (red line) is in the quarkyonic phase 0706.2191 .

However, for finite $N_{c}$, it was argued that the deconfinement transition
would reduce to a crossover 1801.09215  and decrease with the chemical
potential 0701022  implying the shrink of the quarkyonic region 0803.3318 ;
0805.1509 . In addition, it was shown that the quarkyonic transition is close
to the $\chi$SB with the CEP of $\chi$SB at $(\mu_{c},T_{c})$ being not far
from the tri-point 0911.4806 . LQCD also showed that the quarkyonic transition
is in the deconfined phase for small $\mu$ 0905.2949 . Moreover, at low $T$
and very large $\mu$, QCD is conjectured to be in the color superconductor
(CSC) or color-flavor locking (CFL) phase 111We will not consider CSC/CFL
phase in this work.. Combining the above evidences, the conjectured QCD phase
diagram for finite $N_{c}$ is plotted in FIG.1(b). Nevertheless, whether or
not the quarkyonic phase exists for finite $N_{c}$ is still an open question.
Many recent studies suggested that the deconfinement and $\chi$SB coincide for
$\mu\geq\mu_{t}$ implying the vanishing of the quarkyonic region 1104.0873 ;
1602.06699 ; 1705.00718 , while others claimed that there is no direct
relation between the deconfinement and $\chi$SB 1405.1289 ; 1502.07706 ;
1709.05981 .

On the other hand, holographic correspondence 9711200  offers an ideal
framework to explore QCD in the strongly coupled region by studying its dual
gravitational theory, i.e., holographic QCD (hQCD). Several hQCD models have
been constructed to study the phase transitions in QCD 1511.02721 ; 1512.06493
; 1610.09814 ; 1810.07019 ; 1908.02000 ; 1910.02383 . However, different
conclusions have been obtained from various models.

In this work, we use an extensively studied hQCD model based on the Einstein-
Maxwell-Scalar (EMS) system to investigate QCD phase structure. A probe scalar
and open strings are added into the system to realize the $\chi$SB and
deconfinement phase transitions, respectively. Nevertheless, the QCD dual of
the black hole phase transition in the bulk spacetime is still not clear. In
this work, we propose to interpret the black hole phase transition as the
quarkyonic transition in QCD and show that it coincides with the $\chi$SB.
Furthermore, we introduce the bypass mechanism in the holographic viewpoint
and show that the quarkyonic transition also coincides with the deconfinement
transition. We thus conclude that the dynamically stable quarkyonic state does
not exist and explain why the deconfinement and $\chi$SB transitions coincide
despite their different physical origins.

Holographic Model We follow the hQCD model constructed in 1301.0385 ;
1406.1865 ; 1506.05930 ; 1703.09184 ; 1705.07587 ; 2004.01965 ; 2009.05694  by
considering the 5-dimensional EMS system with probe matter fields
$S=S_{b}+S_{\chi}+S_{NG}$. In string frame, labeled by the super-index s,

$\displaystyle S_{b}\sim\int
d^{5}x\sqrt{-g^{s}}e^{-2\phi^{s}}\left[{R^{s}-\frac{f^{s}(\phi^{s}){F}^{2}}{4}}+4|\partial\phi^{s}|^{2}-V^{s}(\phi^{s})\right],$
(1) $\displaystyle S_{\chi}\sim\int
d^{5}x\sqrt{-g^{s}}e^{-\phi^{s}}Tr\\{\nabla_{M}X^{\dagger}\nabla^{M}X+m_{\chi}^{2}X^{\dagger}X\\},$
(2) $\displaystyle S_{NG}\sim\int d^{2}\xi\sqrt{-G},\text{
}G_{ab}=g_{\mu\nu}^{s}\partial_{a}X^{\mu}\partial_{b}X^{\nu}.$ (3)

where the background action $S_{b}$ includes a gravity field $g^{s}_{\mu\nu}$,
a Maxwell field ${F}_{\mu\nu}=\partial_{\mu}{A}_{\nu}-\partial_{\nu}{A}_{\mu}$
and a neutral scalar field $\phi^{s}$. The gauge kinetic function $f^{s}$ and
the scalar potential $V^{s}$ are functions of the scalar field $\phi^{s}$.
$S_{\chi}$ is for a composite massive scalar $X$ describing the quark
condensation $\langle\bar{\psi}\psi\rangle$ in the vacuum and $S_{NG}$ is the
Nambu-Goto action describing open strings which is the test objects for the
deconfinement transition. We treat all matter fields as probes by ignoring
their backreaction.

Background Solutions In practice, it is more convenience to study the
background in Einstein frame by a Weyl transformation:
$\phi^{s}=\sqrt{\frac{3}{8}}\phi$,
$g^{s}_{\mu\nu}=g_{\mu\nu}{e}^{\sqrt{\frac{2}{3}}\phi}$,
$f^{s}=f{e}^{\sqrt{\frac{2}{3}}\phi}$ and
$V^{s}={e}^{-\sqrt{\frac{2}{3}}\phi}V$. The background action Eq.(1) becomes

$S_{b}\sim\int
d^{5}x\sqrt{-g}\left[{R-\frac{f(\phi)}{4}{F}^{2}}-\frac{1}{2}|\partial\phi|^{2}-V(\phi)\right]$
(4)

To study QCD at finite temperature and chemical potential, we consider an
asymptotic AdS black hole ansatz,

$\displaystyle ds^{2}$
$\displaystyle=\frac{e^{2A\left(z\right)}}{z^{2}}\left[-g(z)dt^{2}+d\vec{x}^{2}+\frac{dz^{2}}{g(z)}\right],$
(5) $\displaystyle\phi$
$\displaystyle=\phi\left(z\right),~{}A_{\mu}=\left(A_{t}\left(z\right),\vec{0},0\right),$
(6)

where $z=0$ corresponds to the conformal boundary of the 5-dimensional bulk
spacetime. The equations of motion then can be derived from the action Eq.(4),

$\displaystyle
A_{t}^{\prime\prime}+\left(\frac{f^{\prime}}{f}+A^{\prime}-\dfrac{1}{z}\right)A_{t}^{\prime}$
$\displaystyle=0,$ (7) $\displaystyle\dfrac{\phi^{\prime
2}}{6}+A^{\prime\prime}-A^{\prime 2}+\dfrac{2}{z}A^{\prime}$
$\displaystyle=0,$ (8) $\displaystyle
g^{\prime\prime}+3g^{\prime}\left(A^{\prime}-\dfrac{1}{z}\right)-\frac{fz^{2}A_{t}^{\prime
2}}{e^{2A}}$ $\displaystyle=0.$ (9)

Imposing the asymptotic AdS condition $A(0)=\phi(0)=0$, $g(0)=1$ at $z=0$ and
the regular condition $A_{t}(z_{h})=g(z_{h})=0$ at $z=z_{h}$, the equations of
motion Eqs.(7-9) can be solved analytically,

$\displaystyle\phi\left(z\right)$
$\displaystyle=\int_{0}^{z}dy\sqrt{-6\left(A^{\prime\prime}-A^{\prime
2}+\dfrac{2}{y}A^{\prime}\right)},$ (10) $\displaystyle A_{t}\left(z\right)$
$\displaystyle=\mu\left[1-\frac{I_{2}(z)}{I_{2}(z_{h})}\right]=\mu+\rho
z^{2}+\cdots,$ (11) $\displaystyle g\left({z}\right)$
$\displaystyle=1-\frac{I_{1}(z)}{I_{1}(z_{h})}+\dfrac{\mu^{2}\left|\begin{array}[]{cc}I_{1}(z_{h})&\int_{0}^{z_{h}}I_{1}^{\prime}(y)I_{2}(y)dy\\\
I_{1}(z)&\int_{0}^{z}I_{1}^{\prime}(y)I_{2}(y)dy\end{array}\right|}{I_{1}(z_{h})I_{2}^{2}(z_{h})},$
(14) $\displaystyle I_{1}(z)$
$\displaystyle=\int_{0}^{z}\frac{y^{3}}{e^{3A}}dy\text{,
}I_{2}(z)=\int_{0}^{z}\frac{y}{fe^{A}}dy,$ (15)

where $\mu$ and $\rho$ represent chemical potential and baryon density by
holographic correspondence. Eqs.(10-15) are a family of solutions for
different choices of the functions $A$ and $f$. To realize the linear Regge
trajectories of vector meson spectrum and fit the chiral phase transition
temperature at zero chemical potential $T_{\chi}^{0}\simeq 0.159$ from the
recent lattice QCD simulations 2002.02821 , we fix
$f=e^{-cz^{2}-A-\sqrt{\frac{3}{8}}\phi}$ and $A=-a\ln(bz^{2}+1)$ with the
parameters $a=4.037$, $b=0.0155$ and $c=0.227$ in the following 1703.09184 ;
1812.09676 .

(a) $T-z_{h}$ (b) $F-T$

(c) $\rho-\mu$ (d) BH Phase diagram

Figure 2: (a) Temperature vs. black hole horizon at different chemical
potentials. (b) Free energy vs. temperature at different chemical potentials.
(c) Baryon density vs. chemical potential at different temperature. (d) Black
hole phase diagram in $T-\mu$ plane with the critical endpoint at
$(\mu_{c},T_{c})\simeq(0.216,0.159)~{}GeV$.

Black Hole Phase Transition The black hole temperature
$T=g^{\prime}\left(z_{h}\right)/4\pi$ is plotted in FIG.2(a). The multivalued
behavior at $\mu>\mu_{c}$ implies a phase transition. The transition
temperature can be obtained from the free energy $F=-\int sdT$ at each given
chemical potential with the entropy density $s=e^{3A(z_{h})}/4z_{h}^{3}$. The
free energy at different chemical potentials are plotted in FIG.2(b). At
$\mu>\mu_{c}$, the swallow tail shape implies a first-order phase transition
which reduces to a second-order one at the CEP $(\mu_{c},T_{c})$ and becomes a
crossover for $\mu<\mu_{c}$. The black hole phase diagram is plotted in
FIG.2(c).

How to interpret the black hole phase transition in the dual hQCD is still not
clear. It seems natural to follow the idea of interpreting the Hawking-Page
transition as the deconfinement transition in the dual field theory. However,
since both phases permit a black hole in the current case, this interpretation
cannot describe the confinement phase. We propose to interpret the black hole
phase transition in the bulk spacetime as the quarkyonic transition in the
dual hQCD. To justify our proposal, we calculate the baryon density
$\rho=\mu/2I_{2}(z_{h})$, which is the order parameter of the quarkyonic
transition, and plot it in FIG.2(d). The baryon density jumps from a tiny
value to a huge amount during the black hole phase transition implying that
the system transits to the quarkyonic phase.

(a)  |  (b) $\mu=0$ (c) $\mu=0.5$
---|---
Figure 3: (a) Quark condensation vs. temperature at different chemical
potentials. (b, c) The $\chi$SB transition temperatures at
$\mu=\\{0,0.5\\}~{}GeV$.

Chiral Symmetry Breaking We consider $N_{f}=2$ in this work. The probe scalar
field $X$ in the action Eq.(2) is a $2\times 2$ matrix and can be brought to
the diagonal form $X_{ii}=\chi/2$. By varying the action Eq.(2), we obtain the
equation of motion for the probe scalar field $\chi$,

$\chi^{\prime\prime}(z)+\frac{p(z)}{z}\chi^{\prime}(z)-\frac{q(z)}{z^{2}}\chi(z)=0,$
(16)

where $p=z\left(3A_{s}^{\prime}-\phi_{s}^{\prime}+g^{\prime}/g\right)-3$ and
$q=e^{2A_{s}}m_{\chi}^{2}/g$ are regular at $z=0$. And
$A_{s}=A+\frac{\phi}{\sqrt{6}}$. Near the boundary at $z=0$, Eq.(16) can be
solved by Frobenius method 2009.05694 ,

$\displaystyle\chi_{B}(z)$ $\displaystyle=\alpha\chi_{1}(z)+\beta\chi_{2}(z)$
$\displaystyle=\beta z+\beta C_{1}^{(2)}z^{2}+\alpha z^{3}+\beta Cz^{3}\ln z$
$\displaystyle+(\alpha C_{1}^{(1)}+\beta C_{3}^{(2)})z^{4}+\beta
CC_{1}^{(1)}z^{4}\ln z+\cdots$ (17)

where $\alpha$ and $\beta$ are coefficients of the linear combination and the
coefficients $C,C_{i}^{(1)}$ and $C_{i}^{(2)}$ can be obtained order by order
from Eq.(16). On the other hand, by holographic correspondence, the scalar
field $\chi$ can be expanded at the boundary as

$\chi_{B}(z)=m_{q}\zeta z+\Sigma z^{3}/\zeta+\cdots$ (18)

where $m_{q}\simeq 3MeV$ is the current quark mass, $\zeta=N_{c}/2\pi
N_{f}^{1/2}$ with $N_{c}=3$ and $N_{f}=2$ is a constant, and
$\Sigma=\langle\bar{\psi}\psi\rangle$ is the quark condensate. Comparing
Eq.(17) and Eq.(18), we have $\Sigma=\alpha\zeta$.

Near the horizon at $z=z_{h}$, the scalar field can be expanded as,

$\chi_{H}(z)=\sum\limits_{n=0}^{\infty}D_{n}(z-z_{h})^{n},$ (19)

where the coefficients $D_{n}$ can be calculated order by order from Eq.(16).

To get an analytic solution in the whole regime $[0,z_{h}]$, we smoothly match
the two asymptotic solutions $\chi_{H}$ and $\chi_{B}$ at a medium point
$z_{\epsilon}\in\left(0,z_{h}\right)$,

$\chi_{H}(z_{\epsilon})=\chi_{B}(z_{\epsilon}),~{}\chi_{H}^{\prime}(z_{\epsilon})=\chi_{B}^{\prime}(z_{\epsilon}).$
(20)

For a giving $m_{q}$, the quark condensation can be solved from the above
matching equations,

$\left.\Sigma=-m_{q}\zeta^{2}\frac{(\chi_{2}/\chi_{H})^{\prime}}{(\chi_{1}/\chi_{H})^{\prime}}\right|_{z=z_{\epsilon}}.$

The quark condensation vs. temperature at different chemical potentials are
plotted in FIG.3(a). For the physical quark mass, quark condensation is not an
exact order parameter for $\chi$SB, so that $\Sigma$ smoothly approaches to a
small but nonzero value at high temperature. While at low temperature,
$\Sigma$ becomes large implying the chiral symmetry breaking. For
$\mu<\mu_{c}$, the quark condensation is a monotonous function indicating that
$\chi$SB is a crossover. The transition temperature can be determined from the
point with the maximum changing rate as shown in FIG.3(b). For $\mu>\mu_{c}$,
the quark condensation becomes a multi-valued function, and the phase
transition between black holes induces a transition on quark condensation as
shown in FIG.3(c). Since we have interpreted the black hole phase transition
as the quarkyonic transition in hQCD, this shows that the quarkyonic
transition coincides with $\chi$SB and the point $(\mu_{c},T_{c})$ should be
identified as the CEP of $\chi$SB.

Deconfinement Transition To study the deconfinement transition in hQCD, we add
probing open strings in the black hole background and consider the expectation
value of the Polyakov loop $\left\langle
P\left(\mathcal{C}\right)\right\rangle\sim e^{-V_{q\bar{q}}(r,T)/T}\simeq
e^{-S_{on}}$ as the ordered parameter 9803002 , where $V_{q\bar{q}}$ is the
quark potential and $S_{on}$ is the on-shell Nambu-Goto action on a world-
sheet bounded by a loop $\mathcal{C}$. Choosing the static gauge, $\xi^{0}=t,$
$\xi^{1}=x$ and the boundary conditions $z\left(\pm r/2\right)=0$,
$z(0)=z_{0}$ and $z^{\prime}(0)=0$ with $z_{0}$ the maximum depth that the
string can reach, the shape of an open string can be solved from the equation

$z^{\prime}=\sqrt{g\left(\frac{\sigma^{2}(z)}{\sigma^{2}(z_{0})}-1\right)},$
(21)

where $\sigma(z)=z^{-2}e^{2A_{s}(z)}\sqrt{g(z)}$ is the effective string
tension. The distance and the quark potential are calculated as 1703.09184 ,

$\displaystyle r$
$\displaystyle=2\int_{0}^{z_{0}}dz\left[g(z)\left(\frac{\sigma^{2}(z)}{\sigma^{2}(z_{0})}-1\right)\right]^{\frac{1}{2}},$
(22) $\displaystyle V_{q\bar{q}}$
$\displaystyle=2\int_{0}^{z_{0}}dz\frac{\sigma(z)}{\sqrt{g(z)}}\left[1-\frac{\sigma^{2}(z_{0})}{\sigma^{2}(z)}\right]^{-\frac{1}{2}},$
(23)

which are plotted in FIG.4.

(a)  (b)

Figure 4: (a) The distant between quark and antiquark $r$ vs. $z_{0}$. (b) The
$q\bar{q}$ potential at different temperatures.

At low temperature or small black hole (red line), $r\to\infty$ at
$z_{0}=z_{m}$ where the dynamical wall emerges. The quark potential is linear
at large $r$ and open strings are always connected to form bounded states, so
that QCD is in the confinement phase. While at high temperature or large black
hole (blue line), there exists a maximum distance $r=r_{M}$ beyond that open
strings break and the quark potential becomes constant, so that quarks and
antiquarks become free indicating that QCD is in the deconfinement phase. A
cartoon of the dynamical wall for small and large black holes is shown in
FIG.5. The dynamical wall is reminiscence of the IR cut-off in the hard wall
model, but here the wall emerges dynamically instead of imposing by hand.

Figure 5: The cartoon of the dynamical walls for small and large black holes.
The green circle represents the dynamical wall.

(a)  (b)

Figure 6: (a) The deconfinement (dashed blue line) and $\chi$SB (solid red
line) transitions. We enlarge the region near $\mu\sim 0$ in (b).

The dynamical wall is a crucial concept to understand the deconfinement
transition in the holographic framework. It disappears at the deconfinement
transition temperature. The phase diagram of the deconfinement transition
(dashed blue line) is plotted in FIG.6(a). Since the dynamical wall appears as
the black hole horizon decreasing and passing the critical horizon smoothly,
we deduce that the deconfinement transition is a crossover. In addition, We
notice that the deconfinement transition temperature depends on the chemical
potential, but the dependence is much weaker than the $\chi$SB. To comparing,
we also plot the $\chi$SB line (red curve) in FIG.6(a). The deconfinement and
$\chi$SB lines intersect at a tri-point $(\mu_{t},T_{t})$. The near $\mu=0$
region of the QCD phase diagram is enlarged and plotted in FIG.6(b).

QCD Phase Diagram The phase diagram FIG.6(b) is consistent with the
conjectured QCD phase diagram at finite $N_{c}$ as shown in FIG.1(b). The
deconfinement line is a crossover and slightly bends down with $\mu$. The
quarkyonic line coincides with the $\chi$SB line. There exists the quarkyonic
phase between the deconfinement and $\chi$SB lines. The QGP, hadronic and
quarkyonic phases encounter at a tri-point which is close to the chiral CEP.

(a) $\mu<\mu_{t}$

(b) $\mu>\mu_{t}$

Figure 7: The cartoon for jumping the dynamical wall and the corresponding
temperature at small and large chemical potentials.

(a)  (b)

Figure 8: (a) The conclusive QCD phase diagram. The blue and red lines
represent the deconfinement and $\chi$SB transitions, which intersect at the
reduced tri-point at $(\mu_{t},T_{t})\simeq(0.607,0.150)$. The transition
temperatures for deconfinement and $\chi$SB at $\mu=0$ are $T_{d}\simeq
0.1503$ and $T_{\chi}\simeq 0.1555$. The $\chi$SB CEP is at
$(\mu_{c},T_{c})\simeq(0.234,0.155)$. (b) The chemical freeze-out parameters
from the STAR and ALICE collaborations are plotted in black squares 1212.2431
, green diamonds 1403.4903 , purple circles 1701.07065  and the red square
1710.09425 . The chiral crossover temperature obtained up to the order of
$O\left(\mu^{4}/T^{4}\right)$ from the HotQCD 1812.08235 ; 2002.11957  and
Wuppertal-Budapest collaborations 2002.02821  are marked with the yellow and
blue-net zones.

Nevertheless, this is not the final result for the QCD phase diagram yet. In
fact, the black hole phase transition dramatically affects the deconfinement
transition in an novel way. To explain it, we plot the jumps of the black hole
horizon in FIG.7. For $\mu<\mu_{t}$, the jump of the horizon between
$z=z_{h_{s}}$ and $z=z_{h_{l}}$ (red lines or dots) during the black hole
phase transition does not affect the dynamical wall (blue line or dot) at
$z=z_{m}$. While for $\mu>\mu_{t}$, the jump of the horizon bypasses the
dynamical wall and forces the deconfinement transition to take place at the
same temperature of the black hole transition. Therefore, the deconfinement
line coincides with the quarkyonic phase line (thereby also coincides with the
$\chi$SB line) and the quarkyonic phase disappears. More exactly, the
deconfinement transition occurs in the dynamically unstable region which is
bypassed by the black hole phase transition, as shown in FIG.7(b).

The final QCD phase diagram is plotted in FIG.8(a). At $\mu=0$, both
deconfinement and $\chi$SB are crossovers with $T_{d}\lesssim T_{\chi}$. The
$\chi$SB becomes first order at the CEP, and the deconfinement becomes first
order at the reduced tri-point beyond that the deconfinement transition
coincides with $\chi$SB. The chemical freeze out (CFO) parameters extracted
from the experimental measurements in heavy-ion collisions are plotted in
FIG.8(b) to compare with the phase diagram obtained from our holographic QCD
model. The CFO is due to the rapid decrease of degree of freedom in the
system. At $T=0$, $\mu_{CFO}\simeq 0.938$ represents nuclear mater in the
ground state so that the CFO temperature should be less than the $\chi$SB
temperature for large $\mu$ 1105.3934 . While in the crossover region
$\mu<\mu_{c}$, it was showed that the CFO temperature is close to the $\chi$SB
temperature 1909.02991 ; 2002.11957 . As shown in FIG.8(b), the chemical
freeze-out parameters from the STAR and ALICE collaborations are consistent
with the phase diagram obtained from our model.

Conclusion In this work, we investigated QCD phase diagram in a well studied
hQCD model constructed in the 5-dimensional Einstein-Maxwell-scalar system. To
realize the quarkyonic phase in the holographic correspondence, we proposed to
interpret the black hole phase transition in the bulk spacetime as the
quarkyonic transition in the dual QCD theory and showed that it coincides with
the $\chi$SB. We justified our proposal by showing the significantly jump of
the baryon density, which is conjectured to be the order parameter for the
quarkyonic transition, at the black hole transition temperature. In addition,
we introduced the bypass mechanism to explain why the deconfinement line
coincides with the $\chi$SB line despite that their physical origins are quite
different. The final QCD phase diagram is shown in FIG.8(a). The chemical
freeze-out parameters from the STAR and ALICE collaborations are well
consistent with the phase diagram we obtained from our holographic QCD model
as shown in FIG.8(b).

## Acknowledgements

We would like to thank Song He, Mei Huang, Danning Li, Wen-Yu Wen for useful
discussions. This work of YY is supported by the Ministry of Science and
Technology (MOST 106-2112-M-009 -005 -MY3) and National Center for Theoretical
Science, Taiwan. The work of PHY was supported by the University of Chinese
Academy of Sciences.

## References

  * (1) Cabibbo N and Parisi G 1975 Phys. Lett. B 59 67.
  * (2) Y. Aoki, Z. Fodor, S.D. Katz, K.K. Szabo, ”The QCD transition temperature: results with physical masses in the continuum limit”, arXiv:hep-lat/0609068, Phys.Lett.B643:46-54,2006.
  * (3) Michael Fromm, Jens Langelage, Stefano Lottini, Owe Philipsen, ”The QCD deconfinement transition for heavy quarks and all baryon chemical potentials”, arXiv:1111.4953 [hep-lat], JHEP 01 (2012) 042.
  * (4) Larry McLerran, Robert D. Pisarski, ”Phases of Dense Quarks at Large $N_{c}$”, arXiv:0706.2191, Nucl.Phys.A796:83-100,2007.
  * (5) Zhibin Li, Kun Xu, Xinyang Wang, Mei Huang, ”The kurtosis of net baryon number fluctuations from a realistic Polyakov–Nambu–Jona-Lasinio model along the experimental freeze-out line”, arXiv:1801.09215 [hep-ph], Eur.Phys.J.C 79 (2019) 3, 245.
  * (6) Troels Harmark, Vasilis Niarchos, Niels A. Obers, ”Instabilities of Black Strings and Branes”, arXiv:hep-th/0701022, Class.Quant.Grav.24:R1-R90,2007.
  * (7) Kenji Fukushima, ”Phase diagrams in the three-flavor Nambu–Jona-Lasinio model with the Polyakov loop”, arXiv:0803.3318 [hep-ph], Phys.Rev.D77:114028,2008.
  * (8) H. Abuki, R. Anglani, R. Gatto, G. Nardulli, M. Ruggieri, ”Chiral crossover, deconfinement and quarkyonic matter within a Nambu-Jona Lasinio model with the Polyakov loop”, arXiv:0805.1509 [hep-ph], Phys.Rev.D78:034034,2008.
  * (9) A. Andronic, D. Blaschke, P. Braun-Munzinger, J. Cleymans, K. Fukushima, L.D. McLerran, H. Oeschler, R.D. Pisarski, K. Redlich, C. Sasaki, H. Satz, J. Stachel, ”Hadron Production in Ultra-relativistic Nuclear Collisions: Quarkyonic Matter and a Triple Point in the Phase Diagram of QCD”, arXiv:0911.4806 [hep-ph], Nucl.Phys.A837:65-86,2010.
  * (10) C.E. DeTar, U.M. Heller, ”QCD Thermodynamics from the Lattice”, arXiv:0905.2949 [hep-lat], Eur.Phys.J.A 41 (2009) 405-437.
  * (11) Fukun Xu, Hong Mao, Tamal K. Mukherjee, Mei Huang, ”Dressed Polyakov loop and flavor dependent phase transitions”, arXiv:1104.0873 [hep-ph], Phys.Rev.D 84 (2011) 074009.
  * (12) M. Shifman, ”Dynamically Emergent Flavor in a Confining Theory with Unbroken Chiral Symmetry”, arXiv:1602.06699 [hep-ph], Phys. Rev. D 93, 074035 (2016).
  * (13) Kenji Fukushima, Vladimir Skokov, ”Polyakov loop modeling for hot QCD”, arXiv:1705.00718 [hep-ph], Prog.Part.Nucl.Phys. 96 (2017) 154-199.
  * (14) Takahiro M. Doi, Hideo Suganuma, Takumi Iritani, ”Relation between Confinement and Chiral Symmetry Breaking in Temporally Odd-number Lattice QCD”, arXiv:1405.1289 [hep-lat], Phys. Rev. D 90, 094505 (2014).
  * (15) M. Pak, M. Schröck, ”Overlap Quark Propagator in Coulomb Gauge QCD and the Interrelation of Confinement and Chiral Symmetry Breaking”, arXiv:1502.07706 [hep-lat], Phys. Rev. D 91 (2015) 074515.
  * (16) Hideo Suganuma, Takahiro M. Doi, Krzysztof Redlich, Chihiro Sasaki, ”Relating Quark Confinement and Chiral Symmetry Breaking in QCD”, arXiv:1709.05981 [hep-lat], J.Phys. G44 (2017) no.12, 124001.
  * (17) Juan M. Maldacena, ”The Large N limit of superconformal field theories and supergravity”, Int.J.Theor.Phys. 38 (1999) 1113-1133, hep-th/9711200.
  * (18) Kaddour Chelabi, Zhen Fang, Mei Huang, Danning Li, Yue-Liang Wu, ”Realization of chiral symmetry breaking and restoration in holographic QCD”, arXiv:1511.02721 [hep-ph], Phys. Rev. D 93, 101901 (2016).
  * (19) , Kaddour Chelabi, Zhen Fang, Mei Huang, Danning Li, Yue-Liang Wu, ”Chiral Phase Transition in the Soft-Wall Model of AdS/QCD”, arXiv:1512.06493 [hep-ph], JHEP 04 (2016) 036.
  * (20) Danning Li, Mei Huang, ”Chiral phase transition of QCD with Nf=2+1 flavors from holography”, arXiv:1610.09814 [hep-ph], JHEP 02 (2017) 042.
  * (21) Jianwei Chen, Song He, Mei Huang, Danning Li, ”Critical exponents of finite temperature chiral phase transition in soft-wall AdS/QCD models”, arXiv:1810.07019 [hep-ph], JHEP 01 (2019) 165.
  * (22) Xun Chen, Danning Li, Defu Hou, Mei Huang, ”Quarkyonic phase from quenched dynamical holographic QCD model”, arXiv:1908.02000 [hep-ph], JHEP 03 (2020) 073.
  * (23) Yoon-Seok Choun, Sang-Jin Sin, ”Bridging the Chiral symmetry and Confinement with Singularity”, arXiv:1910.02383 [hep-th], Phys.Lett.B 805 (2020) 135433.
  * (24) Song He, Shang-Yu Wu, Yi Yang, Pei-Hung Yuan, ”Phase Structure in a Dynamical Soft-Wall Holographic QCD Model”, arXiv:1301.0385 [hep-th], JHEP 04 (2013) 093.
  * (25) Yi Yang, Pei-Hung Yuan, ”A Refined Holographic QCD Model and QCD Phase Structure”, arXiv:1406.1865 [hep-th], JHEP 11 (2014) 149.
  * (26) Yi Yang, Pei-Hung Yuan, ”Confinement-Deconfinment Phase Transition for Heavy Quarks”, arXiv:1506.05930 [hep-th], JHEP 12 (2015) 161.
  * (27) Meng-Wei Li, Yi Yang, Pei-Hung Yuan, ”Approaching Confinement Structure for Light Quarks in a Holographic Soft Wall QCD Model”, arXiv:1703.09184 [hep-th], Phys.Rev.D 96 (2017) 6, 066013 .
  * (28) Yi Yang, Pei-Hung Yuan, ”Universal Behaviors of Speed of Sound from Holography”, arXiv:1705.07587 [hep-th], Phys.Rev.D 97 (2018) 12, 126009.
  * (29) Song He, Yi Yang, Pei-Hung Yuan, ”Analytic Study of Magnetic Catalysis in Holographic QCD”, arXiv:2004.01965 [hep-th].
  * (30) Meng-Wei Li, Yi Yang, Pei-Hung Yuan, ”Analytic Study on Chiral Phase Transition in Holographic QCD”, arXiv:2009.05694 [hep-th].
  * (31) Szabolcs Borsanyi, Zoltan Fodor, Jana N. Guenther, Ruben Kara, Sandor D. Katz, Paolo Parotto, Attila Pasztor, Claudia Ratti, Kalman K. Szabo, ”The QCD crossover at finite chemical potential from lattice simulations”, arXiv:2002.02821 [hep-lat], Phys. Rev. Lett. 125, 052001 (2020).
  * (32) Meng-Wei Li, Yi Yang, Pei-Hung Yuan, ”Imprints of Early Universe on Gravitational Waves from First-Order Phase Transition in QCD”, arXiv:1812.09676 [hep-th].
  * (33) Juan M. Maldacena, ”Wilson loops in large N field theories”, arXiv:hep-th/9803002, Phys.Rev.Lett. 80 (1998) 4859-4862.
  * (34) Francesco Becattini, Marcus Bleicher, Thorsten Kollegger, Tim Schuster, Jan Steinheimer, Reinhard Stock, ”Hadron Formation in Relativistic Nuclear Collisions and the QCD Phase Diagram”, arXiv:1212.2431 [nucl-th], Phys. Rev. Lett. 111 (2013) 082302.
  * (35) Paolo Alba, Wanda Alberico, Rene Bellwied, Marcus Bluhm, Valentina Mantovani Sarti, Marlene Nahrgang, Claudia Ratti, ”Freeze-out conditions from net-proton and net-charge fluctuations at RHIC”, arXiv:1403.4903 [hep-ph], Phys. Lett. B 738 (2014) 305-310.
  * (36) STAR Collaboration, ”Bulk Properties of the Medium Produced in Relativistic Heavy-Ion Collisions from the Beam Energy Scan Program”, arXiv:1701.07065 [nucl-ex], Phys. Rev. C 96 (2017) 4, 044904.
  * (37) A. Andronic, P. Braun-Munzinger, K. Redlich, J. Stachel, ”Decoding the phase structure of QCD via particle production at high energy”, arXiv:1710.09425 [nucl-th], Nature 561 (2018) 7723, 321-330.
  * (38) A. Bazavov, H.-T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, N. Karthik, E. Laermann, Anirban Lahiri, R. Larsen, S.-T. Li, Swagato Mukherjee, H. Ohno, P. Petreczky, H. Sandmeyer, C. Schmidt, S. Sharma, P. Steinbrecher, ”Chiral crossover in QCD at zero and non-zero chemical potentials”, arXiv:1812.08235 [hep-lat], Phys.Lett.B 795 (2019) 15-21.
  * (39) Heng-Tong Ding, ”New developments in lattice QCD on equilibrium physics and phase diagram”, arXiv:2002.11957 [hep-lat], Quark Matter 2019.
  * (40) Sourendu Gupta, Xiaofeng Luo, Bedangadas Mohanty, Hans Georg Ritter, Nu Xu, ”Scale for the Phase Diagram of Quantum Chromodynamics”, arXiv:1105.3934 [hep-ph], Science 332 (2011) 1525-1528.
  * (41) Wei-jie Fu, Jan M. Pawlowski, Fabian Rennecke, ”The QCD phase structure at finite temperature and density”, arXiv:1909.02991 [hep-ph], Phys. Rev. D 101, 054032 (2020).