Neural ordinary differential equation and holographic quantum chromodynamics

The holographic principle [11–13], also known as AdS/CFT correspondence, is a profound relation between a d-dimensional quantum field theory (QFT) and a (d + 1)-dimensional gravity theory. It has been successfully applied to a large class of strongly coupled QFTs in high energy theory and condensed matter theory. Despite its success, a constructive way of finding the holographic gravity dual theory for a given QFT is lacking. If we have the experimental response data of a quantum system under external probing fields, can we model it holographically by a classical field theory in some curved geometry? The entanglement feature learning [6, 19] and the AdS/DL correspondence [4, 5] can answer that question in a concrete setup. Here we briefly review the setup of [4], for which we apply the neural ODE method in later sections.

We assume the d + 1-dimensional bulk spacetime coordinated by $(t,\eta, x_1, \ldots, x_{d-1})$ including the time dimension t, the space dimensions x_i and the holographic bulk dimension η. We assume the translation symmetry except for the η direction, and the spatial homogeneity in $(x_1,\ldots,x_{d-1})$, then in the gauge $g_{\eta\eta} = 1$, the holographic bulk spacetime can be described by the following metric (we will consider d = 4 specifically):

$$\begin{equation} \mathrm{d}s^{2} = -f(\eta)\mathrm{d}t^{2}+\mathrm{d}\eta^{2}+g(\eta)(\mathrm{d}x^{2}_{1}+\cdots+\mathrm{d}x^{2}_{d-1}) \, . \end{equation} \tag{ 1 }$$

The dual quantum field theory lives in a d-dimensional flat spacetime spanned by $(t, x_1, \ldots,x_{d-1})$ on the holographic boundary. We call η the radial coordinate and the others are angular directions. The spacetime volume factor is

$$\begin{equation} \sqrt{|g|} = \sqrt{-\det g} = \sqrt{f(\eta)g(\eta)^{d-1}} \, . \end{equation} \tag{ 2 }$$

A scalar field φ in this curved spacetime is described by the action

$$\begin{equation} \begin{split} S[\phi] = \dfrac{1}{2}\int \sqrt{|g|}\left( g^{\mu \nu}\partial_{\mu}\phi \partial_{\nu}\phi+m^{2}\phi^2+\frac{\lambda}{2}\phi^4\right) \, . \end{split} \end{equation} \tag{ 3 }$$

The saddle point equation (the classical equation of motion) δS/δφ = 0 reads

$$\begin{equation} -\dfrac{1}{\sqrt{|g|}}\partial_{\mu}\left(\sqrt{|g|}g^{\mu\nu}\partial_{\nu}\phi \right)+m^{2}\phi+\lambda \phi^3 = 0 \, . \end{equation} \tag{ 4 }$$

Since we are interested in homogeneous static condensate in the dual quantum field theory, we assume that φ is only a function of η. Then equation (4) becomes

$$\begin{equation} \begin{aligned} -\partial_{\eta}^{2}\phi-(\partial_{\eta}\ln \sqrt{|g|})\partial_{\eta}\phi+ m^{2}\phi+\lambda \phi^3 = 0, \end{aligned} \end{equation} \tag{ 5 }$$

or equivalently, we could write it as

$$\begin{equation} \begin{aligned} &\pi = \partial_{\eta}\phi \, ,\\ &\partial_{\eta}\pi+h(\eta)\pi -m^{2}\phi -\lambda \phi^3 = 0 \, , \end{aligned} \end{equation} \tag{ 6 }$$

where the metric function is (with d = 4) [20]

$$\begin{equation} h(\eta) \equiv \partial_{\eta}\ln\sqrt{f(\eta)g(\eta)^{d-1}} \, . \end{equation} \tag{ 7 }$$

The input data is the pair $\phi(\eta \sim \infty), \pi(\eta\sim \infty)$ near the AdS horizon. And the field will propagate following the classical equation of motion equation (6). On the other hand, there is black hole horizon at η ∼ 0. The on-shell static scalar field satisfies the black hole boundary condition:

$$\begin{equation} \left[\dfrac{2}{\eta}\pi -m^{2}\phi -\lambda \phi^3\right]_{\eta\sim 0} = 0 \, , \end{equation} \tag{ 8 }$$

or equivalently, we could require [21]

