Holographic Nuclear Physics with Massive Quarks

We discuss nuclear physics in the Witten-Sakai-Sugimoto model, in the limit of large number $N_c$ of colors and large 't Hooft coupling, with the addition of a finite mass for the quarks. Individual baryons are described by classical solitons whose size is much smaller than the typical distance in nuclear bound states, thus we can use the linear approximation to compute the interaction potential and provide a natural description for lightly bound states. We find the classical geometry of nuclear bound states for baryon numbers up to B=8. The effect of the finite pion mass - induced by the quark mass via the GMOR relation - is to decrease the binding energy of the nuclei with respect to the massless case. We discuss the finite density case with a particular choice of a cubic lattice, for which we find the critical chemical potential, at which the hadronic phase transition occurs.


Introduction
The holographic model of Witten-Sakai-Sugimoto (WSS) [1][2][3] is the top-down holographic theory closest to QCD to date. The model is based on a D4-D8 brane setup in type IIA string theory and the flavor dynamics is encoded in the low-energy action for the gauge fields on the D8 flavor branes in the geometry left by the gauge D4-branes. The model at low-energy reproduces a large-N c SU(N c ) gauge theory with N f massless quarks plus a tower of massive adjoint matter fields. Being a top-down model, it has very few parameters: N c , N f , the 't Hooft coupling λ and the mass scale M KK . In this paper, we will also include a quark mass term. The model shares all the important features with QCD, in particular confinement and chiral symmetry breaking as well as the existence of a low-energy chiral Lagrangian, which is of a Skyrme-type theory coupled to an infinite tower of vector mesons. Baryons in the WSS model are identified with instantons of the gauge theory describing the flavor branes [4][5][6][7], just like baryons can be seen as solitons of the Skyrme model [8,9]. Quantization of low-energy instantonic degrees of freedom provides quantum numbers for the corresponding nucleons. Many techniques developed in the context of quantization of zeromodes of Skyrmions as nuclei have been used in our approach, see e.g. Refs. [10][11][12][13][14][15][16][17].
Composite nuclei are described by multi-instantons configurations. We used this approach to describe composite nuclei from a "solitonic perspective" in a previous work [18], where we restricted to the case of massless quarks. This approach can be viewed as complementary to other approaches to holographic nuclear physics, see e.g. Refs. [19][20][21][22][23][24][25][26]. In the limit which we are considering, the instanton radius scales as λ − 1 2 and the distances between individual nuclei in the bound state configuration scale as λ 0 and so a linear approach can be used for the computation of the dominant two-body potential between the nuclei as an infinite sum of monopole and dipole interactions. In this limit the instantons become point-like but with an SU(2) orientation, which is very similar to the models considered in 3+1 dimensions in Refs. [27,28]. In this work we include the quark masses, which via the Gell-Man-Oakes-Renner (GMOR) relation induces a pion mass, which in turn is felt by the solitonic (pionic) degrees of freedom. In Ref. [18], we found bound states in the large N c and large λ limit and computed their respective nuclear binding energies. In distinction to the previous work (i.e. the massless case), the inclusion of a pion mass makes the nuclear bound states larger and the corresponding binding energies smaller.
The paper is organized as follows. In section 2 we give an overview of the model and discuss the introduction of the quarks mass term in the framework. In section 3 we compute the nucleon-nucleon potential. In section 4 we present the numerical results for the nuclear bound states. In section 5 we discuss the quantization for the deuteron state. In section 6 we discuss the case of finite density and the hadronic phase transition. We conclude in section 7 with a discussion and outlook.

The WSS holographic model with massive quarks
The model encodes color degrees of freedom in a background metric from type IIA string theory [1]: where φ is the dilaton, C n , F n+1 indicate the Ramond-Ramond n-form and its corresponding field strength, Vol 4 and 4 stand for the volume of a unit radius 4-sphere (S 4 ) and its volume form, respectively. The function f (u) makes the geometry terminate at a fixed coordinate u KK , so in order to avoid singularities, the coordinate x 4 has to be periodic with period δx 4 = 4π 3 Throughout this article we will work in units of the intrinsic mass scale of the theory, that is, we set: 3) The inclusion of flavor degrees of freedom is performed via insertion of a couple of N f stacked D8/D8-branes (with N f being the number of light quark flavors), transverse to the color branes in the x 4 direction (that is, localized on the circle) [2,3]. We will work in the setup with antipodal flavor branes, in which case the two stacks merge at the cigar tip labeled by u = u KK : this is in every sense a geometrical realization of the spontaneous breaking of chiral symmetry. A rigorous treatment should include also the backreaction on the geometry due to the presence of these stacks of branes, but since we will only be considering the presence of two light flavors, N f = 2, the effect of the modified geometry can be neglected as a first approximation (see Ref. [29] for a treatment of the backreaction of the flavor branes and the implications).
We will be interested in the theory on the D8/D8-branes, so we employ for the cigar subspace bulk coordinates, defined by with z running on the curve that defines the embedding of the flavor branes, and y being its transverse coordinate. The action on the D8-brane world volume is composed by two terms: a Yang-Mills action in warped spacetime arising from the truncation of the DBI action to quadratic terms, and a Chern-Simons term, originating from the coupling of the D8-branes to the Ramond-Ramond 3-form C 3 : with the warp factors k(z), h(z) and κ given by (2.6) The notation for the indices we use is as follows: In continuity with Ref. [18], it is useful to define a new coupling Λ, and a unique warp factor H(z) as: Moreover, we rescale the action by: so that the full rescaled action in the new notation reads: (2.10)

