The R-matrix bootstrap for the 2d O(N) bosonic model with a boundary

The S-matrix bootstrap is extended to a 1+1d theory with $O(N)$ symmetry and a boundary in what we call the R-matrix bootstrap since the quantity of interest is the reflection matrix (R-matrix). Given a bulk S-matrix, the space of allowed R-matrices is an infinite dimensional convex space from which we plot a two dimensional section given by a convex domain on a 2d plane. In certain cases, at the boundary of the domain, we find vertices corresponding to integrable R-matrices with no free parameters. In other cases, when there is a one-parameter family of integrable R-matrices, the whole boundary represents integrable theories. We also consider R-matrices which are analytic in an extended region beyond the physical cuts, thus forbidding poles (resonances) in that region. In certain models, this drastically reduces the allowed space of R-matrices leading to new vertices that again correspond to integrable theories. We also work out the dual problem, in particular in the case of extended analyticity, the dual function has cuts on the physical line whenever unitarity is saturated. For the periodic Yang-Baxter solution that has zero transmission, we computed the R-matrix initially using the bootstrap and then derived its previously unknown analytic form.


The S-matrix and R-matrix bootstrap programs
New insights were recently found on the old idea [1] of determining the Smatrix directly from its analytic structure, symmetries, crossing, and unitarity. This certainly works in two dimensional integrable theories but only after using the factorization constraint, namely the Yang-Baxter equation. Without that, those constraints are not enough to completely determine the S-matrix. However, recently it was found that maximizing the coupling between particles and their bound states led to well-known theories such as a subsector of the sine-Gordon model. It can be also applied to 3+1 dimensional theories, and multiple amplitudes [2,3,4,5]. The main physical argument is that, when the spectrum of bound states is fixed, there is a limit on the value of the coupling since stronger couplings will lead to more bound states. This is a very powerful idea, namely that certain theories lay at particular points of the space of allowed theories (or S-matrices) and that those particular points can be found by maximizing certain functionals in that space. For this paper, the case of interest is the 2d O(N ) nonlinear sigma model studied in [6], exactly solved in [7] and more recently revisited with the S-matrix bootstrap approach in [8,9,10]. In particular, in [8] it was argued that maximizing a linear functional in a convex space generically leads to vertices of the space. In fact, it was shown that the O(N ) non-linear sigma model (NLSM) lies at one of those vertices, the functional just being a way to find it. Later this was made more manifest in [11] where a section of the space was plotted with a clear vertex at the NLSM. However, it was also found that other theories did not appear to be at vertices. Further work on other models [12] showed that sometimes full regions of the boundary correspond to interesting theories if such theories have free parameters. Various other ideas have been discussed in the context of the S-matrix bootstrap and similar methods applied to gapped theories [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30].
Motivated by this, we consider the two dimensional bosonic O(N ) model with a boundary [31,32,33,34] and compute the reflection matrix (Rmatrix). The procedure is similar. We impose all constraints of analyticity, crossing, and unitarity and map the allowed space of R-matrices looking for special points at the boundary of the space. A difference is that crossing depends on the bulk S-matrix that we have to specify initially. It might be interesting to find the S-matrix and R-matrix simultaneously, but the constraints are non-linear. In fact, it seems more straightforward to compute first the S-matrix and then the R-matrix. In any case, this is the procedure we follow here and find the same variety of phenomena previously discussed for the S-matrix. When an integrable R-matrix exists with no free parameters it usually appears at a vertex of the allowed space. If this is not the case, we apply a new procedure that we call extended analyticity where we extend the analytic properties of the functions beyond the physical cuts. This severely restricts the allowed space by eliminating the R-matrices that have poles in that extended region. This is similar to removing R-matrices with bound states, but in this case, we can say that we remove resonances 1 . We find that R-matrices that were not at vertices now appear at the vertices of the restricted space. In other cases, there is a one parameter family of integrable R-matrices. In that case, we find that all the boundary corresponds to integrable R-matrices. We also find vertices that do not seem to correspond to any known theory.
The paper is organized as follows: In the rest of this section, we describe the properties of the S-matrix and R-matrix and exact results that follow from integrability with several examples. One result is new, we obtain an integrable reflection matrix for the periodic Yang-Baxter solution with no transmission (pYB). In the following section, we forget the requirement of integrability and just map out the allowed space of R-matrices from generic constraints. There we find the aforementioned properties. In the next section, we discuss the dual problem and find some useful properties of the problem with extended analyticity. In particular, we argue that it has to be regularized and that now, unitarity saturation does not follow automatically. In spite of that, all the R-matrices we found numerically actually saturate unitarity. In the last section, we give our conclusions.