$$\begin{equation} \pi(\eta\sim 0) = 0 \end{equation} \tag{ 9 }$$

The mapping between the asymptotic value of the scalar field $\phi(\eta \sim\infty)$ and the data of the dual quantum field theory is given by the AdS/CFT dictionary with the asymptotically AdS spacetime with the AdS radius L [4],

$$\begin{equation} L^{3/2}\phi\sim \alpha e^{-\eta/L} + \beta e^{-3\eta/L}-\frac{\lambda \alpha^3}{2 L^2}\eta e^{-3\eta/L}, \end{equation} \tag{ 10 }$$

for an operator ${\cal O}$ whose dimension is three, corresponding to the bulk scalar field φ with the mass $m^2 = -3/L^2$. The coefficients are related to the condensate as

$$\begin{align} \alpha = \frac{\sqrt{N_c}}{2\pi}m_{\cal O} \, , \quad \beta = \frac{\pi}{\sqrt{N_c}}\langle {\cal O}\rangle L^3 \, . \end{align} \tag{ 11 }$$

Here, $m_{\cal O}$ is the source for the operator ${\cal O}$ of the quantum field theory, and N_c denotes the color number in QFT and hence we set N_c  = 3 as we focus on QCD later. Therefore, the data of one-point function of the quantum field theory $\{m_{\cal O},{\langle\cal O\rangle}\}$ is given, it is mapped to on-shell configuration of φ(η) and π(η) (by taking derivative on both sides of equation (10)) near the holographic boundary $\eta\sim\infty$.

The experimental data pairs ($\phi(\eta\sim \infty), \pi(\eta\sim \infty)$) can be viewed as the positive data. And they will satisfy the black hole boundary condition equation (9) after following the classical equation of motion equation (6). We could also view pairs of data $\phi(\eta\sim \infty), \pi(\eta\sim \infty)$ that does not lie on the experimental curve as negative data. We expect those negative data will not satisfy the black hole boundary condition. Therefore, this becomes a binary classification problem, with the propagation equation (6). Here, for a given data of the condensate, the parameters in the differential equation to be learned are: the continuous metric function h(η), the AdS radius L and interaction coupling λ are in general unknown.

We regard equation (6) as a neural network, and the network weights are the metric function and other parameters. For that purpose, the numerical method known as the neural ODE is a perfect framework to find the optimal estimation for those unknown parameters. In the following, we will briefly review the neural ODE method.

The neural ODE [14] is a novel framework of deep learning. Instead of mapping the input to the output by a set of discrete layers, the neural ODE evolves the input to the output by a differential equation, which is trainable. The general form of the differential equation reads

$$\begin{equation} \dfrac{\mathrm{d}z(t)}{\mathrm{d}t} = f_{\theta}(z(t), t), \end{equation} \tag{ 12 }$$

where the vector z denotes the collection of hidden variables and θ denotes all the trainable parameters (which could also be t-dependent) in the neural network. Without loss of generality, suppose we have observations at the beginning and end of the trajectory: $\{(z_{0}, t_{0}), (z_{1}, t_{1})\}$. One starts the evolution of the system from $(z_{0}, t_{0})$ for time $t_{1}-t_{0}$ with parameterized velocity function $f_{\theta}(z(t), t)$ using any ODE solver. Then the system will end up at a new state $(z_{1}, t_{1})$. Formally, we could consider optimizing the general loss function $\mathcal{L}$, which explicitly depends on the output z₁ as

$$\begin{equation} \mathcal{L}(z_{1}) = \mathcal{L}\left(\int_{t_{0}}^{t_{1}}\mathrm{d}t~f_{\theta}(z(t),t)\right). \end{equation} \tag{ 13 }$$

To back-propagate the gradient with respect to the parameters θ, one introduces the adjoint parameters $a(t) = \frac{\partial \mathcal{L}}{\partial {z}(t)}$ and their corresponding backward dynamics:

$$\begin{equation} \dfrac{\mathrm{d}{a}(t)}{\mathrm{d}t} = -{a}(t)\cdot\dfrac{\partial {f}_\theta}{\partial {z}}. \end{equation} \tag{ 14 }$$

After solving equations (12) and (14) jointly, the parameter gradient can be evaluated from

