Effective field theory bootstrap, large-N $\chi$PT and holographic QCD

We review the effective field theory (EFT) bootstrap by formulating it as an infinite-dimensional semidefinite program (SDP), built from the crossing symmetric sum rules and the S-matrix primal ansatz. We apply the program to study the large-$N$ chiral perturbation theory ($\chi$PT) and observe excellent convergence of EFT bounds between the dual (rule-out) and primal (rule-in) methods. This convergence aligns with the predictions of duality theory in SDP, enabling us to analyze the bound states and resonances in the ultra-violet (UV) spectrum. Furthermore, we incorporate the upper bound of unitarity to uniformly constrain the EFT space from the UV scale $M$ using the primal method, thereby confirming the consistency of the large-$N$ expansion. In the end, we translate the large-$N$ $\chi$PT bounds to constrain the higher derivative corrections of holographic QCD models.


Contents
-ii -For long distances beyond a certain characteristic scale 1/M , low-energy effective field theories (EFTs) are utilized to describe physical processes and predict observables. In these scenarios, ultra-violet (UV) effects manifest as tails shaped by higher-dimensional operators, which are suppressed by 1/M . However, at short distances, UV effects intensify and become nonnegligible, raising a profound question: What is the allowable space of EFTs that ensure a consistent UV completion, such as in quantum gravity? Particularly regarding completion to quantum gravity, there are numerous intriguing conjectures and arguments, known as the swampland program [1][2][3][4]. Examples include the weak gravity conjecture [5][6][7], the distance conjecture [1,8], and others. These are inspired by string theory and studies of black hole physics, and they provide conceptual criteria for gaining insights into this question. Even without the gravitational degree of freedom, this question remains profound and warrants further investigation. It is known that EFTs without gravity can still be pathological. Having oversized Wilson coefficients [9] or possessing the wrong sign for some EFT Wilson coefficients [10] can violate causality.
The EFT bootstrap program has recently been developed to quantitatively and systematically explore this question, assuming the unitarity and causality of the underlying UV theory above M , as well as Regge boundedness [11][12][13][14][15][16]. The strategy involves studying 2to-2 scattering amplitudes in EFTs, denoted as M EFT , and then searching for the allowed space of Wilson coefficients. Causality and Regge boundedness provide a bridge between the EFT amplitudes M EFT and the underlying UV amplitudes, which is known as the dispersive sum rules. Built upon dispersive sum rules, the unitarity can then be used to optimally carve out the EFT space. This whole procedure is known as the dual bootstrap algorithm, as it rigorously rules out disallowed values of Wilson coefficients. The dual EFT bootstrap can also incorporate dynamical gravity [17][18][19][20], thereby providing sharp bounds on some of the swampland conjectures [21,22]. There are many relevant works utilizing this idea to constrain EFTs and their UV completions, see, e.g., .
However, there is a weakness in the current version of the dual EFT bootstrap. Typically, to optimize the EFT bounds, it is essential to ensure that we measure only a finite number of Wilson coefficients in which we are interested. Meanwhile, the null constraints should be employed [12][13][14][15][16]. These null constraints are constructed using crossing symmetry, a crucial ingredient of quantum causality. They complement the original dispersive sum rules because the latter are not fully crossing symmetric and thus lack important information. These requirements are usually met using the improved sum rules [17], which subtract the forward limit expansions from the original sum rules. This subtraction ensures they measure only those Wilson coefficients that saturate specific Regge boundedness. However, this procedure can be vulnerable to loop effects since the forward limit scale might compete with the loop expansions. This competition essentially hinders an efficient generalization that powerfully constrains EFTs at the loop level (for exploration on this subject, see, e.g., [35,49]). Even at tree-level, this forward limit subtraction can sometimes complicate numerical exploration.
In this paper, we will explore the numerical usage of the crossing symmetric dispersive sum rules [23,50,51], which automatically incorporate the crossing symmetry. In other words, the null constraints are encoded in the crossing symmetric sum rules, and these sum rules only retain a finite number of Wilson coefficients. We will show that this type of sum rule is an optimized version of the "improved sum rules", as it fulfills all the requirements that the "improved sum rules" satisfy and it is free of any forward limit subtractions. Although our discussions in this paper are limited to tree-level, we believe that the crossing symmetric sum rules are an excellent playground for understanding the loop effects of EFT bounds.
There is a different approach termed the primal S-matrix bootstrap [52]. The primal bootstrap is a powerful tool to constrain quantum field theories (QFT) non-perturbatively: it is built upon an appropriate ansatz of the S-matrix designed to obey causality and directly searches for optimal couplings under unitarity constraints. The term "primal" indicates that this method is used to rule in allowed values of couplings. Although the S-matrix bootstrap was proposed to constrain the dynamics of non-perturbative QFTs (with a large amount of applications, see, .e.g., [53][54][55][56][57][58][59][60][61][62]), it can be easily adapted for studying EFTs, e.g., [63][64][65]. Some natural questions then arise: how do the resulting EFT bounds from the primal method compare to those from the dual method? In what sense are the primal and dual bootstraps dual to each other in the context of optimization theory? Regarding the first question, intuitively, we expect that when both methods are applied to the same EFT with the same assumptions, their resulting bounds should converge to each other, eventually leading to the optimal constraints of EFTs. This convergence has indeed been observed in relevant investigations for scalar EFTs [64]. For the second question, it has been recognized that the primal bootstrap is the "primal problem" in optimization theory such as semidefinite program (SDP) [66,67], and therefore one can construct its optimization Lagrangian and identify the "dual problem". This "dual problem" should be the method one can construct using the dispersive sum rules as the dual EFT bootstrap [68].
In this paper, we will show that the EFT bootstrap can indeed be formulated as an infinite-dimensional SDP problem. This guarantees that the dual and primal bounds should converge to each other, as long as the strong duality is satisfied [69]. The strong duality condition can be translated to conditions of EFT bootstrap, which then indicates that we can extract UV physics from the primal solutions under the guidance of the dual extremal functionals.
We then apply the EFT bootstrap with the crossing symmetric sum rules, now formulated as an SDP, to the case study of large-N chiral perturbation theory (χPT). χPT is an EFT that describes light meson physics, arising from chiral symmetry breaking in the lowenergy regime of quantum chromodynamics (QCD). Large-N χPT is the meson EFT that emerges from large-N QCD [70], which generalizes the colour group to SU(N ) with N → ∞ [71], and provides a qualitative understanding of many aspects of hadron physics. The dual EFT bootstrap program for large-N χPT was initiated in [44], and it is still under active investigation [45,47,72]. One advantage of studying the large-N limit is that it retains only tree-level physics, where the dual EFT bootstrap is efficient. For finite N , such as in real QCD, χPT may become strongly coupled around the threshold and therefore the loop cannot be neglected. To our knowledge, in this case, only primal studies has been be performed [63,73]. On the other hand, the large-N limit is also expected to have a string description [71], so it might provide rich "experiments" for understanding quantum gravity and holography [74][75][76]. In this paper, we will study large-N χPT using both dual and primal methods, and we observe excellent convergence. We also extract the physical spectrum from our primal solutions and the dual functionals, not only confirming some understandings from [44,45], but also revealing some novel and hidden physics that seems to be accessible only when using the primal method.
Typically, for O(p 4 ) Wilson coefficients, a segment of the boundary corresponds precisely to the Skyrme model [44]. However, as the Wilson coefficients increase beyond the kink [44,45], the Skyrme model is excluded. On the other hand, large-N QCD is known to have string and holographic descriptions, such as the well-known Witten-Sakai-Sugimoto model [77][78][79], which yields precisely the Skyrme model at low energy [78,79]. It is intriguing to translate the constraints on large-N χPT to see if there are any problems with holographic QCD models. Not surprisingly, all known holographic QCD models produce the Skyrme model that exists below the kink, and is thus consistent. However, these models have all been analyzed only at the leading order as EFTs of gauge fields. We will show that including the higher dimensional operators on the gravity side leads to a large-N χPT that deviates from the Skyrme model, controlled by the bulk Wilson coefficients. This allows us to translate the large-N χPT bounds to constrain the bulk EFTs of gauge fields on non-trivial backgrounds.
The rest of the paper is summarized as follows. In section 2, we review the basic ideas of dispersive sum rules and provide the construction of crossing symmetric sum rules. After reviewing the structures of amplitudes and the unitarity constraints, we focus on the positive unitarity condition, explaining why the EFT bootstrap is an infinite-dimensional SDP and how we can construct its Lagrangian formulation and derive physics from SDP duality. In section 3, we review large-N χPT, including its Lagrangian, the flavour structure of the pion amplitudes, and the partial waves and unitarity conditions. We then explicitly construct the associated crossing symmetric dispersive sum rules and the primal S-matrix ansatz as our dual problem setup and the primal problem setup, respectively. In section 4, we obtain the EFT bounds using both the dual and primal algorithms and present the convergence between the two methods. We also display the spectral density and S-matrix, which are numerically solved using the primal methods for saturating certain EFT bounds. Using a simple sample bound, we demonstrate that modifying the Regge behaviour of the primal ansatz does not alter the resulting bounds, as long as it stays below the Regge boundedness assumption. As an ad hoc approach, we incorporate the upper bound of unitarity to uniformly bound the O(p 4 ) Wilson coefficients in terms of the cut-off scale M , which does not contradict the large-N bound, thereby confirming the consistency of the large-N limit. In section 5, we study holographic QCD. We include the higher dimensional operators built from gauge fields and show that holographic QCD with higher derivative terms produces the most general χPT Lagrangian at order O(p 4 ). We verify that the Witten-Sakai-Sugimoto model has no issues with the leading string corrections. Afterwards, we translate the chiral EFT bounds to constrain the 5D EFT with double gauge fields. We conclude the paper in section 6. Appendix A formulates other EFT bootstrap scenarios as SDP; appendix B records the SU(N f ) projectors that we used to organize the pion amplitudes and partial waves; and appendix C provides more details on the bootstrap Lagrangian for fixing one parameter and bounding another. One essential component of EFT bootstrap is the dispersive sum rules. The strategy is to design vanishing integral identities along a large circle at infinity in the complex s plane, e.g.,