The 2d O(N) bosonic model, general properties and exact S-matrices
Consider a two dimensional theory with O(N ) symmetry with N-species of bosonic particles with equal mass m labeled by a = 1 . . . N and two particle scattering (p 1 , a) + (p 2 , b) → (p 3 , c) + (p 4 , d) given by a generic S-matrix of the form where S T (s), S R (s) and S A (s) represent the transmission, reflection and annihilation amplitudes. Equivalently we can write The functions S I (s) and S ± (s) represent the scattering amplitudes in the three isospin channels: isoscalar, symmetric and antisymmetric. It is convenient to use the relative rapidity variable θ defined through If there are no bound states, the functions S a (θ) are analytic in the physical strip defined by 0 ≤ Imθ ≤ π. On the real line, θ ∈ R unitarity implies S I (θ) ≤ 1, S − (θ) ≤ 1, and S + (θ) ≤ 1. Finally crossing implies In this formula, the indices A, B, . . . label isospin channels and the matrix C satisfies C 2 = 1 and is given by (1.7) Up to now we have only described general constraints on the S-matrix due to standard properties of the field theory. If the theory is integrable the S-matrix should also satisfy the Yang-Baxter equation The simplest integrable model is just a free theory whose reflection matrix we briefly discuss in appendix C. In the main body of this paper we consider two non-trivial integrable models. One is the 2d O(N ) non-linear sigma model solved by Zamolodchikov and Zamolodchikov [7] with S-matrices given by and we introduced the function F a (θ) (see [11]): The other is the periodic Yang-Baxter solution (pYB) given by [35,9] S pY B (1.14) . (1.17) The function H ν is studied in more detail in appendix A, including an easier to evaluate definition and various properties needed to check unitarity, crossing, and the YB equation. It will appear again in the computation of the R-matrix. In [8,9,10] the S-matrix bootstrap procedure was applied to the O(N ) symmetric case with no bound states. In [8] it was observed that the NLSM appears at a vertex of the boundary space, a fact made manifest in [11]. In [9] the pYB solution was found and seen to correspond to some earlier work [35].

The 2d O(N) bosonic model, reflection matrices
Quantum field theories on a half-line show up in many areas of theoretical physics like D-branes in string theory and defects in condensed matter systems. The main object of interest is the reflection matrix giving the reflection amplitude for a given initial incoming state to be reflected into an outgoing one (see fig.1). Since the boundary breaks Lorentz invariance, the reflection matrix is a natural function of the incoming energy which, for a single incoming particle can be parameterized in terms of the rapidity θ as The incoming momentum is k = −m sinh θ. If the outgoing state is a single particle state it will have the same energy but opposite momentum. Taking into account that the theory has N species of particles of mass m labeled by indices a = 1 . . . N the 1 → 1 R-matrix can be formally defined as Allowing for particle production, unitarity requires that (1.20) The variable θ can be analytically continued to the strip 0 ≤ Imθ ≤ π 2 (see fig.2) or equivalently Re ε ≥ 0. Poles on the real axis of ε in the segment 0 ≤ Re ε ≤ m correspond to bound states of the particle to the boundary. We are going to assume no such bound states exist. For ε ≥ m there is a cut corresponding to intermediate single particle states. Further along the real axis we can have multi-particle cuts as in the case of the S-matrix. A double Wick rotation relates the imaginary axis of the energy plane to the real axis of a situation where the boundary is at t = 0 and therefore becomes an initial state. The spatial axis would now run from −∞ to +∞ allowing the definition of the usual S-matrix. The 1 → 1 reflection process is replaced by particle pair production from the initial state. Using this idea, Ghoshal and Zamolodchikov [36] showed that the R-matrix satisfies the following crossing equation for an integrable bulk S-matrix, where S ab cd (θ) is the (flavor part) bulk S-matrix in eq.(1.2) for the model under consideration. We are still allowing for particle production in the reflection t x t x bdy Figure 1: Pictorial description of the 1 → 1 R-matrix describing the amplitude for a particle to bounce from the wall possibly changing its identity a → b. A double Wick rotation relates this process to pair production from an initial boundary state B⟩. In the second case we can define a bulk S-matrix using the usual asymptotic states.
It might seem that in the first case analyticity is guaranteed by crossing but this is not so because of the nontrivial crossing equation (1.21).
process. Even when no particles are produced in a reflection, the R-matrix generically breaks integrability. The final condition for integrability is the reflection Yang-Baxter equation [36]: (1.22) This boundary Yang-Baxter equation along with unitarity and crossing constraints can be used to find exact R-matrices [31,33]. As we see below, these solutions lie at distinguished points of the convex space of allowed R-matrices. We emphasize that we do not assume integrability or particle number conservation in the numerical R-matrix bootstrap procedure we use later in the paper. Let us now consider a few important integrable cases.