$$\begin{equation} \dfrac{\partial \mathcal{L}}{\partial \theta} = \int^{t_{0}}_{t_{1}}a(t)\cdot\dfrac{\partial f_\theta}{\partial \theta}\mathrm{d}t \, . \end{equation} \tag{ 15 }$$

The derivation equation (15) can be found in the appendix.

3.1.1. Neural ODE and bulk equation

In the form of the first order differential equation, the equations of motion for the bulk field equation (6) can be translated to the neural ODE equation (12) by the following identifications:

$$\begin{align} (\pi, \phi) \leftrightarrow {z} \, , \quad \eta \leftrightarrow t. \end{align} \tag{ 16 }$$

The bulk metric function h(η) corresponds to the neural network weights θ. To make the network depth finite, we introduce the UV and IR cutoffs for the metric as $\eta_\mathrm{ini} = 1$, and $\eta_\mathrm{fin} = 0.1$ in units of the AdS radius L.

There are two big advantages of using neural ODE. First, the metric function is smooth and we do not need to add penalty terms for smoothness. Therefore, we can largely reduce the number of hyper-parameters needed in the network. Second, our neural ODE uses an adaptive ode solver, called 'dopri5.' This gives us much more accuracy in the integration, and it turns out that the equation of motion in the curved geometry is sensitive to the discretization in some region of η. This adaptive method provides accuracy and efficiency simultaneously.

3.1.2. Bulk metric parameterization