EFT bootstrap and SDP problem
where k 0 is the ceiling of the Regge spin J 0 for UV amplitudes Causality is also assumed, which is thought to imply both the crossing symmetry and the analyticity of the S-matrix in the complex s plane, except for poles and branch cuts in the real axis (see [80,81] for more details on the analyticity of the S-matrix). Analyticity allows us to deform the contour of integral identities (2.1) towards the real axis, with a smaller arc within the regime where low-energy EFTs remain valid (|s| < M 2 ) 1 , as illustrated in Fig 1. This procedure establishes the dispersion relations (or dispersive sum rules) that relate low-energy to high-energy physics This gives where the discontinuity is defined by On the other hand, crossing symmetry allows for the relationship between the u-channel contribution and the s-channel contribution. Crossing symmetry enables us to modify and Figure 1: The contour deformation leads to the sum rules given by eq. (2.3). The red branch cut represents the UV branch cut, which is beyond our knowledge, while the blue branch cut represents the low-energy cut contributed by loop effects in low-energy EFTs. The final contour relates low-energy EFT data along the arcs to the discontinuity along the UV branch cuts.
improve the dispersive sum rules into a more convenient and powerful basis by subtracting the null constraints [12]. For instance, it is demonstrated in [17] that improved spin-k sum rules can be constructed by subtracting the forward-limit expansions of higher spin sum rules (i.e., k ′ > k). This method leaves only a finite number of Wilson coefficients with Regge spin k in the low-energy measurement, which is exceptionally beneficial 2 . For example, the improved spin-2 sum rule for four-dimensional scalar EFT is given by [17] B imp However, it is obvious that the construction of this improvement is heavily dependent on details; see [19,20] for more complicated scattering processes. In addition, the construction relies on the forward-limit expansion, making the loop effects vulnerable. In this note, we will 2 The Regge spin of a Wilson coefficient can be defined by the exponent of s in the fixed-t Regge limit of the associated tree-level amplitude in EFTs. For example, for a higher-dimensional operator giving amplitudes gs k in the fixed-t Regge limit, we say the Regge spin of g is k.
instead use the crossing symmetric sum rules [23,50,51], which we will introduce momentarily, to build causality (i.e., analyticity plus crossing symmetry) directly into the dispersive sum rules.

Crossing symmetric sum rules
We review the crossing symmetric sum rules in this section 4 . The essence lies in the analytic and crossing-symmetric parameterization of the Mandelstam variables in terms of a complex variable z and an auxiliary momentum p s(z, p) = − 3p 2 z 1 + z + z 2 , t(z, p) = s(z ξ, p) , u(z, p) = s(z ξ 2 , p) , (2.8) where ξ = e 2/3iπ and 0 < p 2 ≤ M 2 /3. In terms of the complex z-plane, the Mandelstam variables are geometrically symmetrical, with an angular difference of 2/3π from each other. The Regge limit in this parametrization is associated with special points on the unit circle. For example, the fixed-t Regge limit |s| → ∞ corresponds to z = ξ 2 ; other channels follow similarly. To build the crossing symmetric sum rules, it is necessary to study the full crossing symmetric amplitudes and find the crossing symmetric kernel that can perform the sufficient subtractions in the Regge limit. This kernel is easy to construct, from which we obtain the fixed-p identity Note that the integration contour consists of small circles surrounding the Regge limit point. The transformation of the integration variable to s yields Here, the Mandelstam variables t and u are parameterized in terms of s by In order to deform the contour in (2.9) and obtain the sum rule, we must understand the analyticity of M sym (z, p 2 ) in terms of the complex variable z. We have three pieces of UV branch cuts: s ≥ M 2 , u ≥ M 2 for fixed-t; s ≥ M 2 , t ≥ M 2 for fixed-u; and t ≥ M 2 , u ≥ M 2 for fixed-s. In the complex z plane, these branch cuts all reside on the unit circle and sandwich the Regge limit points in the associated channels. The UV branch cuts are summarized as follows Figure 2: The analytic structures and the contour deformation for crossing symmetric sum rules. The red branch cut represents the UV branch cut, and the blue branch cut represents the low-energy cut contributed by loop effects in low-energy EFTs.
We continue to have (2.3), and more specifically, it now yields It is worth noting that we express the UV part in terms of s, which will be convenient when performing the partial-wave expansion. The sum rule (2.13) is constructed to be crossing symmetric, and it naturally subtracts all null constraints in a nonlinear manner. Indeed, for example, the low-energy part for scalar EFT at tree-level is precisely the same as in (2.7), as observed by [42], but we are not taking any forward-limit! We argue that the crossing symmetric sum rule is a more natural and well-defined approach for addressing low-energy loops.

Manipulate sum rules using functionals
To harness the powerful capabilities of dispersive sum rules, it is instructive to build functionals that manipulate sum rules and measure the interesting couplings at low-energy (2.14) To derive constraints on EFTs, one can the search for functionals that optimize quantities at low-energy, subject to the unitarity that we will introduce later. Generally, we can define the functionals by smearing the sum rules against wave functions where p 2 max = M 2 for (2.4), while p 2 max = M 2 /3 for (2.13).
There are two kinds of functionals in the literature, which we call the forward-limit functional and the impact parameter functional [17]. The forward-limit function is achieved by taking the wave function ψ(p 2 ) = i c i ∂ i p 2 , which then performs the forward-limit expansion; on the other hand, the impact parameter functional measures the sum rules at small impact parameter b ∼ 1/M .
• Impact parameter functional ψ(p 2 ) has finite support in the momentum space and decays fast enough in the impact parameter space ψ(b) := d d−2 p e ib·p ψ(p 2 ); such functionals are usually chosen as 5 as well as its variants for numerical benefits [17,19,20].
The forward-limit functional is much simpler and requires fewer computational resources, therefore it is more often used in EFTs without graviton. However, this functional can be singular at low-energy when dealing with graviton propagation due to the 1/t graviton pole at low-energy. In contrast, the impact parameter functional would suppress the graviton pole, making the gravitational low-energy behaviour well-defined under its action 6 . In addition, the impact parameter functional also provides a bonus for allowing one to weaken the Regge boundedness assumption (2.2). The essential reason behind this bonus is that smearing amplitudes against the fast decay wave function ψ(b) would suppress the higher spin contributions at high energy, effectively enhancing the Regge behaviour under the smearing [19] (see also [82] for a more evident proof).
In this paper, we do not include the graviton, therefore we will use the forward-limit functional for quick convergence of numerics.

Low-energy amplitudes
It is worth noting that at this stage, the dispersive sum rules formally utilize the full amplitudes, even along the low-energy arc. For low-energy part, since this arc is inside the EFT regime, we may expect to replace the amplitudes there with EFT amplitudes. The simplest situation is that the underlying theory is weakly coupled at low energy, where we can replace 5 It is worth noting that one has to pay attention when choosing the starting point of the polynomial p i 0 , which controls the numerics in the large impact parameter regime b → ∞ [17]. 6 In 4D, the graviton pole is not completely resolved. Nevertheless, the divergence can be improved from polynomial divergence to the logarithmic IR divergence log M/mIR using the functionals integrated from m 2 IR [17,19]. This logarithmic IR divergence reflects the behaviour of the classical Newton potential. As a consequence, the functional becomes ineffective beyond bmax ∼ 1/mIR, a region that we should simply discard. low-energy amplitudes along the small arc with tree-level EFT amplitudes. These amplitudes contain only simple poles, allowing us to evaluate the arc integral by picking up the residues of simple poles. This is the simplest case that has been extensively studied. Generally, we have where the full amplitudes M(s, t) are expanded in a Taylor series in terms of 1/M . The EFT amplitudes are computed using the effective Lagrangian, which is dependent on the scale through logarithmic structures such as log(s/µ) and log(m 2 /µ). An additional matching piece often appears because the order of Taylor expansions and the integrals do not commute, i.e., The matching piece contains terms like log(M 2 /µ). For simplicity, we often choose µ = M 2 , which gives us a simpler relation Therefore, the dispersive sum rules measure the Wilson coefficients at scale µ = M 2 ; it is then necessary to apply the renormalization group equation to evolve Wilson coefficients back to other scales. In this paper, we will study the large-N chiral EFT, therefore the tree-level approximation is sufficient.