NLSM, diagonal R-matrix
In the case of the NLSM, there are various forms of the R-matrix that are known to be integrable [33]. Here we consider only two. The first is a diagonal R-matrix of the form where R 0 (θ) is an auxiliary function, F a (θ) was defined in eq.(1.13) and λ in eq.(1.12). The cases k = 0 and k = N − 1 were studied in [31] and correspond to Neumann and Dirichlet boundary conditions for all O(N ) indices. The intermediate values of k were studied in [33] and correspond to mixed boundary conditions.

NLSM, block diagonal R-matrix
When N is even, another possibility is an R-matrix of the form [33] 27) namely N 2 equal anti-symmetric blocks along the diagonal. The symmetry is U ( N 2 ). The functions A(θ) and B(θ) are real analytic. It is not our purpose to consider all integrable cases, see [34,33] for a nice summary, but this particular one has a free parameter α and leads to a case where the boundary of the space corresponds to integrable solutions. The functions A(θ), B(θ) are given by (1.28) (1.29)

pYB, diagonal R-matrix
Now we can consider the pYB solution 2 . To our knowledge, the reflection matrix was not previously worked out in this case. Interestingly, we found it first numerically using the bootstrap and then we used the YB equation to get the analytic answer. We only consider the case of an R-matrix of the form . Using the Yang-Baxter equation and the explicit form of the S-matrix for the pYB solution leads to The constant θ 0 ∈ iR defines a one-parameter family of integrable R-matrices. It is easy to see that, given a point θ 1 on the imaginary axis of the physical region, the ratio R 1 (θ 1 ) R 2 (θ 1 ) takes every possible value as we change θ 0 ∈ iR. To obtain the R 1,2 (θ) functions themselves we use the ratio to write the crossing identity in the form . (1.33) We have to solve equation (1.32) together with the condition that R 1 (θ) is real analytic and saturates unitarity on the real axis R 1 (θ)R 1 (−θ) = 1. It is useful to notice that if, for a moment, we ignore S Y B + (2θ − iπ) the general solution to eq.(1.32), up to CDD factors [38], iŝ . (1.34) In fact S Y B + (2θ − iπ) can be absorbed in h(θ) and then compute R 1 (θ) with the result (up to the CDD factors described below) with χ 12 (θ) as in eq.(1.31), the function H ν as defined in appendix A, and It is also straightforward to use the identities in appendix A to check that this function solves eq.(1.32). One can also check that the reflection matrix has the same periodicity as the S-matrix: θ → θ + 2π 2 ν . However this is not necessarily a minimal solution. We consider a minimal solution one with no poles on the physical region 0 ≤ Imθ ≤ π 2 and also with no common zeros of R 1,2 in that region. Whether there are such zeros or poles depends on the value of the parameter θ 0 . With a bit of effort one can see that, if we define θ 0 = 2πiξ 0 and a periodic CDD factor 3 as then we need to multiply R 1,2 (θ) by one or both of the following CDD factors is the largest integer smaller or equal to x. For a given value of θ 1 we can plot the curve (R 1 (θ 1 , θ 0 = iξ 0 ), R 2 (θ 1 , θ 0 = iξ 0 )) parameterized by ξ 0 ∈ R. That curve determines the boundary of the allowed values of R 1 (θ 1 ), R 2 (θ 1 ) and is depicted in fig.9 for N = 6 and k = 1, 2, 3. It has vertices where a new CDD factor comes in, namely for ξ 0 = 0 and ξ 0 = − α 2ν .