To make the integration variable monotonically increase from the AdS boundary to the black hole horizon, we made a change of variable $\widetilde{\,\eta} = 1-\eta$ for the metric function, and we model the metric function $h(\widetilde{\,\eta}$ using the following two ansatz:

$$\begin{align} \textrm{ansatz 1: } & h(\widetilde{\,\eta}) = \sum_{n = 0}^{8}a_{n}\widetilde{\,\eta}^{n} \, , \end{align} \tag{ 17 }$$

$$\begin{align} \textrm{ansatz 2: }& h(\widetilde{\,\eta}) = \sum_{n = 0}^{8}b_{n}\widetilde{\,\eta}^{n}+\dfrac{1}{1-\widetilde{\,\eta}}. \end{align} \tag{ 18 }$$

The first one is the Taylor series around the AdS boundary. The second choice explicitly encodes the divergent behavior of the metric function near the black hole horizon at η = 0. Any black hole horizon with a nonzero temperature has $f(\eta) \propto \eta^2$ with g(η) being nonzero constant. Hence, equation (7) leads to h(η)∼ 1/η as the generic behavior of h(η) near the horizon η = 0. The second ansatz equation (18) explicitly encodes this prior knowledge.

3.1.3. Lattice QCD data as input

We use the lattice QCD data of RBC-Bielefeld collaboration [23] as our input data. The data is the chiral condensate $\mathcal{O} = \bar{q}q$, as a function of its source, the quark mass m_q . A plot is given in figure 1(left). We take the T = 0.208 (GeV) temperature data (the black line in figure 1(left), and the detail of the data is listed in table 1 [24].

Zoom In Zoom Out Reset image size

Table 1. Chiral condensate as a function of quark mass [23], at the temperature T = 0.208 (GeV), converted to physical units [4].

mq (GeV)	$\langle\bar{\phi}\phi\rangle\,[(\mathrm{GeV})^{3}]$
0.00067	0.0063
0.0013	0.012
0.0027	0.021
0.0054	0.038
0.011	0.068
0.022	0.10

We generate positive data and negative data in such a way that if the data's vertical distance to the experimental curve is less than 0.005, then it is labeled as positive (the label is 0). Otherwise, it is labeled as negative (the label is 1). We collected 10 000 positive data and 10 000 negative data used for training, as shown in figure 1(right). Our goal is to obtain a holographic description of our QCD data using the neural ODE method. The variation parameters are λ, L and h(η).

3.1.4. Loss function

As for the loss function $\mathcal{L}$, we use

$$\begin{equation} \begin{split} \mathcal{L} = \frac{1}{N_{\textrm{data}}} \sum_\mathrm{data} & \left[ \bigm| T(\pi(\eta_{\textrm{fin}});\epsilon,\sigma)-l\bigm|^2 \right. \left.+\,\,\,\,\beta\left(h(\eta_{\textrm{int}})-4\right)^{2} \right], \end{split} \end{equation} \tag{ 19 }$$

where the first term is the mean square error of the classifier loss function for the output data to approach the true result, equation (9). The function T(x; , σ) is a specific differentiable non-linear activation function that maps region [−, ] to 0, and otherwise to 1, in a fuzzy manner:

$$\begin{equation} T(x;\epsilon,\sigma) = 1+0.5\left(\tanh\left(\dfrac{x-\epsilon}{\sigma}\right)-\tanh\left(\dfrac{x+\epsilon}{\sigma}\right)\right). \end{equation} \tag{ 20 }$$

The parameter σ controls the slope of the boundary as shown in figure 2. In the mean square error, l is the label of the data (l = 0 for positive data and l = 1 for negative data). The second term in equation (19), the β penalty term, is to impose the condition that the emergent metric needs to be asymptotically AdS near the boundary $\eta = \eta_\mathrm{ini}$. Due to non-linear nature of the ODE function and sensitivity of neural ODE, one may need to modify the hyperparameters (, σ) to ensure nonzero value of the gradient during the training.

Zoom In Zoom Out Reset image size

With the architecture described above, we perform the training. We first choose equation (17) for the ansatz of the metric function h(η). We randomly initialize the training parameters. The initial configuration of the metric function is given in the subplot (c) of figure 3. As shown in the subplot (a) of figure 3, the machine with the initial metric judges all the orange + green data as positive data.

Zoom In Zoom Out Reset image size

After training with 13 000 epochs, the loss is reduced to 0.02. The result is shown in subplot (b) and (d) of figure 3. As we can see the predicted data agrees well with original positive data. We also observe that the emergent metric is a smooth function. The trained metric function reads

$$\begin{eqnarray} \begin{aligned} h(\eta) = & \, 8.2352\widetilde{\,\eta}^{8} +8.0109\widetilde{\,\eta}^{7}+ 7.6072\widetilde{\,\eta}^{6} \\ & +6.9469\widetilde{\,\eta}^{5} + 150.89\widetilde{\,\eta}^{4} -130.81\widetilde{\,\eta}^{3} \\ & + 55.539\widetilde{\,\eta}^{2}-22.223\widetilde{\,\eta}^{1}+ 3.7720. \end{aligned} \end{eqnarray} \tag{ 21 }$$

The machine also finds the optimal values of the coupling constant and the AdS radius:

$$\begin{align} & \lambda = 0.0004, \end{align} \tag{ 22 }$$

$$\begin{align} & L = 5.1640\ [\mathrm{GeV}^{-1}]. \end{align} \tag{ 23 }$$

As we can see in subplot (d) of figure 3, the metric function h(η) which the neural ODE found has tendency to grow significantly near η ∼ 0. This is indeed the black hole horizon behavior. It is quite intriguing that the machine automatically captures the divergence behavior of the metric function h(η) near the black hole horizon.

As a check, we also perform the training with the second ansatz for the metric function h(η), i.e. equation (18), which encodes the prior knowledge about the black hole horizon. As shown in figure 4, the result looks almost the same as that of the first ansatz that does not use the prior knowledge. Therefore, the regularization to implement the black hole horizon in h(η) is not necessary. This result indicates that neural ODE can automatically discover the black hole geometry in the holographic bulk and recover the near-horizon metric behavior without prior knowledge. For convenience, we use the training results of the second ansatz to calculate a physical observable (Wilson loop) in the next section.

Zoom In Zoom Out Reset image size

We also applied the above method to the multi-temperature QCD data given in table 2. During the training, we require different neural networks to share the same value of AdS radius L, and the training results are summarized in table 3. The model discovers the optimal emergent metric as well as the coupling constant λ at each temperature.

Table 2. Chiral condensate as a function of quark mass [23], at the values of temperature T = 0.188, 0.192, 0.196, 0.200, 0.204, 0.208 (GeV), converted to physical units [4]. The quark mass m_q is in (GeV), and the chiral condensate $\langle\bar{\psi}\psi\rangle$ is in $[(\mathrm{GeV})^{3}]$.

T = 0.188		T = 0.192		T = 0.196		T = 0.200		T = 0.204		T = 0.208
mq	$\langle\bar{\psi}\psi\rangle$	mq	$\langle\bar{\psi}\psi\rangle$	mq	$\langle\bar{\psi}\psi\rangle$	mq	$\langle\bar{\psi}\psi\rangle$	mq	$\langle\bar{\psi}\phi\rangle$	mq	$\langle\bar{\psi}\psi\rangle$
0.00061	0.056	0.00062	0.049	0.00064	0.034	0.00065	0.019	0.00066	0.011	0.00068	0.0064
0.0012	0.058	0.0012	0.053	0.0013	0.042	0.0013	0.027	0.0013	0.018	0.0014	0.012
0.0024	0.064	0.0025	0.059	0.0025	0.052	0.0026	0.040	0.0026	0.029	0.0027	0.022
0.0049	0.07	0.005	0.068	0.0051	0.065	0.0052	0.058	0.0053	0.048	0.0054	0.038
0.0098	0.08	0.010	0.081	0.010	0.081	0.010	0.079	0.011	0.075	0.011	0.068
0.020	0.095	0.020	0.098	0.020	0.10	0.021	0.10	0.021	0.10	0.022	0.10

We have two observations of the trained results shown in table 3. First, the obtained metric h(η) and the AdS radius L do not depend on the temperature T. Second, the only dependence on the temperature is encoded solely in the coupling constant λ of the scalar field theory.

The former sounds counter-intuitive, since normally the metric itself should be highly dependent on the temperature, and the change in the metric will modify the gravitational fluctuation, which corresponds to the gluon physics. It is easy to resolve this issue. The obtained function is h(η) and not the full metric components f(η) and g(η). Even for the case of the AdS Schwarzschild geometry in which the metric is temperature-dependent, we find $h(\eta) = \frac{4}{L} \coth \frac{4\eta}{L}$ which is temperature independent. In the next section, to compute physical quantities from the emergent h(η), we assume some functional form of g(η) and discuss the temperature dependence of the metric components.

What the machine found is that the reproduction of the input data mainly relies on the temperature dependence of the coupling constant λ in the holographic bulk theory. For lower temperature, we find a strong non-linear interaction, i.e. larger λ. The value of λ is directly related to the self-coupling of sigma meson. Although we cannot compare our trained results with experiments since the self-coupling has never been precisely measured due to the broad width of the sigma meson, our result provides a unique view of the QCD phase transition, in particular about the mysterious relation between the chiral transition and the deconfinement transition. In addition, our analyses are only for the temperature values around the QCD critical temperature. The behavior of the chiral condensate at lower or zero temperature is also important in QCD, and the consistency of the bulk scalar model for QCD observables at zero temperature may improve the understanding of a gravity dual of QCD.

Since in our case the machine learns only h(η), to compute physical quantities such as Wilson loop, we need to assume the form of g(η) to get f(η). Here we assume the functional form of the AdS Schwarzschild configuration,

$$\begin{align} g(\eta) = A \left(\cosh\frac{2\eta}{La}\right)^a, \end{align} \tag{ 24 }$$

where A and a are temperature-dependent constant. In particular the constant a encodes the dimensionality of the AdS_d + 1-Schwarzschild as a = d/4, and here we just set it as a free parameter. The ansatz equation (24) also satisfies the criterion that g is a monotonic function of η, which is normally required for spacetimes without a bottle neck. The Hawking temperature T constrains the function f(η) as

$$\begin{align} f(\eta) \sim (2\pi T)^2 \eta^2, \end{align} \tag{ 25 }$$

so, for our calculation we define a new function F(η) as

$$\begin{align} f(\eta) = (2\pi TL)^2 \left(\tanh \eta/L\right)^2 F(\eta), \end{align} \tag{ 26 }$$

which satisfies the boundary condition

$$\begin{align} \lim_{\eta \to 0} F(\eta) = 1 \, . \end{align} \tag{ 27 }$$

Substituting equations (24) and (26) to equation (7), and perform the integration over η with the integration constant fixed by equation (27), we obtain

$$\begin{align} F(\eta) = \exp \int_0^\eta \left( 2 h(\eta) - \frac{4}{L\sinh (2\eta/L)}-\frac{6}{L} \tanh\frac{2\eta}{La} \right) d\eta \, . \end{align} \tag{ 28 }$$

The overall factor A in g(η) in equation (24) can be fixed by the following asymptotically AdS₅ constraint at $\eta \gg L$,

$$\begin{align} f(\eta) \simeq g(\eta) \simeq e^{2\eta/L + \textrm{const.}}, \end{align} \tag{ 29 }$$

which implies h(η)≃ 4/L according to equation (7). To determine this constant which we require temperature independent, we expand equation (28) around $\eta \gg L$ as

$$\begin{align} &\int_0^\eta \left( 2 h(\eta) - \frac{4}{L\sinh (2\eta/L)}-\frac{6}{L} \tanh\frac{2\eta}{La}\right) d\eta \nonumber \\ & \qquad= \frac{2\eta}{L} + c(a) + {\cal O}(1/\eta) \, . \end{align} \tag{ 30 }$$

Using this constant c(a), the constraint equation (29) determines the normalization of g(η) as

$$\begin{align} g(\eta) = (2\pi T L)^2 e^{c(a)} \left(2\cosh\frac{2\eta}{La}\right)^a. \end{align} \tag{ 31 }$$

Now, since we require that the constant in equation (29) is temperature independent, we have a condition:

$$\begin{align} \frac{\partial}{\partial T} \left[ T^2 e^{c(a(T))} \right] = 0 \, . \end{align} \tag{ 32 }$$

Up to an integration constant, we can numerically solve this equation. Assuming that at T = 0.208 (GeV) we have a = 1, we find numerically c(a = 1) = 11.1952. Then the equation above leads to $c(a(T = 0.188 \, \mathrm{(GeV)})) = 11.3984$ and $a(T = 0.208 \, \mathrm{(GeV)}) = 1.098$. We are going to use g(η) given by equation (31) and f(η) given by equation (26) with equation (28) for the calculation of physical quantities below.

Following the standard method [25–27] for calculating the expectation value of the Wilson loop holographically, we evaluate the Wilson loop for a quark and an antiquark separated by the distance d, using our emergent spacetime. The logarithm of the Wilson loop $\langle W \rangle$, which is proportional to the quark potential V, is the area of the Euclidean worldsheet of a string hanging down from the AdS boundary. The string reaches η = η₀ at the deepest, and both the quark potential V(d) and the quark distance d are functions of η₀, as

$$\begin{align} d &= 2\int_{\eta_0}^\infty \frac{1}{\sqrt{g(\eta)}}\sqrt{\frac{f(\eta_0)g(\eta_0)}{f(\eta)g(\eta)-f(\eta_0)g(\eta_0)}}d\eta \, , \end{align} \tag{ 33 }$$

$$\begin{align} 2\pi \alpha^{^{\prime}} V &= 2\int_{\eta_0}^\infty \sqrt{f(\eta)}\sqrt{\frac{f(\eta_0)g(\eta_0)}{f(\eta)g(\eta)-f(\eta_0)g(\eta_0)}}d\eta \,. \end{align} \tag{ 34 }$$

Here $1/(2\pi \alpha^{^{\prime}})$ is the string tension which is undetermined in this work. Eliminating η₀ from these expressions implicitly defines V(d). Note that the integration in V(d) diverges at $\eta = \infty$, and we need to introduce a cut-off for the asymptotic AdS boundary for the calculation.

The quark potential V(d) has another saddle, which is just two straight strings connecting the black hole horizon and the asymptotic boundary,

$$\begin{align} 2\pi \alpha^{^{\prime}} V_\mathrm{Debye} = 2\int^{\eta_0}_0 \sqrt{f(\eta)}d\eta\, . \end{align} \tag{ 35 }$$

We need to adopt V(d) in equation (34) or $V_\mathrm{Debye}$, whichever is smaller.

Using the metric obtained in the previous subsection, we calculate the quark potential for each temperature. In figure 5, we present the quark potential for T = 0.188 (GeV) data and T = 0.208 (GeV) data. They exhibit three phases: at short d, the potential is Coulombic, while at large d, the potential is flat and Debye-screened, and in the middle range of d, the potential is linear, signifying the quark confinement. The set of these features is well-known in lattice QCD simulations (see figure 6), and, interestingly, our holographic results reproduce these features [28].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This reproduction was reported in [4], and here we further investigate the temperature dependence. As we see in figure 5, the two plots are identical with each other except for the height of the Debye screening parts. The higher temperature corresponds to the lower height of the flat potential, which is qualitatively consistent with the lattice QCD result, as shown in figure 6.