Regge theory in a Holographic dual of QCD in the Veneziano Limit

We initiate the study of Regge theory in a bottom-up holographic model for QCD in the Veneziano limit, where the backreaction of the quarks to the gluon dynamics is included. We determine the parameters of the model by carrying out a precise fit to the meson spectrum in QCD. The spectrum for spin-one and pseudoscalar mesons is well reproduced. We then generalise the model to incluce higher spin fields in the bulk trajectories dual to the Pomeron and meson Regge trajectories at the boundary. With this setting, we fit the masses of the mesons with spins $J=2$, $3$, and $4$,as well as the experimental data of the total cross-sections $\sigma(\gamma \gamma \to X)$, $\sigma(\gamma p \to X)$ and $\sigma(p p \to X)$. For the cross sections we obtain a $\chi^2_{\mathrm{d.o.f.}}$ of $0.74$ for a total of 199 experimental points.


Introduction
Since it was conjectured that the QCD Pomeron is dual to the graviton Regge trajectory [1], holographic techniques have been successfully applied to the description of QCD processes where Pomeron exchange dominates . Once the external scattering states and the dynamics of the higher spin J fields of the graviton trajectory are modelled, comparisons can be made with experiment, provided we are in a kinematical window where QCD is dominated by a gluon rich medium. In this regime the Bjorken variable x is small, or the Mandelstam variable s is large, corresponding to high center of mass energies.
In order to describe the total cross-section data of hadronic processes in QCD, one also includes, besides the Pomeron trajectory, a meson trajectory (see, for e.g, [43]). This trajectory can be obtained by a linear fit of the meson spins against the meson squared masses. With the resulting straight line, one extrapolates from t > 0 to t ≤ 0 in order to make predictions in the scattering region. In particular, for the total cross-section we are interested in the value of the trajectory for t = 0, also known as the meson intercept, and the above linear fit yields an intercept of 0.55. There is no theoretical justification to assume that the meson trajectory is linear aside from the fact that it has described successfully scattering data, provided t is not too negative. In this work we use holography to study this issue by first fitting the meson spectrum for J = 0, 1, 2, 3, 4 and check if the resulting holographic intercept is able to describe the experimental data of the total cross-sections of γγ, γp and pp scattering. In previous works [34,38,42] we have studied the dynamics of the higher spin J fields in the graviton Regge trajectory by generalising the bulk graviton equation of motion. This was done using effective field theory inspired by Regge theory of a 5D string theory. In this work we will not only follow the same procedure for the Pomeron trajectory, but also apply it to the meson trajectory by generalising to higher spin the equation of motion of the bulk field dual to the vector mesons.
Before we start such procedure, we need to guarantee that our model is describing with accuracy the spectrum of the vector mesons. This will be done by considering the extension of the Improved Holographic QCD model of [44,45], with a backreacted quark sector [46,47], as presented in [48]. This model (V-QCD) consists of five-dimensional dilaton gravity dilaton coupled to a tachyon described in terms of a generalised Sen-like tachyonic Dirac-Born-Infeld (DBI) action [49]. Our numerical solution includes the full backreaction of the tachyon in the dilaton and metric. The asymptotic behavior of the model at weak and strong coupling is chosen such that various generic features of QCD, such as asymptotic freedom and confinement, are reproduced [45,48]. The remaining free parameters of the model will be determined through an extensive comparison of the spectrum of the quadratic fluctuations against the experimental meson masses. This paper is organized as follows. In section 2 we discuss in detail the holographic model, as well as how to compute the spectrum of the quadratic fluctuations. This section ends with a fit to the meson spectrum, fixing our background fields for the remaining of the paper. In section 3 we derive holographic expressions for the total cross-sections that will be used later to fit data from the Particle Data Group [50]. In section 4 we focus on the holographic duals of the pomeron and meson trajectories, and in particular in constructing the analytic continuation of the spin J equations that govern the dynamics of fields in these trajectories. These equations contain two parameters that will be fixed by the soft-pomeron intercept and the spin J = 2, 3, 4 meson masses in section 5. This fixes the pomeron and meson kernels that are used in the total cross-section fits. We discuss our results and suggest further work in section 6.

Holographic model for QCD in the Veneziano Limit
We consider a slight generalisation of Quantum Chromodynamics which consist of a gauge field in the adjoint representation of SU (N c ) coupled to N f fermions (quarks) in the fundamental representation of SU (N c ). This generalisation has been studied in great depth in the 't Hooft large-N c limit, where N c → ∞ and λ = g 2 YM N c and N f are kept fixed. This limit is also known as the quenched limit since nontrivial quark contributions to observables are suppressed in powers of N f Nc → 0. Another interesting large-N c limit is the Veneziano limit [51] where with x and λ fixed. In this limit the quark contributions are not suppressed, and their backreaction to the gluon dynamics must be taken into account. A holographic dual (V-QCD) that reproduces several expected features of QCD in the Veneziano limit was presented in [48]. It consists of a system of a dilaton and tachyon coupled to five-dimensional gravity. Let us first discuss the field content and the action of the model, and then present the precise structure of the various potentials appearing in the action.

The model
The action in the gravitation sector is the same as in the Improved Holographic QCD (IHQCD) model [44,45]. We work at zero temperature so that Poincaré invariance is intact. The metric Ansatz is therefore where the warp factor A is identified with the logarithm of the energy scale in the field theory at the boundary. The exponential of the dilaton field λ = e Φ is dual to the TrF 2 operator with its background value equal to the 't Hooft coupling (near the boundary where the coupling can be unambigously defined). The action for the metric and dilaton fields is given by five-dimensional Einstein gravity coupled to a scalar field, where M p is the five-dimensional Planck scale. The dilaton potential V g (λ) will be specified below. We choose the potential to be qualitatively similar to that studied in [52,53] -see [54,55] for an alternative approach focused on the fit to QCD thermodynamics.
To add matter we insert space-filling D 4 −D 4 branes that give rise to a tachyon field T and the gauge fields A L , A R living on the branes [46,47]. A similar approach has been considered in the probe limit in [56,57], and also in the Witten-Sakai-Sugimoto model [58][59][60][61]. In the boundary theory the operators with lowest dimension involving fermions areψ i R ψ j L with spin 0 and the two spin 1 conserved currentsψ i L γ µ ψ j L andψ i R γ µ ψ j R , where i, j are the flavour indices. The spin 0 and spin 1 operators are dual to bulk complex scalars T ij and two bulk gauge fields A µ L, ij and A µ R, ij , respectively. The tachyon transforms as (N f ,N f ) of the flavour symmetry U (N f ) R × U (N f ) L , while the fields A µ L,ij and A µ R,ij transform in the adjoint representations of U (N f ) L and U (N f ) R , respectively. In string theory the three bulk fields can be modelled by considering N f flavour branes (R) and N f flavour antibranes (L). In this configuration the complex scalar fields T ij are the lowest modes of open strings with one end in a D-brane and another in the anti-D-brane, while the bulk gauge fields are the lowest open string modes with both ends on a D-brane or on an anti-D-brane. The system obeys a tachyonic Dirac-Born-Infeld (DBI) action [47,48,62] where Str is the symmetric trace over the (hidden) flavour indices and the determinant is taken with respect to the five dimensional space-time indices a, b (since we are going to work up to quadratic order we can actually replace the symmetric trace by the usual trace of matrices for the purposes of this article). The functions V f , κ, and w will be given explicitly below. The normalisation convention for the symmetrisation of indices is The covariant derivative terms are given by 5) and the field strengths by In this work we are assuming that the light quark masses are equal and under this assumption the tachyon is just T = τ I N f . Furthermore, for the QCD vacuum we have A R a = 0 = A L a . Using these conditions in the action (2.4) we obtain the flavour action The brane action also contains a Wess-Zumino term [47] which we will not need as it does not contribute to the background solutions or mass spectra considered here. There is also an additional pseudo-scalar axion field a, which is dual to the operator TrF ∧ F and therefore sources the θ angle in QCD [47]. Its action takes the form [44,45,62,63] will be specified below, and we allowed the tachyon to have an overall phase ρ, i.e. we took T = τ e iρ I N f .