Numerics
For the numerical bootstrap, we need to parameterize analytic functions on the strip b 1 ≤ Imθ ≤ b 2 for some real b 1,2 . In principle, we only need the physical region corresponding to b 1 = 0, b 2 = π 2 but it also turns out to be useful to consider other values of b 1,2 . In addition, we impose an extra periodicity along the real axis θ ≡ θ + 2ω for some ω that we choose. This periodicity is introduced to facilitate the numerics, however, some R-matrices (pYB) are actually periodic. Therefore, we are considering analytic functions on a cylinder that can be parameterized in terms of Fourier coefficients or in terms of the values of the real part at the boundary on a set of equally spaced points. We describe more details in appendix B. The numerical procedure gives excellent results on the imaginary axis and, on the real axis on a region roughly −ω 2 ≲ Reθ ≲ ω 2 since beyond that, boundary effects can be important. In the plots, we plot this central region where the solution is smooth. Although initially, we make all the plots for functions analytic in the physical region b 1 = 0, b 2 = π 2 , we found that the allowed space of R-matrices drastically reduces if we extend the region of analyticity. Also, an integrable theory which was originally at a regular point of the boundary, now appears at a vertex of the reduced region. This also shows that all the excluded points corresponded to R-matrices with at least one pole in the extended region. This is an interesting way to determine where different R-matrices have poles.

NLSM, diagonal R-matrix
Here we consider the bulk S-matrix to be given by the NLSM and assume the ansatz in eq.(1.23): (2.1) We then choose a point θ 1 = iξ 1 , 0 < ξ 1 < π 2 and plot the boundary of the two dimensional region of allowed values of (R 1 (θ 1 ), R 2 (θ 1 )). Notice that they are real by real analyticity. The results are plotted in fig.4 for the case N = 6, k = 1. The largest region corresponds to the allowed values of (R 1 (θ 1 ), R 2 (θ 1 )) when imposing analyticity in the physical region. An integrable solution appears at a vertex indicated by the red circle corresponding to Dirichlet boundary conditions for one of the fields. The analytic form of the functions R 1,2 (θ) is given in eq.(1.26) and is shown in fig.5 to agree perfectly with the numerics. Another integrable solution exists for which R 1 (θ) = R 2 (θ), the one corresponding to Neumann boundary conditions and given also by eq.(1.26) with k = 0. It is at the center of the purple circle but it is not at a vertex of the curve. There is a clear vertex on the upper left part of the curve although we were not able to identify it with a known theory. To make the other integrable vertex manifest we find now the allowed region for R-matrices analytic up to 0.9π. This drastically reduces the allowed space to the region surrounded by the orange curve. The integrable solution now is clearly at a vertex!. Further increasing the region of analyticity to 1.1π gives the green curve that does not seem to display any new vertices. As can be seen in fig.6 the same phenomenon happens when taking N = 6 and k = 2 and k = 3.

NLSM, block diagonal R-matrix
Here we consider the bulk S-matrix to be given by the NLSM and assume the ansatz in eq.(1.27) We then choose a point θ 1 = iξ 1 , 0 < ξ 1 < π 2 and plot the boundary of the two dimensional region of allowed (real) values of (A(θ 1 ), B(θ 1 )). Numerically we can verify that the relation B(θ) = −iαθA(θ) is satisfied all around the boundary where the constant α depends on the boundary point. In fact we can plot the boundary curve using the analytic solution in perfect agreement with the numerical results, see fig.7. Moreover, we can check various points at the boundary and see that the R-matrix itself agrees as exemplified in fig.8 Figure 4: Plot of the (R 1 (θ 1 ), R 2 (θ 1 )), θ 1 = 0.2i, allowed region for the NLSM with N = 6, k = 1 and various analytic regions b 1 = 0, b 2 = π 2 , 0.9π, and b 1 = −0.1π, b 2 = 1.1π. The vertex at the center of the red circle (see rotated and rescaled inset) is the R-matrix in eq.(1.26) with N = 6, k = 1. A new vertex appears in the smaller regions, namely the one at the center of the purple circle that corresponds to a solution with R 1 = R 2 , or N = 6, k = 0.   Figure 6: Plot of the (R 1 (θ 1 ), R 2 (θ 1 )) allowed region for the NLSM and θ 1 = 0.9879i. We consider the case N = 6, and k = 2, 3 and analyticity regions b 1 = 0, b 2 = π 2 , π.