Quark mass
The addition of a quark mass to the model is performed by the insertion of a Wilson line operator in the dual QCD-like theory: this is the best we can do given the geometry, since flavor and anti-flavor degrees of freedom are not localizable at the same position (they will always remain separated along the x 4 direction). Hence, a quark mass term will be nonlocal and take the form: where OW j i is the open Wilson line operator: To obtain this object, we insert an open string stretching between the flavor branes, and provide the action with the Aharony-Kutasov term [30]: with S str being the action of the stretched string. The action S str is composed of two terms: the Nambu-Goto action for the free string (S NG ) and an interaction term of the string endpoints with the flavor branes to which they are attached: where the subtraction of the identity matrix accounts for the subtraction of the vacuum. After factorizing away the contribution from the Nambu-Goto term in S str , we are left with an action that includes the interaction of the endpoint of the string with the D8 world volume: the Nambu-Goto part will produce the quark mass and chiral condensate, as well as prefactors. We pack all this information into the constant c and in the quark mass matrix M . Finally, since the endpoints are forced to move on the y = 0 curve in the cigar subspace, the action in terms of the gauge fields is: At this point it is useful to evaluate c in terms of Λ: 16) and to perform the rescaling (2.9): (2.17)

Deformation of the instanton
The additional term in the Lagrangian induces two effects: it obviously changes the mass of the 1-instanton configuration, but could also potentially modify the size of the instanton. The classical YM instanton possesses a modulus ρ, that is interpreted as the size of the instanton and does not affect the energy. As explored in Ref. [5], the combined effect of the Chern-Simons term (that tends to dilate the instanton) and of the gravitational field (that tends to shrink it) is reflected by the fact that ρ ceases to be a modulus, and the energy gains a ρ dependence. The size of the instanton is then determined by minimizing the energy with respect to ρ. Here we will determine the contribution of the mass term to the total mass, and size of the instanton. We will follow the computations made in Ref. [31], adapting them to the gauge that was used in Ref. [7].
The energy of the static model is obtained by computing the opposite of the static action. This action is obtained by considering A I andÂ 0 as the only nonvanishing fields. The energy then reads (2.18) We will use the same field configuration as was used in Ref. [7]. In fact, as argued in Ref. [31], modifications to the field profile are subleading in Λ, and the first-order effects only affect the energy. 1 The field configuration we use is where ρ = √ x I x I and the profiles a and b are given by where r = √ x i x i . With this configuration, we can perform the integrals. The integral of the first two lines of Eq. (2.18) is done in great detail in Ref. [7]: resulting in: (2.21) For the third line of Eq. (2.18), using the arguments of Ref. [31] to get the leading order contribution in Λ, we can write the integral as To evaluate this integral in a divergenceless fashion, we must perform a gauge transformation to avoid singularities: we will thus work in the gauge (2.23) Performing the integration, the matrix exponentiation and the trace, we get tr M (U + U † − 21) = −4m 1 + cos π r 2 r 2 + ρ 2 , (2.24) where we assumed the quark masses to be degenerate, and thus M = m1, with m u = m d = m. This term cannot be integrated in closed form. The integral can be written as The total energy of the system is then As expected, the mass term works as an inertia factor, favoring configurations of smaller radius.
Minimization of Eq. (2.27) cannot be done analytically. In order to obtain some insight of the deformed quantities, we proceed by writing ρ = ρ 0 + x, where and making a power series expansion of the derivative of Eq. (2.27) around x = 0, truncating at first order and then finding the zero of the derivative (ramification point). Then, the result can also be expanded in a power series using m or 1 Nc . The deformation is (2.29) The mass of the system is equal to the energy at this value of ρ: we get (after rescaling the energy) We see that the corrections brought by the mass term can be written as a power series in m Nc (rescaled by an initial power of N c ), so the first correction given by the mass is subleading in N c .

