Integrable Boundary for Quad-Graph Systems: Three-Dimensional Boundary Consistency

We propose the notion of integrable boundary in the context of discrete integrable systems on quad-graphs. The equation characterizing the boundary must satisfy a compatibility equation with the one characterizing the bulk that we called the three-dimensional (3D) boundary consistency. In comparison to the usual 3D consistency condition which is linked to a cube, our 3D boundary consistency condition lives on a half of a rhombic dodecahedron. The We provide a list of integrable boundaries associated to each quad-graph equation of the classification obtained by Adler, Bobenko and Suris. Then, the use of the term"integrable boundary"is justified by the facts that there are B\"acklund transformations and a zero curvature representation for systems with boundary satisfying our condition. We discuss the three-leg form of boundary equations, obtain associated discrete Toda-type models with boundary and recover previous results as particular cases. Finally, the connection between the 3D boundary consistency and the set-theoretical reflection equation is established.


Introduction
Discrete integrable systems arise from various motivations in applied or pure mathematics like the need to preserve integrability of certains continuous equations when performing numerical (and hence discrete) simulations or the theory of discrete differential geometry. In previous work concerning an important class of such systems known as integrable quad-graph equations [11], one of the original motivations was to study discrete differential geometry of surfaces. In this context, a general construction allows one to always obtain a discretization of the surface in terms of quadrilaterals [11], at least as long as one is only concerned with the bulk of the surface and does not worry about its boundary (if it has one). The vertices of the graph thus obtained can be seen as discrete space-time points where the field is attached. The dynamics of the field is then specified by an equation of motion involving the values of the field at the four vertices forming an elementary quadrilateral. Typically, this is of the form Q(u 00 , u 10 , u 01 , u 11 ; a, b) = 0, where u 00 , u 10 , u 01 , u 11 are the values of the field at the vertices of the quadrilateral and a, b are parameters (see l.h.s. of the Fig. 3).
There exist several integrability criteria that different authors use to characterize the notion of discrete integrability. Let us mention for instance algebraic entropy [10], singularity confinement [19] or the 3D consistency/consistency around the cube condition [11,26]. The latter is deeply related to the notion of discrete Lax pair, discrete Bäcklund transformations and other classical notions of integrability and, combined with a few other assumptions, led to the important ABS classification of quad-graph equations [2]. This fact and the similarity of the above structures with those existing in the case of continuous integrable systems make it a popular criterion.
In this paper, we want to propose a way of defining an integrable discrete system on a quadgraph with boundary as arising from the discretization of a surface while taking into account its boundary. From the geometric point of view, this is a natural generalization of the discrete differential geometry motivation to the case of a surface with boundary. From the point of view of discrete integrable quad-graph equations, this then allows us to tackle the problem of formulating the analog of the 3D consistency condition together with its consequences, i.e. Bäcklund transformations and zero curvature formulation. We also introduce Toda-type systems with boundary through the three-leg form of integrable equations on quad-graphs and we recover the previous approach to boundary conditions for discrete integrable systems presented in [20].
We give a precise definition of the discretization procedure and of a quad-graph with boundary in Section 2.1. Then we show how to define a discrete system on it. In addition to the elementary quadrilateral and the corresponding equation of motion (1), the new crucial ingredient is an elementary triangle together with the corresponding boundary equation of the form q(x, y, z; a) = 0, where x, y, z are the values of the field at the vertices of the triangle, and a function σ determining the effect of the boundary on the labelling (see r.h.s. of the Fig. 3). We provide our definition of integrability in this context by defining the 3D boundary consistency condition involving Q, q and σ in Section 2.3 and then go on to present a method that allows us to find solutions for q and σ for a given Q in Section 3. In Section 4, we justify further our introduction of the 3D boundary consistency condition by discussing Bäcklund transformations and the zero curvature representation in the presence of a boundary. Section 5 introduces the three-leg form for boundary equations and we use this to define systems of Toda-type with boundary. As an example, we recover as particular case the approach to boundary conditions of [20] in (fully) discrete integrable systems. Finally, in Section 6 we establish the connection between the 3D boundary consistency condition introduced in this paper and the set-theoretical reflection equation [14,15] in the same spirit [28] as the 3D consistency condition is related to the set-theoretical Yang-Baxter equation [17]. Conclusions and outlooks are collected in the last section.

Integrable quad-graph systems with boundary
In this section, we first define the notion of quad-graph with boundary and then use it to define the elementary blocks needed to study integrable equations on quad-graphs with boundary.