Choice of potentials
Let us now discuss the choices for the various potentials (V g , V f , κ, w, V (a) and Z). The generic picture is that the leading IR asymptotics of the various functions is chosen to agree with known features of QCD (which will be soon specified). The UV asymptotics is set by rough agreement with perturbative QCD, in particular by the perturbative UV dimensions of the QCD operators. For intermediate scales the functions need to be determined by and extensive comparison to experimental and lattice QCD data; in this article the functions will be fitted to the meson spectrum in QCD. In this section we will present the Ansätze for these functions which will obey the asymptotics at small and large λ determined by qualitative comparison to QCD, and the fit to data will be carried out in section 2.4.
We start with the action for the gluon sector, i.e. the action of improved holographic QCD. In IHQCD the dynamics of the dilaton is set by its potential V g (λ), which is constrained [44,45] in the UV to reproduce the YM β-function and in the IR to yield confinement and a "good" singularity in the classification of [64]. In this work we will take [65,66] V g (λ) = 12 The values of V 1 and V 2 are fixed by the gluon sector contribution to the QCD β-function, while the parameters α λ and V IR will be fitted by the spectrum, in particular by computing meson mass ratios and comparing them to experimental results, following [67]. The spectrum of the theory can be found by analysing the action of the quadratic fluctuations around the background solution followed by the reduction to the four-dimensional dynamics. For example, in the case of pure Yang-Mills (i.e. x = 0), the quadratic fluctuations are dual to glueballs with quantum numbers J P C = 0 ++ , 2 ++ and with J P C = 0 −+ by considering an axion term. For pure Yang-Mills we have found that α λ ≈ 2.504 and V IR ≈ 3.478 reproduce 1 the lattice ratios of 1.46 and 1.87 of m 2 ++ /m 0 ++ and m 0 * ++ /m 0 ++ respectively. In the IR, V g ∼ λ 4 3 (log λ) 1 2 , which gives linear asymptotic trajectories for glueballs. The DBI action that we described above is analogous to the flat space Sen action for the D −D system [49]. Since we are in the presence of a curved space-time and other non-trivial background fields which fully backreact to the metric, we correct it by including the general potentials V f (λ, τ ), κ(λ) and w(λ). However these potentials must satisfy some properties. The tachyon potential V f is expected to have a regular series expansion in λ and τ near the boundary (i.e. λ → 0, 11) and to vanish exponentially in the IR when τ → ∞ [63]. In particular in the flat space string theory where V τ (τ ) = 1 + a 1 τ 2 e −a 2 τ 2 , (2.13) (2.14) This Ansatz is identical to that considered in [66] except for the introduction of the coefficients a i in V τ (τ ). These are motivated by the observation that the correct mass gap of the mesons at large quark mass can only be reproduced if the coefficient of τ 2 in the exponent determining the large-τ asymptotics of V τ (i.e. a 2 above) differs from the second order series coefficient of V τ at small τ (here a 2 − a 1 ) [68]. Notice that one of these coefficients can be eliminated by adjusting the normalisation of the τ field.
As both κ(λ) and w(λ) are coupling functions under the square root of the DBI action we expect them to have similar qualitative behaviour. On the other hand, in order to have the correct UV dimension of theqq operator we need to impose The IR asymptotics of the potentials κ and w directly affect the meson spectrum. We are interested in the case where the mesons have an asymptotic linear spectrum and the meson towers have the same asymptotics. This can be achieved with the following IR asymptotics κ(λ) ∼ λ −4/3 (log λ) 1/2 and w(λ) ∼ λ −4/3 log λ [62,69] (see also [70]). Taking into account these considerations, we adopt the following Ansätze for κ and w, In order to compute the profiles of the background fields we need to specify the values of the parameters that appear in the definition of the potentials presented above. These parameters will be fitted to the ratios between the low-spin meson masses and the ρ 0 meson mass as predicted by the model. Finally, we need to specify the potentials in the CP-odd action S a (2.8). In flat-space tachyon condensation V (a) (λ, τ ) is independent of λ and is the same as the tachyon potential that appears in the DBI action. However in principle it may be different, so we will take V (a) to be V f defined above without the V 0f term. This form guarantees that it becomes a field-independent constant at τ = 0 and that it vanishes exponentially at τ = ∞. The Z(λ) function is defined by The definition is constrained by Yang-Mills theory [45,63,71]. In this work the parameters Z a and c a will be determined by fitting the spectrum of singlet axial vector mesons.

Evaluation of the meson spectrum
The quadratic fluctuations around the background fields can be mapped to the spectrum of mesons and glueballs. The normalisable fluctuations of dilaton Φ, QCD axion a, and the (traceless part of the) metric g ab , correspond to glueballs with J P C = 0 ++ , 0 −+ , and 2 ++ , respectively 2 . Here J is the spin, P refers to parity, and C refers to charge conjugation. The meson sector comes from the normalisable fluctuations of the tachyon T and of the gauge fields A L/R a . They correspond to mesons with J P C = 0 ++ , 0 −+ , 1 ++ , and 1 −− .
The fluctuations can be further classified according to how they transform under the vectorial SU (N f ). They can be grouped in flavour singlet and flavour non-singlet modes, i.e. mesons transforming in the adjoint of SU (N f ). The fluctuations that come from S f and S a are only flavour singlet, while the ones coming from S g include singlet and non-singlet terms. The singlet terms from S f will mix with with the singlet terms coming from S g and S a .
The masses of the different glueballs and mesons can be obtained after expanding the action S = S g + S f + S a to quadratic order of the fluctuations of the background fields. Due to flavour and rotational covariance the fluctuations decouple in separate sectors [62], apart from the mixing of the flavor singlet sectors mentioned above. In summary, there are flavour singlet rank-two tensor fluctuations (J P C = 2 ++ ), flavour singlet and non-singlet vector mesons (J P C = 1 −− ), flavour singlet and non-singlet axial vector mesons (J P C = 1 ++ ), flavour singlet and non-singlet scalars (J P C = 0 ++ ) and flavour singlet and non-singlet pseudoscalars (J P C = 0 −+ ). These fluctuations generate towers of 2 ++ glueballs, singlet and non-singlet vector mesons, singlet and non-singlet axial vector mesons, non-singlet scalar mesons and mixtures between 0 ++ glueballs and σ mesons, and non-singlet pseudoscalar mesons and mixtures between 0 −+ glueballs and η mesons, respectively. These towers of mesons and glueballs come as solutions of a Schrödinger problem associated with the equation of motion of the associated fluctuation. The eigenvalues correspond to the square of the mass of the glueball or meson and their holographic wave functions are the associated eigenfunctions. In this work we will not consider the flavour singlet states of J P C = 0 ++ , 0 −+ , as they involve mixing of the 0 ++ glueball with the flavour singlet σ meson and mixing between the 0 −+ glueball with the η meson, respectively. Therefore their analysis is considerably more challenging than that of the flavor nonsinglet states (see [62,63,72]) and would slow down the computer code for the spectrum significantly. Notice that these states are not central for the Regge analysis which is the main application of this work.
A detailed derivation of the equations of motion and Schrödinger problems associated with each fluctuation has been done in [62] and hence we will just summarise the main results relevant for the present work (see Appendix B for the analysis of the spin 1 fluctuations). The singlet and non-singlet vector mesons have the same equation of motion where ψ V = ψ V (z) is their wavefunction. By performing the change of variable defined by and rescaling one can rewrite the equation of motion in the Schrödinger form The singlet and non-singlet axial vector mesons have Schrödinger potentials differing by a term coming from the action S a . The potentials of the non-singlet axial vector mesons and of the singlet axial vector mesons are, respectively, The non-singlet scalar mesons have the potential where the expression for H S differs from that of [62] because our Ansatz for V τ is different. Finally the equation of motion of the non-singlet pseudoscalar fluctuations is given by 3 with associated Schrödinger potential The numerical determination of the spectrum proceeds by first finding the background solution (the metric and the scalar fields λ and τ ) of the equations of motion defined by the action S = S g + S f . Details of the numerical procedure can be found in appendix A. We then solve the fluctuation equations on top of the numerical background. For the cases of singlet and non-singlet vector and axial vector fluctuations, and for non-singlet scalar fluctuations, we compute the Schrödinger potential and use a pseudospectral method based on Chebyschev polynomials to compute the predicted masses of this model. The number of Chebyschev points used was 1000 and we checked the results were stable by computing the masses with a higher number of points. The reliability of the results was also studied by considering different IR and UV cutoffs on the background fields used to solve the Schrödinger problems. The masses of the pseudoscalars were computed using the shooting method. These methods were implemented in C++ and all results were also cross-checked against the (significantly slower) Mathematica code used in [62].