Equations of motion
We will eventually be interested in the configuration with A M = 0, since those components arise as N c −1 corrections (induced via the Chern-Simons term by the presence of a slow motion, be it rotational or translational), so we neglect the factor exp − i dz A z . We also still impose the condition that the masses of the two light flavors are degenerate. The Aharony-Kutasov action contributes to the equations of motion with: Keeping only the first order of the sine power series, we obtain the new equation of motion for the component A z (i.e. it is a set of three equations, since A z = A a z T a ): while the other equations are unaltered (to be precise, the equation of motion for A z would receive a correction if we go beyond the static approximation, but that is beyond our goal to include time derivatives only). Since the holonomy of the A z field is dual to the pion field, and since the first-order term in the equations of motion arise from a Lagrangian term quadratic in it, the modification we performed to the equations of motion corresponds completely with the relaxation of the massless pion regime, in favor of a finite pion mass. This is fully consistent with the quark description: by giving finite mass to single quarks, pions are now pseudo-Goldstone bosons, hence massive. However, since the quark masses are degenerate, we still have a residual symmetry, reflected in the degenerate masses of the pion triplet. 2 For completeness and later convenience, we provide once again the full set of equations 2 It could look like the mass of the η meson is also bound to be the same as that of the pion triplet (in this limit with only two flavors), since the corresponding mass term has to originate from this very same Aharony-Kutasov action, and the fields are equally normalized. However, the Chern-Simons term for the Ramond-Ramond form, C 7 , induces the holographic realization of the Witten-Veneziano mechanism, thus removing this degeneracy. See Ref. [2], and Refs. [32,33] for more detailed discussions of motion restoring the parity indices and the sources: We define the fields in terms of Green's functions G(x, z, x , z ) and L(x, z, x , z ) as: and we take the functions G(x, z, x , z ) and L(x, z, x , z ) to obey with G(x, z, x , z ) and L(x, z, x , z ) given as in Refs. [6,18] by Our goal is to derive the values of k n dual to the meson masses again, to check how the presence of finite quark masses influences the masses of the bound states. To start off, we define a scalar product as in Ref. [18]: The Ansätze (2.42)-(2.43) and Eq. (2.39) give us the usual condition for the profile functions ψ n (z): so that ψ n (z) are found as solutions to an eigenvalue problem of an Hermitian (with respect to Eq. (2.44)) operator. They must then obey a completeness relation given by so that at this stage there are no differences compared to the massless quarks problem. The pion profile function corresponds to φ 0 so let us turn our attention to that: Eq. (2.41) is solved by imposing, for every n, the conditions Eq. (2.47) is solved by choosing φ n (z) = ∂ z ψ n (z) and d n = k 2 n c n : note however that d 0 is not determined by this condition, since there is no normalizable mode ψ 0 . Finally, Eq. (2.48) imposes the shape of the pion wave function φ 0 (z) = H − 3 2 (z). Now we substitute into Eq. (2.35) obtaining Here all contributions to the sum in the second term vanish due to the relation φ n (z) = ∂ z ψ n (z) and the normalizability of the eigenfunctions ψ n (z): so that only the n = 0 contribution survives. In the first term we can make use of a completeness relation for the functions φ n (z): We define a new inner product which satisfies φ n , φ n = d n = k 2 n c n for all n but n = 0. Extending the notation to φ 0 and performing the integration we identify d 0 = π and obtain the completeness relation Exploiting this new relation, we can check that the first term of Eq. (2.49) cancels out with the deltas on the right hand side. Hence all we are left with is the following equation: that we can further simplify by recalling that φ 0 (z) = H − 3 2 (z). After doing so, it is possible to identify the pion mass by recalling that the k 2 n are dual to the masses of the mesons, so we obtain An alternative, perhaps more intuitive way to obtain the same result is to obtain the effective Lagrangian for the mesons by expanding the gauge fields: doing so for the Aharony-Kutasov action, and remembering that the holonomy of the A z field is related to the pseudoscalars as we can easily derive the Gell-Man-Oakes-Renner relation for this model: it is easy to verify that m 2 π indeed coincides with k 2 0 of Eq. (2.54). Indeed, moving to higher orders in derivatives, further terms are induced in the action by the presence of the quark mass: in chiral perturbation theory a term can be written that has the form which is the leading correction to the tr (M U ) term, and generates mass terms for (axial) vector mesons. However, as argued in Ref. [31], this term is holographically realized as a correction of relative order Λ −2 , hence we neglect this contribution in the usual 1/Λ expansion.