Unitarity constraints
How do we use (2.3) to constrain EFTs in terms of the UV amplitudes when the details of the UV theory are absent? It is instructive to study amplitudes at high energy using the partial wave expansion where ρ labels the irreducible representation of SO(d) and π ρ is the associated partial waves [20], and we slip off the indices of possible global symmetry 7 . The unitarity then implies a strong constraint on the partial wave coefficients This is the scenario that gives rise to the positivity bounds considered in the literature [11]. If the couplings of the underlying theory are weak enough that the quadratic terms |a ρ (s)| 2 ≪ Disc a ρ (s) can be ignored, then the positivity condition robustly constrains the low-energy EFTs. However, it is important to understand that a weakly coupled EFT does not necessarily mean its UV completion will also be weakly coupled. The essential condition is the existence of a parametrically small but positive parameter 0 < g ≪ 1 in low-energy EFTs that can be measured by dispersive sum rules. This leads us to where Y ρ (m 2 ) is any appropriate function which is generated by functionals acting on partial waves. This obviously shows that Disc a ρ has to be parametrically small as g. In this case, the optimal bounds on Wilson coefficients will be scaling with g. Gravitational EFTs studied in [19,20] (and their couplings to scalar and photon [21,22,84]) fall into this category, where the parametrically small parameter is the Newton Large-N chiral EFT that will be studied in the following sections also falls into this category, where the small parameter is the inverse of the pion decay constant 1/f 2 π ∼ 1/N ≪ 1 [70] (see section 3 for more details).
In other cases, the full unitarity (2.22) will be providing more stringent constraints on low-energy EFTs, such as scalar EFTs studied in [64]. It turns out that the full unitarity (2.22) can be linearized by formalizing it as a positive matrix [52] where Re a ρ (s) = 1/2 a ρ (s+i0)+a ρ (s−i0) . It is then easy to see that for small Disc a ρ (s) we only need to consider the first diagonal element of this matrix; on the other hand, if Disc a ρ (s) is not necessary small but Re a ρ (s) is small, we can also ignore the off-diagonal terms and impose the linear constraint 0 ≤ Disc a ρ (s) ≤ 2 (which is considered in, e.g., [48,64,85,86]). These discussions allow us to classify the scenarios of EFT bootstrap, which invoke different unitarity conditions according to the basic assumptions (we follow the terminology invented in [64])

III. Nonlinear unitarity
From now on, for simplicity, we will denote the discontinuity as the imaginary part, using the notation Disc → Im, when there is no confusion.

Semidefinite programming
Since the unitarity of the S-matrix can be expressed as a semidefinite matrix, carving out the allowed space of EFTs can then be transformed into the semi-definite program (SDP), subject to those unitarity constraints. In this section, we provide a crash course on SDP. We will then show that EFT bootstrap is an infinite-dimensional SDP.
The SDP can be formulated as the following primal optimization procedure [69,87] • Primal problem In this algorithm, S K is the space of K × K symmetric real matrices. The primal S-matrix bootstrap [52], as applied to EFTs, falls into this problem with infinite dimensions. For simplicity, we only consider the positivity constraint. We can approximate the full amplitude (which is valid for both UV and low-energy EFT) by using an infinite number of analytic but simpler functions M i with the assumed analyticity and Regge behaviour (2.26) The positivity of the unitary condition now becomes where a i ρ is the partial wave coefficients contributed by M i . If we imagine that we put Im a ρ (s) at all values of s ≥ M 2 for all spins into an infinite-dimensional diagonal matrix, the primal EFT bootstrap is, in principle, an infinite-dimensional primal problem with K, N = ∞ when taking C ≡ 0. 8 Here, B T x = b is simply the normalization condition for fixing a particular Wilson coefficient, and we can minimize the target Wilson coefficient by choosing appropriate c. This is because any Wilson coefficient can be represented in terms of a linear combination of x i by taking the low-energy limit of the full amplitudes. However, in practice, we cannot reach N, K = ∞ in numerics. Instead, one takes a maximal value of N max , and we also consider up to a certain J max , imposing unitarity for a finite but large number of s-grids [52]. This truncation procedure gives a well-defined optimization, after which one must extrapolate the bound [58,59,64].
Duality plays an important role in SDP. Typically, the primal problem has a dual formulation, known as the dual problem We can easily construct the dual version for primal EFT bootstrap with positivity that we described previously. The maximization target is b · y since we choose C ≡ 0, we then have where we have already explicitly evaluated the trace by integrating over s ≥ M 2 and summing over all spins in the irreducible representation. Note that in this language, we use Y ρ (s) to denote an infinite-dimensional matrix in which the diagonal elements take the values of all s ≥ M 2 and spins in ρ. What does (2.28) mean? It becomes clear if we dot (2.28) into x, we find In other words, the dual problem is to find a positive function Y ρ (s) such that its average against the spectral density can represent low-energy Wilson coefficients c · x, and we bound c · x by maximizing b · y according to the positivity. The dual problem (2.29) can obviously be achieved by using the dispersive sum rules, which is precisely the dual bootstrap algorithm studied in [12]. In this algorithm, the function Y ρ (s) can be constructed by acting the functionals on the UV part of sum rules Other scenarios can also be formulated as SDPs, see [67,68] and appendix A for further discussions. It is also natural to ask whether the primal and dual problems yield the same optimal value for our target. This question can be answered by the duality theory in SDP, which we will review in the following subsections.

The Lagrangian formulation, SDP duality and physical implications
SDP has the Lagrangian formulation, which manifests the logics of optimization. A SDP is described by the following Lagrangian where the first line is intended to manifest the primal problem, while the second line targets the dual problem. Equality is achieved using the expression X = N i=1 A i x i − C along with some basic algebra. It's worth noting that we don't use the subject identity, because we want to emphasize that every component of SDP can be seen from the Lagrangian. Using this Lagrangian, the primal problem can be formulated as [69] The dual problem can then be constructed by interchanging the ordering of "minimize" and "maximize", namely This indeed gives rise to the standard dual problem by noting According to our previous discussions, we can then immediately traslate (2.31) to the Lagrangian of positivity EFT bootstrap, built from any reasonable functionals acting on the dispersive sum rules (2.36) The first term is the dual objective, and the second term represents functionals acting on sum rules for our targets, giving rise to (2.14). This Lagrangian is compact and is the guide throughout this paper. To be more concrete, let's say we want to minimize a Wilson coefficient g F in terms of g 0 > 0. We can simply expand F • B(p 2 ) We then have We can now easily read off either the primal or dual algorithm. The primal problem is straightforward to read; we fix g 0 = 1 and then minimize g F subject to the positivity Im a ρ (s) ≥ 0. It is worth noting that setting g 0 = 1 is not physical; this is a scaling trick to deal with the numerics, and one can always set other values of g 0 . The physical result is the lower bound of g F /g 0 ; similarly, the solved Im a ρ is also in the unit of g 10 . For the dual problem, we This is precisely the dual algorithm proposed by [12]. A natural question arises: For the same bootstrap problem, do the primal and dual methods yield the same constraints? This question can be answered using the duality theory of SDP. In general, if we find feasible solutions for both the primal and dual problems, we have weak duality, which states that the primal bound is always greater than the dual bound, leading to a duality gap as follows This is the reason the primal bound is always referred to as the rule-in bound, while the dual bound is termed as the rule-out bound. To ensure the duality gap vanishes, we clearly need XY ≡ 0 and the Slater's condition X ≻ 0 or Y ≻ 0 [69]. This implies that a solution to the problem must satisfy Physically, this condition is powerful. For example, numerical conformal bootstrap employs Y ≡ 0 to identify the physical spectrum, a method known as the extremal functional method [88,89]. In terms of positivity EFT bootstrap, the second possible condition is trivial, it simply provides a free field theory solution. The nontrivial physical implication is then clearly The first condition makes physical sense, as the spectral density has to be nonzero for nontrivial physics. We can also use the second condition to locate the physical bound state or resonance above the cut-off [17,44]. It is worth noting, however, that this is the most ideal situation. As previously mentioned, in practice, it is not possible to treat EFT bootstrap as an infinite dimensional SDP. Therefore, practically, we expect the duality gap is not zero but will converge to zero as we increase the dimension of the problem. In this situation, the weak duality can serve as a double check criteria, helping us diagnose any numerical mistakes. Practical implementation of EFT bootstrap also complicates the task of solely using Y ρ (s) ≡ 0 to determine the physical points [44]. Nonetheless, we will demonstrate later that the first condition from (2.42) with sufficiently small Y ρ (s) can still be insightful for extracting physical information.
Other scenarios of EFT bootstrap can be formulated similarly, we keep the discussions in appendix A.

Large-N chiral perturbation theory
Starting from this section and in all subsequent sections, we will apply the EFT bootstrap to large-N χPT. We are adopting the set-up presented in [44], where the dual algorithm was employed.
Our innovations compared to [44,45] are twofold. For the dual algorithm, we use the crossing symmetric sum rules. As we previously indicated, these rules automatically incorporate all null constraints, making them more efficient and allowing us to easily explore higher dimensional operator; Additionally, we will establish the primal method and demonstrate the convergence between the primal and dual approaches.