Fitting the spectrum
We now proceed to fix the parameters that appear in the potentials by comparing the predictions of our model with the experimental values of the meson masses quoted by the Particle Data Group [50].
The overall energy units in the model is also a free parameter. Its effect on the background and spectrum is trivial due to a scaling symmetry of the holographic model [48] which reflects the scale independence of the QCD Lagrangian. The scaling 3 Notice that the UV boundary condition for the pseudoscalar fluctuations is nontrivial and also depends on whether the quark mass is finite or not [62]. A consistent way which leads to UV finiteness of the fluctuated action in all cases is to require that the factor in square brackets (rather than the wave function ψ P ) in (2.28) vanishes in the UV. For the pseudoscalars we actually solved the differential equation by using a different method than in the other sectors (i.e. by shooting) because of the complication with the boundary condition.
1.300 0 −+ 1 π 0 (1800) 1.812 Table 1. Light unflavoured mesons from [50] used in the spectrum fit. The quantum number I = 1 means the meson is a flavour non-singlet state while I = 0 means the meson is a flavour singlet state.
symmetry implies, in particular, that the equations of motion are unchanged under the transformation By applying this transformation to the spectrum we see that all masses are scaled by the factor Λ. We will in effect choose Λ such that the numerical mass of the ρ meson matches the experimental result in GeV units. This is equivalent to fitting the parameters of the to the numerical values of ratios of masses (with respect to the ρ meson mass) instead of numerical values of masses. We will only consider mesons made of light up and down quarks. This sets the x parameter coming from the flavour sector to be 2/3. In table 1 we show all the mesons listed in [50] under light unflavoured mesons with the values of J P C mentioned before. The exceptions are the flavour singlet scalars and pseudoscalars and the a 0 (980). Whether the latter is a quark-antiquark state or a four-quark state is still debatable, although the literature favours more the four-quark state hypothesis. For this reason we did not include it in this work. In table 2 we have the mesons listed in [50] under other light unflavoured mesons, which are still not well established. We also included theses masses in our fit, therefore in case some of these states are not confirmed this work should be updated.
As our goal is to include the Regge behavior of vector and axial vector mesons, the most important criterion for the fit will be the deviation of the vector meson masses of the model from the experimental results. We have explored different fitting strategies. Table 2. Other light mesons from [50] used in the spectrum fit. The quantum number I = 1 means the meson is a flavour non-singlet state while I = 0 means the meson is a flavour singlet state.
We tested fits where the parameters α λ and V IR are fitted either independently to Yang-Mills data or together with the other parameters to "final" meson mass data. We also tried including lattice data for glueball masses. The result of these tests was that the optimal method, which lead to a physically sound solution for the metric and a good fit of the spin 1 states, was to do a global simulatenous fit of all parameters, excluding the glueball masses, and also imposing specific constraints to the fit parameters. We will explain the details below. The profile of the background fields and the mesons masses (excluding axial vector singlet states which will be discussed below) are determined by 17 parameters. 16 of these parameters (α λ , α κ , α w , W 0 , w 0 , κ 1 , w 1 , V IR , W IR ,κ 0 ,w 0 , W 1 ,κ 1 ,w 1 , a 1 , and a 2 ) are parameters of the potentials appearing in the action and τ 0 is a parameter that characterises the IR asymptotics of the tachyon field. In our fits a 1 and a 2 are fixed by imposing the following constraints: we choose a 2 − a 1 = 1 by rescaling the τ field, and set a 2 = 2κ(0) in order for the mass gap of the mesons to be correct at large quark mass [68]. This reduces the number of free parameters to 15. It turns out that it is useful to set extra constraints for the behavior of the tachyon which guarantee that the fit parameters remain in the domain of physically reasonable solutions. The first is related to chiral symmetry breaking. The chirally symmetric vacuum solution of the model flows to an IR fixed point [48]. We require that there is an instability towards forming a tachyon condensate in the IR around this fixed point, which will imply chiral symmetry breaking on the field theory side.
The presence of the instability, and therefore chiral symmetry breaking, is guaranteed if the Breitenlohner-Freedman (BF) bound [73] of the tachyon is violated at the fixed point. This means that −m 2 Here the location of the fixed point is the maximum of the effective potential, V eff (λ * ) = 0, and * is the IR AdS radius. Actually, while violation of the BF bound guarantees tachyon condensation and chiral symmetry breaking, it turns out that, in practice, values close to the bound are enough to trigger condensation. Therefore we will in fact require −m 2 τ 2 *

