Bootstrapping boundary-localized interactions

We study conformal boundary conditions for the theory of a single real scalar to investigate whether the known Dirichlet and Neumann conditions are the only possibilities. For this free bulk theory there are strong restrictions on the possible boundary dynamics. In particular, we find that the bulk-to-boundary operator expansion of the bulk field involves at most a `shadow pair' of boundary fields, irrespective of the conformal boundary condition. We numerically analyze the four-point crossing equations for this shadow pair in the case of a three-dimensional boundary (so a four-dimensional scalar field) and find that large ranges of parameter space are excluded. However a `kink' in the numerical bounds obeys all our consistency checks and might be an indication of a new conformal boundary condition.


Introduction and summary
The classification of conformal boundary conditions for a bulk CFT is a difficult problem. Besides well-known results for rational boundary conditions in rational CFTs (reviewed in [1]), very little is known even for relatively simple theories. It is natural to ask whether a systematic approach is feasible -one which does not rely on explicit constructions but leverages instead the modern conformal bootstrap methods [2] (see [3] for a review and [4] for a first application to BCFT which relied on results from [5,6]). A promising methodology is to start from theories that are as simple as possible in the bulk. In this work we pursue precisely such a direction in the case where the bulk theory is a single real free scalar field.
In any spacetime dimension a free scalar can certainly have Dirichlet or Neumann conformal boundary conditions. The question we try to answer here is whether more general conformal boundary conditions are possible, for example by coupling the bulk scalar to new boundary degrees of freedom and flowing to the infrared. These putative boundary conditions should modify the behavior of the scalar near the boundary and produce non-trivial boundary correlators, analogous to those of an interacting one lower-dimensional CFT. We find numerical evidence for at least one such 'exotic' boundary condition in four dimensions, and more generally very strong constraints on the space of potential conformal boundary conditions. 1 In exploring consistent boundary conditions for a free scalar theory we obtained a very special set of 'shadow-related' crossing symmetry equations, as follows. First of all, the φ = 0 equation of motion implies that the bulk-boundary expansion of the bulk field φ is restricted to contain at most two operators that we denote as O 1 and O 2 ; their dimensions are ∆ φ and ∆ φ +1, respectively. At most one of these two operators can vanish, and if so then we are in the Dirichlet or Neumann case and the two operators are immediately recognizable as the restriction of φ or ∂ ⊥ φ to the boundary. If they are both non-vanishing then the operators can be thought of as a 'shadow pair' in the sense of Ferrara, Gatto, Grillo and Parisi [11][12][13][14]. Their dimensions match this observation since 2∆ φ + 1 = d − 1, the dimension of the boundary, but the picture extends to their three-point functions: for a generic third defect operator O with dimension ∆ and spin l we find the relations (1.1) 1 The existence of such strong constraints is remarkable. In some cases, for example Maxwell theory in four dimensions, it is known that the space of conformal boundary conditions is vast, since it includes the space of all CFTs with a U (1) symmetry in three dimensions [7][8][9][10].
and a φ 2 is the one-point function of the operator φ 2 in the presence of the boundary. This relation between OPE coefficients agrees with the result obtained from applying a shadow transformation to the relevant three-point functions. In section 2 we derive this equation by demanding the absence of unphysical singularities in a three-point functions involving two bulk operators φ. 2 Note that the relations as written are still valid when the dimension of O is such that the gamma functions have poles; this is precisely when the operators are of 'double-twist' type and the shadow transformation is singular. The properties of the previous paragraph already lead to a remarkable bootstrap problem. Indeed, up to the special 'double-twist' operators there is one spectrum and set of OPE coefficients that needs to solve the five different crossing symmetry equations corresponding to the possible four-point functions of O 1 and O 2 . 3 This is intriguing in itself, but in the numerical analysis we can actually impose three more constraints. The first one is related to the Ward identity for the displacement operator. The second one is that of locality of the BCFT setup, which translates to the absence of any vector operators of dimension d in the O 1 × O 2 OPE. Both of these are described in section 2.2.3. The third one is imposed to separate local three-dimensional CFTs, which do not interest us here, from boundary conditions: this requires the scaling dimension of the first spin 2 operator to be strictly greater than 3.
We have numerically explored the system of crossing equations originating from the 1111 , 1122 and 2222 four-point functions in four bulk dimensions subject to all the above conditions. It might be tempting to conjecture that no non-trivial conformal boundary conditions exist that meet such stringent criteria, but surprisingly this is not quite what we find. On the one hand, there does exist a large range of possible values of a φ 2 where the first spin 2 operator must have a dimension less than about 3.1. Since this value is likely to decrease even further when increasing computational power, it is natural to conjecture that it must converge to 3, and then no conformal boundary conditions would be possible in this range. On the other hand, for a φ 2 near its Neumann value we suddenly find room for interesting physics: as indicated in 2 In [15] the same analysis was carried out for defects with a higher co-dimension in the free scalar theory, leading to similar shadow relations and a proof of triviality in many cases. For non-integer dimensions this setup can also be used to describe the long-range Ising model, where the relations can be found from the non-local equation of motion [16]. More details can be found in [17,18] and a first numerical analysis in this context was done in [19]. 3 The corresponding four-point functions should be related by the integral transformation that implements the shadow symmetry. However it is not clear to us whether the conformal block decompositions of such shadowtransformed four-point functions are automatically consistent. For example, the integral transformation is sensitive to contact terms and it seems unlikely that it can be swapped with the sum over conformal blocks. figure 2, the first spin 2 operator can have a dimension of nearly four without violating any other constraint. A corresponding kink in all the numerical plots points towards at least one possible exotic boundary condition for a free scalar in this neighborhood. We subject this point to a detailed analysis in section 5 and show that it passes all the obvious consistency checks for a proper local boundary condition. Using the 'extremal functional' methods of [20,21] we also estimate the dimensions of several low-lying operators and in particular conclude that the higher-spin symmetries of the bulk theory are broken. We leave for the future the interesting question of identifying a microscopic candidate for such a conformal boundary condition. The rest of the paper is organized as follows. In section 2 we study the constraints that any conformal boundary condition of a free scalar field must obey, and in particular derive the shadow relations (1.1). In section 3 we provide some examples, starting by reviewing the free boundary conditions and then using conformal perturbation theory to construct interacting ones (under the assumption that we can appropriately tune the parameters of the local degrees of freedom on the boundary). In section 4 we derive the set of crossing equations for the four-point functions involving O 1 and O 2 , we organize them in a way that takes advantage of the exact relations, and we explain the approximations of the resulting 'superblocks' that we use in our numerical implementation. In section 5 we present the numerical results in the case of d = 4, showing plots that involve several different characteristic observables, and in particular we show the kink that we mentioned above. We finally discuss possible future directions in section 6. A summary of the conventions and various technical results that we used along the way are relegated to the appendices.
2 Analytic constraints on the free scalar BCFT Consider a free massless scalar field φ in d > 2 dimensions with a planar boundary. We use the coordinate y ≥ 0 for the direction orthogonal to the boundary, and x for the directions parallel to the boundary. We denote the components of x = ( x, y) ∈ R d−1 × R + as x µ , µ = 1, . . . , d with x d = y, and those of x ∈ R d−1 as x a , a = 1, . . . , d − 1. We are interested in unitary boundary conditions that preserve the boundary conformal symmetry SO(d, 1).

Two-point functions with the scalar field
In this section we discuss the bulk-to-boundary OPE (bOPE) of the scalar field for a generic conformal boundary condition, and the constraints imposed by crossing symmetry on the bulk two-point function of the scalar field.