Chiral Lagrangian and low-energy amplitudes
We consider the chiral limit of large-N QCD (with SU(N ) gauge group) [71], where the fermionic sector possesses U L (N f ) × U R (N f ) chiral symmetry. Usually, there is an axial anomaly which breaks the global symmetry by . However, the axial anomaly is suppressed by the large-N limit [90][91][92]. At low-energy, we then have the spontaneous symmetry breaking pattern This results in N 2 f − 1 pseudo-Goldstone bosons (which are massless in the chiral limit) in the adjoint representation of SU(N f ). There is also a singlet meson that can mix with the gluon. However, due to large-N understanding of the OZI rule, this mixing is suppressed by 1/N , making it the trivial U(1) part of U(N f ). For N f = 2, the adjoint representation contains pions π; for N f = 3, the adjoint representation contains pions π, kaons K and the eta η; while the singlet is referred to as the eta prime η ′ . Nevertheless, we will show momentarily (see also [44]) that the unitarity constraints are independent of N f at the strict large-N limit. We therefore follow [44] to refer to the Goldstone bosons to as the large-N pion.
We can formulate the large-N pion physics using the coset construction of and then construct a non-linear sigma model by parameterizing the symmetry breaking in terms of low-energy field where f π is the pion decay constant that scales as √ N in the large-N limit. Here Π a denotes the large-N pion. For example, for N f = 3 we have where I is a 3 × 3 identity matrix. Then the chiral Lagrangian describing the EFT can be constructed [93]. Up to p 4 at large-N limit, it is generally given by [94] L It is worth noting that we drop all sub-leading terms that are not single-trace because a flavor trace comes from a quark loop and thus also acquires a color trace [94,95]. This leaves us only two independent Wilson coefficients. Interestingly, for finite N but N f = 2, there are also just two independent Wilson coefficients up to p 4 ; while for finite N but N f = 3, there are three independent Wilson coefficients at this order. Since we only consider 2-to-2 pion scattering, we then simply drop all background gauge fields (see [47] for dual bootstrap including background gauge fields). At higher orders, it becomes quite challenging to enumerate a complete set of higher-dimensional operators without redundancy, especially when identities exist that allow trading one operator for another in N f = 2, 3. However, using global symmetry and Bose symmetry, one can easily write down tree-level amplitudes to any order in p.
It turns out to be useful to parametrize the amplitudes using the generators of U(N f ) , making the Bose symmetry manifest. At large-N limit, M(s, t) is trivially zero because it can only be contributed by nonplanar diagrams and is therefore suppressed (see, e.g., [96] for an explicit one-loop result).
Using the permutation symmetry, one can easily construct tree-level amplitudes at low-energy up to any orders in p where we have removed p 0 term as it is forbidden by the Adler's zero [97]. The low-lying identification with the Lagrangian is [44] (3.8)

Partial waves and unitarity
We now turn our attention to reviewing the partial waves of large-N pion scattering. This can be understood by considering the amplitude M ab , cd as a sum over the U(N f ) irreducible representations that label the Π a Π b → X three-point vertices. Here, X represents the intermediate states in the Π a Π b → Π c Π d scattering. Following the notation in [98], the representation theory behind this physical picture is given by 9 The multiplicity of 2 for the adjoint representation arises because the vertices can be either symmetric or anti-symmetric in two legs. The resulting adjoint representation labels meson states consisting of bilinear quarks (i.e., qq states) in the intermediate channel, while the other representations label different exotic resonances. In addition to these global symmetries, Π a behaves like a scalar, making the partial waves trivially correspond to the Legendre polynomials, as in scalar scattering. Therefore, one has the following s-channel partial wave decomposition [44] M ab cd (s|t, u) = Here, R denotes the irreducible representations in (3.9), and P R is the projector associated with R. The crucial difference from the pure scalar case is that the Bose symmetry also permutes the projector P R . The unitarity condition is then also straightforward The projector associated with the adjoint representation can be easily constructed because there are only two simple Π a Π b X c vertices: either the structure constant f abc or d abc := 2Tr T a T b , T c . The two associated projectors are therefore proportional to f abe f ecd and d abe d ecd . Other projectors are more intricate but can be constructed using the Casimir operators and the relevant eigenvalues [98]. We document all projectors in appendix B. From the permutation symmetry of the projectors, we can easily write down the spin selection rules for (3.12) 0, adj S , ss, aa : even J , adj A , as ⊕ sa : odd J . (3.13) To perform the EFT bootstrap, it is beneficial to explicitly know the relations between the generator basis (3.6) and the basis from irreducible representations (3.10). This can be readily achieved if we are aware of all the projectors, and we have (3.14) After we take M = 0 due to the large-N limit, we can reproduce the relations outlined in [44].
Using these relations, we can easily translate the unitarity condition where we have c su J = (−1) J c st J . This aids in constructing both the dual and primal problems. In the large-N limit, it suffices to use the positivity bootstrap (this will be justified in the next subsection). Additionally, all exotic mesons are suppressed, and therefore Im c ss ≡ 0. We then have [44] Im a st J (s) = Im a su J (s) = a 0 (s) , for even J , It is obvious that, in the large-N limit, the bootstrap constraints would be independent of the number of flavours. However, other bootstrap scenarios depend on N f .

Dual problem set-up
Let's set up the dual problem for large-N pion scattering. The most crucial component is the crossing symmetric dispersive sum rules. As we previously described, the construction of these sum rules is universal. The theory-dependent inputs involve constructing the crossing symmetric amplitudes and making assumptions about their Regge behaviour.
In QCD, the Regge intercept is at J 0 ≃ 0.52 for both M(s, t) and M(s, u) [99,100], therefore k 0 = 1. We follow [44] to assume that we can still take k 0 = 1 at large-N limit.
It turns out that one can construct three independent crossing symmetric amplitudes [30] M (1) = M(s, t) + M(t, u) + M(s, u) , To construct well-defined crossing symmetric sum rules from these amplitudes, we need to determine their Regge behaviors using the Regge boundedness of the building blocks M(s, t).
We find Thus, in terms of these symmetric amplitudes, the sum rules for M (2,3) are super-convergence sum rules. The complete set of sum rules is therefore where k = 1, 3, 5 · · · , denoting the Regge spin of the sum rules with respect to M(s, t). The low-lying low-energy contributions from these sum rules are As we noted, using these sum rules, we don't need to construct the null constraints as in [44,45,47]. At high energy, we have where x = (m 2 − 3p 2 )/(m 2 + p 2 ) 1/2 . The average is defined by We can now decide what bootstrap scenario that we should use. Let's simply look at B 1 , at leading order in p 2 → 0, we have This does not only prove the positivity of g 10 [44], but it also satisfies the condition of using only positivity, because g 10 ∼ 1/f 2 π ∼ 1/N ≪ 1. Therefore, we will focus on the positivity bootstrap, with the exception of subsection 4.4. In subsection 4.4, we will employ the linear unitarity bootstrap to verify that the large-N expansion is meaningful at the EFT level. We use SDPB [87,101] to implement the algorithm.

Primal problem set-up
Let's now focus on the setup of the primal problem. The fundamental building block of the primal problem is the S-matrix ansatz, which approximates the S-matrix [52]. To ensure that we are indeed constructing a primal problem that is dual to the previously described dual problem, this ansatz must satisfy the assumptions of analyticity and Regge boundedness. Consequently, it has to validate all sum rules.
To ensure analyticity, we follow the approach in [52] to define a function using Mandelstam variables This function obviously has branch cut starting at s = M 2 . Then the ansatz can be built by polynomials in (ρ s , ρ t , ρ u ) under the restrictions of Bose symmetry and momentum conservation s + t + u = 0. Let's now use ρ s,t,u to construct the ansatz for M(s, t), which is symmetric in s and t. It is crucial not to overcount or miss any terms in the ansatz. Since the primary problem is to determine the optimal coefficients of the ansatz terms, any redundancy or omission could either produce unfaithful bounds or simply disrupt the numerical calculations. We start with listing the generators of our ansatz that are symmetric in (s, t) (3.25) We do not need to consider ρ a u because the large-N limit suppresses the u-channel cut of M(s, t), as previously reviewed. Typically, the relation s + t + u ≡ 0 constrains the number of independent polynomials at each order, starting from order 5, which need to be subtracted [52]. In our case, since there is no ρ u available in M(s, t) in the large-N limit, the polynomials of the form (s + t + u) × (· · · ) do not exist. Hence, polynomials constructed using the above generators are independent. We can easily write down the ansatz 10 where the lowest order is 1 rather than 0 due to the Adler's zero when expanding in s, t ≪ M 2 . The overall function R(s, t, u) is for controlling the Regge behaviour of the amplitudes. Because k 0 = 1, for simplicity, we choose We can modify R(s, t) to exhibit a more refined Regge behaviour by specifying J 0 . However, since the only ingredient needed for constructing the solutions of the positivity EFT bootstrap comes from the sum rules, and these sum rules are sensitive to k 0 rather than J 0 , it's natural to speculate that modifying R(s, t) won't alter the resulting bounds as long as it grows at a rate below s k 0 =1 at high energy. We will verify this point in subsection 4.3. By expanding the ansatz (3.26) in the low-energy limit where s, t, u ≪ M 2 , we can derive the low-energy tree-level amplitudes and establish a dictionary that translates Wilson coefficients to α ab . For example, the dictionary for low-lying coefficients is g 10 = 1 4 α 10 , g 20 = 1 16 (2α 01 + α 02 ) , g 30 = 1 64 (5α 01 + 4α 02 + α 03 ) , g 21 = 1 32 (2α 02 + α 10 ) .