3.5.
The other condition is to require that the tachyon diverges fast enough in the IR to set all potential IR boundary terms arising from the flavour action to zero. This is required, among other things, for the correct implementation of the flavour anomalies [47,63]. For our choice of potentials the asymptotics of the tachyon in the IR is and the flavor action vanishes fast enough in the IR if τ c > 1, which is roughly what we require below. We then fit these parameters to ratios between the meson masses in tables 1 and 2 and the ρ mass, by minimising the function excluding the singlet axial vector mesons. This gives a total of 22 data points. The sum term is the absolute relative difference between the predictions of our model and the ones obtained by using experimental data, while the other two terms are the constraints we want our background to satisfy. We repeated the fit with different values of the W τ parameter in order to balance the ability of the model to reproduce the observed ratios and still be consistent and stable. We have found that W τ = 0.1 is a good choice and the results that we present below were obtained with such value.
Having fixed the background, we fit the parameters Z a and c a of equation ( Several remarks are in order. Firstly, the fit is stiff: while the number of parameters is large, the dependence of the results on their values is relatively mild. This is because the fit parameters appear through only a few functions of λ, the asymptotics Parameter value Parameter value Parameter value Table 3. Best fit parameters of the background potentials to the mass ratios between the mesons listed on tables 1 and 2 and the ρ meson.  at large and small coupling of which have already been determined by qualitative arguments and comparison to perturbation theory. Therefore the fit parameters essentially only affect the functions in the middle, at λ ∼ 1. Also there is limited parameter space where the functions are simple, monotonic functions, and one can check from the fit result that it indeed lies within this regime of the parameter space. Given the stiffness of the fit, the results for the spin-one mesons are really good. There are a few isolated states for which the deviation is 10%, but in general the deviations are in the ballpark of 1% or even less than that. The masses of the pseudoscalar mesons are also reproduced at a very good precision. There are, however, significant deviations in the scalar sector. While the scalar sector is challenging to explain in any model among other things due to the presence of significant four-quark contribution [50], the predicted nonsinglet scalar masses are still clearly too low, unlike in the probe limit study of [56,57] which used similar flavor action as the current article with the fixed background of [74]. While exploring different fit procedures, we noticed that there are parameter values for which the scalar masses are reproduced to a much better precision, but such parameter values are not favored by the overall fit which stresses the masses of the spin 1 mesons. Understanding this shortcoming requires further study. Notice that the scalar states are not needed for the analysis of the Regge trajectories which is the topic of discussion in the remainder of this paper. The fitted value of Z a = Z(0) in Table 3 is negative. Due to the positive value of c a , however, the function Z(λ) is mostly positive so that Z(λ) has a node at small λ. This behavior is unexpected and does not agree with phenomenology, i.e. the physics of the θ-angle in QCD and in particular the value of the topological susceptiblity [45,63]. Apparently this issue arises because the singlet axial meson masses are relatively unsensitive to the shape of the function Z(λ), and would be cured if additional observables (e.g. the topological susceptiblity) would be included in the fit. The precise functional form of Z(λ) is again irrelevant for the Regge analysis of the following sections.
It is also interesting to compare the fit results to those obtained in the same holographic model [66] by fitting the potentials independently to lattice data for QCD thermodynamics. Namely, most of the values are very close to those obtained in that study, deviations are typically in the ballbark of 10%. In particular, the scale parameters α λ and α κ are close to the value 1 used in this reference, whereas α w is a bit higher than what was obtained through the fit to thermodynamics (α w equals to 3w s of [66]) but still smaller than the other scale parameters, in agreement with the earlier fit. We also note that there is rough agreement with [75,76] where the model was compared to lattice results at finite magnetic field and temperature by using a sligthly different Ansatz for the potentials of the model. The parameter c of these references, which controls the scale in the λ dependence of the w(λ) function, maps roughly to the ratios α w /α λ or α w /α κ in this article. For the ratios we obtain numbers close to 0.5 whereas c ≈ 0.25 was preferred by the thermodynamics at finite magnetic field. That is, the numerical values are different, but clearly smaller than one in both cases. Moreover the value of W 0 (which was a free parameter in [66]) is determined by the spectrum fit to be near 2.5. This results is therefore an important constraint with respect to the earlier fit. The most significant difference between the fits is the value ofw 0 , which here is smaller by a factor of about 5 to 10 with respect to the various fits of [66]. This is apparently connected to the change in the value of α w .

γγ, γp and pp total cross-sections in holographic QCD
In this section we present the necessary ingredients to compute the total crosssections of γγ, γp and pp scattering in holographic models of QCD in the Veneziano limit. First we will discuss the kinematics of each process. Then we will present generic holographic expressions of the forward scattering amplitude, in the Regge limit, via the exchange of higher spin J fields. We conclude by deriving the holographic expression of the total cross-sections by taking the imaginary part of the amplitudes and using the optical theorem.

Kinematics
For all processes we will use light-cone coordinates (+, −, ⊥), with flat space metric For the γ * p → γ * p process the incoming and outgoing off-shell photons are the following while the incoming and outgoing protons with mass M have momenta The momentum transfer q ⊥ is a R 2 vector and is related to the Mandelstam variable t through t = −q 2 ⊥ . We work in the Regge limit of large Mandelstam variable s. For the forward scattering amplitude the momentum transfer q ⊥ = 0 and the photon virtualities satisfy Q 3 = Q 1 = Q, since the outgoing off-shell photon has k 3 = −k 1 and the outgoing proton k 4 = −k 2 . The incoming and outgoing photon polarizations are the same. The possible polarization vectors are where λ is just the usual transverse polarization vector.
For the γ * γ * → γ * γ * process the incoming photons have the four momenta while the outgoing photons have are the corresponding virtualities. As in the case of γ * p scattering, for the forward scattering amplitude the momentum transfer is null. The possible polarization vectors are, respectively, since in the forward amplitude the incoming and outgoing off-shell photons have the same polarizations. Notice that the transverse photons (λ = 1, 2) are normalized such that n 2 = 1, while for the longitudinal photons (λ = 3) n 2 = −1.
Finally, the large s kinematics of pp scattering is given by where k 1 and k 2 are the incoming proton momenta and k 3 and k 4 are the outgoing proton momenta. As in the other processes we will only compute the forward scattering amplitude for which q ⊥ = 0.

Holographic scattering amplitudes
Before we start the computation of the forward scattering amplitudes we need to define the external states as well as the interaction between them and the spin J fields that are exchanged in the Witten diagrams of figure 1. An external photon is a source of the conserved currentψγ µ ψ, where the quark field ψ comes from the open string sector. According to the gauge/gravity duality this field is dual to the nonnormalizable mode of a vector field in the bulk. In the context of this model the natural candidate is the linear combination of the A L and A R gauge fields In the string frame the action of this field is where g s is the determinant of the metric in the string frame,g ab is the inverse of g ab = g sab + κ s (λ)∂ a τ ∂ b τ , κ s (λ) = λ 4/3 κ(λ) and w s (λ) = λ 4/3 w(λ). Working in the gauge V z = 0 and ∂ µ V µ = 0, the vector field components describing a boundary plane wave solution with polarization n µ take the form where f Q satisfies the differential equation subject to the boundary conditions f Q (0) = 1 and ∂ z f Q (z → ∞) = 0. For the computation of the Witten diagrams that have photons as external states, in particular to compute the bulk interaction vertex, it is convenient to know the field strength of a given mode where we use the notationḟ Q = ∂ z f Q . For the proton external state we consider that it is dual to the normalizable mode of a bulk scalar field Υ (x, z) = e iP ·x υ (z) that represents an unpolarised proton. We will see that the contribution of the proton wavefunction to the scattering amplitudes is inside an integral that will be absorbed in the coupling constants and hence the precise details will not be important.
The last ingredient of our model are the higher spin fields h a 1 ···a J that will mediate the interaction between the external states in the considered scattering states. In this work we will consider bulk spin J fields that are dual to the spin J twist two operators made of the gluon field, as well as bulk spin J fields dual to the spin J twist two operators made of the quark bilinears. This extends the previous works [32,34,38], where only bulk fields dual to the gluon operators were considered. As discussed in appendix B, we will consider a coupling between the U (1) gauge field and these spin J fields given by while for the scalar field Υ dual to the proton state the coupling is In equations (3.14) and (3.15) ∇ a is the covariant derivative, while k J andk J are the couplings constants between the the U (1) gauge field and the bulk scalar with the spin J field, respectively. The higher spin J field h a 1 ···a J is totally symmetric, traceless and satisfies the transversality property ∇ a 1 h a 1 ···a J = 0. This implies that, in the Regge limit, it is not important in which external fields the covariant derivatives in (3.14) and (3.15) act. Below we assume that the spin J field has a propagator, without specifying its form. In the next section we focus on the dynamics of this field in detail for the case of pomeron and meson trajectories. Now we will show how to compute the forward scattering amplitude for the case of γp scattering. The calculations for γγ and pp follow the same pattern and bring no additional dificulty. For those cases we will simply present the results. In the Regge limit, the amplitude describing the spin J exchange between the incoming gauge field V (1) a ∼ e ik 1 ·x and scalar field Υ (2) ∼ e ik 2 ·x , and outgoing gauge field V (3) a ∼ e ik 3 ·x and scalar field Υ (4) ∼ e ik 4 ·x , can be written as where bars denote quantities evaluated atx in the Witten diagrams. The tensor Π a 1 ···a J ,b 1 ···b J (x,x) is the propagator of the spin J field. Using the kinematics of equations (3.1), (3.2) and (3.3), the expressions (3.13) and summing over the photon polarisations, the amplitude takes the form To make progress we make the change of variable x −x = (w + , w − , l ⊥ ) ≡ w and define the transverse propagator in transverse space G J (z,z, t) through that is valid both for spin J fields of the graviton Regge and meson trajectories.
In the next section we will propose phenomenological equations of motion for the higher spin fields h a 1 ···a J of both trajectories. In particular, it will be shown that the function G J (z,z, t) for the pomeron and meson trajectories admits a spectral decomposition associated to a Schrödinger potential that describes spin J glueballs or spin J mesons. This function can be written in terms of the eigenfunctions ψ n (J, z) and eigenvalues t n (J) of this Schrödinger potential in the following way (3.20) The function B(z) depends on the holographic QCD model as well on the trajectory that the higher spin field belongs to. We will determine them in the next section. Finally, in order to get the total amplitude we need to sum over the spin J fields with J ≥ J min , where J min is the minimal spin in the corresponding Regge trajectory. Then we can apply a Sommerfeld-Watson transform which requires analytic continuation of the amplitude for the spin J exchange to the complex J-plane. The contour on the complex plane consists of circles around simple poles at integer values of J. Then, we assume that the J-plane integral can be deformed from the poles at physical values of J to the poles J = j n (t) defined by t n (J) = t. The scattering domain of negative t contains these poles along the real axis for J < J min . The scattering amplitude for t = 0 is then By analysing (the regular solution to) the equation (3.12) we see that It therefore follows from the optical theorem that where we have made the change of variable du = Gdz.
As mentioned before the procedure to compute the holographic total crosssections for the other processes is similar to the one just presented. These calculations are done in appendices C and D. Hence, we finish this section by presenting the final results. For γγ → X we have we have The first set of operators are gluonic while the second set are quark operators. These twist 2 operators are dual to bulk spin J fields whose dynamics will be specified below. We shall then follow an effective field theory approach, by proposing a general form of the equations of motion with phenomenological parameters that can be fixed by data. We will propose two different equations of motion, one that describes the gluon sector, which includes the energy-momentum tensor T αβ , and the other the quark sector, which includes the quark bilinear current J α . These equations will satisfy two requirements: i) compatibility with the graviton's equation of motion for the case J = 2 for gluonic operators, and compability for J = 1 with the equation of motion of the U (1) current dual to the operator J α ; ii) reduction to the conformal limit case (pure AdS space and constant dilation and tachyon).
In pure AdS the equation of motion of a spin J field is where L is the AdS length scale and ∆ is the dimension of the dual operator O J . These fields are symmetric, traceless and transverse (TT). The independent components are the ones along the boundary direction (i.e. h α 1 ···α J ) due to the transversality condition. Decomposing these into irreducible representations of the Lorentz group SO(1, 3), the TT components h TT α 1 ···α J decouple from the others and they describe the operator O J in the dual theory. The UV asymptotics of these fields are with the free theory value ∆ = J + 2, where J is the source of O J and O J its vacuum expectation value.
We will now motivate two equations of motion for spin J, one for the fields in the Pomeron trajectory and another for the fields in the meson trajectory.