Bulk-boundary two-point functions
In this section we review the existence of two operators with protected scaling dimension in the bOPE of the free scalar field [4,[22][23][24][25]. In a BCFT the bulk-boundary two-point function of a scalar bulk operator O of scaling dimension ∆ O and a scalar boundary operator O of scaling dimension ∆ O is [5,6,4] The bOPE coefficient b O O , which is real for Hermitian operators, is not determined by symmetry.
Specializing O to be a free scalar φ of scaling dimension ∆ φ = d 2 − 1, the equation of motion φ = 0 gives Therefore the possible scaling dimensions for boundary primaries with b φÔ = 0 are Without loss of generality, we can assume there is at most one boundary operator of dimension ∆ 1 with b φÔ = 0, that we denote as O 1 , and similarly for ∆ 2 , the corresponding operator being denoted as O 2 . As observed in [25], the scaling dimensions of these operators add up to d − 1, which suggests that the two operators might be thought of as a 'shadow pair'. In the next subsection we will show that also their three-point functions are compatible with such a 'shadow relation'. The bOPE of the free scalar is [5,6,4] φ( x, y) The explicit form of the differential operator C ∆ i [y, ∇ 2 ] is given in appendix A. The bOPE can be used to reconstruct bulk correlators starting from the boundary ones.

Bulk two-point function
Next, we consider the two-point function This correlator is not completely fixed by the symmetry as it depends on the cross-ratio We can compute (2.5) by plugging twice the bOPE (2.4), using the diagonal and unit-normalized boundary two-punt functions and resumming the contributions from the descendants. The resulting boundary channel decomposition of (2.5) is An alternative way of computing the two-point function (2.5) is to invoke the bulk OPE φ×φ, namely where T µν is the stress tensor, and the operators J µ 1 ...µ with ≥ 4 are the tower of higher-spin conserved currents present in the free scalar CFT. The OPE data involving the stress tensor are [26] (2.10) Plugging in (2.5), we write the two-point function as a sum over bulk one-point functions and their derivatives. Boundary conformal invariance allows only for scalar bulk one-point functions [5,6,4], hence from the φ × φ bulk OPE the only non-trivial contributions are due to Resumming the contribution from bulk descendants we obtain the bulk channel decomposition of (2.5): (2.12) Equating the two different decompositions (2.8) and (2.12) gives the bulk-to-boundary crossing equation [4]. Since everything else in the equation is fixed, the only dynamical data are the onepoint function coefficient a φ 2 on the l.h.s. and the bulk-to-boundary couplings (b 1 , b 2 ) on the r.h.s. The solution is [4,24] This result tells us that in any boundary condition for a free scalar the parameter a φ 2 is constrained by unitarity to lie in an interval (2.14) As we indicated above, the boundaries of the interval correspond to the Dirichlet (b 1 = 0) and Neumann (b 2 = 0) boundary condition. These elementary boundary conditions will be discussed in detail in section 3, but in the remainder of this section we will assume that b 1 b 2 = 0 because we would like to explore the possibility of more exotic boundary conditions.

Three-point functions with the scalar field
In this section we consider three-point functions with two insertions of the free scalar φ and a generic boundary operator O. Note that, by Lorentz invariance, these correlators can be nonvanishing only if the third operator transforms as a symmetric and traceless tensor of SO(d − 1). Without loss of generality we can place the boundary operator at infinity and consider Following the standard procedure [27], in the above expression we contracted the tensor indices with a boundary polarization vector θ a as follows We will show that the boundary channel expansion of this correlation function exhibits unphysical singularities, which can be removed only if special conditions are met. Therefore these conditions have to be satisfied in any conformal boundary condition of the free scalar.

Boundary channel computation
The bOPE (2.4) allows to completely determine the correlator (2.15) in terms of the three-point functions between the operators O i , i = 1, 2, and O (l) . Conformal invariance fixes the latter three-point functions to take the form [26,27] where ∆ denotes the scaling dimension of the operator O (l) which carries SO(d − 1) spin l. The dependence on the polarization vector is through the following polynomial 19) which implies that only even spins l are allowed in (2.17) To compute (2.15), we act twice with the bOPE on the boundary three-point functions (2.17). After some algebra to resum the contributions from the descendants, we obtain the following boundary channel expansion (2.20) This expression depends on two cross-ratios w ± , which we take as follows: The functionsF ij ∆,l (w + , w − ) are computed in Appendix B.1 and their explicit expressions are given in (B.7). In the next section we will study the analyticity properties of the correlator (2.20).

Constraints from analyticity of the bulk OPE
Next, we study the same three-point function using the bulk φ × φ OPE. Since the only singular term in this OPE is given by the identity operator, which does not contribute to the three-point function, we conclude that the three-point function must be free of singularities when the two bulk points coincide. In terms of the cross-ratios w ± , this limit corresponds to w + → ∞ with any fixed w − .
As we show in the appendix C, for generic values of their parameters the boundary blocks on the r.h.s of (2.20) become singular in this limit. These unphysical singularities can be removed if the OPE coefficients are related in the following waŷ For certain special values of the parameters ( ∆, l) some of the blocks on the r.h.s of (2.20) are themselves regular. These special values correspond to the poles of the gamma functions in (2.22) and read (see also table 1) We denote these operators as • ∆ = d+l+2n−2 with n a positive integer. We have κ 1,2 ( ∆, l) → ∞ and sof 12 O (l) = 0 whilê Given that generically they appear in both OPEs, we could denote these operators both as The special cases listed above are related to the higher-spin symmetry of the bulk theory, as we will now explain. We recall that the φ × φ OPE (2.9) contains infinitely many higher-spin conserved currents J , with even spin ≥ 2 and scaling dimensions ∆ = d + − 2. The conservation of these currents is generically violated by terms localized on the boundary, leading to the following Ward identities In this formula any subset among the − 1 symmetric traceless indices {µ 1 . . . µ −1 } can be taken to be parallel to the boundary, with the remaining indices being orthogonal, i.e. in the y direction. Therefore, the BCFT generically contains boundary operators D (l) and V (l+1) of spin l and l + 1, respectively, and protected dimensions ∆ = d + − 2, where l is an even integer ranging from 0 to − 2. By 'generically' we mean that some of these operators might actually be absent from the spectrum in special cases. The equations (2.23) can be equivalently rephrased in terms of the bOPE, namely the operators D (l) , V (l+1) have the correct dimensions and spins to appear in the bOPE of the bulk higher-spin current J in a way that is compatible with its conservation in the bulk. Furthermore, spin selection rules and bulk conservation imply that V (l+1) cannot appear in the bOPE of any J with = , while the only other bulk current besides J that can contain D (l) in its bOPE is J l .
The relation to the special cases of (2.22) now stems from the observation that when −l = 2n, with n non-negative integer, D (l) has the right dimension to be the special operator

Displacement operator, flux operator and bulk-to-boundary crossing
The case = 2 deserves special attention because it corresponds to the bulk stress tensor T µν . Then the scalar operator D 2 ≡ D is the so-called displacement operator, and we will refer to the spin 1 operator V (1) 2 ≡ V (1) as the flux operator. Their general importance stems from the conservation of momentum P µ along a time coordinate chosen parallel to the boundary. If we split x µ = (τ, z, y) then where in the second equality we used the conservation to trade the time derivative with a spatial one, and then rewrote the integral of the spatial derivative as a boundary term. Choosing µ = y orthogonal to the boundary we find that the limit y → 0 gives the displacement operator D, which therefore measures the breaking of translations orthogonal to the boundary and must be non-zero for any sensible boundary condition. Choosing µ parallel to the boundary, on the other hand, we find the vector operator V (1) and so we conclude that it measures the flux of energy into the boundary. Theories with a non-trivial flux operator V (1) = 0 may still have a conserved boundary-translation charge, if there is an additional boundary contribution to the charge P a tot = P a + P a satisfying d dτ However the flux operator must be absent in any local unitary BCFT setup. To see why, note that the locality condition on the boundary is that P a , if non-trivial, should be expressible as the integral of a local boundary operator with two indices t ba . The condition (2.25) locally takes the form Moreover, by repeating the argument for the other generators of the conformal group on the boundary, one can easily show that the operator t ba has spin 2, i.e. it is symmetric and traceless.
Recalling that V (1)a has scaling dimension d and therefore t ba has scaling dimension d − 1, we see that eq. (2.27) with V (1) = 0 is incompatible with the unitarity bound of a spin 2 operator in d − 1 dimensions. We conclude that indeed in any unitary BCFT locality implies that in which case P a is trivial. In practice this means that if we couple a bulk CFT (not necessarily our free scalar theory) to some local boundary degrees of freedom, perhaps triggering an RG flow to a new conformal boundary condition, then the flux operator must never appear. 4 This is because local boundary degrees of freedom should not be able to hold a macroscopic amount of energy. It might be instructive to consider some non-local setups that do feature a flux operator. The first is a conformal interface, where there is an entire new CFT living on the half space with y < 0. In that case the stress tensor for each side 'sees' a flux operator, but if the interface setup is local then these two flux operators are in fact the same operator and the interface cannot act as a simultaneous energy sink for both sides. Such a setup can be generalized to the case where our d − 1 dimensional boundary is at the same time a conformal defect in some d -dimensional auxiliary space in which it is coupled to an arbitrary d -dimensional bulk CFT, perhaps even triggering a boundary/defect RG flow to some new conformal configuration. According to the general structure of the operator expansion near a defect of dimension d − 1 in a d -dimensional CFT, the d -dimensional stress tensor can always provide a vector operator of precisely the requisite dimension 5 (which is d) and unless the two sides decouple we will observe this as a flux operator in the d-dimensional BCFT. 6 Lastly one could also try to create a non-local setup by adding a GFF on the boundary and coupling it, perhaps with other degrees of freedom, to the bulk field. But this scenario is captured by the previous one, because GFFs are just regular local fields in an auxiliary higher-dimensional space (albeit with non-integral d ).
The previous discussion applied to any BCFT, but for the free scalar theory there are a few additional results that we can derive. To this end we return to the φφ O three-point function of equation (2.20) and take the third operator to be either the flux operator V (1) or the displacement operator D. Let us thererefore first consider: (2.29) The relations (2.22) in this case give thatf 11V (1) =f 22V (1) = 0 and leavef 12V (1) undetermined.
On the other hand, the bulk channel decomposition of (2.29) receives contribution only from the bulk stress tensor T µν , because -as we explained in the previous subsection for the general case -spin selection rules and bulk conservation imply that that V (1) cannot appear in the bOPE of any other operator in the OPE (2.9). We find where the ellipses represent contributions of bulk descendants, and the OPE data of the stress tensor are given in (2.10). The two-point function on the r.h.s. is completely fixed by the boundary conformal symmetry up to a single bOPE coefficient b T V (1) [6,4]. We can further relate b T V (1) to the two-point function coefficient From the exact correlator (2.20) we find (details in appendix C.2) and we conclude that the flux operator appears in the O 1 × O 2 OPE if it also appears in the bOPE of the stress tensor. 7 Next we consider the three-point function involving the displacement operator In this case, the relations (2.22) give thatf 12D = 0 while leavingf 11D andf 22D undetermined. Using again the general argument from the previous subsection about spin selection rules and conservation, we have that the only operators in the OPE (2.9) that contribute to the bulk channel decomposition of (2.33) are φ 2 and the bulk stress tensor. Therefore we have where again the ellipses denote contributions of bulk descendants. Using the Ward identities for the displacement operator [6], the bulk-boundary two-point functions in the r.h.s. are determined in terms of the parameter a φ 2 in (2.14), as well as the coefficient C D in the two-point function of the displacement operator Comparing with the boundary channel correlator (2.20), after some algebra which we relegate to appendix C.3, we find (2.36) The unitarity requirement C D ≥ 0 implies that: (2.37)