(3.28)
We can then easily verify that every term in (3.26) satisfies the crossing symmetric sum rules. The primal problem involves imposing the positivity condition (3.16) on the ansatz (3.26) and solving for the coefficients α ab by optimizing the targeted Wilson coefficients using the dictionary like (3.28). The last technical question is how we read off the partial wave coefficients from our ansatz (3.26)? We use the standard inversion formula 11 It is crucial to note that, although we only consider the imaginary part in the positivity SDP Lagrangian, we can still solve the full S-matrix from the optimal primal solutions. While 10 It is worthing noting that this ansatz actually has maximal analyticity, which is stronger than the analytic assumptions of the dispersive sum rules. Nevertheless, as we show below, the dual and primal bounds converge, seemingly suggesting that the SDP set-up of the positivity EFT bootstrap does not use the maximal analyticity.
We are grateful to Miguel Correia for the discussions on this point. 11 This formula may explain why the primal bootstrap typically does not employ maximal analyticity, thus converging to the dual one with weaker analyticity. For s ≥ M 2 , this formula solely relies on the physical regime −M 2 ≤ t < 0, making the resulting data sensitive only to the analyticity for t < 0. We, therefore, propose examining the subtlety of 'maximal analyticity vs. partial analyticity' using the gravitational EFT. In this context, unitarity must be demanded beyond integer spin, e.g., for J ∼ b √ s at high energy [17,19].
it may seem that the primal method always provides more information than the dual, this perception is mistaken. On the dual side, one can also rely on the extreme functional to employ the analytic "rule-in" method, which enables the construction of a relevant UV theory [12,16,44]. We use SDPB [87,101] to implement the algorithm.
Using the crossing symmetric sum rules, it is easy to obtain the dual bounds, we have (3.23) as well as The primal bounds for g 10 , g 20 > 0 are also trivial to obtain, where the solutions are all α ab ≡ 0, since g 10 = 0 or g 20 = 0 would a trivial free theory. This is consistent with the Slater 's condition and the complementary condition previously reviewed: the dual functional Y J (s) is strictly positive, therefore the strong duality gives Im a J (s) ≡ 0.
The first nontrivial example is g 21 , since its dual functional Y J (s) can be zero for even spins, suggesting that nontrivial UV amplitudes with only even spin particles exist. It turns out that g 21 ≥ 0 converges trivially for a low N max = 5 and a low J max , which we choose to be J max = 60. A nontrivial S-matrix profile with J = 0 that saturates g 21 = 0 can be illustrated, as shown in Fig 3. We observed that the spectral density for all higher spins is zero. This preliminary study thus confirms the statement from [44] that the UV theory atg 21 = 0 is a scalar theory. However, we observe from Fig 3 that the UV scalar spectral density doesn't show an extreme peak at certain points; instead, it presents a continuum. This suggests that the UV theory is not a single scalar but a scalar theory with all possible mass values where m ≥ M . To bootstrap the upper bounds, we simply flip the overall sign of g 20 , g 21 in the Lagrangians (4.1). The upper bound of g 20 in the unit of g 10 is also trivialized by the dual method The upper bound ofg ′ 2 = 2g 21 /g 10 M 2 , although it is not straightforward, it can still be easily solved from the dual algorithm using SDPB [87]. Using the crossing symmetric sum rules, we reproduced the result of [44] These two bounds are nontrivial from the primal side, since low spin sampling and low N max would give us trash, which does not extrapolate well to an infinite dimensional SDP. For a more involved S-matrix bootstrap, which involves either linear unitarity or even complete unitarity, the strategy is to fix N max and then increase J max so that one can extrapolate the bounds to be valid for all J; subsequently, one should vary N max and extrapolate the bounds to N max = ∞ [58,59]. However, for the positivity primal bootstrap, we find that we can simply fix J max to a large value without doing the extrapolation. We choose J max = 60, and we can see the nice convergence of bounds by varying N max from 5 to 25, as shown in Fig 4. When N max takes a small value, the approximation of the positivity EFT SDP is not good. However, we still expect the weak duality to be valid. This is precisely why we see that the primal upper bounds are always smaller than the dual upper bounds. Ultimately, we find that N max = 25 is enough to conclude the strong duality, as the relative error of primal bounds from the dual bounds is roughly ∼ 1%.  Figure 4: We determined the upper primal bounds ofg 2 andg ′ 2 by varying N max . These bounds converge quickly to the dual bounds, which are represented by red lines, approaching from below as guaranteed by the strong duality of SDP. Now we can combine the dual and primal functionals to analyze the physical spectrum that saturates bounds. For primal side, we simply use the solutions from N max = 25. The crucial point to understand is that since we are still far from the actual infinite-dimensional SDP, we cannot rely exclusively on either the dual functional or the primal solution to extract physical information. The strategy is as follows: initially, examine the dual functional. If the dual functional is precisely zero at a particular point, then we should trust the spectral density from the primal solution at that point, irrespective of its magnitude size. Conversely, if the dual functional is strictly positive and large, we would expect the corresponding primal "spectrum" to be small. Ideally, this primal spectrum should be vanishingly small. If it's not, it should be small enough to be interpreted as a numerical artifact, and we should simply discard it. The most subtle situation arises when the dual functional is strictly positive and small enough to be approximated as zero. In this case, we should estimate the gap Y J (s)Im a J (s) at that point. If the gap is sufficiently small, we can trust the primal spectrum; otherwise, we discard the data. However, this is also difficult to implement. It is worth noting that the dual functional is usually a polynomial of M 2 /s, and is small for sufficiently high s numerically. It is challenging to numerically detect that a small number is a zero or it is simply suppressed by 1/s. This suggests that the EFT bootstrap does not have sharp implications in the deep UV, but a vague picture of the physics there can still be captured by the primal solutions: for large s, we believe that it is reasonable to trust the primal spectrum density, because we can always treat Y ρ (s) there as zero with small errors.
From the dual functional (4.3), we see that Y J (s) can be zero only when s = M 2 and it is strictly positive for s > M 2 , which is robust against adding more functionals. We indeed observe from the primal solution that there is a single peak around s = M 2 for J = 0, see Fig 5; while Im a J for J ≥ 1 is vanishing. Besides, we also checked that the other Wilson coefficients at this point areg ′ 2 = 0 andg 3 = g 30 /g 10 M 4 ≃ 0.976 ∼ 1. This analysis confirms that the UV theory withg 2 = 1 corresponds to a single scalar theory with mass m = M , as first pointed out by [44]. The relevant scalar mode with s 0 Regge behaviour is A comparison of this scalar amplitudes with our numerical solution is illustrated in Fig 7a. Nevertheless, it is important to note that the plot in Fig. 7a is drawn for the global region, which significantly suppresses the differences between the analytic and numerical amplitudes. The largest difference between amplitudes obtained by the two methods occurs at |s| → ∞ and t → 0, and is approximately 0.018.  Figure 5: The spectral density at J = 0 that saturatesg 2 = 0. To generate this plot, we used a polynomial order in S-matrix of N max = 25.
The associated dual functional is complicated, but we can nevertheless easily observe that it is strictly positive for J = 0 but it is zero for J > 0, s = M 2 . From the primal side, we indeed observe that Im a J=0 , although not exactly zero, is parametrically small as of order 10 −5 ; in addition, J > 0 spectral density exhibits a sharp pump around s = M 2 , which is, however, getting smaller and smaller for larger J. See Fig 6 for an illustration with J = 0, 1, 2, 3. In addition, we find that this point givesg 2 ≃ 0.99 ∼ 1 andg 3 ≃ 0.97 ∼ 1. This analysis confirms the statement of [44] that the UV theory withg ′ 2 ≃ 3.26 is a theory with J ≥ 1 and m = M . However, it is important to note that our "numerical theory" is radically different from the su-model which also saturatesg ′ 2 ∼ 3.26 [44], because the su-model has Regge behaviour s −1 while our ansatz grows like s 0 . We can modify the su-model to describe  Figure 6: The spectral density from J = 0 to J = 3 is solved using N max = 25 atg ′ 2 ≃ 3.26. The J = 0 spectral density is nonphysical, as it is of the order 10 −5 , which is parametrically small compared to others. For higher J values, we have smaller spectral densities, serving as a reminder of the low-spin dominance. a spin-1 theory with Regge behaviour s 0 (4.6) We can then compare our numerical rule-in with this analytic rule-in in Fig. 7b, where the maximal difference is around 0.017