pYB, diagonal R-matrix
Here we consider the bulk S-matrix to be given by the pYB and assume the ansatz in eq.(1.30) We then choose a point θ 1 = iξ 1 , 0 < ξ 1 < π 2 and plot the boundary of the two dimensional region of allowed (real) values of (R 1 (θ 1 ), R 2 (θ 1 )). Again, there is a one-parameter family of analytic R-matrices and we can plot (see fig.9) the allowed region analytically in perfect agreement with the numerical one. For various points at the boundary there is also good agreement between the numerical results and the exact ones eqs.(1.35), (1.40) as seen in fig.10. To find agreement it is important to include the CDD factor of eq.(1.41). The required CDD factor depends of the point at the bounday. In this case, the vertices of the shape correspond to the points where the CDD factor changes, namely ξ 0 = 0, − α 2ν . In the case N = 6, k = 3 we have α = 0 and therefore there is only one vertex. That the CDD factor changes means that, at a vertex, a pole (or a zero) is moving into the physical region.   Figure 9: Plot of the (R 1 (θ 1 ), R 2 (θ 1 )) allowed region for the pYB model for θ 1 = 0.9879i. We consider the case N = 6, and k = 1, 2, 3. In this case, the whole boundary is integrable and we can compare the numerics (dots) with the analytical result (solid lines). The curves intersect at R 1 (θ) = R 2 (θ) common to all k in the ansatz (2.3). The vertices (encircled) appear when a pole or zero enters that physical region and has to be cancelled by changing the CDD factor.  There is perfect agreement between the R-matrix bootstrap and the analytic solutions. We take N = 6, k = 1, 2, 3, and θ 0 = i π 4 , −0.329i, 0. 20 In this section, we discuss the dual problem. For the physical region, it works along the lines found for the S-matrix. When considering a region of extended analyticity some new results are found that are also applicable to the S-matrix. We do not provide numerical results since they all agree with the primal problem.

Dual problem in the physical region
The allowed space of R-matrices can be identified by maximizing various linear functionals while imposing crossing and unitarity constraints over linear combinations of a basis of functions. As we increase the basis, we cover a larger region and approach the boundary from the inside. We can also formulate a dual problem where we exclude points that do not satisfy the constraints by using a set of testing functions that partially impose the constraints. As we increase the set of functions, more constraints are imposed and we contract the region ending up with the same allowed region [11]. Let the functional we wish to maximize in the primal problem be for some constant real coefficients n b a and a given point θ 1 usually taken on the imaginary axis in which case taking real part is redundant due to real analyticity. Consider now a set of dual analytic functions K b a (θ) with a simple pole at θ = θ 1 and residue n b a and perform the following integral along the real axis: Now we move the contour of integration from the real axis to the line Imθ = π 2 . Taking into account that we cross the pole we get If we now impose anti-crossing condition for K j i (θ) then we can change variable θ → −θ in the integral and use crossing (1.21) to get Adding eqs.(3.3) and (3.5) we get F = F P , namely, Putting a bound on the right hand side of eq.(3.6) is easy if R b a and then K b a are diagonal as in eq.(1.23) since we just use that the diagonal elements of R have modulus less than one. If we want to be more generic we can use that 4 and that the unitarity constraint can be written as Taking into account that, if A and B are positive semi-defininte then TrAB ≥ 0 we get where k a (θ), a = 1 . . . N are the positive square root of the eigenvalues of K ab (θ) = K c a (θ)K c b (θ). This defines the dual functional F D and the dual problem is the minimization of F D subject to the anti-crossing constraint (3.4). As mentioned, the simplest case is when R b a (θ) is diagonal in which case we take K b a to also be diagonal and then k a is given by the diagonal elements as k a (θ) = K a a (θ) (no sum over a). If it is not diagonal then it is easier to introduce two hermitian matrices Y 1,2 (θ) and define the dual problem as subject to the constraint for all θ ∈ R. Here K denotes the matrix K b a . When the duality gap closes we have (at each θ ∈ R) From here we get 14) The last equality implies that unitarity is saturated RR † = 1 or, otherwise, K is not invertible, i.e. KK † has at least one zero eigenvalue. If we assume that K is invertible then unitarity is saturated and we can compute (using that Y 1 is hermitian) The inverse square root is well defined since K † K is a positive definite matrix under the assumption that it is invertible. We verified the previously obtained numerical results using the dual formulation and the duality gap is indeed zero as expected.