The three-point function with the boundary modes of φ
Another interesting special case of (2.20) arises when the boundary operator is one of the boundary modes of φ, i.e.
On general grounds, due to Bose symmetry, there are four independent boundary OPE coefficients that enter these correlators:f 111 ,f 112 ,f 221 ,f 222 . The latter are further related to each other, by means of the three independent constraints provided by regularity of the φ × φ OPE (2.22): Hence, the bulk three-point function φφφ is completely controlled by a single boundary OPE coefficient, e.g.f 111 . The latter can be non-zero only if the boundary condition breaks the Z 2 global symmetry φ → −φ , under which both O i are odd.

Examples
In this section we explore some examples of conformal boundary conditions for a free scalar. We start by reviewing the free boundary conditions, i.e. Neumann and Dirichlet, and then construct examples of interacting boundary conditions using conformal perturbation theory around the free ones. As we will see, these constructions rely on some ad-hoc assumptions on the spectrum of an additional local 3d sector living on the boundary, which we couple to the bulk, and therefore they do not prove rigorously the existence of interacting boundary conditions. On the other hand, they will provide useful benchmarks to compare to our numerical results in section 5.

Free boundary conditions
Suppose the theory is fully described by the free bulk action 8 without any boundary-localized interaction. In order to have a stationary action, besides the bulk equation of motion φ = 0 we need to set to zero the boundary term The two solutions to this condition that preserve boundary conformal invariance are We can rephrase these conditions in terms of the bOPE of the scalar field. Namely, in (2.4) we have b 2 = 0 and φ| y=0 ∝ O 1 in the case of Neumann boundary condition, and b 1 = 0 and ∂ y φ| y=0 ∝ O 2 in the case of Dirichlet boundary condition. In either case, there is only one boundary operator in the bOPE of φ, and the full set of boundary correlators can be simply characterized as the mean-field theory of this operator. This implies that all the correlation functions of these BCFTs are completely disconnected, i.e. they are computed by Wick contractions as products of twopoint functions. For this reason we call these boundary conditions 'free boundary conditions'. We can also consider additional free boundary conditions that are not Neumann and Dirichlet. Such a boundary condition is obtained requiring both O 1 and O 2 to appear in the bOPE (2.4), i.e. 8 Note that using this canonical normalization of the action the operator φ has a different normalization compared to the one in equation (2.9), namely φ| (3.1) = φ| (2.9) . We will specify which normalization we are using whenever important. b 1 b 2 = 0, and postulating that these operators are two decoupled generalized free fields. 9 However this implies that there is a spin 1 operator of dimension 4 in the spectrum of the boundary theory, namely the vector 'double-trace' operator in the OPE of It is easy to check that this operator also appears in the bOPE of the bulk stress tensor, hence for these boundary conditions we have a non-vanishing flux operator V (1) = 0. Therefore, following the discussion in the previous section, these are non-local boundary conditions. We conclude that the only local free boundary conditions are the familiar Neumann and Dirichlet boundary conditions reviewed above.

Interacting boundary conditions in perturbation theory
In order to look for examples of interacting boundary conditions, a natural strategy is to couple the bulk scalar to a CFT d−1 living on the boundary. We turn on some relevant interaction between the two sectors and then flow to the IR, hoping to reach a non-trivial BCFT fixed point. Concretely, we add to the free bulk action (3.1) a boundary action of the form where σ I are some scalar composites made of φ| y=0 or ∂ y φ| y=0 , depending on whether we start with Neumann or Dirichlet boundary condition, as well as of local operators of the CFT d−1 . In order to have perturbative control over the resulting RG flow, we will assume that the operators σ I have scaling dimensions i.e. the deformations are weakly relevant. Then one can systematically expand observables of the BCFT at the putative IR fixed point as a series in I . We will further assume that the boundary degrees of freedom are local. Technically, this means that in the absence of bulk-boundary couplings, i.e. for g I = 0, the spectrum of the CFT d−1 contains a stress tensor t ab , which is a conserved, spin 2 primary operator, of protected dimension ∆ t = 3. At the perturbative BCFT fixed point this operator gets a small anomalous dimension, which must be non-negative by unitarity, and actually strictly positive if the bulk and the boundary are not decoupled. We refer to this spin 2 operator at the interacting fixed point as 'pseudo stress tensor'. In the next subsection we show how to compute the leading order contribution to this anomalous dimension for a rather generic interaction of the form (3.4), using 9 This includes the case of 'no boundary', or more precisely the 'trivial interface', where the theory on the full R d is re-interpreted as a BCFT. In that case a φ 2 = 0 so according to (2.13) this corresponds to b 2 1 = 1 and b 2 2 = (d − 2). multiplet recombination. We will then consider more specific examples for the perturbation, and compute the leading order corrections to the observables a φ 2 and C D , defined in eqns. (2.11) and (2.35), respectively. Typically, when computing (B)CFT observables in perturbation theory, one first computes the corrections as a function of the coupling constants, and then plugs the value of the coupling constants at the fixed point, obtained by solving for the zeroes of the beta functions. However, by restricting to the case with a single bulk-boundary coupling, we can also avoid the computation of the beta function and simply assume that a perturbative fixed point exists. This is sufficient because we can consider ratios of the leading order corrections to the observables mentioned above, in such a way that the coupling cancels out from the ratios. It would be interesting, but much more laborious, to actually compute the beta functions in terms of the data of the CFT d−1 . This would actually be necessary if one wants to verify the existence of the fixed point, or consider higher order corrections/multiple bulk-boundary couplings. The beta function needed in this setup starts at cubic order in the coupling, and the coefficient of the cubic term is given by a regularized integral of the four-point function of the deformation, see e.g. [29,17] and also [30] for the case of 1d CFTs. 10