Nucleon-Nucleon potential
We now turn to computing the interaction potential between two nucleons: this is done along the lines of the previous work (i.e. Ref. [18]). As a first step we build a two-instanton configuration by positioning two instantons at a distance far larger than their size and giving them arbitrary orientations: . In general, this field configuration is not a solution to the equations of motion, but it is an approximation at leading order in Λ in the limit R 2ρ with R being the distance between the instanton centers. This is a viable approximation scheme since the holographic model extrapolates from the large λ regime, and the size of the single instantons is of order ρ − 1 2 while R λ 0 . Space is then divided into three regions: two balls of radius ρ centered at X = X (p) and X = X (q) , where the corresponding field is strong and the other is in its linear approximation, and the rest of space where both fields can be approximated by their linear form.
To compute the (static) energy, we employ the definition S = − dtE so that is the rescaled energy, while E is the energy in physical units (L is the Lagrangian density in physical units). Part of the integral in Eq. (3.2) will account for the self energies of the solitons (the masses of the single baryons), so to compute the interaction potential, we have to subtract off these self energy terms: the Ansatz (3.1) allows us to easily do so by keeping just the cross-terms that involve both fields A p,q . The full rescaled energy is obtained by the on-shell action as:  [18], is the presence of the last term in Eq. (3.3), and that now k 0 = 0. We introduce the symmetric tensor: as well as the rotation matrix which gives the spatial rotation corresponding to an SU(2) rotation implemented via the matrix G ∈ SU(2). The full potential at the end of a somewhat lengthy computation (but analogous to that of section 3.1 of Ref. [18]) is given in physical units by: In this formula, r i is the relative position of the single soliton cores, while r 2 = r i r i . We see immediately that the difference with respect to the potential obtained in Ref. [18] is entirely contained in the last term, corresponding to the contribution of the n = 0 mode, that is, the pion: here we have simply another short-range interaction (actually absorbable in the previous term by extending the sum to n = 0 and remembering that d 0 = π), mediated by a particle of mass k 0 , as can be argued by the exponential decay. In Ref. [18] we had a long-range interaction, as appropriate for a massless mediating boson.
Here we assumed fixed a position in the z coordinate and fixed the instanton sizes: It is possible to consider these quantities as massive moduli as long as we consider small oscillations around their equilibrium value, so that they will actively enter in the expression of the potential to give the more general formula (3.7) As we can see now, every sum runs over all values of n, not alternating between even and odd values: It is easy to verify that this generalization is due to the presence of Z i = 0, and that in the case Z i = 0 the appropriate terms vanish because of the parity properties of the functions ψ n (z), and φ n (z).

Numerical evaluations 4.1 Fit of the free parameters
We now proceed to the numerical evaluation of the various physical quantities that we have found. In our model, we have three free parameters 3 , M KK , Λ and m, of which M KK is the only dimensional quantity. We calibrate the parameters to the following values: • For the mass parameter M KK , we choose the value 949 MeV. This is chosen to fit the mass of the ρ meson, that is the second-lightest massive meson in our model. Even if the usual practice is to fit the mass scale using the lightest massive particle in the model, we choose to use the ρ meson for the fit, following Ref. [18]. This way the introduction of the pion mass will appear numerically in the model as a perturbation of the massless model.
• For the coupling Λ, we choose the value 1.569. This is done to fix the pion decay constant in the massless model.
• The new parameter m is interpreted as the quark mass. We have chosen the value 3.066 × 10 −3 for this parameter. This yields a quark mass of 2.910 MeV, which is between the masses of the u and d quarks.
With the above parameter values, we get the following values for the observables: • Pion mass: Evaluating Eq. (2.54), we get a pion mass of 134.8 MeV. This is a very good result, as the experimental pion masses are 135.0 MeV for the neutral pion (π 0 ) and 139.6 MeV for the charged pions (π ± ). This result confirms that the model works well at low energies.
• Baryon mass: Evaluating Eq. (2.30), we get a baryon mass of 1.628 GeV for the light baryons (neutrons and protons). The baryons in our model are significantly heavier than the physical baryons, whose masses sit around 938 MeV for the proton and neutron. This is to be expected, as the model is expected to break down at energy scales larger around 949 MeV.