Pomeron in Holographic QCD in the Veneziano limit
The Pomeron is dual to the graviton Regge trajectory [1]. To derive the equation of motion of the spin J fields in the graviton's Regge trajectory we follow the same approach as in [32,34], that is we derive first the equation of motion of the graviton and then generalize it based on the two requirements mentioned above. Only the gluon part of the action S g , contributes to the fluctuations in the spin-2 sector. Therefore the equation of motion of the TT components of the graviton in this class of models is, in the Einstein frame, the same as in IHQCD [45] We propose that the equation of motion of the spin J fields in the graviton's Regge trajectory is 6) where e g is a constant that will be fixed later by setting the soft pomeron intercept to 1.08. We note that: i) for J = 2 we get the the graviton's equation (4.5); ii) for the conformal case, i.e. A = − log (z/L) and Φ and τ constant, the equation reduces to (4.2); iii) the second term comes from the tree level coupling of a closed string, as appropriate for the graviton Regge trajectory in a large N approximation; iv) following an effective field theory rational we could have included other terms proportional to derivatives of A s , Φ and τ , that is terms proportional to e −2As Ä s −Ȧ 2 s , e −2AsΦ2 , e −2AsΦ , e −2AsȦ sΦ . (4.7) All these terms, of dimension inverse squared length, are compatible with the constraint i) and also with constraint ii) provided they multiply J − 2. Like the term proportional toτ 2 , these terms are all subleading in the UV. However in the IR, where the wavefunctions of the associated Schrödinger problem are localised, the termτ 2 dominates and for this reason we will only consider this term. We have also not included the termsτ ,τȦ s andτΦ because the background is symmetric under τ → −τ . The third term in (4.6) is a mass term obtained by the analytic continuation of the dimension of the exchanged operators ∆ = ∆ (J). We shall set In the boundary theory the dimension of the operator O J can be written as ∆ = 2 + J + γ J , where γ J is the anomalous dimension. In free theory γ J = 0. The term we added in the right-hand side of (4.8) ensures the correct UV asymptotic behaviour of (4.6) leading to the free theory result h TT +···+ ∼ z 2 . Beyond perturbation theory, the curve must pass through the point J = 2 and ∆ = 4 as it represents the energy-momentum tensor which is protected. Equation (4.8) guarantees such properties.
The propagator for the spin J fields in the graviton's Regge trajectory is the solution of where the notation g a 1 (b 1 · · · g |a J |b J ) means that symmetrisation is applied only among the indices b i , i = 1, · · · J. As we have seen in the previous section, in the Regge limit we are only interested in the component Π +···+,−···− . By using the identity (3.18) one can show that where ψ satisfies the Schrödinger equation (4.14) The spectrum for each integer J discretises t = t n (J) and the corresponding eigenfunctions satisfy the identity n ψ n (z)ψ * n (z) = δ(z −z). Hence, the solution to equation (4.10) is given by

Meson Trajectory
Next we discuss how to describe the dynamics of spin J fields in the mesons trajectory. The Chew-Frautschi plot of figure 2 shows two important properties of the mesons trajectories. The first is linearity which is usually used to extrapolate from the tchannel physical region to the s-channel scattering region where t < 0. That is, we find the best line to the points and use it to find j(t) for t < m 2 ρ . In particular, for total cross-sections we are interested in the intercept value at t = 0. Another important property is the near degeneracy of the four trajectories {f 2 (1270), f 4 (2050), · · · }, {a 2 (1320), a 4 (2040), · · · }, {ω(780), ω 3 (1670), · · · } and {ρ(770), ρ 3 (1690), · · · }.
Using the fact that the four trajectories are nearly degenerate and that we had a very good description of the spectrum for non-singlet and singlet vector mesons (ρ and ω) we will construct the holographic meson trajectory by generalizing the equation of motion of vector mesons to any spin J field in the same trajectory. We will follow the same procedure of last section. Later we will validade our approach by showing how, with just one single parameter, one can simultaneously have a good description of scattering data, the spectrum of the meson trajectory, and an approximately linear Regge trajectory for the mesons.
From the action (3.10) for the vector field in the string frame, one gets the equation of motion ∇ a e − 10 3 Φ V f w 2 s Gg acgbd F cd = 0 . In order to simplify the notation, from now on we assume all the warp factors and effective metrics are in the string frame. As a first step to generalize (4.16) we write it asg whereṼ = e − 10 3 Φ V f w 2 s . In the conformal limit (pure AdS with constant dilaton and tachyon) the above equation reduces to with (LM ) 2 = −4, as expected for a bulk field dual to the current operatorψγ µ ψ with spin J = 1 and protected dimension ∆ = 3.
We now propose an equation of motion for the symmetric, traceless and transverse spin J field h a 1 ···a J in the meson Regge trajectory. As in the graviton case, we are interested in the TT part h TT α 1 ···α J which are the propagating degrees of freedom that decouple from the other components after decomposing the field in SO(1, 3) irreducible representations. We propose the equation We note that: i) this equation reduces to the vector meson equation of motion for J = 1; ii) in the AdS case the second, third and fifth terms vanish, reducing this equation to the equation of motion of the TT components of spin J fields in AdS; iii) Following the same logic as for the fields in the Pomeron trajectory we only included the two-derivative term inτ 2 that is compatible with i) and ii). Also, as in the case of the Pomeron, the ∆(J) curve follows from the analytic continuation of the dimension of the exchanged operators and imposing the correct UV asymptotic behaviour of the spin J fields: (4.20) The last term guarantees the correct UV behaviour of the spin J fields given in (4.3), while the first one ensures that the curve passes through the protected point ∆ = 3, J = 1. Equation (4.19) can be brought to Schrödinger form by first rewriting it in terms of the u variable defined previously. After that we can write the spin J field as where ψ satisfies the Schrödinger equation The Schrödinger potential is given by where V V is the Schrödinger potential of the vector mesons given in (2.23). Here the dots mean derivatives with respect to the u variable. We will denote the eigenvalues of this Schrödinger potential as t n (J) and the corresponding eigenfunctions as ψ n (J, u).
The propagator for this spin J field obeys an equation of the type where D is a differential operator defined by (4.19). In the Regge limit we will be interested in the components Π +···+,−···− and for this particular case Consider now the integral where D 3 is the differential operator we conclude that G J satisfies where y andȳ are coordinates in transverse space, as defined after (4.10).
To solve (4.29) we consider its homogeneous version and the following Ansatz It follows that ψ is a solution of the Schrödinger problem (4.22). Since the eigenfunctions of the Schrödinger potential satisfy n ψ n (u)ψ n (ū) * = δ(u −ū) the solution to (4.29) is Thus, the implications of this result for (3.25), (3.26) and (3.28) is that when considering exchange of reggeons in the meson trajectory.