Anomalous dimension of the pseudo stress tensor
We now consider a slightly more specific bulk-to-boundary interaction, with a single coupling, of the form In the expression above, χ denotes an operator in the CFT d−1 and Ω is any local boundary operator built out of φ| y=0 or ∂ y φ| y=0 , depending on whether we are perturbing a Neumann or Dirichlet free boundary condition. The assumption (3.5) in this case takes the form In the presence of the interaction (3.6) the conservation and the tracelessness of the stress tensor t ab of the CFT d−1 is violated as follows The computation of the beta functions for bulk-boundary couplings in terms of the data of the CFT d−1 was performed in [31] for some examples of perturbations around Dirichlet and Neumann. Some perturbative constructions of interacting boundary conditions for free theories can also be found in [32,33,25].
Assuming a nearby fixed point with g 2 ∝ , 11 we have two seemingly problematic features in the above equations, namely the divergence is not expressed in terms of a primary operator of the undeformed theory, and the operator does not have only a spin 2 component because the trace is non-zero. Both these issues are solved by defining the 'corrected' operator Taking the divergence we then obtain The new operator τ ab is a symmetric traceless tensor, and its divergence (3.10) is a primary spin 1 operator of the undeformed theory, making the recombination of the multiplets manifest. Note that (3.10) is a manifestation in perturbation theory of the locality condition that we discussed in 2.2.3. If the boundary degrees of freedom were non-local they would not have the operator t ab and then the right hand side of (3.10) would be a primary operator of spin 1 and protected dimension d (it is easily checked that indeed this operator would appear in the bulk-to-boundary OPE of the bulk stress tensor). We can exploit the recombination to compute the leading order anomalous dimension of τ ab at the interacting fixed point. Let us consider computing the two-point function On the one hand, we can take derivatives of the two-point function of τ ab , which is fixed by boundary conformal invariance to be [26] τ ab ( x) τ cd (0) = C τ (g) The definition of I ab was given in (2.18), and we introduced The constant C (0) τ is the 'central charge' of the CFT d−1 that the bulk scalar couples to, i.e. the coefficient appearing the two-point function of the stress tensor t ab before we turn on the interaction. On the other hand we can compute (3.11) at the leading order in g by directly using the r.h.s. of (3.10). By comparing the two results, we find O denotes the coefficient of the two-point function of the boundary operator O computed at g = 0. With the canonical normalization (3.1) of the bulk action we have (3.14) We note that the leading order anomalous dimension is essentially controlled by the central charge C

Modified Dirichlet boundary conditions and perturbation theory
We now further specialize to the case in which the free boundary condition is Dirichlet, and the operator Ω is ∂ y φ| y=0 , namely we take a deformation of the form The interaction term leads to the the following modified Dirichlet boundary condition 12 In this case the condition (3.5) gives ∆ χ = d 2 − 1 − , with 0 < 1. As we discussed above, we assume the existence of a perturbative fixed point with g 2 ∝ . Plugging in eq. (3.13) we obtain (3.17) Let us now consider the leading order correction to the one-point function coefficient a φ 2 (g) The coefficient δa (1) φ 2 must be non-negative as a consequence of the unitarity bound (2.14). To compute its value, note that the modified Dirichlet boundary condition (3.16) determines the bOPE coefficient b 1 to be 13 Plugging this result in the crossing relations (2.13), we find Having obtained two observables we can form a ratio that does not depend on the value of the coupling at the putative fixed point, namely ( 3.21) This quantity depends on the central charge C (0) τ of the CFT d−1 that the bulk scalar couples to. Next, we consider the leading order correction to the coefficient C D in the two-point function of the displacement operator D denotes the value at the free Dirichlet boundary condition. The displacement operator in this theory is [5,6] Note that this formula makes sense even for an interacting boundary condition if we interpret the composite operators on the right hand side as products of φ| y=0 and ∂ y φ| y=0 , made finite by subtracting all the singular terms in the OPE. This is because D( x) = lim y→0 T yy ( x, y) and the bulk operator T yy is always equal to (3.23). (The scaling dimension of D is guarantueed to come out correctly because φ| y=0 and ∂ y φ| y=0 have protected dimensions.) This observation allows us to easily compute the two-point function of D in conformal perturbation theory for the modified Dirichlet condition (3.16). We find that the contributions from φ| y=0 are O(g 4 ) whereas the two-point function of ∂ y φ| y=0 is corrected already at O(g 2 ) by the interaction term (3.15) and is given by: (3.24) Note that the integrals have a power-law UV divergence for u ∼ x and u ∼ x that we subtracted. As a check, the result (3.24) implies which is in perfect agreement with the correction (3.20) that we computed for a φ 2 and the crossing relations (2.13). Using (3.24) to compute the two-point function of 1 2 (∂ y φ) 2 y=0 and therefore of D, we obtain (3.26) We can then form another ratio of observables that is independent of the coupling at the putative perturbative fixed point Note that this ratio does not depend on any data of the CFT d−1 and therefore it is a universal result for deformations of the form (3.15) of the Dirichlet boundary condition.

Modified Neumann boundary conditions and perturbation theory
As a final example, we consider deformations of the Neumann free boundary condition by the following interaction The interaction gives rise to the following modified Neumann boundary condition The condition (3.5) now gives ∆ χ = d 2 − , with 0 < 1, and again we will assume the existence of a perturbative fixed point with g 2 ∝ . Plugging in eq. (3.13) we obtain (3.30) To compute the variation of the parameter a φ 2 we use the same strategy as in the previous example, namely it follows from the modified Neumann condition that and using the crossing relations (2.13) this gives (3.32) Note that this has an opposite sign compared to eq. (3.20), in agreement with the unitarity bounds (2.14). The coupling-independent ratio then is (3.33) Like in the previous example we now compute the correction to C D , again using the definition (3.23) as the starting point. The main difference is that in this case the leading-order correction comes from the second and third terms in eq. (3.23), namely those involving φ| y=0 , while the first term involving ∂ y φ| y=0 only starts contributing at subleading order O(g 4 ). We will then only need the two-point function of φ| y=0 up to O(g 2 ) corrections, that is (3.34) The integrals have a power-law UV divergence for u ∼ u that we subtracted. As a check, from (3.24) we obtain which, upon substitution in the crossing equations (2.13), gives a correction to a φ 2 in agreement with (3.32). Using (3.34) we obtain (3.36) Comparing with (3.26) we see that the value at the free boundary condition is the same for Neumann and Dirichlet, while the leading correction differs by a factor of d − 2. Taking the ratio with δa 2 φ we get which notably is the same as the one obtained for the deformation of Dirichlet in eq. (3.27). Like in that example, this ratio is universal for deformations of the form (3.28) of the Neumann boundary condition, because it does not depend on data of the CFT d−1 .

Crossing equations
In this section we present the crossing equation for the mixed system of four-point functions of the boundary modes of φ, namely The crossing equations for a generic mixed system of scalars, labelled by indices i, j, m, n, were derived in [34] and read where u = . The functions F ij,mn ±, ∆,l are the following combinations of the ordinary s-channel conformal blocks g Note that not all equations in this system (4.2) are independent, since 14 14 Recall that [35,36] g ∆12,∆34 ∆, If we specialize all these ingredients to our problem where i, j, m, n ∈ {1, 2} then we find 7 independent crossing equations: (4.6) In our case the operators must also obey the OPE relations (2.22). Imposing those, the system of equations can be rewritten as follows (details can be found in appendix D) . (4.7) The