Quad-graph with boundary
Our starting point is the definition of a quad-graph from a cellular decomposition of an oriented surface S containing only quadrilateral faces. As explained in [11], a quad-graph can always be obtained from an arbitrary cellular decomposition G by forming the double D (in the sense of [11], otherwise called a diamond in [25]) of G and its dual cellular decomposition G * . So far, these notions apply to surfaces without a boundary. For the case of a surface with boundary, Figure 1. Example of a cellular decomposition G and of the associated G * . The "horizontal" curved line is the boundary of the underlying surface S. The black dots and the straight lines connecting them are the vertices and edges of the initial cellular decomposition G. The white dots and the straight lines connecting them are the vertices and edges of G * . Edges of G on the boundary and the boundary itself are identified in this picture. the notion of dual cellular decomposition does not exist in general. In [25], Definition 1 gives the definition of an object Γ * associated to a cellular decomposition Γ of a compact surface with boundary and it is noted that Γ * is not a cellular decomposition of the surface. Nevertheless, in Section C.3 of [24], a notion of double is given for a surface with boundary (and not necessarily compact, like a half-plane) and we adapt it here for our purposes. In particular, we will see that the generic structure that arises is what we call a quad-graph with boundary with its faces being either quadrilaterals or "half quadrilaterals", i.e. triangles.
So let us consider a cellular decomposition G of our surface S with a boundary. We denote, respectively, by F , E and V the set of faces, edges and vertices of this cellular decomposition. Following [24] and adapting slightly, we define G * as the following collection of cells with F * , E * and V * respectively the sets of faces, edges and vertices.
• To each face in F , we associate a vertex v * in V * (called white), placed inside the face.
• To each edge e ∈ E not on the boundary, we associate the dual edge e * ∈ E * which cuts it tranversally and forms a path between the two white vertices contained in the two adjacents faces in F separated by e.
• To each vertex v ∈ V (called black) not on the boundary, we associate the face in F * that contains it, i.e. the face whose boundary is made of the edges in E * that cross the edges in E having the vertex v under consideration as one of their ends.
Compared to Definition 6 of [24], in our definition of G * , we include neither the dual edge corresponding to an edge on the boundary nor the additional white vertex in the middle of an edge in E belonging to the boundary. A typical configuration of G and G * is shown in Fig. 1. We now form the structure that we will call a quad-graph with boundary below. Let us denote it D. The vertices of D are all the black and white vertices, i.e. V D = V ∪ V * . The edges of D are all the edges in E that lie on the boundary of S together with those edges (w, b) obtained by connecting a white vertex to each of the black vertices sitting on the face that contains the white vertex 1 . The faces are then taken to be the interiors of the "polygons" thus obtained. The result of this procedure for G and G * as in Fig. 1 is shown in Fig. 2. The faces are therefore of two types: • Quadrilaterals, with two black and two white vertices, in the bulk.
• Triangles, with two black vertices (on the boundary) and one white vertex (inside S), alongside the boundary.

Definition 1.
A quad-graph with boundary is the collection of vertices and edges of D obtained as described above from a cellular decomposition of a surface with boundary.
As an important by-product, this procedure to get a quad-graph with boundary allows one to obtain a bipartite graph 2 . In addition, all the vertices on the boundary are of the same type. This property will be important in the construction of the Toda models in Section 5.

Discrete equations on quad-graph with boundary
We are now ready to define discrete equations on a quad-graph with boundary. As usual, we associate a field to this quad-graph (i.e. a function from V ∪ V * , the set of vertices, to C) and a constraint between the values of the field around a face. This constraint can be seen as the equation of motion for the field. For each quadrilateral face, the constraint is, as usual, defined by Q(u 00 , u 10 , u 01 , u 11 ; a, b) = 0, where u 00 , u 10 , u 01 , u 11 ∈ C are the values of the field at each vertex around the face and a, b ∈ C are parameters associated to opposite edges. One usually represents this equation as on the l.h.s. of Fig. 3. We will refer to equation (2) as describing the bulk dynamics. Sometimes, one demands additional properties [11] for the function Q(u 00 , u 10 , u 01 , u 11 ; a, b) such as linearity on each variable u 00 , u 10 , u 01 , u 11 (affine-linearity), symmetry by exchange of these variables (D 4 -symmetry), the tetrahedron property or the existence of the three-leg form (see Section 5).
The new elementary building block needed to define a discrete system on a quad-graph with boundary is an equation of the following type defined on each triangular face q(x, y, z; a) = 0,  where x, z are values of the field on the boundary, y a value inside the surface and a is a parameter associated to one edge (the other edge is associated to a function σ(a) of the parameter a). We represent this equation as on the r.h.s. of Fig. 3 where the dashed line represents the edge on the boundary of the surface. We will refer to this equation as describing the boundary dynamics. In general, there is no special requirement on q or σ but we will see that our methods to construct expressions for q result in q having certain properties: • Linearity: the function q(x, y, z; a) is linear in the variables x and z. Let us emphasize that it may not be linear in y.
• In some cases, a three-leg form for q inherited from the three-leg form of the corresponding bulk equation Q = 0. We go back to this point in Section 5.
Let us remark that for a given quad-graph with boundary (and for σ being not the identity), it is not always possible to label the edges with parameters as prescribed above. We show on the l.h.s. of Fig. 4 a quad-graph with boundary for which we cannot find a suitable labelling. On the r.h.s., we show a case where this is possible. In the following, we consider only those quad-graph with boundary that can be labelled. It would be very interesting to study this combinatorial problem in general but this goes beyond the scope of this paper.