Fit of γγ, γp and pp total cross-sections in holographic model
In this section we will test the presented phenomenological model for γγ, γp and pp total cross-sections against the hadronic cross-section data files from the Particle Data Group [50]. These data sets are formed from experimental results obtained by different groups over the last decades. The datasets of σ(γp → X) and σ(pp → X) have cross-section values as a function of the laboratory momentum of an incoming on-shell photon or proton, respectively. A calculation of the respective center of mass energy √ s was performed before starting the fits. We also considered only subsets of data with √ s > 4 GeV for σ(γγ → X), σ(γp → X) and σ(pp → X), yielding 39, 45 and 115 experimental points, respectively.
We find the best set of parameter values α i by minimising the χ 2 quantity that is, the sum of the weighted difference squared between experimental data and model predicted values where the weight is the inverse of the experimental uncertainty. In our fits the parameters α i are the couplings k g/m jn andk g/m jn defined in equations (3.23), (3.27) and (3.29), where the superscript refers the coupling to the pomeron or meson trajectory. Usually a fit is deemed of good quality if the quantity where N par is the number of parameters α i to be fitted. In (5.1) O k represents a generic data point of one or several observables mentioned  previously and, usually, σ k is the experimental uncertainty associated with the measurement. Some data points have uncertainties in the values of s (e.g. in γγ → X that is always the case because it is measured experimentally). To account for this we calculate the total cross-section for s + ∆s and s − ∆s, and evaluate For these cases σ k = (σ exp. ) 2 + (σ eff. ) 2 where σ exp. is the experimental error. By minimising (5.1) with all the data mentioned above we can find the best values for the potentials parameters e g and e m , as well as for the coupling values k jn and k jn with n = 1, 2, · · · for each trajectory. We will follow another approach. As shown in [34], the first two trajectories of the gluon kernel can be identified to the hard and soft pomeron trajectories. Here we will fix e g by demanding that the intercept of the soft pomeron trajectory is 1.08. We also identify the meson trajectory with the first trajectory of the meson kernel and we will fix e m such that we have agreement with the meson masses of (f 2 , a 2 ), (ρ 3 , ω 3 ) and (f 4 , a 4 ) for J = 2, 3, 4, respectively. The numerical values of these parameters were found to be 0.246 and 1.712 for e g and e m , respectively. The intercept values obtained from these parameter values can be found in table 5 and the leading trajectories are shown in figure 3. The corresponding masses of the mesons for J = 2, 3, 4 are present in table 6. The intercept of the hard pomeron is close to 1.17 which was the one found in the IHQCD model [34,38]. We also note that the intercept of the 4th Pomeron trajectory and of the meson trajectory are very close to 0.55 found in the hadronic cross-section fits of [43].
After fixing the kernel parameters e g and e m , we fit the total cross-section data by minimising (5.1) with respect to the couplings k  γ coupling value proton coupling value

12.8753
in figures 4, 5 and 6. These results show that one does not need to assume a linear meson trajectory after fixing it at t > 0 in order to describe total cross-section data, as it is often assumed in the literature, leading to an the intercept value of 0.55. The Regge trajectory is convex line, which yields a slightly higher intercept.

Conclusions
In this article, we studied Regge theory in a full-fledged holographic model (V-QCD), which includes backreaction of quark degrees of freedom to the gluon dynamics in QCD. The main new results can be divided into two categories: Firstly, we made progress with the comparison of the model with QCD data by carrying out a detailed fit of the model paratemeters to the meson spectrum. Secondly, we developed a scheme to describe higher spin mesons and Regge trajectories in this model, and applied it to analyse the total QCD cross sections of scattering precosses having protons and photons in the inital state. As explained in section 2, the holographic model is strongly constrained by the requirement that it agrees with known features of QCD such as confinement, chiral symmetry breaking, asymptotic linearity of meson trajectories, qualitatively correct dependence of the spectrum on the quark mass, correct response to small chemical potential at small temperatures, and asymptotic freedom with correct dimensions (and anomalous dimensions) of the most important operators of QCD at weak coupling. Most of the remaining parameters, which are not determined by such qualitative considerations, amount to tuning of the various potentials of the V-QCD action at intermediate values of the coupling. In this article we have chosen to tune these parameters such that the meson spectrum of the model agrees well with experimental QCD data. The number of fitted parameters is large, because we want to make sure that our Ansatz for the potentials covers essentially all of the parameter space left free by the constraints listed above. However since the effect of all these parameters on the potentials is relatively small, the dependence of the result for the meson masses on the parameters is weak. In other words, one obtains a rather good description of the QCD spectrum for any reasonable values of the parameters, and the task carried out in this article is to tune the masses to agree as well with experimental values as possible.
Because we were mostly interested in the Regge physics in this article, we chose a strategy where we only fitted the meson masses with spins J = 0 and J = 1, but did not consider other data such as decay constants or thermodynamic potentials. Notice also that we included radial excitations with high masses and these state were fitted with the same weight as the "important" low-lying states such as the pions and the ρ-meson. That is, the fit was tailored for the purpose of studying the Regge physics where reproducing the correct asymptotics of the trajectories is important. It is anyhow interesting that the results of the fit are in good agreement with those obtained in [66] where essentially the same Ansatz was compared to the lattice data for the thermodynamics of Yang-Mills theory and QCD at finite temperature, and in rough agreement with the value of the parameter c of [75,76] as we explained in section 2.4. We also remark that the fit carried out here appears to be more constraining than those carried out in the earlier references: our fit favors W 0 ≈ 2.5 whereas this parameter was left unconstrained by the comparison to lattice thermodynamics. We obtained a very good fit for the spin 1 and pseudoscalar meson masses. This fit was more extensive than that carried out in the probe limit in the closely related model of [56,57]. It also compares favorably to work in simpler holographic models and in models inspired by gauge/gravity duality, such as the hard [77,78] and soft wall models [79], light front holography [80,81], and the holography inspired stringy hardon model [82,83].
There are several ways to further develop the holographic model and the fitting procedure in the future. A simple project would be to redo the fit to meson spectrum with the aim of producing a model for all purposes, which would mean to weight more the mesons having low masses and also consider other experimental data relevant for the zero temperature vacuum such as decay constants, the flavor singlet scalar and pseudoscalar states, topological susceptibility, the S-parameter and so on. The potentially challenging issue, which we noticed while doing the fit of this article, is that the scalar meson masses agree poorly with the experimental values. This issue would need to be solved. A more ambitious project would be to carry out a simulatenous comparison of the zero temperature and finite temperature data, i.e. to also include the data for lattice thermodynamics at finite temperature and potentially also at finite magnetic field. Good agreement of the model parameters obtained by fitting the zero and finite temperature data independently suggests that such a project is feasible. As a part of the project, it might make sense to also generalise the model to include flavor dependent quark masses and flavor dependent coupling of the magnetic field to the quarks.
Having fixed the action for the geometry and the actions for scalar and vector mesons through the fit, we proposed the dynamics for the spin J fields dual to the gluon and quark twist two operators. This dynamics is controlled by two parameters that were fixed by making the second Pomeron trajectory intercept to be 1.08 (the known soft-pomeron intercept) and by reproducing quite accurately the masses of the mesons with J = 2, 3, 4. We have found that the hard pomeron intercept is close to 1.17 in IHQCD. We also note that the fourth Pomeron trajectory and the leading meson trajectory that we obtained have intercepts close to the meson intercept value of 0.55 commonly found in the literature. Using the first three pomeron trajectories and the first meson trajectory we had a very good fit of the total cross-sections of γγ, γp and pp scattering. We note here that although our meson intercept is close to 0.55, the corresponding trajectory is non-linear in the t > 0 region (see figure 3), as it is often assumed. This suggests that the linearity of the trajectory might be at best a very good approximation to obtain an intercept that explains total cross-section data, and not its true shape in this region.
This work can now be extended to other processes like the proton structure functions F p 2 and F p L and the photon structure function F γ 2 , as done in [42]. In this model their holographic expressions are given respectively by equations (E.3), (E.4) and (D.3) from the appendix E. From these expressions one can see that: i) the Bjorken x behaviour depends solely on the value of the intercept; and ii) the Q 2 dependence comes from an integral involving the non-normalizable mode of the U (1) gauge field f Q , functions of the background fields and the normalisable wave functions ψ n of the pomeron and meson kernels. Having the background fixed, the first two classes of functions are uniquely determined while the ψ n 's might be controlled by phenomenological parameters (similar to e g and e f ) of the equations of motion of the spin J fields. In this work, with the values of e g and e f fixed by the spectrum of higher spin mesons, we may have obtained the correct Bjorken x behaviour, but it is not guaranteed that the wave functions yield the observed dependence of the structure functions on Q 2 . Hence, one needs to consider the neglected terms of equation (4.7) in order to get good agreement between our model and data. Of course this involves adding more parameters to these fits, making it harder to find a minimum for the χ 2 function that describes satisfactory all the data being considered. If successful, this would also extend [42] by including also pp total cross-section data in a consistent holographic Regge analysis.
Since this model reproduces well the masses of the towers of ρ and π mesons, one could also test this approach against the Vector Meson Production (VMP) data from HERA with a ρ meson in the final state and to include π 0 p total cross-section data. This will introduce extra parameters in this holographic Regge model. In the π 0 p total cross-section we could introduce couplings between the pions and the spin J fields of the different trajectories. On the other hand, the VMP data would promote the coupling constants k J to functions of t since in differential cross-section data we work with more than one value of t. In order to control the amount of parameters to be introduced in this model, it would be interesting to fix the couplings, or the functional form of k jn(t) , to other observables where these couplings might be important. This could be done by first using a theoretical well motivated ansatz for k jn(t) with some free parameters. These free parameters could then be fixed by the experimental values of decay rates. As an example of this, k m and Γ(f 2 → γγ) are related through equation (3.14) for the J = 2 case.
To describe successfully VMP in this holographic setting we need a better approximation to the proton state than assuming it to be a scalar field. In this work this is sufficient, since the integrals that involve the proton wavefunction are absorbed in the fitting parameters. However, in the VMP case the proton state is in an integrand multiplying wavefunctions of the Pomeron and meson kernel that depend on t. Moreover, having a good model for the proton state could also allow to extend this work directly by including available data of pp differential cross-section data.
In holography, baryons are dual to solitons in the bulk and hence the problem is reduced to the calculations of these solutions [84]. Such solitons have been studied in the literature in the Witten-Sakai-Sugimoto model [85,86] and in hard-wall models [87,88]. In the V-QCD model, baryons have so far only been considered by employing an approximation scheme [70], and the construction of the soliton solutions is work in progress. Such a solution could be a good starting point to model the proton state in holographic Regge theory. with the fellowship PD/BD/114158/2016. The research of MJ was supported by an appointment to the JRG Program at the APCTP through the Science and Technology Promotion Fund and Lottery Fund of the Korean Government. MJ was also supported by the Korean Local Governments -Gyeongsangbuk-do Province and Pohang City.