Implementing the exact relations
We will numerically 'bootstrap' a set of crossing equation in the sense of [2]. For most problems, the fastest program available for this task is the semidefinite program solver SDPB [37]. We have used the recent version [38] which supports the 'hotstarting' algorithm suggested in [39]. In order to use SDPB efficiently, the ingredients of the crossing equations must be approximated as rational functions of the scaling dimension such that all poles of odd order lie at or below the unitarity bound. This is always possible when the basis functions are conformal blocks and we refer the reader to [34] for details. The complication in this work is that the vector V +, ∆,l in (4.7) has several occurences of ∆ which are not in conformal blocks. As shown in Appendix D, we must consider linear combinations in which the coefficients are various products of κ 1 ( ∆, l) and κ 2 ( ∆, l) -the functions from (2.22). While these types of blocks were first introduced for studying the long-range Ising model, in the numerical analysis of [19] only the last two components of (D.6) were used. In addition, κ 1 ( ∆, l)κ 2 ( ∆, l) was treated as a vector with > 1000 discrete evaluations, thereby eschewing some of the benefits of semidefinite programming. In this work, we do not need to limit ourselves to those crossing equations that involve only the product κ 1 ( ∆, l)κ 2 ( ∆, l) where b 1 and b 2 cancel out. Thanks to (2.13), the other products only introduce one new parameter and it has a clear physical meaning in the BCFT context. To remedy the second problem, we need to discuss the approximation theory of gamma functions.
To start, it is easily verified that As such, a rational approximation 15 for κ( ∆, l) 2 will cover the cases of κ 1 ( ∆, l) 2 , κ 2 ( ∆, l) 2 and 15 With infinitely many poles above the unitarity bound it is clear that any rational approximation for κ( ∆, l) 2 is going to have significant errors in the semi-infinite range of allowed values for ∆. Extremely large values of ∆ should however be unimportant for the numerical results, and for a finite window of values a rational approximation is perfectly feasible. We have attempted to account for this in practice by using the simpler crossing equations (4.6) to cover the 'tail' of the more constraining crossing equations. Each time we approximate κ( ∆, l) 2 with ∆ ∈ [∆ 0 , ∞), we allow these extra conformal blocks, multiplying independent f ij O coefficients, to have an exchanged scaling dimension in [∆ 0 + 20, ∞). κ 1 ( ∆, l)κ 2 ( ∆, l). The most expensive step for our purposes will be the Weierstrass formula which introduces a new series of poles for each gamma function. Unlike [19], which suggested using (4.10) on the full function, we will only apply it to the − ∆ part of κ( ∆, l) 2 .
The + ∆ part, since it is regular, should be approximated with one of the many expressions for the Wallis ratio. This is a quantity which has attracted interest for hundreds of years due to the application of calculating π. In particular, we note the asymptotic formula It was found in [40] that (4.11) is the n = 1 case in a sequence of approximants that have the schematic form (z n + . . . ) 1 2n . We cannot use these higher radicals due to the requirement that κ( ∆, l) 2 be a rational function but it is still possible to make (4.11) arbitrarily accurate. One simply applies the functional equation n times to arrive at We have not found it necessary to choose a large value of n. For example, even when z = 1 4 , Weierstrass does not become better than (4.11) until k = 36.
We will now quote expressions for the two factors of κ( ∆, l) 2 . Using (4.12) with n = 1, (4.13) The singular part requires a cutoff which we call k max . (4.14) While it is optional to resum the exponent in the second step, the logarithmic behaviour of the harmonic series makes it convenient.
The expressions (4.13) and (4.14) have been implemented as a patch for the helper program PyCFTBoot [41]. For the poles exhibited here, which are two units apart, we have taken k max = 20. The poles coming from conformal blocks [34] are only one unit apart so we take k max = 40 for those. The standard way to account for poles in ∆ is to absorb them into the OPE coefficients of (4.7) so that crossing symmetry becomes a statement about polynomials. The most desirable type of problem for SDPB is one in which these polynomials can be expressed in terms of an orthonormal basis [37]. Recently, [42] gave an example of a problem which cannot be optimized in this way. In our case, this privileged basis of polynomials is again unavailable due to the 20 double poles of (4.14) that are above the unitarity bound. For this reason, we have opted to still use the simple crossing equations (4.6) for spins above a certain cutoff l 0 . For most of the bounds in the next section, this is l 0 = 4 while some of them have been redone with l 0 = 6. Seeing almost no difference, we conclude that the exact relations for l = 0 and l = 2 are doing most of the work. The other limitation of our approach is that the square root in (4.12) can only be eliminated when the κ i ( ∆, l) appear quadratically. This forces us to drop . According to standard lore, the bounds should be unaffected as these two correlators exchange the same operator families as the other three.

Numerical results for 4d/3d systems
Let us collect the constraints used for the numerical bootstrap analysis. First, we take the crossing equations for (4.6). These apply to any (possibly non-local) CFT containing the operators O 1 and O 2 . Second, we have the exact OPE relations (4.13) and (4.14) which are necessary for a solution of these crossing equations to be an admissible boundary condition for a massless free scalar in d dimensions. These conditions reduce the crossing equations to equation (4.7) at the cost of introducing a new parameter, a φ 2 , that our bounds will depend on. Notice that this in particular implies that the odd-spin operators can only have the scaling dimensions of the generalized free theory. Third, the boundary spectrum cannot have a stress tensor, so it is natural to demand that the first spin 2 operator has a dimension ∆ τ strictly larger than d − 1. Fourth, we should demand that the flux operator V (1) of dimension d is absent to avoid interfaces and other possible sources of non-locality on the boundary, see the discussion in 2.2.3. Fifth, we have the Ward identities for the displacement operator (2.36) which restrictf 11D andf 22D to a curve parametrized by C D .
We will set d = 4 throughout in order to work with a correlator system that involves 3d conformal blocks. As discussed in [43,44], similar problems with 2d blocks often require more experimentation with the gaps being imposed. These works are concerned with maximizing the gap in the scalar sector, and indeed, we can provide a nice preview of our results by doing the same. Figure 1 bounds the dimension of the lightest exchanged scalar, which we call ∆ ε , as a function of the pseudo stress tensor dimension ∆ τ . To obtain this plot we scanned over all the allowed values of a φ 2 . The so obtained blue region is clearly smaller than the pink region, obtained without imposing the OPE relations, or the single correlator region delineated by the upper black line. Three further comments are worthwhile. Figure 1: A plot showing the upper bound on the dimension of ε, the first scalar, other than the identity, seen by any of the OPEs in our correlator system. The unshaded region is the one that follows from a single correlator O 1 O 1 O 1 O 1 . The pink region, which is more restrictive, uses the multi-correlator system but the only inputs it uses from the exact relations are the odd-spin operator dimensions given in table 1. The blue region, more restrictive again, follows from a genuine use of the exact relations. Since these depend on a φ 2 , we have extremized ∆ ε over a third axis which is not shown.
First, we observed that much of the constraining power came from our fourth constraint, i.e. the exclusion of the dimension d vector V (1) = [ O 1 O 2 ] 0,1 from the spectrum. In fact, if we were to reinstate just this vector then the blue region would expand to almost the same size as the pink region. We emphasize that the OPE relations are essential to meaningfully impose this constraint: they prevent the appearance of vector operators of dimensions very close to d that would numerically be indistinguishable from V (1) . Furthermore, because of the fake primary effect [45,46] the block for V (1) can be mimicked in our numerical analysis by a spin 2 operator of dimension 3, and therefore the constraint that ∆ τ > 3 (strictly) is also essential to ensure that it is really absent. This latter argument relies on the observation that, for a spin 2 operator whose dimension ∆ → 3, the corresponding combination of blocks that enters in the crossing equation (4.7) is: by virtue of the OPE relations discussed in section 4, whose notation we follow here. Therefore we can recover a vector operator if we assume that its overall coefficientf 12 Second, for the Dirichlet or Neumann boundary conditions either O 1 or O 2 vanishes so in some sense they are not within the reach of our numerical analysis. On the other hand we can find a more general solution including both O 1 and O 2 as independent GFFs, with any value of a φ 2 . This solution does contain the vector V (1) , and because of the fake primary effect we just discussed it corresponds to the point with ∆ ε = 2 and ∆ τ = 3, which is well within the allowed region.
Third, one may check that these bounds are saturated by the following two extremal solutions. The point ( ∆ τ , ∆ ε ) = (4, 2) represents the 'single GFF' solution where O 1 is a GFF and ε = O 2 = O 2 1 . This satisfies our crossing equations because it consists entirely of protected operators in (4.7). Also, the aforementioned vector is indeed absent from the spectrum of primaries because 1 in this theory. Since the bound in figure 1 can only decrease as a function of ∆ τ , we can be confident that it will stop changing once it hits ∆ ε = 2. In the blue plot, this turns out to happen well before ∆ τ = 4. We can also understand the point ( ∆ τ , ∆ ε ) ≈ (3, 2.95): here the four-point function of O 1 can be the extremal solution for a local three-dimensional CFT with ∆ = 1, which according to [47] has ∆ ε ≈ 2.95, and then O 2 can be a disconnected GFF. This setup satisfies all of the constraints we have imposed (we are of course not imposing ∆ τ > 3 here) except for the absence of V (1) , which again manifests itself as a spin 2 operator of dimension 3.
We find it plausible that, with infinite computational power, the drop from ∆ ε ≈ 2.95 becomes infinitely sharp leading to a value of ∆ ε = 2 almost everywhere. The remainder of this section is about what lies below ∆ ε = 2. 16 16 The bulk has a global reflection symmetry φ → −φ under which O 1 and O 2 are odd but ε is even. Since ∆ ε < 3 always, and ∆ ε ≤ 2 seems likely, any non-trivial boundary condition must be (strongly) unstable even for RG flows that preserve the Z 2 symmetry. The cases of Dirichlet and Neumann are not included in this discussion because as we explained in these cases one should remove many operators from the spectrum. The Neumann condition is also (strongly) unstable due to the operator φ 2 , while in the Dirichlet case the leading Z 2 even deformation is (∂ y φ) 2 so this boundary condition is stable. with our best estimate for the allowed region shaded in blue. The curves have n max = 5, 6, 7, 8 in the notation of [48,49]. As for the number of derivative components being kept in each crossing equation, these correspond to 21, 28, 36, 45 respectively. The dotted line shows the maximum possible value for ∆ τ from leading order conformal perturbation theory under the assumption that the Ising model is the 3d CFT with the lowest central charge.