Integrability: the 3D boundary consistency condition
As explained in the introduction, we adopt the 3D consistency approach to integrable quadgraph equations. Let us first recall what it means in the bulk case [11,26] before we introduce We propose, in the following, the main new equation of this article which is a similar consistency condition for the function q (and σ) that we call the 3D boundary consistency condition. This condition is in fact a compatibility condition between the bulk equation Q = 0 and the boundary equation q = 0. Instead of the cube for the 3D consistency, the 3D boundary consistency lies on a half of a rhombic dodecahedron as displayed in r.h.s. of Fig. 5. Let us remark that this half of a rhombic dodecahedron has the same combinatorial structure than the r.h.s. of the Fig. 4. On each face (resp. 4 quadrilaterals and 4 triangles), we attach its corresponding equation of motion (resp. Q = 0 and q = 0). Then, the 3D boundary consistency is the statement that, given the values x, x 1 and x 2 of the field, the different schemes to compute w (see Fig. 5) lead to the same results. More precisely, following the notations in Fig. 5, we get that, given x, x 1 and x 2 : • the three equations gives respectively the values of y 1 , y 2 and y 3 .
• Then, both equations gives respectively the values of z 1 and z 2 .
Definition 2. We say that we solve the 3D boundary consistency for q if, given a 3D consistent Q, we find a function q of x, y, z, a and a function σ of a such that the scheme (4)- (6) gives a unique value for w given values for x, x 1 and x 2 . In this case, q is called a solution of the 3D boundary consistency condition (we omit explicit reference to σ which is taken as part of the solution). We also say that q is compatible with Q.
Note that this should not be confused with the notion of a solution of the actual dynamics described by Q and q. Such a solution would consists in finding an expression for the field u at each vertex of the quad-graph that satisfies the bulk and boundary dynamics. This is an exciting open question which deserves separate attention and is beyond the scope of the present paper.
In the bulk case, when one finds a Q which satisfies the 3D consistency (see Fig. 5), the associated system is called an integrable equation on quad-graph [2,11,26]. To introduce a boundary which preserves the integrability, for a given Q, we must solve the 3D boundary consistency condition, i.e. find compatible functions q and σ. Definition 3. We call integrable equation on a quad-graph with boundary the data of a quadgraph with boundary with compatible labelling as well as functions Q, q and σ which satisfy the 3D consistency and the 3D boundary consistency conditions. We justify the adjective "integrable" in Sections 4.1 and 4.2 by showing the presence of Bäcklund transformations and a zero curvature representation.
Let us emphasize that similar approaches have already appeared in the literature to introduce integrable boundaries in different contexts. Indeed similar figures to those in Fig. 5 appeared in [4,9,18] as the face representation of the reflection equation [16,30]. The right-hand-side of Fig. 5 is also half of the figure representing the tetrahedron equation [7,8]. There exists also a close connection between the set-theoretical equation introduced recently in [14,15] and the 3D boundary consistency (see Section 6 for more details). Similar connections have been studied previously in the bulk case where the set-theoretical Yang-Baxter equation is linked to the 3D consistency condition [2,28].

Solutions of the 3D boundary consistency
In this section, we provide a list of solutions of the 3D boundary consistency condition associated to the bulk equations Q classified in [2]. In other words, given Q, we want to find solutions q and σ which satisfy the 3D boundary consistency conditions. The underlying idea is that this should provide integrable boundary conditions for integrable discrete equations characterized by Q.