A Solving the Equatios of Motion
In this appendix we give details on how the equations of motion are solved. It turns out to be convenient to write the resulting equations in terms of A instead of the radial coordinate z. This change of coordinates stretches distances close to the boundary which eases the numerical UV analysis. To implement this one introduces the new variable q(A) = dz dA e A , (A.1) so that the metric reads From now on all the background fields will be functions of A and the differential equations we present are all with respect to A being the independent variable.

Yang-Mills Equations of Motion
If we set x = 0 we obtain the action of the IHQCD model and hence we are dealing with a pure Yang-Mills theory in the large N c limit. The equations of motion are The solution of the pure Yang-Mills background is needed in our method for constructing the full V-QCD solution, as we explain below. We will therefore solve this coupled set of differential equations by shooting from the IR to the UV. The coordinate A runs from −∞ to +∞, but naturally we will need to intriduce cutoffs for the numerical solution. We set A IR = −150 and A UVYM = 50 as the lower and upper bound on A. In the IR the geometry ends in an IR singularity of the "good" kind according to the classification of [64]. At the singularity the warp factor A and the dilaton Φ have the asymptotic forms [48] A as z → ∞. This allows us to determine z IRYM such that A IR (z IRYM ) = −150 and use its numerical value to compute Φ (A IR ) = Φ IR (z IRYM ) and q (A IR ) = e A IR / dA IR dz and hence defining boundary conditions to solve equations (A.3) and (A.4). Note that α λ and V IR determine uniquely both the initial conditions as well as the evolution of the background fields in pure YM.
After we solve the YM equations we can determine z (A) by using equation (A.1). We solve this numerically by taking z(A IRYM ) = z IRYM as initial condition and at the end we perform the shift z (A) → z (A) − z (A UVYM ) in order to have the UV singularity at z = 0.
We finish this section by specifying which algorithms we have used to find the numerical YM background. To compute z IRYM we have used the root finding Van Wijngaarden-Dekker-Brent method described in [89]. The differential equations (A.3) and (A.4) were solved using integrate_const with a Runge-Kutta-Dormand-Prince stepper from the Boost C++ library.