A universal bound
The next parameter to introduce is a φ 2 , which through (2.13) determines b 1 and b 2 , in order to more fully exploit the exact relations. When scanning over a φ 2 , it is instructive to first determine its value for the two extremal solutions at ( ∆ τ , ∆ ε ) = (4, 2) and at ( ∆ τ , ∆ ε ) ≈ (3, 2.95) discussed above. In the first we found a spin 2 boundary operator with ∆ τ = 4 corresponding to an unprotected block in the first line of (4.7). Since the overall coefficient of this combination iŝ to maximize ∆ τ since this can be interpreted as a measure of how non-local a CFT is. Figure 2 presents this result. The four different lines correspond to four different search spaces, giving a sense of how close we are to having an optimal bound. Every other plot in this section uses the number of components corresponding to the second most restrictive region in figure 2. Figure 2 enables a comparison with the results of conformal perturbation theory, particularly around a φ 2 = − 1 4 where the bound is very strong and we can make a meaningful comparison with the deformation of the Dirichlet boundary conditions by a putative 3d CFT with a scalar operator with ∆ χ = 1 − , that we studied in section 3.2.2. Recall that, according to equation (3.17), the anomalous dimension of the spin 2 operator depends on the unperturbed central charge C (0) τ . Even though we do not know a theory with an operator that can play the role of χ, it is clear that the central charge of a unitary 3d CFT cannot be arbitrarily small. In fact, we believe that it is not unreasonable to assume that the 3d Ising CFT with 17 is the theory with the lowest possible central charge. An early indication for this conjecture was the local minimum corresponding to the Ising CFT found in [50,51], and recently in [52] a rigorous lower bound was found that, with sufficient numerical precision, is likely to lie between about 0.6C free τ and 0.95C free τ . For us this implies that as a bound on the anomalous dimension of the first spin 2 operator. In figure 2 it follows that every such example must lie below the dotted line. The possibilities obeying (5.3) are all within the allowed region for now, but we will see that many of them are ruled out when we add more constraints.

The kink and the extremal spectrum
We now come to the most striking feature of figure 2 which is the jump near the right hand side. Since the convergence appears to be rapid in this vicinity, we can be confident that the coordinate at which the curve flattens again is not tending towards a φ 2 = 1 4 . In other words, there is a kink at (a φ 2 , ∆ τ ) ≈ (0.215, 3.966) which obeys the exact relations and cannot be one of the free boundary conditions. If it is truly a new boundary condition then it must obey a further constraint that we have not yet imposed: the Ward identity (2.36) for the displacement operator. This is one of the reasons we would like to investigate the spectrum at the kink in more detail. 17 We are using conventions such that C free  Table 2: The low-lying spectrum at two points in (a φ 2 , ∆ τ , C D ) space. The point associated with the left table is still visible after projecting down to just (a φ 2 , ∆ τ ) -it is the kink in figure  2. Due to our maximization choice, we see every possible operator with odd spin. Protected operators (the ones with integer scaling dimensions) of even spin have vanishing mixed OPE coefficients and they start above the leading twist. Unprotected operators have their OPE coefficients related by (2.22) and thus we have shown calculated values in red.
The extremal spectrum The extremal functional method [20,21] allows for the extraction of an approximate spectrum and OPE coefficients for any point on the boundary of an allowed region. We have done so at two points: the first is the kink in figure 2 and the second involves a tuning of the displacement central charge C D using the procedure that will be explained in the next subsection. The CFT data for operators with ∆ < 6.5 are listed in table 2. The black numbers were obtained from the output of the script in [53] and the red numbers were computed using the OPE relations. Our OPE coefficients are defined such that for the standard cross-ratios z andz approaching zero along the diagonal. Notice that the most stable result is obtained after maximizing an OPE coefficient near the boundary of the plot to make the functional as close to extremal as possible. We have chosen to optimize the coefficient of the V 4 operator, which is a spin 3 operator of dimension 6. Our reason for doing so is to avoid another interference from the fake primary effect: much like V would imply that the higher-spin charge corresponding to the bulk spin 4 current is preserved by the boundary. We do not expect such an 'integrable' boundary, and our optimization minimizes the chances of an unwanted spin 4 operator taking the place of V After going through this maximization, a reassuring feature we observe is that there is no even spin l operator with ∆ = d + l − 2, as anticipated in the range of the sum in (4.7). Such a block, if present, would have to be treated with V 0,d+l−2,l because the exact relations degenerate at this point. But indeed, bulk spin currents only have boundary modes up to l = − 1 in the bOPE and there is no reason for l = to be present as a protected operator.
The displacement Ward identity With the spectrum in hand we can investigate whether the Ward identity (2.36) for the displacement operator is satisfied. However this is again a rather subtle business, this time because there might be other scalar operators of dimension 4. The correct procedure is as follows.
Consider all the scalar operators of dimension 4 in the putative extremal solution at the kink. One of these operators is the displacement operator, and its coefficientsf 11D andf 22D must obey (2.36), i.e. they are constrained to lie on a curve parameterized by C D . Every other operator is then not a displacement operator, meaning that it must be absent from the bulk-to-boundary OPE of the stress tensor: b T D = 0 for any 'non-displacement' D . Repeating the arguments in appendix C.3 this leads to the condition that: This is a matrix of rank two and we need more than one operator. It is natural to try to see if we can fit it with one displacement and one non-displacement operator. Notice that the matrix has three independent entries but for the two operators we only have the two parameters C D and b φ 2 D . Two simple approaches can be taken at this point. In the first approach, we demand that a non-displacement is exactly present and extract b φ 2 D : To fit the rest of the matrix we need the given displacement OPE coefficients, which lead to C D = 0.0053 fromλ 11D and C D = 0.0049 fromλ 22D . Alternatively, in the second approach, we demand that one of the outer products is an exact displacement and extract C D and OPE coefficients for the non-displacement: This leads to b φ 2 D = 1.914 fromλ 11D and b φ 2 D = 1.767 fromλ 22D . Although the small mismatches in both approaches might be due to numerical errors, it seems reasonable to conclude that this solution is not as physical as the solution on the left hand side of table 2.