Composite nuclei
We will now turn to the minimization of the 2-body potential acting on B nuclei as where r ij ≡ |x i −x j | and B i is the SU(2) rotation matrix of the i-th nucleon (instanton). The numerical methods we will use are a simple random walk algorithm (only moving forward when the potential energy is lowered) and a Metropolis algorithm allowing for statistical random movement (with some probability of moving "uphill") as well as a simple gradient flow method.
In Fig. 1 we present the numerical results for classical composite nuclei, stable and meta-stable, up to baryon number B = 8. The stable configurations (ground states) are labeled with their corresponding baryon number B, while metastable states have a Latin letter added to the label. With increasing letter in the alphabet, the higher is the energy. This bond between two nuclei symbolizes the 2-body potential, which is used to find these configurations and the bond is displayed only for distances shorter than 1.5R 0 with R 0 being the optimal distance of two nucleons in deuterium (B = 2) for m π = 0.142. The color scheme utilized is a map from the unit 2-sphere to the Runge's color sphere, which is defined by white at the north pole, black at the south pole and the hue of the colors around the equator, going through red, yellow, green, cyan, blue, and magenta as one goes around the equator. The color scheme illustrates the orientation of the point-like instantons in SU(2) space, with the Runge's color sphere being the standard orientation corresponding to the unit matrix B = 1 2 . All points on the color sphere for a nonstandard orientation of the instanton are then rotated by the rotation matrix M ij (B) of Eq. (3.5).
We will now calculate the binding ratios of the B nuclear bound states, defined by The result is shown in Fig. 2. Qualitatively, the bound states show the same spatial/geome- trical distributions in the case of a finite pion mass as compared to the massless case, see Fig. 1 and Fig. 5 of Ref. [18]. The optimal distance between two nucleons in the 2-body potential grows with an increasing pion mass and this effectively enlarges all the sizes of the nuclear bound states, but does not affect the orientations of the instantons in SU(2) space. In particular, the optimal distance of the 2-body potential for m π = 0 in dimensionless units is R 0 = 2.06 and for m π = 0.142 it is R 0 = 2.23. The binding energies and thus the binding ratios, similarly decrease with an increasing pion mass.
On the other hand, the presence of the pion mass removes the long-range force of the potential, so even though the minimum is pushed a bit away in the 2-body potential, the attraction is also shortened somewhat. This has the consequence that some states are different in the massive case, with respect to the massless case. In particular, a metastable state in the B = 5 sector made of a tetrahedron with a satellite exists in the massless case [18], but has disappeared in the massive case. Similarly, the ground state of the B = 6 state has changed by moving the leftmost and the topmost nucleon in Fig. 1 closer  to the hexagon in the massive case, than in the massless case where it was more compact.
The details of the solutions in the massive case with m π = 0.142 are shown in Tab