Dual to the extended analyticity problem
The dual problem for extended analyticity is quite interesting, and it allows for unitarity non saturation. Extending the analyticity region for the primal problem leads to more freedom in the dual functions and thus a smaller allowed region. Once again, we rewrite the primal functional F P = Re [n b a R b a (θ 1 )] as an integral along the real axis imposing analyticity except for a pole at θ = θ 1 and also impose anti-crossing for K b a (θ), Now, we introduce a new functionK b a (θ) that is analytic in a region bounded by the real axis and a curve C below it, see fig.11. If we require the R-matrix R b a (θ) to also be analytic in that region we have the identity allowing us to rewrite the primal functional as where, similar as before, ∆k a are the positive square root of the eigenvalues of ∆K∆K † with ∆K b a = K b a −K b a . Unfortunately we do not have any bound on the last term in eq. (3.18). In fact, to have a well defined dual problem we need a regularization, namely we have to bound the values of R b a (θ) when θ is on the curve C. We can do that by using a bound similar to the unitarity bound RR † ⪯ M , (3.19) where an identity matrix is assumed on the right hand side. With this in mind we write the dual functional as where θ(s) parameterizes the curve C andk a (θ) are the positive square root of the eigenvalues ofK(θ(s))K † (θ(s)). Again K(θ(s)) is the matrixK b a (θ(s)) on the curve C. The dual problem can be conveniently described by taking K andK as a single analytic function with a cut on the real axis. Then the first term in eq. (3.20) is an integral that measures the jump across the cut. When the duality gap closes we have the same analysis as before except that we should replace K → ∆K. If ∆K is invertible then unitarity is saturated. This � � � 0 0 C physical region extended region Figure 11: Extended region of analyticity. We generically extend the region below the real axis (across the physical cut) to a curve C. A well defined dual problem requires a regularization, namely a bound on the analytic function on C. An example of such curve is in fig.3 in which case the curve C is Imθ = b 1 . There we also extended the region to Im θ = b 2 which can be done in the same way.
is more easily expressed for a diagonal R-matrix and corresponding diagonal ∆K. In that case, a diagonal element R a a of the R-matrix does not have to saturate unitarity if the corresponding ∆K a a = 0. Then, in that region of the real axis, the functions K andK are analytic continuations of each other.

Conclusions
We presented an extension of the S-matrix bootstrap to the case of reflection matrices. We successfully reproduced known integrable R-matrices and found a new one which we then derived analytically. This shows that there is a new rich and interesting playground to test and develop S-matrix bootstrap ideas. In particular, a novel idea was the extended analyticity constraint. By requiring analyticity beyond the physical region, we contracted the allowed space of R-matrices, leading to a new vertex corresponding to an integrable model. This is a promising way to identify theories that can be readily applied to the S-matrix (work in progress). As regards the R-matrix bootstrap, we only explored a few possibilities and many other theories can be investigated including supersymmetric ones, theories with bound states, etc. 25 We are very grateful to Lucía Gómez Córdova, Yifei He, Nima Lashkari, and Pedro Vieira, for comments and discussions. We are also very grateful to the DOE that supported in part this work through grants DE-SC0007884, DE-SC0019202 and the QuantiSED Fermilab consortium, as well as to the Keck Foundation that also provided partial support for this work.