Solving the equations of motion in the holographic model
The solution of the equations of motion in V-QCD at zero temperature is done in four stages. This is because, as it turns out, the equations of motion are stiff in particular close to the IR singularity. Therefore for the numerical code to be stable, it is better to divide the range of A into region where different kind of approximations can be used that reduce the stiffness problem. In the first two regions we are deep in the IR and the tachyon is decoupled from the other background fields and q (A) so that Φ (A) obey equations (A.3) and (A.4). In these stages we therefore use the YM background constructed as explained above for the metric and for the dilaton. In the third stage the tachyon will couple to the other background fields until it reaches the UV where it decouples again starting the fourth and last stage of the solution of the EOMs.
As in the YM case we start by first specifying the IR boundary conditions. We first define A IR = −150, A UVYM = 50, A UVc = 100 and A UVf = 1000. Again we compute z IRYM as in the YM case and we use it to compute q (A IR ) and Φ (A IR ) as in YM and τ (A IR ) using the tachyon IR asymptotics where τ 0 is a parameter that is going to be fitted to the spectrum. The constants τ cut = 1000 and V f cut = 10 −8 are defined and they mark the end of the first and second stages of the construction of the numerical background, respectively. The difference between the first and the second stage is that in the first stage, the tachyon is so large that nonlinear corrections ∼ 1/τ 2 can be ignored in the tachyon equation of motion. In the second stage we need to use the full tachyon equation of motion, while the tachyon remains decoupled (as signaled by the smallness of the tachyon potential, We start the construction of the background by computing the YM profile of q (A) and Φ (A) from A IR to A UVYM as explained above. If τ IR > τ cut the tachyon profile will be given by the solution of the differential equation (linearized at large τ ) until τ < τ cut , marking the end of the first stage. The value of A such that this condition is met is called A UV1 The values of q and Φ used are the ones given by YM.
In the second stage we solve the (full) tachyon differential equation In the stage where the tachyon is coupled the dynamics of the background fields obeys This is a stiff system of differential equations and for this reason we had used the function integrate_adaptive with the rosenbrock4 stepper from the Boost C++ library. This system is solved from A UV2 to A UVc . The value A UVc = 50 was chosen such that at this point we are in the UV and the tachyon decouples again from q and Φ: this is guaranteed as the tachyon is suppressed exponentially in the UV (τ ∼ m q e −A ).
The UV equations of motion are 16) where τ n = e A τ . This system is solved from A UVc to A UVf . With the profile of τ n and λ in this region we can compute an estimate the quark mass m q through the UV asymptotic formula where UV is the UV AdS radius. We then obtain the final estimate for the quark mass by evaluating this estimate at two large values, i.e. at A = A UVf and at B EOM and couplings of the U (1) gauge field In this appendix we start deriving the action of the vector meson sector since the nonnormalisable solution of the corresponding equations of motion is dual to the external photon in the boundary. Then, by linearising the action we find the coupling of the U (1) gauge field to the graviton of the bulk theory. All of this is made in the Einstein frame. Since we will be dealing later with calculations on the string frame we will explain how to translate these results to the string frame. Finally we generalise the coupling to the graviton to any spin J field in the graviton's Regge trajectory.
The action of the vector U (1) gauge field As mentioned previously, for the QCD vacuum we set A R a = 0 = A L a . When we turn on the gauge fields, the flavour action becomes Here we remember that Tr is a trace over the flavour indices while the determinant is computed with respect to the space-time indices. Writing the determinants inside the square roots as defining X L/R = g −1 eff. T L/R and using the identity The expressions for the matrix elements of g −1 eff. and T L/R are where the vector and axial gauge fields V a and A a are linear combinations of the left and right gauge fields From these identities and after some calculations one gets The first term is the background term, while the the two other terms can be packed into the actions which are the actions for the Non-singlet Vector and Axial-Vector mesons presented in [62]. The equation (B.14) can still be rewritten as where here the notation F ab F ab means i.e. S f in the string frame is obtained by simply substituting the metric by the string frame metric g s , substituting κ(λ) and w(λ) by κ s (λ) = e 4Φ/3 κ(λ) and w s (λ) = e 4Φ/3 w(λ) and multiplying V f by the factor e −10Φ/3 . The derivation of S V in the string frame will be formally the same and hence the results equal to the Einstein frame ones but with w s , κ s and e −10Φ/3 V f in place of w, κ and V f respectively. The action (B.16) in the string frame takes the form Couplings with the spin J fields We will now determine the gravitational coupling between the U (1) gauge field and the spin J fields in the graviton's Regge trajectory. We first compute the coupling with the graviton and generalise to any even spin J field. All of this is done in the Einstein frame. To find the coupling in the string frame we just substitute the functions w, κ and V f by w s , κ s and e −10Φ/3 V f respectively, as discussed previously. Again, we start by writing the square roots of the determinants as The coupling with the graviton is found by linearising equation (B.21) around the background metric, i.e. g ab =ḡ ab +h ab . To study the graviton Regge trajectory in our background we need to decompose the metric in SO(1, 3) irreducible representations. We will be only interested in the graviton TT components h αβ , satisfying ∂ α h αβ = 0 and h α α = 0, and also set h zα = h αz = h zz = 0. For our purposes we can ignore the perturbation of √ − det g eff. because it involves only a term proportional to h = h a a = 0. We will also neglect the terms involving the axial vector mesons A µ since we are only interested in the coupling of V µ with the graviton for now. We wish then to compute δtr g −1 eff. T L/R and δtr g −1 eff. T L/R g −1 eff. T L/R , where by δ we mean a perturbation relative to the background metric. Using the identity δg ab eff. = −g am eff. g bn eff. h mn one can show that Using the expressions for the matrix elements of g −1 eff. and T L/R we get δtr g −1 eff. T L/R = 0 , (B.24) Hence the coupling between the vector U (1) gauge field and the graviton is given by In the string frame this coupling takes the form We now generalise this coupling to the case of an interaction between the gauge field and a symmetric, transverse and traceless spin J field, h a 1 ···a J . The pomeron trajectory includes such higher spin fields of even J. Again there are several possibilities, but we shall focus on the simplest extension of the graviton coupling considered above. For a spin J field we take the coupling We note that the transverse condition of the spin J field h a 1 ···a J guarantees that this term is unique up to dilaton and tachyon derivatives. We now consider the coupling between the external photon states with the bulk spin J fields dual to the spin J twist two operators made of quark bilinears. To determine the coupling to any spin J in this trajectory we could proceed analogously with the case of the graviton's Regge trajectory. In the case of the coupling with the meson trajectory, one could first determine the coupling between the non-normalizable mode dual to the photon with the the ρ meson states and generalise the result to higher spin J fields. To do this we attempted to expand the DBI action to cubic order in the fluctuations and keep only the terms with three vector gauge fields V a V b V c . We start by writing where we have dropped terms involving products of TrX L/R because they contribute only to axial gauge fields A a . The quadratic term leads to the action of the quadratic fluctuations of the vector gauge field and hence the coupling V a V b V c if exists must be contained on Tr X L/R 3 . That is, we evaluated the expression Another approach is to find the coupling between the vector gauge field with the bulk field dual to the f 2 meson and extrapolate to all the other spin J fields in the meson trajectory. In [90] the tensor meson f 2 state is the first Kaluza-Klein mode of a bulk spin-2 h ab field that has the same equation of motion as the graviton in AdS 5 .
The coupling of f 2 to the photon is also the same as the one between a graviton and a bulk gauge field in AdS 5 . The geometry is basically AdS 5 with a wall whose position is fixed by the mass of the ρ meson. After this they are able to predict not only the mass of f 2 but also the decay width Γ(f 2 → γγ). f 2 also has the same quantum numbers of the tensor glueballs J P C = 2 ++ which are the normalisable modes associated with the graviton's equation of motion. For these reasons, in this work, we will assume that the coupling of the U (1) gauge field with the f 2 meson is the same as the coupling with the graviton and hence, in general, the coupling of any bulk spin J field in the meson trajectory is also given by equation (B.29).

C pp scattering
In this appendix we will present the computation for the total cross-section of pp scattering. The steps of the computation are the same as in the case of γp show in the main text.
The scattering amplitude for spin J exchange between two incoming scalar fields Υ (1) ∼ e ik 1 ·x and Υ (2) ∼ e ik 2 ·x is As in the γ * p case, in order to get the total amplitude we need to sum over the spin J fields with J ≥ J min , where J min is the minimal spin in the corresponding Regge trajectory. Then we can apply a Sommerfeld-Watson transform In the scattering domain of t < 0 these poles are in the real axis for J < J min . This procedure yields equations (3.28) and (3.29).
D γ * γ processes In this section we derive the holographic expressions for F γ 2 and σ(γγ → X) in the context of Holographic QCD in the Veneziano limit. We will consider the photon structure function F γ 2 and the total cross-section σ(γγ → X). Like in the case of the proton structure function F p 2 the photon strucutre function F γ 2 is related to the transverse and longitudinal total cross-sections of γ * γ scattering by The calculation of the forward scattering amplitude for γ * γ is the same as in γ * p scattering except that the external state in the Witten diagram of figure 1 is an on-shell photon. This means that the definition of Im γγ n should be proportional to

E Holographic structure functions
The forward scattering amplitude for γ * p scattering is given by equation (3.22). That equation was obtained summing the contributions of the transverse and longitudinal polarisations of the off-shell photon. In particular the term with f 2 Q is the contribution from transverse polarisations while the termḟ 2 Q is the contribution for the longitudinal polarisation.
The proton structure functions F p 2 and F p L are related to the transverse and longitudinal total cross-sections of the process γ * p by Using the optical theorem and the last relations one finds the contribution of the holographic expressions to the structure functions are where the definition of g γp n is the one of equation (3.23). The function B will depend on whether the spin J fields belong to the pomeron or meson trajectory.