The deuteron state
The quantum description of the deuteron is obtained by quantization of the collective modes starting from the B = 2 classical system that we examined in section 3. The procedure has already been described in detail in Ref. [18]: Here we will review the analysis and obtain the binding energy of the deuteron in the system at hand.
A generic field configuration in the two-instanton sector is written by considering only the degrees of freedom that do not change the potential computed in section 3, and fixing all other degrees of freedom to have the potential as attractive as possible. The unfixed degrees of freedom form the zeromode manifold, and can be interpreted in the following way: • Overall position of the center of mass of the system, that we denote as x = (x 1 , x 2 , x 3 ).
Due to the gravitational field in the z-direction, we do not allow motion along the z-axis. Quantizing those degrees of freedom gives the system an overall momentum.
• Isospin of the system: This is represented by an SU(2) matrix called U .
• Spin of the system: This is represented by an SU(2) matrix called E and the corresponding SO(3) matrix, M ij (E) is given in Eq. (3.5).
• Parity: We can always switch the two components of the system. As this is a discrete symmetry, it will not induce a momentum. We will indicate parity with a binary variable P = 0, 1.
The relative position and phase are minimized by assuming maximal attraction: This is attained at the separation R = (R 0 , 0, 0) (with both nuclei centered at the origin of the holographic coordinate) with phase opposition, B † C = iσ 3 . The numerical value of R 0 is the value at which the potential (3.7) is minimized.
A configuration A I in the two-instanton sector belongs to the zeromode manifold, if it can be written as The kinetic energy is computed by starting from the kinetic energy of two baryons, computed in Ref. [2], and then comparing configuration (5.1) with configuration (3.1) to understand how the collective coordinates are related to the coordinates describing the single baryons. After performing the transformation and freezing the coordinates that are not zeromodes, we obtain the kinetic energy for the zeromode manifold. Since the derivation is the same as in Ref. [18], we will just cite the result here. In terms of the angular left-invariant velocities Ω i = −i tr(U †U σ i ) and ω i = −i tr(E †Ė σ i ), we have the kinetic energy The mass M and radius ρ have been computed in subsection 2.2, and the parameter R 0 is obtained from the potential (3.7) in the attractive channel, by finding the minimum. The Hamiltonian for the deuteron is computed straightforwardly. First, the momenta are defined as In terms of those momenta, the Hamiltonian reads where we have included rest mass and defined the inertia tensors X, Y, Z as Quantum states are expressed as where k and l are eigenvalues for K 2 and L 2 , of values k 2 ( k 2 + 1) and l 2 ( l 2 + 1), respectively, k 3 and l 3 are eigenvalues for K 3 and L 3 , and i 3 and j 3 are eigenvalues for the right-invariant momenta I 3 and J 3 , that commute with the (left-invariant) momenta K 3 and L 3 and have the properties I 2 = K 2 and L 2 = J 2 . The right-invariant momenta can be obtained from the left-invariant momenta as J i = −M (E) ij L j and I i = −M (U ) ij K j (where E and U are promoted to position operators).
Due to the fact that the configuration space has discrete symmetries (as an example, multiplying U and E by iσ 1 on the right leaves the configuration invariant) we have to impose Finkelstein-Rubenstein constraints. The analysis in Ref. [18] is still valid: The only states that are compatible with the quantization of the single baryons as fermions and that include the constraints are |D = |0, 0, 0, 1, 0, j 3 , |I 0 = |1, 0, i 3 , 0, 0, 0 , Of those states, |D is the only one having the quantum numbers of the deuteron (isospin 0 and spin 1). The energies of the states are computed by acting on them with the Hamiltonian: the final result is It is evident that E D is always the smallest energy of the three ones above, so the deuteron state is effectively the ground state of the system. Although the result is identical to the result in Ref. [18] in form, we emphasize that in this setting the values of M B , ρ, V min and R 0 are different, so the final numerical result will be different.
With E D obtained from the first line of (5.10), the most relevant quantity we can compute is the binding energy E D − 2M B . Various parameters are needed to obtain this value. The new parameters R 0 and V min are not free, but they are fixed by looking for the minimum of the classical potential (3.7) in the attractive channel (B † C = σ 3 , r = (R 0 , 0, 0), ρ i = ρ, Z i = 0). In the massive model, we have found R 0 = 2.23 and V min = −0.220 by including 40 massive mesons in the potential computation. Furthermore, the instanton size can be computed from Eq. (2.29), obtaining a value of ρ = 2.307. We can combine these parameters to obtain the binding energy. The classical binding energy in MeV is given by −V min = 208.5 MeV, and the rotational corrections have a small influence on this value: the quantum binding energy is given by E D −2M B = 208.4 MeV. This is to be compared against a binding energy of 2.224 MeV found in experiments. As in our previous work, the binding energy is two orders of magnitude larger than the expected binding energy. The pion mass improves the prediction, as the binding energy of the massless model is 275.5 MeV. A study of the 1 Nc and 1 Λ corrections is required to understand if the prediction for the binding energy can be improved, as the parameters N c and Λ are not that large.