A A useful function
Given four complex numbers α, β, γ, δ with α + β = γ + δ and a positive real number ν ∈ R >0 we define 5 where (a, q) ∞ = ∏ ∞ j=0 (1 − aq j ) is the q-Pochhammer symbol. The function H ν can also be written as with the understanding that the infinite products are computed as in the previous equation. This is an analytic function of α, β, γ, δ on the whole complex plane except for poles whenever γ = inπ − 4νj or δ = inπ − 4νj for some n ∈ Z, j ∈ Z ≥0 . It also has zeros whenever α = inπ − 4νj or β = inπ − 4νj, n ∈ Z, j ∈ Z ≥0 . For values of α, β, δ, γ other than the poles, convergence 6 is 5 The equivalence of both definitions can be shown by checking that both have the same zeros and poles and their (quasi)-periodicity properties listed below, e.g. behavior under α → α + 4ν, etc. 6 This can be used to accelerate the convergence of the product of Γ functions, namely expanding the log in inverse powers of n and summing ∑ manifest if we take logs and expand the logs for large n (in the case of Γ using B Discussion of numerical approach Consider functions Φ a (θ) analytic in the strip 7 0 ≤ Imθ ≤ b and periodic along the real axis Φ a (θ + 2ω) = Φ a (θ). In this case, they are the nonzero components of the R-matrix. Generically, the R-matrix is not periodic, the periodicity is imposed to facilitate the numerics by cutting off the energy range. For numerical purposes, we parameterized (a subset) of such functions in two different, equivalent ways -(i) by the real parts of the boundary values of the function, and (ii) by the Fourier coefficients. We now discuss them both in brief.

Boundary values
Let v a j ,ṽ a j be the real part of the function Φ a (θ) at the points θ j and ib + θ j where We also define the functions Now the functions Φ a (θ) are taken to be It is easy to check that v a j = ReΦ a (θ j ) andṽ a j = ReΦ a (ib + θ j ). Notice that, for consistency, we need to impose In the main text we consider a strip b 1 ≤ Imθ ≤ b 2 that can be obtained from here by a simple translation along the imaginary axis.
which is equivalent to stating that the contour integral of the function around the domain is zero since there are no poles.

Fourier coefficients
We can expand the function Φ a (θ) as where M ∈ N is a high frequency cut-off. Since we are working with real analytic functions, the coefficients α a n are real. Now, consider In both cases, M is typically taken from 40 to 400 depending on the accuracy needed. The period ω is taken large enough to determine the salient features of the functions taking into account that beyond Reθ ∼ ω 2 the boundary effects due to the imposed periodicity are notable. On the imaginary axis, the functions are determined very accurately.
In terms of either set of variables (v a j ,ṽ a j ) or (α a 0 , α a n , andα a n ), the crossing constraints are linear (they can be imposed on the line Imθ = π 2 or between the upper line Im θ = b and the line Im θ = π − b). The unitarity constraints are imposed at the points θ j and are quadratic in the variables. The functional to maximize is taken as a linear function of the variables. Therefore the maximization problem becomes a standard conic convex optimization problem that can be solved by standard methods [39]. For example we can take a point θ 1 on the imaginary axis and impose R 1 (θ 1 ) = t cos ξ, R 2 (θ 1 ) = t sin ξ , (B.11) and maximize t. By sweeping the values of ξ ∈ [0, 2π] and plotting the resulting R 1 (θ 1 ), R 2 (θ 1 ) we find the boundary of the allowed region as depicted in figs.4, 6, 7 and 9. Moreover, for any point at the boundary of the region, the functions R 1,2 (θ) are determined with good accuracy except for boundary effects due to the imposed periodicity.

C Free bulk theory
When the theory in the bulk is free, i.e. S cd ab = δ c a δ d b , the boundary Yang-Baxter eq.(1.22) is trivially satisfied. In this appendix, we analytically obtain the allowed space of R-matrices with a free bulk and compare with numerics. The crossing equation for this case reads Once again, we consider diagonal (1.23) and block diagonal (1.27) reflections. For the diagonal case, we have two self-crossing functions that are bounded in the physical strip by modulus less or equal to one (by unitarity and crossing). So the allowed region for (R 1 (θ 1 ), R 2 (θ 1 )) should be contained in a square with vertices (±1, ±1). These vertices satisfy all the constraints and so, the allowed region should contain them and by convexity also their convex hull, namely the said square. Therefore, the allowed region is the square with vertices at (±1, ±1) corresponding to the usual Dirichlet/Neumann boundary conditions in different directions.

(C.2)
This problem can be solved analytically by mapping the strip 0 ≤ Im θ ≤ π to a unit disk (z = i−e θ i+e θ ) and using the results from [40]. In the language of chapter 8 of [40], the problem corresponds to a rational kernel with two poles (n = 2) located at z 1 and −z 1 where z 1 = i−e θ 1 i+e θ 1 and we are maximizing within the space H p=∞ of bounded analytic functions. Using the result below eq.(12), page 138 of [40]   where we used that n = 2 and a is real by real analyticity. Here −1 ≤ a ≤ 1 can be used as an arbitrary parameter that parameterizes the boundary of the allowed region which is then given by A(θ) = 2a sinh(θ) (a 2 + 1) sinh(θ) − i(a 2 − 1) , (C.4) B(θ) = i(a 2 − 1) cosh(θ) (a 2 + 1) sinh(θ) − i(a 2 − 1) . (C.5) As a check, one can match these analytic results with the ones from the numerical bootstrap (see fig. 12). Although very simple, we are not aware of work where this off-diagonal free theory reflection matrix was further studied.