The Skyrme bound and a mysterious Regge trajectory
There is an interesting linear bound, giving rise to precisely the Skyrme model [102,103], as first noticed in [44]. This bound involves g 21 and g 20 , and we can formulate it as The bound is again trivial on the dual side, we have According to the previous experience, N max = 25 is good enough for us to perform a nice primal algorithm. We then find that the primal bound is The error from the dual rigorous bound is around 1.54%. Let's now turn to analyze the physical spectrum of the Skyrme model. We observe that the simple dual functional Y J (s) can be zero only for odd J, which seems to suggest that we should discard all data with even spin in the primal solution. However, this conclusion is not robust against expanding the space of the dual functionals. When using 126 functionals, we can find that the dual functional is small at J = 0, s = M 2 (around 10 −3 ), while it is strictly positive and not small close to s = M 2 for other J. With this in mind, we can examine the primal solution, and we find that the contributions from J ̸ = 1 are relatively smaller than those from J = 1. This behaviour can be described as the vector meson dominance [104]. Specifically, for J = 1 we identify a sharp peak around m = M , which can be interpreted as a vector ρ meson. There is also significant physics at higher energies: a continuum with a resonant bump around m ≈ 7.4M and a width of roughly 12.6M , which might be a numerical artifact when compared to the m = M peak. For other spins, we should only trust the behaviour at sufficiently high s = m 2 , and we indeed observe dominate resonances at energies m > 7M with a relatively wide width. It is important to note that these resonances are stable upon increasing N max . Interestingly, if we consider M as the mass of the ρ meson, approximately 770MeV, then the mass of all those heavy resonances exceeds 4000MeV and thus would contain, e.g., a bottom quark. See Fig 8 for an explicit illustration for J = 0, 1, 2, 3. We can extend our primal analysis up to J = 10 and find subsequent resonances. More surprisingly, these resonances can be organized as an approximately linear Regge trajectory 12 , as shown in Fig 9. For higher J, we observe a significant deviation from the fitted Regge trajectory. This deviation is likely due to poor numerical shooting. Unfortunately, we have no clear explanation for this trajectory. It is possible that we should not take it from such high energy behaviour of primal solutions, and as inferred by [44], those resonances are actually pushed to infinity in the large-N limit. Refining the numerics to better understand this Regge trajectory would be interesting in the future. Now we can compare the numerical amplitude to the model proposed in [44] M (UV) where it reduces to a single ρ meson model in the limit m ∞ → ∞. Our strategy of comparison is to first solve m ∞ for requiring the Regge limit of this analytic model equals the numerical amplitude at fixed t and then compare the amplitudes with other energy. We take t = −1/10 and find that m ∞ t=−1/10 ≃ 14.7M , the comparison is displayed in Fig 10. We observe that although the extremely high energy limit is required to be the same, two amplitudes become clearly distinguishable at s ∼ 15M 2 . The difference at large s but below |s| → ∞ is anticipated, because (4.10) is just a toy model with a single resonance m ∞ to adjust, while the primal solutions

O(p 4 )
So far, we have only dealt with simple linear bounds. The power of the EFT bootstrap lies in its ability to search for allowed spaces that involve multiple Wilson coefficients. These represent nonlinear bounds, implying that the boundary is not a linear function. Let's focus on the space spanned byg ′ 2 andg 2 and reproduce the exclusion plot made using dual methods in [44].
Our strategy is to search in different directions in theg ′ 2 −g 2 plane. The corresponding Lagrangian is where c ∈ [0, 1). On the dual side, this corresponds to fixing the objective and varying the normalization condition; while on the primal side, this precisely means fixing the normalization to g 10 = 1 and bounding different linear combinations of g 20 and g 21 . There is another approach, different from these angle-searching methods. In this approach, one can fix a particular value of one parameter, for example, g 20 /g 10 , and search for the upper and lower bounds of another parameter g 21 /g 10 [12]. While this method is typically as efficient as the previous one on the dual side, it complicates the search on the primal side. This is because it fixes two parameters, g 20 and g 10 , making the positive matrices degenerate. Consequently, an additional transformation is needed to generate a valid positivity input. Nevertheless, near the corner, it is aways better to adopt the "fixing-parameter" method rather than the "angle-searching". We leave the details of "fixing-parameter" methods in appendix C. By sampling a sufficient number of c values (roughly 100 points) and few "fixing-parameter" points near the corner, we can make a sufficiently nice exclusion plot using both dual and primal methods (where we choose N max = 25). The plot is displayed in Fig 11, where the dashed black boundary is drawn using the dual method, coinciding with the findings in [44]; on the other hand, the solid boundary is derived from the primal method. We find that the primal bounds efficiently converge to the dual bounds. To generate the dual bounds, we use 126 functionals that are constructed from the crossing symmetric sum rules B 1 , B 3 , B 5 13 .
We observe from Fig 11 that the Skyrme line represents only a small segment of the entire boundary. The point where the bounds begin to deviate from the Skyrme model is referred to as a "kink" in [44]. Nontrivial physics is anticipated at this kink [44,106]. Further studies suggest that the kink is located at the point (g 2 ,g ′ 2 ) = (1/3, 4/3) [45], which is ruled-in by the analytic model (4.10) with m ∞ → ∞. We present several points on the boundary in Fig 11, retaining three digits after the decimal, as shown in Table 1. We note that aroundg 2 ≃ 0.26, the value ofg ′ 2 is already slightly smaller than what the Skyrme bound predicts, which is, Figure 11: The exclusion plot involvesg 2 andg ′ 2 , and compares the primal bounds to the dual bounds. The red line represents the linear Skyrme bound and the black dot is the position of (1/3, 4/3) kink.
however, likely to be numerical error. Therefore, using our dual numerical results does not provide sensible way to clearly pinpoint the position of the kink.  Table 1: Few boundary points with smallg 2 , obtained using the dual methods.