Local boundary conditions
As with numerical bounds on the gap, the exact relations also lead to significant improvements for bounds on OPE coefficients. Consider again the (unit-normalized) displacement operator which appears with the coefficientsλ 11D andλ 22D . To constrain them, we set as in [54], then apply standard methods for bounding the magnitude of an OPE-space vector [55]. Figure 3 shows the results of this exercise for different values of a φ 2 .
A first thing to note is once more the importance of the exact OPE relations in (2.22). Without them the allowed region would certainly be the union of all the regions in figure 3. However for a φ 2 → 1 4 we observe an unbounded growth in the vertical direction (note the different vertical scales), and thereforeλ 22D is really only bounded by virtue of the OPE relations.
The dotted line represents the combinations of OPE coefficients that obey the Ward identity (2.36), parameterized by the displacement central charge C D . As discussed above, the non-trivial fact about the kink in figure 2 was that it happened to sit on this line. The intersection of the dotted line with the allowed region also translates into a lower and upper bound for C D for each value of a φ 2 . This is shown as the pink region in figure 4. This bound is certainly valid but rather crude: it does not take into account the restriction (5.5) on additional scalar operators of dimension 4 that are not the displacement. To do better we can assume a fixed displacement operator with a certain C D in the crossing equations by replacing and removing the scalar of dimension 4 from the special operators in the crossing equation (4.7). We then bisect in C D to find the allowed region, and this leads to the much improved blue region in figure 4. One may wonder if the blue region allows for other scalar operators of dimension 4 that are not the displacement operator. The answer is that it does, because such operators lie in the allowed continuum of operators. Furthermore, 14) which implies that the limit of a scalar operator as∆ → 4 in the continuum is actually exactly a D operator whose OPE coefficients automatically obey (5.5). So fixing C D not only allows one to single out a displacement operator for which the Ward identities (2.36) are obeyed, it also ensures that (5.5) holds for every other scalar of dimension 4.
We can compare figure 4 to the conformal perturbation theory results (3.27) and (3.37). The lines corresponding to potential perturbative fixed points saturate the lower bound on C D once the proper constraints on dimension 4 scalars are imposed. These emanate from the points (a φ 2 , C D ) = ± 1 4 , 1 2π 4 at which the upper and lower bounds are forced to meet by the Ward identity. The other point that can be explained analytically is the origin which is associated with no boundary at all. This point has to be allowed by the pink region since it corresponds to adding zero in (5.13). Once we classify dimension 4 scalars into displacements and non-displacements, it appears that there are no longer any nearby solutions that would allow us to see this point in the blue region. To see that they cannot arise from a GFF construction, consider the explicit displacement operator (3.23). We may rewrite it as by using the bulk two-point function (2.8) to relate the double-traces of φ and ∂ y φ to those involving O 1 and O 2 . The rules of GFF then allow us to go from (5.15) tô For a generic a φ 2 , there is no C D which can make both of these coefficients satisfy the Ward identity. After producing universal bounds in the (a φ 2 , ∆ τ ) and (a φ 2 , C D ) planes, it is natural to try scanning in all three parameters. This means choosing a grid of points in the allowed blue region of figure 4 and maximing the spin 2 gap at each one. The best feature of this plot is that every point with ∆ τ > 3 is guaranteed to obey all the constraints given above: exact OPE relations, no stress tensor, no flux operator, and the Ward identity for the displacement operator. For a φ 2 ≥ 0.20, which is the vicinity of the kink, the results are shown in figure 5.
Before we discuss this figure, let us comment first on the analysis for more general a φ 2 and for which the data is not shown. This analysis indicated that the maximum spin 2 gap, so the points on the boundary of figure 2, correspond to the largest possible values of C D , so the points near the upper boundary of figure 4. On the other hand, near the lower boundary of figure 4 the spin 2 gap remains very close to 3. Since the perturbative line in figure 4 is near this lower boundary, it indicates that the corresponding line in figure 2 must actually be quite a bit flatter than the slope determined by the Ising model central charge (5.3). In short, the (non-rigorous!) extrapolation of the one-loop analysis to small but finite values of δa φ 2 indicates that the Dirichlet boundary condition can only be driven to a weakly coupled fixed point if the 3d CFT that triggers the RG flow has a large central charge.
Let us return to figure 5. Our best candidate for a non-trivial boundary condition, the left of table 2, may be found by hugging the upper edge of figure 4 and looking for where ∆ τ jumps. As suggested by the extremal spectrum, this happens at (a φ 2 , ∆ τ , C D ) = (0.215, 3.966, 0.0050) and corresponds to the red point in the figure. We see that the two-dimensional kinks in previous figures have become a three-dimensional feature: a cliff appears to develop around this point. The spectrum on the right of table 2 corresponds to (a φ 2 , ∆ τ , C D ) = (0.218, 3.970, 0.0050) which appears to be a more generic point in this three parameter space. We originally chose this point by hugging the lower edge of pink region in figure 4, i.e. by bisecting in ∆ τ without the constraint (5.5) for non-displacement scalars of dimension 4. This produces a jump at (a φ 2 , C D ) = (0.218, 0.0044) in that plot. However, re-interpreting the extra dimension 4 scalars found in that solution as discussed above shifted C D from 0.0044 to 0.0050. (As discussed above, this is under the assumptions that the Ward identities hold for this point.) Notice also that ∆ τ → 3 rather smoothly as C D approaches its lowest possible value. According to the dashed line in figure 4 this is where we could find potential weakly coupled fixed points from the Neumann end. As for the Dirichlet end discussed above, one might take this as an indication that the anomalous dimension of the three-dimensional stress tensor cannot be too big.

Outlook
We set out to investigate whether a free real scalar field could have conformal boundary conditions other than Dirichlet or Neumann. The bulk equation of motion restricted the two-and threepoint functions of φ so strongly that we found that all non-trivial boundary conditions must support a shadow pair of boundary operators of dimensions ∆ φ and ∆ φ + 1. The numerical analysis in four bulk dimensions (so three boundary dimensions) of correlation functions of this shadow pair yielded interesting results. On the one hand, for a large range of values of a φ 2 (the one-point function of the bulk φ 2 operator) there must be a spin 2 operator relatively close to the unitarity bound, providing some evidence for the absence of interesting boundary conditions. On the other hand, for a φ 2 near its upper bound of 1/4 this maximal value shoots up and we observed an interesting kink in the data at about a φ 2 = 0.215 with a spin 2 operator of dimension 3.966 and C D approximately equal to 0.0050. More numerical data is provided in table 2. This could be a new conformal boundary condition for the free scalar field.
It the future it would be interesting to see whether the shadow relations can be explored analytically rather than numerically. Indeed, one could ask whether the shadow transform can be applied to four-point functions and conformal blocks? As we have seen, the shadow transformation is singular for three-point functions when the scaling dimension of the third operator is of double-twist type, so it is not entirely obvious that shadow transforming one or more operators in a consistent four-point function automatically leads to another consistent fourpoint function. Our expectation is instead that contact terms will become important because they get magnified to non-trivial functions by the shadow transformation. It would of course also be interesting to understand the possible new conformal boundary condition that corresponds to the kink in our numerical plots. One approach would be to try to 'move' the kink by changing the problem. For example, we could try different spacetime dimensions d or generalize the problem to N > 1 free scalar fields. 18 These would of course be interesting studies in their own right, but if we can dial a parameter like d or N to a value where the kink merges with a free boundary condition then that would provide strong indications for a possible perturbative approach to the problem. Another approach would be to put the free bulk theory in AdS: then we can add a mass term to the bulk field which would continuously change the scaling dimensions but which is not expected to spoil the conformality of the boundary and a conformal bootstrap analysis should always be possible [56].
An obvious direction for the future is to try to extend the analysis of this paper to other examples of free theories in the bulk, such as the free scalar in other spacetime dimensions, the free fermion in any d or the free vector in d = 4. In the latter case it would be extremely interesting to see if there is any signature of the continuous family of boundary conditions [10] in the bootstrap approach, perhaps along the lines of the previous bootstrap study of conformal manifolds in [57].
More generally, the 'landscape' of boundary conditions for a given CFT d is a wide open problem. It therefore remains an interesting target for further explorations. The subject is even richer because, as this paper exemplifies, we need to modify the usual crossing symmetry equations in surprising ways when defects, boundaries, or interfaces are taken into account.
(Simons Collaboration on the non-perturbative bootstrap). LD is partially supported by INFN Iniziativa Specifica ST&FI. LD also acknowledges support by the program "Rita Levi Montalcini" for young researchers. CB has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement #787185). The numerical calculations were done on the University of Oxford Advanced Research Computing (ARC) facility [58] and the SeaWulf high-performance computing system maintained by Stony Brook Research Computing and Cyberinfrastructure and the Institute for Advanced Computational Science at Stony Brook.

A.1 bOPE
Consider a scalar bulk operator O, not necessarily free. The bOPE of O is completely determined by SO(d, 1) symmetry, up to a certain collection of CFT data [5,6] O One can check that the expression above reproduces the bulk-to-boundary correlators (2.1), once applied to the boundary two-point functions We will take unit-normalized boundary two-point functions, except for the protected operators that can appear in the bOPE of the bulk conserved currents J . Such operators, collectively denoted by J (l) (with l = 0, . . . − 1) have their normalization fixed by the Ward identities (2.23), and therefore the coefficients in their two-point functions are physical Conventionally, we choose unit normalization for the boundary modes of φ

A.2 Boundary OPE and physical OPE coefficients
We denote generic boundary operators as O k , where the label k indicates collectively all possible indices of the local operator. The OPE between two boundary operators O i is (up to boundary descendants) The boundary two-point functions are normalized as in (A.2). We use the Zamolodchikov metric C O O to raise and lower indices off ij k s. Concretely (sum over repeated indices) With these conventions, we have that the displacement operator, whose normalization is taken as in (A.3), enters the generic boundary OPE (A.6) as and a generic boundary four-point function as In the equation above we introduced the conformal blocks, which are normalized as (5.4). Alternatively we can think of D to be unit-normalized, such that the physical boundary OPE coefficient is Similar remarks apply for other protected operators that can appear in the bOPE of the bulk conserved currents J .

B Three-point function conformal blocks
In this appendix we derive the conformal block decomposition of the free field φ three-point function with a generic boundary operator O (l) We will obtain closed-form expressions for all the conformal blocks exchanged in the boundary channel of this three-point function. We also compute some bulk channel blocks, while leaving a more complete study for the future.

B.1 Blocks in the boundary channel
We start from the blocks in the boundary channel. As we explained in the main text (see section 2.2), the expansion of the correlator (B.1) in a basis of boundary conformal blocks can be obtained by acting twice with the bOPE on the generic three-point functions and then resumming the contributions from boundary descendants. The polynomials P (l) are defined in (2.18). Applying the bOPE (2.4) and using the identity ∇ 2n x 12 we can rewrite (B.1) as In the above formula we introduced The infinite sum in the r.h.s. of (B.4) can be explicitly performed, and the result takes the form The quantitiesF ij ∆,l are hypergeometric functions of the cross-ratios w ± (defined in (2.21)) Note that in terms of two cross-ratios the cross-ratios w ± can be rewritten as (B.9) We have checked that (B.6) satisfies the Klein-Gordon equation with the correct conditions. As a further consistency check, we note that the defect channel blocks for the two-point function (2.8) can be recovered from (B.7) by setting ∆ = l = 0 andf ij1 = δ ij .