An infinite lattice of baryons
The attractive channel of the potential between two nucleons is obtained by employing a relative rotation between their SU(2) orientations: The iso-rotation must act by an angle π with an axis orthogonal to the relative position vector. Using the axis-angle notation, the matrix M ij (B † C) that describes the relative orientation is given by It is possible to arrange nucleons in a cubic lattice in such a way that the interaction between all nearest neighbors is in the attractive channel, see fig. 3. The Ansatz for the lattice we are using is equal to that of Ref. [34]. The relation between the instanton lattice studied here and the Skyrmion lattice of Ref. [34] is given by the holonomy of the instanton gauge field [35]. We make the following choices: The relative iso-orientation B † C between every soliton and its nearest neighbors is given by • ±iσ 2 for r = (±R, 0, 0), corresponding to α = π,û =x 2 , • ±iσ 3 for r = (0, ±R, 0), corresponding to α = π,û =x 3 , • ±iσ 1 for r = (0, 0, ±R), corresponding to α = π,û =x 1 .
This choice for the nearest neighbors fixes completely the lattice in a self-consistent way (since σ 1 σ 2 = iσ 3 ), so it is now possible to compute the energy density associated to it.
Before moving to the computation, let us make some useful considerations: The relative orientation enters the potential formula via the combination M ij P ij (r, k), which reads so that we have two classes of iso-orientations: α = 0 and α = π: To compute the full energy density it is sufficient to evaluate the potential between a single chosen nucleon and all the others: Then translational symmetry ensures that the full energy is given by multiplying this contribution by the number of nucleons (which is infinite) with a factor of one half to avoid double counting. The nucleon number can then be traded for the number of cells, as a function of the lattice volume V L , so the full energy will be computed as with p indicating a lattice site.
To make the calculation more straightforward, we choose our reference frame in such a way that r p = (0, 0, 0) and B p = 1, so we are left with the following formula for the energy density: where p runs over all lattice sites with the exception of (0, 0, 0). The lattice sites are uniquely parametrized as The position r not only enters the expression of M ij P ij , but it also determines whether α = 0 or α = π: We must therefore be careful and classify iso-orientations with respect to (n 1 , n 2 , n 3 ). There are four possible scenarios, depending on how many even or odd instances of n i are present in the position vector: • all n i are EVEN ⇒ C p = 1, • all n i are ODD ⇒ C p = 1, where in the last two lines σ j indicates a rotation around the axis determined by our rule for the rotation of nearest neighbors (that is, j = 2 for i = 1, j = 3 for i = 2 and j = 1 for i = 3). For example, a nucleon sitting at R(2b 1 − 1, 2b 2 , 2b 3 ) with b i ∈ Z will have an orientation of ±iσ 2 , i.e. the same as a nucleon sitting at R(2b 1 , 2b 2 − 1, 2b 3 − 1). The last two scenarios fall into the class of α = π, and the axis of rotation enters Eq. (6.2) only via (û · r), which selects the component of r along the axisû. But sinceû is a unit vector corresponding to one of the coordinate axes, this simply selects the j-th component Rn j .
Now we will compute the contributions to the energy density: Let us start with the first two possibilities in the classification, corresponding to all even or all odd n i , for which we obtain with r 2 = 4R 2 (b 2 1 +b 2 2 +b 2 3 ) for the even case (E 1 ) and for the odd case (E 2 ). We now turn to the more difficult scenarios: We start by considering r = (2b 1 − 1, 2b 2 , 2b 3 )R with b i ∈ Z. In this case the energy density becomes where we have defined e −k 0 r r 3 1 + rk 0 − bR 2 r 2 3 + 3k 0 r + k 2 0 r 2 , (6.9) with r 2 = R 2 [(2b 1 − 1) 2 + 4b 2 2 + 4b 2 3 ]. We can see that the other terms labeled by E (2,3) 3 will have the same form except for the substitution of every instance of b 2 with respectively b 3 , b 1 : Since the sum runs over all three parameters, then each of these terms will give contribution equal to the result We are thus only left with the terms of the case in which one coordinate is an even multiple of the lattice spacing, while the others are odd. As before we begin by analyzing the situation in which this coordinate is r 1 , so that r = R(2b 1 , 2b 2 − 1, 2b 3 − 1): The result is just as the previous one, with the exception that the projection on the j-th component of r due to the rotation axis will now pick up an odd value instead of an even one, so that and r is now given by . Again, as in the previous case, the sum of all three possible combinations with only one even coordinate will just result in a factor of three, so we get The total energy density is then given by the sum of the contributions from all four possibilities in the classification, giving  The total lattice potential can then be computed numerically and plotted as a function of the lattice spacing R: the sum over lattice sites converges with good precision (∼ 10 −5 ) for b j ∈ [−L max , L max ] with L max = 15, ∀ j = 1, 2, 3. Notice that the cutoff in terms of lattice sides is 2L max in each direction (i.e. both in the x+ direction and in the x− direction) and the factor of two is due to b j being an integer parametrizing even or odd n j 's. The resulting potential density has a shallow minimum at R = 4.86 (see Fig. 4).
It is not straightforwardly clear that the quadruple sum in Eqs. (6.7) and (6.9) are convergent in the limit of L max → ∞, i.e. in the limit of summing over all lattice sites b j and all massive vectors n (being the index of k n ). The convergence of the quadruple sum would prove that the density in the given form is finite. We give a mathematical proof of the convergence of each term in the sums separately in appendix A.

Hadronic phase transition
We will now use the baryon lattice to study the presence of a hadronic phase transition: To do so we need to compute the free energy of the configuration and look for a critical density at which it becomes negative.
We can do this without relying on holography, simply by calculating the Legendre transform of the energy density. First of all, we recover the total energy density from the interaction potential density by adding the mass density terms: the cubic lattice cell has unit net baryon number, so the mass density as a function of the lattice spacing R is simply with M B given by Eq. (2.30).
The chemical potential is defined as the derivative of the energy density with respect to the baryon number density (which we can trade for a derivative with respect to R, since the baryon number density is d B = R −3 ): with EV −1 L given by Eq. (6.13).
We want to compare the free energy density of the infinite lattice of baryons constructed above, with that of the vacuum: the embedding of the flavor branes is the same for both phases, so we neglect its contribution to both phases, effectively setting to zero the free energy of the vacuum. All we need to compute is then the change in free energy due to the presence of baryons: the lattice configuration will become favored over the vacuum (that is, a purely mesonic phase) for negative values of this free energy. The free energy density F is then obtained from the energy density as:  Figure 5: The free energy density as a function of the lattice spacing R. As the density increases (i.e. R decreases), the free energy becomes negative, signaling a phase transition. L max is the cutoff of lattice sites: |b j | ≤ L max .
In Fig. 5, we plot the free energy density as a function of R, showing that it becomes negative at a finite density, signaling the presence of a first-order hadronic phase transition. Note that with this analysis we cannot conclude that the favored of all phases is indeed given by this lattice of baryons, but only that it is favored over the vacuum. Different configurations of baryonic matter can potentially give rise to even lower free energies, and become favored over this particular lattice.