O(p 6 )
For O(p 6 ), we consider a similar Lagrangian but for g 30 , g 31 L p 6 = cos(2πc)g 30 + sin(2πc)g 31 + λ(g 10 Following the same strategy as noted previously, we can create the exclusion plot forg 3 = g 30 M 4 /g 10 andg ′ 3 = g 31 M 4 /g 10 . There is a good convergence between the primal and dual methods, as seen in Fig 12. Some comments are in order. Interestingly, we find that the linear upper bound ofg 3 andg ′ 3 is the same as ing 2 andg ′ 2 . This is because that the corner can still be analytically ruled in by (4.6). and we also have the Skyrme-like linear boundg ′ 3 ≤ 4g 3 . However, it is obvious from Fig 12 that the whole allowed region is much smaller than the region enclosed by the linear bounds, even though there is no clear kink. In this case, the nonlinearity largely shrinks the allowed space of EFT. The same phenomenon was also observed in gravitational EFT [19,26,86].

Is the positivity primal bootstrap sensitive to the Regge behaviour?
In this subsection, we aim to address the question of whether the Regge behaviour in the primal ansatz affects the primal bounds. The answer should be "No" for the positivity bootstrap, as long as the Regge behaviour is below k 0 = 1 ensuring that our dual set-up remains valid. This is confirmed by the duality between the dual and primal methods: the dual bounds are only sensitive to k 0 , which provides all sum rules. Therefore, the primal bounds should not depend on the specifics of the Regge behaviour of the amplitudes, as long as their growth rate is below k 0 . However, if the Regge behaviour of the primal ansatz reaches k 0 or exceeds it, it cannot measure and provide bounds on the Wilson coefficients with Regge spin k 0 .
Our strategy is to choose R(s, t) in the ansatz (4.1) as (4.13) For J 0 = 0, we recover the case that we have been studying. In general, this factor modifies the dictionary that relates the ansatz parameters α ab to the low-energy Wilson coefficients, but it preserves the analyticity and grows as s J 0 in the Regge limit. We will study the primal upper bound ofg ′ 2 , which can be measured by B 1 sum rules and is sufficiently nontrivial. We will vary J 0 by taking several values: J 0 = (0, 0.5, 0.75, 1, 1.1) and observe how the primal bounds change accordingly. This exploration is displayed in Fig 13 below. In Fig 13, we present only the results for J 0 = (0, 0.5, 0.75), and it's clear that they all converge to the dual rigorous bound at sufficiently large N max . For J 0 = 1, 1.1, as expected, we do not obtain any bounds.

An ad hoc: primally confirming the large-N assumption
So far, we have been considering the constraints of the large-N χPT, as suggested in [44]. As previously noted, we assume the large-N limit from IR to UV; therefore, the positivity bootstrap is sufficient and strongly constraining. In this way, we can bound the Wilson coefficients in terms of g 10 = 1/(2f 2 π ) ∼ 1/N . Nevertheless, it is essential to question whether this assumption is justified from a low-energy point of view. In other words, do the large-N bounds fall into the allowed regime of a more complete unitary region like linear unitarity?
This question may seem trivial, since there is no doubt that the spectral density scaling as 1/N falls into the linear unitarity Im a J ∼ 1/N ≪ 2. Indeed, this simple argument trivializes the linear bound: the large-N limit yields bounds that scale in g 10 , i.e., g ≤ O(1)g 10 /M dim−2 , which are significantly stronger than the bounds provided by linear unitarity, g ≤ O(1)/M dim , because g 10 ≪ 1/M 2 . However, it is important to emphasize that this straightforward argument and power counting do not obviously work for nonlinear bounds. The dimensional analysis can only infer that the large-N exclusion plot resides in the near-zero corner in the linear unitarity plot, however, its boundary may bend outside of the region of the linear unitarity exclusion plot. This concern arises due to results of [64,85], where it turns out that the boundary of the linear unitarity bounds is more curved and is sandwiched between the linear bounds.
Our strategy to justify the large-N bound involves using the linear unitarity bootstrap for pion amplitudes, whose structures remain constrained by the large-N limit. We employ only the primal method. For the dual method that incorporates the upper bound of unitarity, see [12,48,85,86]. For a recent systematic numerical algorithm on the dual side, refer to [48] 14 . For simplicity, we take N f = 2, then (3.16) indicates that the unitarity constraints are (4.14) Follow appendix A, we then consider the following SDP Lagrangian In this case, we physically "normalize" the free theory part S = 1 + iT . Although in this case, the numerical convergence is more slow, we find that N max = 25 and J max = 60 still suffice for our purpose. By searching different values of c, we found Fig 14, where we define g ′ 2 = 2g 21 M 4 and g 2 = g 20 M 4 . From Fig 14, we can confirm that the linear bounds are indeed consistent with the large-N bounds as long as N ∼ O (10). Nevertheless, we do observe dangerous regions. The first dangerous region is around g 2 ∼ 2.7, where there is a sharp kink, and the boundary shrinks away from the vertical line. In contrast, the large-N bounds in Fig 11 do not exhibit this shrinking. To resolve this danger, we need to require The second danger is the Skyrme line. In this plot, the lower boundary behaves similarly to the large-N plot: it first coincides with the Skyrme bound and then bends inward. Therefore, the position of the "Skyrme kink" in Fig 14 must be larger than the one in the large-N case. Fortunately, the kink in Fig 11 is around g ′ 2 ∼ 11.6, for which the condition (4.16) easily resolves the danger. Since (4.16) can be easily satisfied in the large-N limit where f 2 π ∼ N → ∞, we thus confirm that assuming the large-N limit at low energy is valid without paradox. Besides, we emphasize that this exercise also shows that the unitarity can be used to bound the decay constant in terms of the EFT scale M . Indeed, we can set an EFT bootstrap to directly constrain the decay constant f π which deserves further exploration for finite N χPT

Constrain holographic QCD models
Large-N QCD enjoys holographic descriptions, which utilize the dual gravity theory to capture the salient properties of QCD, such as hadrons, in the strong coupling limit. One famous example of these models is known as the Witten-Sakai-Sugimoto model [77][78][79], which can be constructed from string theory. Such models can also be built from low-energy EFTs of gauge fields with gravitational couplings [107][108][109][110][111][112][113], where the fits to the experimental data of hadrons, glueballs, and so on have been extensively studied (see e.g., [114,115] for brief reviews). Typically, at low energy, the holographic models can also give rise to the χPT Lagrangian with Wilson coefficients mapping to parameters on the gravity side [78,79,[108][109][110][111][112][113]. Therefore, we expect that the bounds on large-N χPT can be translated into constraints on those holographic QCD models, carving out the allowed space of EFTs for gauge theories that can be consistently UV completed with gravity. However, to our knowledge, all holographic QCD models so far only include Tr F 2 term when deriving the χPT, giving rise to the Skryme model at order O(p 4 ). The essential reason is that higher derivative terms are relatively small, like suppressed by the string scale [78,79]. Therefore, one can always adjust the fundamental parameters so that the Wilson coefficients live on the boundary of the exclusion plot 11 below the kink. Nevertheless, although the higher derivative terms only give small corrections, these small corrections may still deform the Wilson coefficients outside of the allowed region in Fig 11. In this section, we will show that higher derivative terms Tr F 3 and Tr F 4 cause the low-energy theory deviate from the Skyrme model. Therefore, requiring the consistency with Fig 11 puts constraints on the higher derivative couplings.

Bulk theory and the power counting
We now move to derive the chiral Lagrangian from 5D EFT of SU L (N f ) × SU R (N f ) gauge fields. We consider the following effective action with the background field g where All A, B, · · · refer to the five dimensional indices, and all indices so far are contracted by the background field g AB and a background dilaton φ where we put an UV Randall-Sundrum (RS) bane [116,117] at z UV and an IR RS brane at z IR . For Anti de-Sitter (AdS) space, the UV brane is served as the boundary of AdS. Nevertheless, throughout this subsection, we consider a(z), b(z) and φ to be arbitrary. Their precise forms satisfy the equations of motion for both the gravity sector and the matter sectors (such as dilaton associated with φ [78,79,118]). However, we treat the flavour gauge fields as probes [78,79]. Appropriate boundary conditions are imposed on both the UV brane at z UV → 0 and the IR brane z IR to obtain the background solution (5.3). We ignore fluctuations from the graviton and other matters, as they are irrelevant for our purposes.
How should we think about the EFT power counting of (5.1) without a top-down picture like string theory? The bulk gauge fields are expected to correspond to conserved currents in QCD. Hence, we then have J µa L and J µa R , which can be constructed by quark bilinear operators J µa = qγ µ T a q; their conservation precisely reflects the global chiral symmetry. From large-N counting, we know that we usually scale such quark bilinear operators by 1/ √ N to normalize the two-point function [71]; this corresponds to scaling A by g YM to normalize the kinematic term, suggesting g YM ∼ 1/ √ N in (5.1). Moreover, the Wilson coefficients g H , α i should be suppressed by some EFT cut-off Λ, which could be the mass of higher spin particles or the string states (it's important to distinguish it from the χPT EFT cut-off M for the moment). However, due to the presence of g YM , there are different possible consistent schemes of power counting for g H and α i , as long as their sizes don't grow beyond ⟨F ⟩ # . For instance, we can have g H ∼ 1/Λ 2 , α i ∼ 1/Λ 4 ; or we can have g H ∼ g 2 YM /Λ, α i ∼ g 2 YM /Λ 2 , both of which are then suppressed by the large-N limit. The idea is then to use the EFT constraints to decide the correct size of those Wilson coefficients. It's important to note that we actually ignore some double trace terms like Tr F 2 2 , because Tr maps to the trace in χPT, and coefficients of double-trace operators in large-N χPT are 1/N suppressed, for general N f . Therefore, we conclude for N f > 3, without doing anything, that the Wilson coefficient of 1/g 2 YM Tr F 2 2 must scale at least as g 2 YM /Λ 2 ! Indeed, in the Witten-Sakai-Sugimoto model from the D 4 -D 8 -D 8 intersection, the effective action is the dilaton-Born-Infeld (DBI) action [78], which still has only a single trace when generalized to the non-Abelian case [119].

Routine to obtain the chiral Lagrangian
We follow the strategy of [112] that uses the IR boundary condition to break the chiral symmetry This implies that the chiral symmetry breaks to its vector subgroup at IR. For convenience, we regroup the gauge group by its vector component and axial component It turns out that one can define a Wilson line stretch from the IR point into the bulk which has the property of the pion field under the gauge transformation U → g R U g −1 L , and thus it corresponds to U in χPT. For simplicity, we choose the gauge A Rz ≡ 0, and use U to gauge away A Lz using the gauge fixing prescription of [78,111,112]. This procedure allows us to define the following boundary condition To obtain an effective action, one can solve the bulk gauge fields with respect to these boundary conditions, and then substitute the solutions back to have the on-shell action. For simplicity, in this paper, we only focus on those terms up to O(p 4 ), therefore we can simply solve the equation of motion at leading order [112] where · · · refer to those terms contributing to higher orders like O(p 6 ), and one has It is then straightforward to obtain the chiral Lagrangian (3.5), where the pion decay constant f π and other Wilson coefficients are given by 3 , . (5.10)

Witten-Sakai-Sugimoto model is healthy
The Witten-Sakai-Sugimoto model, constructed from the D 4 -D 8 -D 8 brane configuration in type IIA string theory, is considered a top-down model and is expected to be robust. Therefore, before imposing constraints on more general bottom-up models like (5.1), we aim to verify the health of the Witten-Sakai-Sugimoto model from a low-energy perspective. The essential idea is to start with the D 8 brane embedded in the D 4 configuration Since we have the complete picture from the string theory, we can easily keep track of the power counting (where we keep the leading KK tower mass to be M KK = 1 for simplicity) [78,79] α ′ = ℓ 2 s = 9 2λ , g s = λ λ is the 't Hooft coupling, where g YM,c is the Yang-Mills coupling for the colour sector on D 4 brane. The weakly-coupled supergravity regime is only valid for λ → ∞, N → ∞. To map this model to our bottom-up EFT (5.1), we cut the brane in half for z ∈ (−∞, 0) and introduce an additional gauge field to compensate for the contribution from the other half, where z UV = −∞ and z IR = 0. We find (1 + z 2 ) 5/12 , g 2 YM = 486π 3 λN . (5.14) 15 Our convention of gauge field is different from [78,79]: Besides, expanding in α ′ yields [119] g H = 0 , α 1 = −π 2 (α ′ ) 2 , α 2 = 4π 2 (α ′ ) 2 .
(5.15) Using the dictionary (5.10), we find At the leading order when λ → ∞, we reproduce the results of [78,79]. We can immediately see that the string correction causes the Wilson coefficients to deviate from the Skyrme model. To compare with the large-N χPT bound, we should examineg ′ 2 andg 2 At leading order, this is constrained to be below the kink for M 2 ≤ 0.68M 2 KK . Indeed, the KK spectroscopy analysis suggests that the ρ mass is M 2 ρ = 0.67M 2 KK ! The string correction then pushesg ′ 2 andg 2 upwards from the boundary, ensuring they still fall within the allowed region of Fig 11. We thus conclude that, even when including the leading string correction, we do not identify problems with the Witten-Sakai-Sugimoto model.

Flat and AdS hard wall models
Now we move to constrain two known holographic QCD models with φ = 1. We focus on two models, one is constructed in the flat space as RS scenario [107], another is the hard wall model [111] constructed in AdS. They are both clearly explained in [111]. We believe that our discussions can be generalized to other holographic QCD models (with possible modifications on how the chiral symmetry is breaking), like the soft wall models, where the dilaton φ is nontrivially turned out [118,[120][121][122].
• Flat space scenario [107] For this model, we consider where we fix the UV brane to be z UV = 0, and z IR ∼ 1/Λ QCD . We have where we have used M 2 ∼ M 2 ρ = π 2 /(4z 2 IR ) [107]. Besides, we denote 2α 1 + α 2 = α, which is the unique combination appears. We can easily observe some simple linear boundsg A more complete exclusion plot is depicted in Fig 15a, which looks like a thin river.
• Hard wall in AdS [111] For this model, we consider We obtain where we have used M 2 ∼ M 2 ρ ∼ 5.78/z 2 IR [111]. Interestingly, in general we have two parameters z IR and R AdS , but the resulting bounds suggest that 1/R AdS rather than 1/z IR is the cut-off for bulk EFT. We have similar simple bounds The exclusion plot Fig 15b also shows a thin river.  Figure 15: (a) The bounds on holographic QCD supported by the model in flat space [107] (b) The bounds on holographic QCD from hard wall model in AdS [111]. All plots are still extending like a thin river.
It is important to note that the dictionary provides bounds uniformly scaled by z IR or R AdS , both of which are the IR scales in QCD. From naive dimensional analysis, we expect them to be bounded by the bulk"string" scale Λ. This either means that this method yields bounds that are too weak, or that the IR RS branes disrupt the naive dimensional analysis. On the other hand, we found that even though we identify R AdS with z IR , we can't reproduce the flat-space scenario from the AdS one, which suggests a breakdown of the flat-space limit of causality bounds [18]. It would be interesting to explore all these points in the future for a better understanding.

Summary
We reviewed the EFT bootstrap, especially for the positivity scenario. We demonstrated that the EFT bootstrap is essentially an infinite-dimensional SDP, where the optimization Lagrangian can be formulated. We built the dual problem of the EFT bootstrap using the crossing symmetric dispersive sum rules, which embrace the crossing symmetry without the IR danger, and thus it serves as a better version of the improved sum rules. For the primal problem of the EFT bootstrap, we adapted the S-matrix primal ansatz and optimized the target Wilson coefficients.
We then applied the EFT bootstrap program to large-N χPT, which is a low-energy pion EFT from the chiral symmetry breaking of large-N QCD. Due to the large-N limit, the positivity EFT bootstrap is sufficiently strong to carve out the allowed EFT space. Our dual bounds match with earlier literature [44,45], and we demonstrated that the primal bounds are also converging to the dual rigorous bounds. This is consistent with the strong duality of SDP. We then focused on some converged bounds and used the primal solutions to extract the physical spectrum and S-matrix that saturate those bounds. By doing this, we confirmed some of the analytic rule-in amplitudes studied in [44,45]. Interestingly, for the Skyrme bound, we also observed a mysterious heavy Regge trajectory, which seems to suggest meta-stable exotic states with heavy quarks. In addition, we showed that the Regge behaviour of the primal ansatz does not affect the bounds, if it stays below the assumed Regge boundedness. This is consistent with SDP, as the dual problem is only sensitive to the Regge boundedness rather than the explicit Regge behaviour. Eventually, we incorporated the upper bound of the unitarity, i.e., the linear unitarity EFT bootstrap to confirm that the large-N limit is consistent.
In the end, we considered the holographic QCD models, which are EFTs of gauge fields in 5D and correspond to large-N QCD in 4D. Typically, we included the higher derivative terms in the bulk and showed that they give rise to the general chiral Lagrangian up to O(p 4 ). We demonstrated that the Witten-Sakai-Sugimoto model with string corrections gives rise to a large-N χPT within the allowed EFT region. Besides, for bottom-up models like the flatspace RS model and AdS hard-wall model, we translated the large-N χPT to constrain the higher derivative couplings of Tr F 3 and Tr F 4 terms.
There are several aspects that deserve further investigations. From formal aspect, it would be interesting to build more precise relation between different methods for EFT and S-matrix bootstrap, include the SDP we reviewed, the moment problem [16,85,86], geometric function [28], iterative algorithm [123,124] and machine learning approach [125]. Focusing on the SDP perspective, typically, we state that the crossing symmetric sum rules are the more natural tools to understand the loop effects on the EFT bounds, since they are free of forward-limit issues. It is then interesting to make this statement concrete by using the crossing symmetric sum rules to study the scalar EFT, χPT and other EFTs, with one and even two loop effects, trying to make the results of [35,49,126] sharp. For this exploration, it is also important to understand the primal-dual convergence when there are loops and the nonlinear unitarity is utilized. For example, we can apply the nonlinear unitarity bootstrap to real χPT like [63,73], and if the primal bounds and dual bounds converge, we can then probably extract the real QCD physics in the UV. Such analysis may also be extended to other important EFTs, like the standard model EFT [48,[127][128][129][130], gravitational EFT [19,20,43,86,131,132], QCD string EFT [55], etc., helping us gain more information about their low-energy space as well as their possible UV completions.
Particularly for large-N χPT, it remains puzzling to us that the Skyrme model is problematic above the kink, since the Skyrme model is a good phenomenological model for understanding many aspects of nuclear physics, e.g., [133,134]. It would be interesting to understand, microscopically, how the Skyrme model goes wrong above the kink. A possible route is to study the pion-nuclei scattering, which can be described as pion fluctuations around the Skyrmion [135,136] (which are solitons of the Skyrme model and serve as the baryon [137,138]), and to detect if there are any causality violations like time advance [9]. Besides, it is also interesting to understand where our constraints on holographic QCD models come from in the bulk. The constraints may again arise from classical causality, and techniques from [131,132,[139][140][141] would then be useful. This investigation may also be generalized to other RS scenarios, which provide the standard model EFTs.
In this appendix, we formulate the linear unitarity bootstrap and the nonlinear unitarity bootstrap as SDP.

A.1 Linear unitarity
For linear unitarity bootstrap, we write the Lagrangian as follows

A.2 Nonlinear unitarity
The nonlinear unitarity bootstrap is more subtle. From primal side, we have further requirement for real part of a ρ (s), however, such object does not appear in our dispersive sum rule. Nevertheless, one can shift the dispersion relation to finite |s 0 | > M , we then have, for example for spin-2 sum rule B 2 (s 0 , p 2 ) = ds 4πi 2s + t (s 0 − s) 2 (s 0 + s + t) 2 = 0 → Dual: Unfortunately, a more concrete dual example of this type of bootstrap is beyond the scope of this paper. We refer the readers to relevant discussions in [29]. It would be interesting to explicitly realize the dual algorithm we propose here and compare with other dual algorithm in the future [29].

B Projectors of irreducible representation in SU(N f )
In this appendix, we record all the projectors of irreducible representation in SU(N f ) [98] that we used to organize the pion amplitudes. (B.1)

C Fixing-parameter method
In this appendix, we explain the "fixing-parameter" method for nonlinearly bounding two Wilson coefficients, as a complementary of "angle-searching" method that was described in the main text. The Lagrangian of "fixing-parameter method" is The interpretation is that we fix g 2 = g * 2 and bound g 1 in the unit of g 0 . From dual method, this corresponds to having which is precisely the fixing-parameter dual method used in [36]. On the primal side, however, the implementation is a bit subtle. g 0 = 1 is the normalization condition in the primal algorithm, then how should we address g 2 = g * 2 ? Recall that any Wilson coefficients are linear combinations of the coefficients in primal ansatz that we aim to solve, this indicates that g 2 = g * 2 put more constraints on the primal ansatz. Effectively, the primal ansatz is then degenerate and the coefficients there are no longer all independent. The strategy is to make an matrix R to reduce the anstaz to independent subspace so that R · g 2 = g * 2 R · g 0 , (C.3) where we understand g 2 and g 0 as vectors spanned by the primal ansatz. Thus we can play with the effective Lagrangian