The ABS classif ication
For completeness, we list the solutions of the 3D consistency equations (2) obtained in [2] (Q1) : a(u 00 − u 01 )(u 10 : sn(a)(u 00 u 10 + u 01 u 11 ) − sn(b)(u 00 u 01 + u 10 u 11 ) − sn(a − b)(u 00 u 11 + u 10 u 01 ) We use the same labels (Q, H, A families) and forms that were used in [2], except for equation (Q4) which is in an equivalent form introduced in [5,21] where sn(a) = sn(a; K) is the Jacobi elliptic function with modulus K. It is worth noting that the 3D consistency condition as well as the affine-linearity, D 4 -symmetry and tetrahedron properties are preserved for each of the equations up to common Möbius transformations on the variables u 00 , u 10 , u 01 , u 11 and point transformations on the parameters a, b.

Method and solutions
Instead of performing brute force calculations where one could make assumption on the form of q (like multilinearity in the variables) and insert in the 3D boundary consistency condition for a given Q, below we describe a simple method that amounts to "fold", in a certain sense, Q to obtain compatible q's. Although it may seem ad hoc and arbitrary, this method has several motivations. First, it gives a simple way to obtain three-leg forms for q knowing the ones for Q and hence a way to introduce discrete Toda-type systems with boundary. This is explained in Section 5. Second, the method is a simple adaptation of the folding method that was used in [14] to obtain reflection maps, i.e. solutions of the set-theoretical reflection equation. This last point is discussed in more detail in Section 6. There, we present an alternative method to find admissible q's starting from reflection maps. This alternative method produces some of the solutions that are also obtained with the method that we explain here. But more importantly, this alternative method establishes a deep connection between solutions of the 3D boundary consistency condition and reflection maps. This is reminiscent of the deep connection between solutions Q of the 3D consistency condition (in particular the ones of the ABS classification) and quadrirational Yang-Baxter maps [3]. Now, given Q satisfying the 3D consistency condition, we look for q(u 00 , u 10 , u 11 ; a) of the form q(u 00 , u 10 , u 11 ; a) = Q(u 00 , u 10 , k(u 00 , u 10 ; a), u 11 ; a, σ(a)), satisfying the 3D boundary consistency condition, where k and σ are the functions to be determined. Equation (7) may be seen as the folding along the diagonal (u 00 , u 11 ) of the quadrilateral in Fig. 3 to get the triangle (u 00 , u 10 , u 11 ). Obviously, one may fold along the other diagonal but, due to the D 4 -symmetry of Q, we get the same results.
To find k and σ (and hence q), for each given Q, we insert our ansatz (7) in the scheme (4)-(6) and try to find functions k and σ that fulfill the resulting constraints. We note that the following "trivial" choice σ(a) = a, k(u 00 , u 10 ; a) = −u 10 , is always a solution of the problem for any Q. It yields q(u 00 , u 10 , u 11 ; a) = au 10 (u 00 − u 11 ). We report in Tables 1 and 2 the nontrivial functions k, σ and q we found for the different Q's of the Q, H and A families of the ABS classification. Note that we make no claim of completeness. A method for a systematic classification is in fact an interesting open problem. Table 1. Results for boundary equations (Q family). µ is a free parameter. The asterisk denotes solutions that are also obtained with the method of Section 6.
ABS σ(a) q(x, y, z; a) k(x, u; a) Other aspects of integrable equations on quad-graphs with boundary In this section, we present results on important traditional aspects of integrability obtained from the 3D boundary consistency equation. They justify a posteriori our definition of integrability for quad-graphs with boundary.

Bäcklund transformations
In this subsection, we prove that the 3D boundary consistency condition proposed in the previous section leads naturally to Bäcklund transformations which are a basic tool in the context of classical integrability and soliton theory. This result is similar to the one without boundary [2] and is summarized in the following proposition: Proposition 1. Let us suppose that we have an integrable equation on a quad-graph with boundary (with the set of all its vertices denoted as V ) as well as a solution, g : V → C. There exist a two-parameter solution g + of the same integrable quad-graph equation and a function f from V Table 2. Results for boundary equations (H and A families). µ is a free parameter. The asterisk denotes solutions that are also obtained with the method of Section 6.
ABS σ(a) q(x, y, z; a) k(x, u; a) for all edges (v, v 1 ) of the quad-graph not on the boundary of the surface (a is the parameter associated to this edge), and for all vertices v 2 on the boundary of the quad-graph. We call the solution g + the Bäcklund transform of g.
Proof . The proof follows the same lines as in the case without boundary. We start with the quad-graph with boundary, called in this proof the ground floor. We consider also two other copies of the surface, called first and second floor, of the ground floor. Then, we construct a 3D graph by the following procedure (see also Fig. 6): • There is a one-to-one correspondence between the vertices of the ground floor and the ones of the first floor but one moved the vertices of the first floor such that no vertex of the first floor lies on the boundary of the surface (see l.h.s. of the Fig. 6). The vertices on the second floor is an exact copy of the ones on the ground floor.
• We copy the edges in the bulk of the ground floor on the first and second floor. The copies carry the same label. We copy the edges on the boundary of the ground floor only on the second floor.
• We add the edges (thin lines on the Fig. 6) linking all the vertices of the ground floor with the corresponding vertices of the first floor and similarly from the first to the second floor. The edges between the ground and first floor carry the label λ whereas the edges between the first and second floor carry σ(λ). We add also the edges between the vertices on the boundary of the ground floor to the corresponding ones on the second floor.
• the set of faces is the union of the following sets: (i) the triangular and quadrilateral faces of the ground and second floor; (ii) the quadrilateral faces of the first floor; (iii) the "vertical" quadrilateral faces made from the edges of the ground floor, the corresponding ones of the first floor and the vertical edges linking the vertices of these edges, and similarly between the first and second floors; (iv) the "vertical" triangular faces made of the edges linking the vertices on the boundary of the ground and second floors and the corresponding vertex of the first floor (which is not on the boundary).
We impose now that the values of the field living on this 3D graph be constrained by Q = 0 on each quadrilateral face and by q = 0 on each triangular faces.
As in the case without boundary, due to the 3D consistency condition, given a function g satisfying the constraints on the ground floor, one can get a function f satisfying the constraint on the first floor depending on λ and on the value of f at one vertex of the first floor. This function f satisfies, in particular, Q(g(v), g(v 1 ), f (v), f (v 1 ); a, λ) = 0.
Knowing the value of the field at one vertex on the boundary of the ground floor (say g(v 1 ) on the figure) and the corresponding value on the first floor (f (v 1 )), we get the value g + (v 1 ) of the field g + at the corresponding vertex of the second floor using the equation q(g(v 1 ), f (v 1 ), g + (v 1 ), λ) = 0 which connects the white, grey and black copies of the vertex v 1 .
We obtain a function g + satisfying the constraints of the second floor using Q on the "vertical" quadrilateral faces between the first and second floors. The important point is to remark that due to the 3D consistency and the 3D boundary consistency conditions, all the different ways to obtain the values of g + give the same result. We may see this geometrically since the only "elementary blocks" of the 3D graph are cubes or half of rhombic dodecahedron and since the functions Q and q are chosen such that they satisfy the 3D consistency and the 3D boundary consistency conditions, the different ways to compute the values of the field g + are consistent.
Finally, using the fact that the second floor is an exact copy of the ground floor, g + satisfies the constraints of the original quad-graph equations.
Let us remark that the main difference between the cases with and without boundary lies in the necessity of an additional, intermediate floor in the case with boundary. This feature appeared previously in the context of asymmetric quad-graph equations [13]. Let us emphasize that the function f defined on the intermediate floor is not in general a solution of the same equations on a quad-graph with boundary since the values of the field on the vertices of the "would-be" boundary of the first floor do not necessarily satisfy an equation of the type (3).

Zero curvature representation
It is well established that a 3D consistent system Q(u 00 , u 10 , u 01 , u 11 ; a, b) admits a zero curvature representation [11,12,26], i.e. there exists a matrix L depending on values of the field on the same edge, the parameter associated to this edge and a spectral parameter λ such that the following equation holds 3 L(u 11 , u 10 , b; λ) L(u 10 , u 00 , a; λ) = L(u 11 , u 01 , a; λ) L(u 01 , u 00 , b; λ).
There exists a constructive way to get L from Q [11,26]: due to the linearity of the function Q, we can rewrite equivalently equation (2) as follows where L is a 2 by 2 matrix describing, with usual notations, a Möbius transformation Geometrically, using Fig. 5, it is easy to show that the matrix (10) satisfies the zero curvature equation (9) if Q satisfies the 3D consistency [11,12,26]. Similarly, we want to show that the 3D boundary consistent system admits a zero curvature representation. For the boundary equation q(x, y, z; a) = 0, we propose the following zero curvature representation K(z; c)L(z, y, σ(a); c)L(y, x, a; c) = L(z, y, σ(a); σ(c))L(y, x, a; σ(c))K(x; c), where K is also a 2 by 2 matrix. We can now state the following results justifying the previous definition: Proof . We give the details of the proof for the case given in the first row of the Table 1, i.e. we deal with the bulk equation Q1 δ=0 given by Q(u 00 , u 10 , u 01 , u 11 ; a, b) = a(u 00 − u 01 )(u 10 − u 11 ) − b(u 00 − u 10 )(u 01 − u 11 ), and the boundary equation characterized by The matrix L associated to this Q is given by It is a known result that equation (9) with this choice for L is satisfied if Q(u 00 , u 10 , u 01 , u 11 ; a, b) = 0 but it is easily verified. Let us mention that the parameters entering in the square roots of the normalisation of L may be negative. Therefore, one must choose a branch cut for the square root appearing in the normalisation: we choose the half-line {ix | x < 0}. The matrix K associated to the function k given by (12) is By algebraic computation, one gets K(z; c) L(z, y, σ(a); c) L(y, x, a; c) − L(z, y, σ(a); σ(c)) L(y, x, a; σ(c)) K(x; c) Therefore, if q(x, y, z; a) = 0, relation (11) holds. That proves the proposition for this case. All the other cases are treated similarly which finishes the proof of the proposition.
Note that equation (11) provides a rather general framework for the representation of an integrable boundary in quad-graph models. In the next section, we show how it contains as a particular case a previous approach to boundary conditions in fully discrete systems.

.1 Three-leg form
It is known that any quad-graph equation Q(x, u, v, y; a, b) = 0 of the ABS classification can be written equivalently in the so-called three-leg form [11], either in an additive form, or a multiplicative form, As demonstrated in [11,12], the existence of a three-leg form leads to discrete systems of Todatype [1]. Indeed, let x be a common vertex of the n faces (x, y k , x k , y k+1 ) (with k = 1, 2, . . . , n and y n+1 = y 1 ) with the parameters a k assigned to the edge (x, y k ). On each face, there is the equation Q(x, y k , y k+1 , x k ; a k , a k+1 ) = 0 written in the presentation (13) or (14). Then, summing the corresponding n equations of type (13), one gets n k=1 φ(x, x k ; a k , a k+1 ) = 0, where a n+1 = a 1 . Similarly, multiplying n equations of the type (14) leads to n k=1 φ(x, x k ; a k , a k+1 ) = 1, where a n+1 = a 1 . When the graph is bi-partite (for example with black and white vertices), we can reproduce this procedure by taking as common vertices all the black vertices to get a Toda-type model on the black subgraph. Now assume that the boundary equation q(x, y, z; a) = 0 can be written as ψ(y, x; a) − ψ(y, z; σ(a)) = ϕ(y; a) or ψ(y, x; a)/ψ(y, z; σ(a)) = ϕ(y; a), where the function ψ is the same as in the bulk case and the new function ϕ depends only on the central vertex of the triangle representing the boundary and on the parameter a. In this case, one can obtain systems of Toda-type with boundary. Indeed, let x be a vertex close to the boundary (i.e. belonging to a triangle but not sitting on the boundary) and common to n − 1 quadrilateral faces (x, y k , x k , y k+1 ) (with k = 1, 2, . . . , n − 1). As in the bulk case, on each quadrilateral face, there is the equation Q(x, y k , y k+1 , x k ; a k , a k+1 ) = 0 written in the presentation (13) or (14). On the triangular face, it holds that q(y n , x, y 1 ; a n ) = 0 (with a 1 = σ(a n )). Then, summing (or multiplying) the corresponding n equations, we get in the additive case ϕ(x; a n ) + n−1 k=1 φ(x, x k ; a k , a k+1 ) = 0, and in the multiplicative case ϕ(x; a n ) n−1 k=1 φ(x, x k ; a k , a k+1 ) = 1.
We illustrate this procedure schematically in Fig. 7 for n = 4. As in the bulk case, if the graph is bi-partite (see footnote 2) and moreover if the vertices on the boundary are of the same type, one can reproduce the above procedure on the whole graph to get a Toda-type model on a subgraph and with a boundary determined by ϕ. The conditions on the graph seem, at first glance, very restricting. Nevertheless, the procedure explained in Section 2.1 provides, starting from any graph, examples satisfying such constraints.

Boundary conditions for Toda-type systems
In general, it seems that not all the solutions for q found previously can be written as in (15). However, it turns out that for each solution q(x, y, z; a) for which the corresponding function k(x, y; a) depends only on y or the function q(x, y, z; a) factorizes as f (y, a)g(x, z, a), there is a way to obtain ϕ from φ. This is reminiscent, at the three-leg form level, of the simple folding procedure (7) used to obtain q from Q. In these cases, the function k(x, y; a) plays the role of the cut-off constraint (or boundary condition) for the Toda-type model in the sense of [20] (see Example 3 below for a precise connection).
For the first case where the function k(x, u; a) = k(u; a) does not depend on x, using equation (7) in (13) or (14) together with the D 4 -symmetry property for Q, one can see that equation (15) with the following function ϕ ϕ(y; a) = φ(y, k(y; a); a, σ(a)), is equivalent to the corresponding boundary equation q(x, y, z; a) = 0. Therefore, by using the explicit forms of φ given in [2] associated to each Q of the ABS classification and the results of Tables 1 and 2, we can derive ϕ and hence three-leg forms of the boundary equation in the form (15). In turn, this allows us to define Toda-type systems with a boundary.
Example 1. The trivial solution (8) always corresponds to ϕ(y; a) = 0 for the additive case or ϕ(y; a) = 1 for the multiplicative case. This boundary condition may be interpreted as a free boundary for the corresponding Toda-type system.
Example 2. We recall that ψ(x, u; a) = a x−u and φ(x, y; a, b) = ψ(x, y; a−b) for Q1 δ=0 . Looking at the fourth solution for Q1 δ=0 of Table 1, i.e. σ(a) = −a + 2µ and k(x, y; a) = −y, we obtain ϕ using (16) and the following discrete Toda-type system with a boundary term: For the second case where the function q(x, y, z; a) factorises, the construction is a bit more involved. The equation q(x, y, z; a) = 0 constrains x and z independently of the values of y. Therefore, we get equations involving only the values of the field on the boundary. Let us suppose we solve these constraints on the boundary and denote byx the corresponding solution for the field x. These values on the boundary play the role of parameters in the function k(x, y; a) appearing in the boundary conditions. Also, one can see that equation (15) involves the following function ϕx ϕx(y; a) = φ y, k(x, y; a); a, σ(a) , and is equivalent to the corresponding boundary equation q(x, y,z; a) = 0. It also appears in the corresponding Toda-type system with boundary. We now illustrate this case.
Example 3. Let us restrict our general framework to a Z 2 lattice system: we consider the quad-graph with boundary represented on Fig. 8 with the bulk equations given by Q1 δ=0 with the label a on the lines m − n = const and the label b = a on the lines m + n = const. The corresponding Toda-type model reads, for n ≥ 2, We use the third solution for q for Q1 δ=0 in Table 1 to generate the boundary condition on the Toda-type system from the following boundary equation 4 q(x m , q m,1 ,x m+1 ; a) = 0 ⇔x m +x m+1 = 0, on the quad-graph with boundary. It is obvious that the general solution for the boundary values of the field isx m = (−1) mx 0 (for anyx 0 ). Then, following the procedure of folding explained in Section 3.2, we obtain where k(x m , q m,1 ; a) = µ(−1) mx 0 q m,1 + (a − µ)x 2 0 (a − µ)q m,1 + µ(−1) mx 0 .

Lax presentation of boundary conditions for a Toda-type model
In the previous Example 3, we showed that we can recover the boundary conditions introduced in [20]. In this subsection, we show that the correspondence goes beyond this and is in fact also valid at the level of the zero curvature representation. Indeed, our zero curvature representation (11) allows us to recover as a particular case the main equation of [20], which we reproduced as (20) below, and which encodes the symmetry approach to integrable boundary conditions applied to integrable discrete chains. For clarity, we restrict our discussion of this connection to the case treated already in Example 3 above. But we believe that the argument is easily generalizable to most Toda-type models. Figure 8. A quad-graph with boundary with an underlying Z 2 lattice structure. The additional white vertices will play the role of boundary vertices for the Toda-type system with boundary living on the white sublattice. Figure 9. The (white) Z 2 (sub)lattice supporting the Toda-type system with a boundary. Our boundary is inserted between the boundary for the Toda-type system at site n = 0 and the first site of the bulk system at n = 1.
Let us first recall that the boundary condition (18) was obtained in [20] by analysing the matrix equation where A(m, n, λ) is the "discrete time" part of the discrete Lax pair (evaluated at the site n = 0 of the boundary), h is some function acting on the parameter λ and H(m, λ) is a matrix encoding an extra linear symmetry at the site of the boundary (see equations (14) and (15) in [20]) and effectively producing allowed integrable boundary conditions. The main result here is that our zero curvature representation boils down to (20) thanks to a remarkable "fusion" property of the two matrices L which then become A and the fact that our matrix K becomes the matrix H. This goes as follows. The zero curvature representation based on Fig. 8 reads µx .
It remains to perform the changes of parameters as in (19), together with to obtain L(x m+1 , q m,1 , σ(a); c)L(q m,1 , x m , a; c) = 2b (2λ where σ 3 is the usual 2 × 2 Pauli matrix. This also yields This is completely equivalent to the results obtained in [20] up to an irrelevant c → 1/c substitution in H. Note that the map λ → h(λ) = −λ − 1 is also correctly reproduced with our choice σ(c) = −c + 2µ and under the identifications (19) and (21).
6 Connection between 3D boundary consistency and the set-theoretical ref lection equation 6

.1 General approach
In [28], a nice approach was described to obtain a relation between a 3D consistent quadgraph equation and a Yang-Baxter map thus yielding a connection between 3D consistency and the set-theoretical Yang-Baxter equation [17]. It is based on the use of symmetries of the equation Q = 0 and the identification of invariants under these symmetries. Our idea is to extend this connection at the level of reflection maps and integrable boundary equations thus yielding a connection between the 3D boundary consistency introduced here and the settheoretical reflection equation introduced in [14,15]. First, let us recall the method of [28]. Given a quad-graph equation Q(u 00 , u 10 , u 01 , u 11 ; a, b) = 0, let G be a connected one-parameter group of transformations acting on the variables u ij G : (u 00 , u 10 , u 01 , u 11 ) → (û 00 ,û 10 ,û 01 ,û 11 ).
with η being the characteristic of G specified by Methods to obtain the characteristic η can be seen, e.g. in [22,23]. Knowing v, the idea is then to define a lattice invariant I of the transformation group G which satisfies vI = 0, and to use it to define the Yang-Baxter (or edge) variables (X, Y, U, V ) from the vertex variables (u 00 , u 10 , u 01 , u 11 ) by X = I(u 00 , u 10 ), Y = I(u 10 , u 11 ), U = I(u 01 , u 11 ), V = I(u 00 , u 01 ), and assign them to the edges of an elementary quadrilateral as shown in Fig. 10. Once the infinitesimal generator v is known, one can solve for I. An important result of [28] is that the variables X, Y , U , V are related by a map (U, V ) = R(X, Y ) which is a Yang-Baxter map, provided the quad-graph equation determined by Q satisfies the 3D-consistency property. To construct boundary equations that satisfy the 3Dboundary consistency property, we propose to use this method "backwards" in connection with our classification of reflection maps associated to quadrirational Yang-Baxter maps. Choosing the invariant properly, the Yang-Baxter map R can be recognized as one canonical form belonging to the classification of quadrirational Yang-Baxter maps [3,27]. Then we can use the corresponding reflection maps h a and σ to construct q according to the following prescription q(u 10 , u 00 , u 01 ; a) = I(u 00 , u 01 ) − h a (I(u 00 , u 10 )).
The origin of this prescription comes from the folding method explained and used in [14] to construct reflection maps. When translated in terms of vertex variables u 00 , u 10 , u 01 , it gives (23). Indeed, the folding method produces reflection maps B acting on the Yang-Baxter variables and the parameters, of the form (V, b) = B(X, a) = (h a (X), σ(a)). Hence, V = h a (X) and recalling that V = I(u 00 , u 01 ) and X = I(u 00 , u 10 ), this becomes equivalent to q(u 10 , u 00 , u 01 ; a) = 0 with q as defined in (23). In particular, our construction ensures that q and Q satisfy the 3Dboundary consistency property since the corresponding Yang-Baxter and reflection maps satisfy the set-theoretical reflection equation. We see that to carry out this programme, the knowledge of the invariants and hence v is crucial. We use the classification for v obtained in [29] of the so-called five-point symmetries, a subset of which is the set of one-point symmetries which are the ones involved in the above method. Note that although the results of [29] were obtained in the context of Z 2 lattices, we can easily adapt them to the present more general context of quadgraphs. For each v of each quad-graph equation, we provide a solution I 0 for the invariant I which is of the simplest form, the latter meaning that other solutions are obtained from our I 0 in the form f (I 0 ) with f being a differentiable bijection depending possibly on the parameters a or b. To the best of our knowledge, only sparse examples of invariants and corresponding YB maps have been given in the literature so far. In Table 3 below, we give a list of invariants satisfying the above criteria and the corresponding family of YB maps. Then, one only has to use formula (23) to find q. In Tables 1 and 2, the solutions for q (and σ) that we have found using this method (on top of the folding method) are shown with an asterisk.
For this family, the classification in [29] gives three one-point symmetry generators 5 The corresponding simplest invariants I 1 , I 2 and I 3 read So, for the first invariant, the Yang-Baxter variables read which satisfy This yields the relations To recognize the family to which this map belongs, we perform the following transformation on the variables We then obtain the F III quadrirational Yang-Baxter map For this family, we have the following reflections maps where µ is a free parameter. Performing the inverse transformation of (25), we obtain the reflection maps that we can use in formula (23) to obtain q. They read h a (X) = µX a or h a (X) = − µX a , and we obtain q(x, y, z; a) = µ(x + y) − a(y + z) or q(x, y, z; a) = µ(x + y) + a(y + z), both valid with σ(a) = µ 2 a . Note that in this case, the two possibilities for the boundary equations q(u 10 , u 00 , u 01 ; a) = 0 are related by the transformation µ → −µ which leaves σ(a) invariant so that, in fact, we only have one boundary equation here.
Let us perform the same analysis for I 2 . The Yang-Baxter variables are X = u 00 /u 10 , Y = u 10 /u 11 , U = u 01 /u 11 , V = u 00 /u 01 , and they satisfy Performing the following transformation we obtain the H II quadrirational Yang-Baxter map For this family, we have the following reflections maps σ(a) = µ 2 a and h a (X) = a + µ − Xµ a or h a (X) = aX aX + µ − Xµ , or, σ(a) = −a + 2µ and h a (X) = −X or h a (X) = a + (X − 1)µ aX + µ − Xµ .
where µ is a free parameter. Performing the inverse transformation of (26) and using (23) we find q(x, y, z; a) = ax(y + z) + (x + y)zµ or q(x, y, z; a) = y(a(y + z) − (x + y)µ), both valid with σ(a) = µ 2 a , and, q(x, y, z; a) = y(x + z) or q(x, y, z; a) = a y 2 − xz − (x + y)(y − z)µ, both valid for σ(a) = −a + 2µ.  η 1 = (−1) k+l I(s, t) = s + t F III η 2 = u 00 I(s, t) = s/t H II η 3 = (−1) k+l u 2 00 I(s, t) = 1/s + 1/t F III A1 δ=1 η 1 = (−1) k+l I(s, t) = s + t F II A2 η 1 = (−1) k+l u 00 I(s, t) = st F I Finally, using I 3 yields the relations (24) so that we obtain the F III quadrirational Yang-Baxter map again. Therefore, we can use the same reflection maps but because the invariant is different we may obtain different expressions for q. A direct calculation gives q(x, y, z; a) = ax(y + z) + (x + y)zµ or q(x, y, z; a) = ax(y + z) − (x + y)zµ, both valid with σ(a) = µ 2 a . These are in fact the same solution under the transformation µ → −µ and it coincides with the first solution already obtained from the H II family. Let us remark that in some cases, point transformations on the parameters are needed on top of Möbius transformations to recognize the canonical Yang-Baxter map. This may affect the form of the map σ to be used in the 3D boundary consistency condition. This happens for Q3 δ=0 , H3 and A2. Finally, let us also mention that this method has some inherent limitations due to the fact that for families Q2, Q3 δ=1 and Q4, there are no one-point symmetry generators in the classification of [29]. Our other method described in Section 3 does also work for these families, as can be seen in the tables. Whether or not the corresponding solutions for q for these "missing" families can be mapped back to reflection maps is an interesting open question.

Conclusions and outlooks
We proposed a definition for integrable equations on a quad-graph with boundary and introduced the notion of 3D boundary consistency as a complement of the 3D consistency condition that is used as an integrability criterion for bulk quad-graph systems. Just like quadrilaterals are fundamental structures when one discretizes an arbitrary surface without boundary, we argued that triangles appear naturally when one considers surfaces with boundary. Therefore, it is natural to associate a 3-point boundary equation to describe the boundary locally just like one associates a 4-point bulk equation to describe the bulk locally. We presented two different methods to find solutions of the 3D boundary consistency condition given an integrable quadgraph equation of the ABS classification. The terminology "integrable boundary" is also backed up by the discussion of other traditional integrable structures like Bäcklund transformations and zero curvature representation for systems with boundary. This is also supported by the connection that we unraveled between 3D boundary consistency and set-theoretical reflection equation. As a by-product of our study, we were also able to introduce three-leg forms of boundary equations and hence to introduce Toda-type systems with boundary.
The present work lays some foundations for what we hope could be a new exciting area of research in discrete integrable systems. Among the many open questions one can think of, we would like to mention a few that we believe are important: finding a method of classification of boundary equations given a bulk quad-graph equation, tackling the problem of posing the initial-boundary value problem for quad-graph equations with a boundary, understanding the connection of our approach with the singular-boundary reduction approach of [6], implementing the discrete inverse scattering method with a boundary with the hope of finding soliton solutions, etc.