Conclusion
In this paper, we have included the quark mass term using the Aharony-Kutasov action in our solitonic approach to holographic nuclear physics. We work in the limit Λ → ∞, where the size of the instantons is much smaller than the typical separation distance between nuclei in a nuclear bound state. This allows us to calculate multi-instanton solutions by gluing together two overlapping instantons in the linear regime; in particular, we compute the 2-body potential which now has acquired a mass term for the pions. This induction of the pion mass from the quark mass term is in line with the GMOR relation. Using numerical methods we find nuclear bound states -stable and metastable -with baryon numbers B = 2 through B = 8.
We find that the main difference by having massive pions in the 2-body potential is that the nuclear bound states grow slightly in size and correspondingly reduce their binding energy. A similar effect applies also to the deuteron bound state and to the other light nuclei.
Using the 2-body potential, we consider the case of an infinite crystal of nucleons in a cubic lattice with each nucleon oriented in SU(2) space so as to put it in the attractive channel with respect to nearest neighbors. We calculate the potential energy density, prove in appendix A that it is finite, and finally use it to calculate the free energy. The free energy becomes negative at a critical lattice spacing R crit ∼ 4.9, which signals a hadronic phase transition of first order. The presence of the hadronic phase transition was also studied in the same model in Refs. [36][37][38][39]: The difference with our approach lies both in the description of the solitons and in the setup of the flavor branes. In the cited works the branes' position at infinity are taken to be close enough to allow the use of the deconfined geometry up to arbitrarily low temperatures, thus being far from the antipodal branes' regime: the baryonic matter is described with various levels of accuracy, starting from exactly pointlike instantons up to an infinite lattice of finite size instantons interacting with the nearest neighbors via the ADHM construction. We work instead with antipodal branes at zero temperature, using the confined geometry, and we employ the flat space instanton approximation, then take into account every single interaction (up to a cutoff in lattice size) with every soliton seeing each other's core as pointlike. The presence of the first order phase transition remains a feature of the model in both regimes, and is becoming of second order if interactions between instantons are ignored [36].
The work carried out in this paper has been done in the large Λ and large N c limit, which is of course none other than approximations to real world physics, as neither parameters take (that) large phenomenological values (i.e. Λ SS = 1.569 and N c = 3). It would be interesting to calculate, for instance, the leading Λ −1 correction to the 2-body potential and see if this could improve some phenomenological properties of our results, for example the large binding energies. Of course, this would turn the problem into a nonlinear one and make it severely more difficult than the one we have solved here.
The large-N c approximation may be a poor approximation for the physics at finite density, such as the baryon crystal studied in the last part of this paper, and subleading 1/N c corrections may alter our results and even the presence or absence of phase transitions, see e.g. Refs. [40,41]. This can be attributed to the fact that the number of nearest neighbors in a cubic lattice is 6 and N c = 3 in Nature, but in the large-N c limit, the number of colors is much larger than the number of nearest neighbors -this has consequences for, e.g. the van Der Waals forces, etc. , m > 0 , R > 0 , n ∈ Z ≥0 , (A.4) and we have defined x > 0 . (A.5) Next we need the Yukawa-like sum where k n ≥ m − + η − n and k n ≤ m + + η + n and the derivative with respect to the argument of G is (A.12)
Proof : Following the lines of the proof of Lemmata 1 and 2, we have where k n ≥ m − + η − n and k n ≤ m + + η + n and the derivative with respect to the argument of G is given by Eq. (A.12) and the double derivative yields Observation 1 The infinite sum ∞ n=0 a n with a n ≥ 0 remains finite by restricting to odd or even n.
Corollary 1 By means of the lemmata 1, 2 and 3 and the observation 1, all terms in the energies E 1,2 and E 3,4 of the cubic lattice in Eqs. (6.7) and (6.9), are separately finite.