B.2 Scalar blocks in the bulk channel
Next, we will be interested in the bulk conformal block expansion of (B.1). For simplicity we will consider only the case where the third operator is a boundary scalar, while leaving the general case for future work. In the bulk channel we plug the φ × φ OPE (2.9) to convert (B.1) into an infinite sum over bulk-to-boundary two-point functions with the ellipsis denoting contributions from bulk descendants, which are fixed by SO(d + 1, 1) conformal symmetry. As discussed in the main text (see subsection 2.2.2), spin selection rules and current conservation imply that the bulk operator φ 2 is the only contribution to the r.h.s of (B.10) for generic ∆ not equal to the scaling dimension of J . In this more generic case we have by plugging the bulk OPE into (2.15) and resumming the bulk descendants. Using the explicit form of the differential operator that controls the scalar exchange (see e.g. [59]), we find the following series expansion where the cross-ratio ξ is defined in (2.6). There are various interesting special situations in which the result (B.12) produces simple closed-form formulae. In the 'cylindrical' configuration y 1 = y 2 = y the infinite sum gives a simple hypergeometric function where we introduced a cross-ratioχχ = | x 12 | 2 y 2 , (B.14) which is nothing but the restriction of ξ defined in (2.6) to the 'cylindrical' configuration. As explained in appendix C.2, this result can be also derived by 'inverting' the boundary channel expansion (B.6). As a final comment, we note that the series representation (B.12) yields simple closed-form expressions, some of which will be presented in appendix C.3, when the third operator is of D (0) type.

C OPE relations and bulk-to-boundary crossing
In this appendix we discuss in detail the derivation of the main results presented in section 2.

C.1 Derivation of the OPE relations
In this appendix we derive the OPE relations (2.22). The starting point is the boundary channel expansion for the correlator (B.1), given in (2.20). Away from other operator insertions, the φ×φ OPE requires this three-point function to be analytic around x µ 1 = x µ 2 (recall that the identity in (2.9) decouples). In order to study this limit, it is convenient to place the two bulk operators at the same transverse distance i.e. y 1 = y 2 = y, such that the expression (2.20) simplifies as follows: whereχ is the cross-ratio defined in (B.14). In this configuration with y 1 = y 2 Bose symmetry (2.19) implies that this expression vanishes when l is an odd integer, so we first consider even l. We then require (C.1) to be analytic around x 1 = x 2 , for finite y. For generic values of d, l, ∆, the r.h.s. of (C.1) contains unphysical singularities, since up to higher powers ofχ. Such unphysical singularities cancel from the r.h.s. of (C.1) precisely when the OPE relations (2.22) are satisfied, such that the analytic result at y 1 = y 2 is When ∆ approaches some special integer dimensions some of the boundary blocks in the r.h.s. of (C.1) are themselves regular and (C.2) is not valid. This can happen when: We then analyse these special cases separately. Requiring the cancellation of any residual singularity on the r.h.s. of (C.1), will again impose certain relations between the boundary OPE coefficients. It is reassuring to see that these special cases are captured by the appropriate limits of the general result (2.22). We now discuss the case when l is an odd integer. Starting from (2.20), we need to setf 11 O (l) = f 22 O (l) = 0 (as dictated by Bose symmetry) so that the three-point function is proportional tô f 12 O (l) . We then study analyticity around x 1 = x 2 for finite y 12 ≡ y 1 − y 2 . For generic values of d, l, ∆ this correlator features unphysical singularities, since up to subleading terms. Because of the first term in the above expression, which is singular for d ≥ 3 (for d = 3 and l = 0 the singularity is logarithmic), for generic ∆ we must setf 12 O (l) = 0.
On the other hand, the boundary blocks are themselves regular and the (C.7) is not valid when the dimension of O equals Again we see that the relations (2.22), together with the constraints from Bose symmetry (2.19), promptly capture these special cases. The analytic correlator (2.20) on the special dimensions (C.8) then reads

C.2 Matching with the bulk
Owing the results from the previous subsection, we are now ready to discuss the consequences of the bulk-boundary crossing symmetry for the three-point functions (B.1). The first step is to derive the leading terms in the bulk channel expansion of the correlator (B.1). To this end, recall that the φ × φ OPE (2.9) contains a scalar φ 2 and infinitely many conserved currents J , with ∈ 2N and ∆ = d + − 2. We now plug the φ × φ OPE into (B.1), impose the selection rules from conservation in order to figure out which bulk primary can couple to O (l) and finally compare to the boundary channel expansion. We conclude that: • When ∆ = d + − 2 and l is odd, O (l) cannot couple to any bulk operator in the φ × φ, and the three-point function must vanish. This perfectly matches with the expectations from the boundary channel.
• When ∆ = d + − 2 and l is even, O (l) can only couple to a spin l bulk current J l (or to φ 2 when l = 0). This is consistent with what we expect from the the boundary channel, where we are left with only one unknown OPE coefficientf 12 O (l) . In either case, from the leading bulk OPE we have and, after comparing to (C.3) we find The result for a scalar (l = 0) operator O is simply obtained from the former by setting We can use the result above in order to re-interpret the expression for (B.1) obtained using the boundary OPE in terms of the bulk OPE, and compute the corresponding bulk block. This procedure is unambiguous, since in both channels there is just one undetermined OPE coefficient. In practice, we solve (C.12) forf 12 O (l) and plug the result into (C.3). We find As a consistency check, note that for l = 0 the above expression reproduces the block W φφ O φ 2 ( x 12 , y, y), which was computed explicitly in (B.13). The same logic can be applied to compute the bulk blocks at generic transverse positions y 1 , y 2 , starting from the boundary channel decomposition (2.20).
• When ∆ = d+ −2 and l is even there are two cases. For > l, the primary O (l) can couple to both J l and J . The number of undetermined bulk OPE coefficients then matches that of the boundary ones (f 11 O (l) andf 22 O (l) ). As an example, in section C.3 we explicitly solve the bulk-to-boundary bootstrap for the case of = 2 with l = 0, but similar results can be obtained for the more general case of D (l) . When = l the operator O (l) can only couple to J l , and this matches with the number of undetermined boundary OPE coefficients (f 11 O (l) ).
• When ∆ = d + − 2 and l is odd the only possible bulk contribution comes from the spin currents J . From the leading bulk OPE limit at | x 12 | = 0 we have φ( x 1 , y 1 )φ( x 2 , y 2 )V (l) (θ, ∞) ∼ c φφ C J b J V (l) (y 12 ) −l ( x 12 · θ) l + . . . , (C.14) where y 12 = y 1 − y 2 . So, after comparing to (C.9) (note that − l = 2n + 1) we find (C. 15) The Ward identity (2.23) further relates the coefficient b J V (l) to the coefficient in the twopoint function of V (l) , e.g. for the flux operator V 2 ≡ V (1) the coefficient in eq. (2.31). In the case of V (1) , upon plugging the value of c φφT and C T given in eq. (2.10), the result (C.15) gives precisely the first equality of eq. (2.32). The second equality is obtained upon using that b T V (1) = 2 C V (1) , as dictated by the Ward identity (2.23).

(C.18)
The second term is the contribution from the bulk stress tensor and reads [5,6] T µν (x)D(∞) =b T D δ µy δ νy − 1 Note bulk descendant operators of T µν do not enter into (C.17), since (C. 19) is a constant. One can further use the Ward identities for the displacement operator [5,6] to relate the bOPE coefficient b φ 2 D to the one-point function of φ 2 : (C.20) (C.21) The final formula (2.36) is obtained by plugging in the above expression the values (2.10) of c φφT and C T corresponding to a d-dimensional free scalar field. It is pleasant to see that the final result (2.36) is consistent with the Ward identity [5] d d−1 x φ(x 1 )φ(x 2 )D( x) = (∂ y 1 + ∂ y 2 ) φ(x 1 )φ(x 2 ) . For numerical purposes it is convenient to isolate, in the crossing equations above, the contributions from the primaries with fixed dimensions from those in the continuum. Such special primaries are the identity 1 (for which (2.22) impliesf 121 = 0 and we choose the normalization f ii 1 ≡f ii = 1), as well as the boundary modes of the bulk higher-spin currents, D (l) and V . (D.7)