Two-dimensional translation-invariant probability distributions: approximations, characterizations and no-go theorems

We study the properties of the set of marginal distributions of infinite translation-invariant systems in the two-dimensional square lattice. In cases where the local variables can only take a small number d of possible values, we completely solve the marginal or membership problem for nearest-neighbours distributions (d = 2, 3) and nearest and next-to-nearest neighbours distributions (d = 2). Remarkably, all these sets form convex polytopes in probability space. This allows us to devise an algorithm to compute the minimum energy per site of any TI Hamiltonian in these scenarios exactly. We also devise a simple algorithm to approximate the minimum energy per site up to arbitrary accuracy for the cases not covered above. For variables of a higher (but finite) dimensionality, we prove two no-go results. To begin, the exact computation of the energy per site of arbitrary TI Hamiltonians with only nearest-neighbour interactions is an undecidable problem. In addition, in scenarios with d≥2947, the boundary of the set of nearest-neighbour marginal distributions contains both flat and smoothly curved surfaces and the set itself is not semi-algebraic. This implies, in particular, that it cannot be characterized via semidefinite programming, even if we allow the input of the programme to include polynomials of nearest-neighbour probabilities.


Introduction
The distribution of stars at large scales, the vacuum state of a quantum field theory, the thermal state of any solvable spin model: these are examples of systems with infinitely many sites where the description of a bounded environment does not depend on its location within the lattice. We call this property translation invariance.
An agent exploring an infinite translation-invariant (TI) word would find that the statistics of the local random variables which he can access are constrained by the requirement of infinite TI. However, and despite a long history of research on TI systems, driven by the needs of statistical physics (e.g. [1]), it is far from clear what those constraints exactly are.
This conundrum is at the essence of the TI MARGINAL problem, where a number of probability distributions of finitely-many variables are provided and the task is to certify if they correspond to the marginals of a TI system. MARGINAL arises naturally at the intersection of quantum information science and condensed matter physics, when we try to determine whether the dynamical structure factors of a large spin system are compatible with an underlying Bell local quantum state [2].
ENERGY, the dual problem of computing the minimum energy per site of a translationinvariant Hamiltonian, is also of special concern for the mathematical physics community. While ENERGY is efficiently solvable in one-dimensional (1D) systems [3][4][5], in two dimensions very few solvable instances are documented.
What do we know about MARGINAL or ENERGY in two dimensions? Not much. We know that the set of TI marginals is closed and convex [5,6]. We also know that the problem of determining the existence of near-neighbour marginals of a TI distribution with a given (finite) support is undecidable [4,7,8]. In [9], the case where the random variables are dichotomic is studied and necessary and sufficient conditions are provided for the existence of a TI extension, given the infinitely many marginal distributions of any n lattice sites. In [5], explicit examples of probability distributions for lattice sites within a 2 × 2 square are given which, despite satisfying all local symmetries associated with translation invariance, do not admit two-dimensional (2D) TI extensions. Similar instances of 2D non-extendible distributions constructed from 3 × 3 squares satisfying additional rotation and reflection symmetries are also provided. In [6], the approximability of the MARGINAL problem via convex polytopes is considered.
It is still open under which scenarios one can solve MARGINAL and ENERGY exactly, or how to attack these problems numerically. 1 Even the geometry of the set of all marginal distributions of 2D TI systems remains a mystery: in one dimension, the marginal distributions of the variables of finitely many sites of TI systems form a convex polytope in the space of probabilities. 2 Should we expect this to hold in 2D systems as well?
Motivated by the desire of understanding the nature of quantum non-locality in 2D materials [2], in this paper we will considerably advance the above fundamental questions. We will show that, in scenarios where the random variables take a small number d of possible values, the MARGINAL problem for nearest-neighbour (for d = 2, 3) and next-to-nearest neighbour distributions (for d = 2) is exactly solvable. We also prove solvable a natural variant of MARGINAL where the TI extension is also required to satisfy invariance under reflection of the horizontal and vertical axes. An immediate consequence of these results is that, in any of the above cases, it is possible to solve ENERGY exactly. For the scenarios not covered above, we provide an algorithm to approximate the set of TI marginals or solve the ENERGY problem up to arbitrary precision.
On the negative side, we show that, for arbitrary TI Hamiltonians, ENERGY is undecidable. We also prove that, contrary to the solvable cases, for local random variables of high cardinality, the set of nearest-neighbour distributions admitting a TI extension is no longer a convex polytope in the space of probabilities or even a semi-algebraic set. 3 This implies, in particular, that standard tools from convex optimization, such as linear programming [13] or semidefinite programming [14] cannot be used to fully characterize Bell non-locality in large 2D condensed matter systems.
The structure of this paper is as follows: in §2, we will define and list known properties of the set of marginal 2D TI distributions. We will introduce the problems MARGINAL and ENERGY, and show how to solve their 1D versions. We will also describe a simple symmetrization process that will play an important role in the mathematical proofs to come. The remaining sections present our original contributions. In §3, we will present a practical algorithm to characterize the set of 2D marginals up to arbitrary accuracy. Later, in §4, we will provide a few instances of the marginal problem which are exactly solvable. In §5, we will prove two no-go theorems about the set of 2D TI marginals, namely: (a) the problem of computing ENERGY exactly is undecidable; (b) for random variables with support d greater than or equal to 2947, the set of nearest-neighbour distributions is not a semi-algebraic set.

The set of marginal 2D TI distributions: definition and known properties
Consider an infinite 2D square lattice, and suppose that each site (x, y) ∈ Z 2 has access to a local random variable a x,y that can take d < ∞ possible values. We will call d the local dimension of the lattice. For any finite subset K of Z 2 , we will assume that the variables a K ≡ {a z : z ∈ K} follow a probability distribution P K (a K ). We say that the system is translation-invariant (TI) if, for any finite set K and any vector c ∈ Z 2 , TI is a very common property in nature. It is satisfied approximately in large crystal structures, and exactly by any stationary state in quantum field theory. In addition, most exactly solvable models in statistical physics comply with this symmetry in the thermodynamic limit of infinitely many sites.
In this paper, we consider a scenario where an agent conducts local observations within an infinite 2D TI system. Intuitively, the condition of TI means that the statistics of the variables within the immediate vicinity of the agent do not give any clue about his position in the lattice. Our goal is to determine how the statistics of such a sample are constrained by the requirement that these variables arise from a 2D TI system (figure 1).
With a slight abuse of notation, any function P assigning probability distributions P K (a K ) to any finite set K ⊂ Z 2 satisfying equation (2.1) and the consistency conditions a L P K∪L (a K , a L ) = P K (a K ) for K ∩ L = ∅ will be called a 2D TI probability distribution. Each of the distributions P K (a K ) will be denoted the marginal of a 2D TI distribution, or simply a TI marginal. A few special finite subsets of Z 2 will appear frequently throughout this article, so we will need special names for them. v, h, +, −, will denote, respectively, the sets )) hence describe the statistics between nearest neighbours (next-to-nearest neighbours). We will also use m × n to describe the rectangle {(x, y) : 0 ≤ x < m, 0 ≤ y < n}; and |m × n|, for {(x, y) : |x| < m, |y| < n}. Finally, 1, . . . , n will denote the set ∪ n k=1 {(k, 0)}. We will now introduce the two problems which will concern us through the rest of the article.  The 2D TI marginal problem. In the example, there is only one region, K = 3 × 3. A probability distribution Q 3×3 for the nine variables {a i,j : i, j = 0, 1, 2} is given (the green square), and the question is whether there exists a 2D TI system with a marginal probability distribution P 3×3 (the yellow squares) equal to Q 3×3 . (Online version in colour.)

Definition 2.1 (The TI MARGINAL problem). Given a finite number of finite subsets of
In the affirmative case, the distribution P will be called a TI extension of {Q i (a K i )} i .
In the following, each configuration of cardinality d < ∞ and finite subsets of Z 2 {K i } n i=1 will be called a scenario. In any scenario, the set of TI marginals is convex and closed [5] (figure 2).
Sometimes, e.g. in statistical physics, we are not interested in solving the marginal problem as much as in optimizing linear functionals over the set of marginals of 2D TI distributions. This motivates the following problem:

Definition 2.2 (The ENERGY problem). Given
where the minimization is carried out over all TI distributions P.
The ENERGY problem is the dual of the MARGINAL problem: given an oracle to solve one of them approximately, one can devise an algorithm that invokes the oracle a polynomial number of times in order to solve the other problem with a similar accuracy [15].
, consider the family of (classical) 2D Hamiltonians of the form The reader can find a proof at the end of this section. As we will see, the TI marginal problem is very hard. However, there exists an interesting variant whose solution turns out to be trivial. The solution of this problem is folklore among the community of condensed matter physicists: a distribution Q 1,...,s (a 1 , . . . , a s ) admits a TI extension iff Q 1,...,s−1 = Q 2,...,s [3][4][5].
Given an arbitrary distribution Q K (not necessarily a TI marginal) of local variables over a large rectangle K = m × n, one can derive a 2D TI distribution P such that P s×t equals a spatial average of all the s × t possible rectangles inside K with a correction of the order O(max(s, t)/min(m, n)).
Definition 2.4 (The symmetrization procedure). Given Q K , define a distributionP over the whole plane by tiling it with copies of the distribution Q K . That is, for K(i, j) ≡ K + (mi, n, j), P ∪ ij K(i,j) (a ∪ ij K(i,j) ) ≡ i,j Q K (a K(i,j) ). Then, the TI distribution P given by will be called the symmetrization of Q.
That P is TI invariant can be seen by noting thatP is invariant under translations of m (n) sites in the horizontal (vertical) axis. Random translations of {1, . . . , m} and {1, . . . , n} sites in those axes thus turnP into a 2D TI distribution.
It can be easily checked that In general, Taking the limit N → ∞, we conclude that the energy value of any configuration can be matched by a TI marginal. Conversely, given any 2D TI

Approximations of the set of TI marginals
In this section, we will prove that, for any given scenario, the set of TI marginals admits an approximate characterization up to arbitrary accuracy. The symmetrization protocol suggests a simple (but expensive) converging sequence of relaxations. Given Q K (a K ), with K = s × t, a necessary condition for Q K (a K ) to be a TI marginal is that Q K (a K ) is the marginal of a distribution P (n) over the square n × n, subject to the rules: Intuitively, P (n) is modelling the marginal for the region n × n of an overall 2D TI distribution containing Q K . The verification can be carried out via linear programming [13]. Linear programming is a branch of convex optimization concerned with the resolution of problems of the form maxc ·x, Ax ≥b,x ≥ 0. (3.2) wherec ∈ R p , the q × p matrix A andb ∈ R q are the inputs of the problem;x ∈ R p are the problem variables; ands ≥ 0 is used to denote that all the components of the vectors are non-negative. For each (primal) linear programme, there exists a dual problem minb ·x, Remarkably, the solutions of both primal and dual problems coincide. Numerical algorithms aimed at solving one problem hence run optimizations over the primal and dual problems. This allows the solver to give rigorous upper and lower bounds on the optimal solution. At present, there exist numerous free software implementations of interior-point methods for linear programmes [13,16]. In addition, all linear programmes can be solved exactly via the costly Fourier-Motzkin elimination method [17]. In our case, we need the solver to verify that there exists a probability distribution P (n) n×n for n 2 variables satisfying the linear conditions (3.1), together with P (n) K = Q K . This can be formulated as a linear program by regarding each probability P (n) n×n (a K ) as a free variable, and choosing A, b so that P (n) n×n (a K ) satisfies the corresponding linear constraints. As for the objective function to optimize, we can takec = 0, i.e. the solution of the primal linear programme will be zero provided that there exists a feasible point. If there exists a distribution P (n) n×n (a K ) compatible with Q K , the solver will find it. Conversely, if there is no such distribution, then the solver will return a solution for the dual problem with an objective value smaller than 0. Such a solution is, in effect, a certificate of infeasibility, that is, a computer-generated proof that Q K is not a TI marginal.
The method described above is a relaxation of the property of being a TI marginal: if Q K does not pass the nth test, then it clearly is not a TI marginal. If, on the other hand, Q K passes the nth test, then we can apply the symmetrization protocol over the distribution P n×n and obtain a 2D TI distributionP that, due to equation (2.6) and the condition P K = Q K , satisfiesP This algorithm is highly inefficient, though, since the time and space complexity of the computations scales as O(e αn 2 ). Actually, there is a much more practical relaxation of the set of TI marginals achieving the same accuracy that just involves O(e βn ) operations.

Definition 3.1 (Approximate solution to the marginal problem). Given
This can again be formulated as a linear programme, and is obviously a relaxation of the property of being a TI marginal. The condition (3.4) will appear often in the rest of the article, so we will give it a name. Any distribution P n×t in the rectangle n × t satisfying the above conditions will be called locally translation invariant (LTI).
To see that the relaxation above is O(s/n)-close to the actual set of TI marginals, suppose that Q K is the marginal of P (n) n×t . We will next extend P (n) n×t to a distribution P n×n with the property that the marginal of any s × t rectangle equals Q K .
In order to derive P n×n , we regard P (n) n×t as the t-site marginal of a TI 1D system with local variables of dimension d n . We can do so because P (n) n×t satisfies the second line of (3.4). From the triviality of the marginal problem for TI 1D systems, we know that there must exist a distribution P n×n satisfying the afore-mentioned properties. Finally, it is easy to see, from equation (2.5), that the symmetrization of P n×n will be a 2D TI distributionP with marginalP To conclude, we would like to remark that the above sequence of relaxations of the set of marginals also allows the user to approximately solve ENERGY. Indeed, given {(F i , K i )} i , one can carry out, via linear programming, the optimization: From what we reasoned above, it follows that

The exact marginal problem in two dimensions: characterizations
In this section, we will identify certain variants or scenarios where the TI marginal problem is exactly solvable. That is, where there exists an algorithm that will determine with certainty if the given distributions are TI marginals or not.

(a) The exact marginal problem with reflection symmetry
Let us start with a relevant variant of the marginal problem: the 2D TI marginal problem with reflection symmetry. This is the case where, in addition to demanding the existence of a TI extension, we require this to be invariant under reflections on both axes. Since Nature is approximately invariant under parity reflection, this condition holds for the thermal states of many physically relevant Hamiltonians, such as the isotropic Ising model and the Potts model [18,19]. and P n×n (a n×n ) = P n×n (a H n×n ) = P n×n (a V n×n ), for all n. Here, a H n×n , a V n×n denote, respectively, the permutation of the variables (a {(x,y)} )| 0≤x≤s−1,0≤y≤1 associated with a reflection over the horizontal and vertical axis, respectively.
Interestingly, the s × 2 TI marginal problem with reflection symmetry can be solved exactly.

Proposition 4.2. The existence of a 2D TI distribution with reflection symmetry is equivalent to the conditions:
Proof. The second condition implies that, viewed as a 1D system of s sites with local variables of dimension d 2 , Q s×2 is the marginal of a TI system. In particular, for any n one can find P 1 n×2 with the property that P 1 This property is kept if we make the distribution invariant under reflection on both axes.
Reflection under the horizontal axis implies, in particular, that P 2 n×1 = P 2 n×1+(0,1) . Therefore, we can view P 2 as the 2-site marginal of a 1D TI system with variables of local dimension d n and extend it to an n × n square. From that point on, we can consider the symmetrization of this distribution P 3 n×n , which will return a 2D TI reflection-symmetric distribution with Q K + O(1/n) as a marginal. Invoking the closure of the set of marginals of 2D TI distributions [5], we conclude that Q K admits a 2D TI extension P. Finally, since Q K satisfies reflection symmetry, one can choose P to be symmetric as well.
The characterization (4.1) of the set of s × 2 TI marginals with reflection symmetry implies that, for is exactly solvable. Indeed, letP K be any TI marginal and call E the value of the functional in equation (2.2) evaluated inP K . Then it is easy to see that P 2 s×2 , as defined by equation (4.2) (replacing n by s), admits a TI, reflectionsymmetric extension. Moreover, the value of the functional is also E. It follows that, for , one can assume that the minimizer of ENERGY(F K ) is a TI, reflection-symmetric marginal. Hence one can use linear programming to solve the problem.

(b) Binary local random variables
Take the local dimension of the random variable at each site to be d = 2 (bits), and suppose that we just want to characterize the bivariate distributions P h (a, b), P v (a, b) between nearest-neighbours. The next proposition states that the problem is solvable even for lattices of spatial dimensions higher than 2.
For d = 2 and k = 2, 3, the nearest-neighbour marginal problems are thus trivial. This perhaps explains why the only known solvable classical models in two dimensions are bit models.
Proof. For simplicity, consider the 2D case k = 2. Then, for d = 2, LTI on either P h or P v (i.e. a)) implies that both distributions P h and P v are symmetric. Now, suppose that the pair (P h , P v ) admits a 2D TI extension. If we make this extension reflection-invariant-via equation (4.2), the distributions P h , P v will not change. In other words, (P h , P v ) is a TI marginal iff it is a TI reflection-symmetric marginal. From proposition 4.2, we know that the marginals of 2D TI reflection-symmetric distributions correspond to the marginals of symmetric, LTI squares. This allows us to completely characterize the set of TI marginals (P h , P v ) via linear programming. The above symmetrization argument holds not only in two dimensions, but in any spatial dimension, and so does the proof of proposition 4.2.
It follows that the set of nearest-neighbour marginals in any dimension is described by a convex polytope, i.e. a set defined by a finite number of linear inequalities or facets . . , are linear functionals on the probabilities P h (a, b), P v (a, b), etc. Using standard combinatorial software [20], we managed to derive the facets which define the 2D set. We verified, using linear programming, that the set of We now move to the problem of characterizing the distributions (P h , P v , P + , P − ) corresponding to horizontal and vertical nearest-neighbours and north-east (+) and south-east (−) next-tonearest neighbour distributions in two dimensions. This problem is not trivial (i.e. it does not reduce to verifying LTI), since the distribution a) for s, t = h, v, +, − and yet it does not admit a TI extension. 4 And yet, as the next result shows, the set of nearest and next-to-nearest neighbour marginals can also be characterized via linear programming.

Proposition 4.4.
For d = 2, the nearest and next-to-nearest neighbours marginals P h , P v , P + , P − admit a TI extension iff they are marginals of a distribution P 2×2 satisfying LTI.
Proof. If P h , P v , P + , P − are TI marginals, then they must constitute an approximate solution of the marginal problem, in the sense explained in definition 3.1. In particular, they must belong to the polytope of distributions P h , P v , P + , P − admitting a 2 × 2 LTI extension. We find, using the combinatorial software Panda [20], that this polytope has 13 extreme points, each of which belongs to six classes modulo rotations, reflections and relabellings of the variables: (C1) The two points in this class are deterministic: 1 . They obviously admit a deterministic 2D TI extension. (C2) This class has two elements of the form and and It is generated by the 2 × 2 deterministic distribution 1 0 0 1 . (4.10) (C4) The class has two elements of the form with s = 0, 1. They are, respectively, generated by symmetrizing the 2 × 2 deterministic distributions: 1 0 0 0 , 0 1 1 1 . (4.12) (C5) This class has four elements, generated via applying reflections to the generator: and 14) The latter, in turn, can be generated by the 3 × 3 deterministic distribution Since all these points are TI marginals, any convex combination thereof will also be a TI marginal. It follows that, for d = 2, the first level of the hierarchy presented in definition 3.1 already characterizes all nearest and next-to-nearest neighbour TI marginals.
From proposition 4.4, one can infer the exact solvability of the ENERGY problem for TI Hamiltonians with nearest and next-to-nearest neighbour interactions. Indeed, since any TI marginal is a convex combination of the above 13 points, one just needs to evaluate the energy per site of each of them and take the smallest result. It is worth noting that, for each extreme point P listed above, there exists a unique (modulo translations) deterministic distributionP for the whole lattice such that supp(P g+z ) ⊂ suppP g , with g = h, v, +, − and z ∈ Z 2 . Now, let the value of ENERGY(F h , F v , F + , F − ) be achieved by just one extreme point P of the set of TI marginals. If P is generated by tiling the plane with an irreducible rectangle of size m × n, then the number of lattice configurations achieving the minimum energy per site of the corresponding Hamiltonian is mn. Assuming that, at zero temperature, the system is in an equal mixture of all such configurations, its entropy will be S 0 = ln(mn). As we will see in the next section, the uniqueness of the extended distribution on the whole lattice is lost when each site can take more than two values. This is in contrast to the situation in one dimension, where maximal entropy extensions of LTI distributions exist and are unique [5].

(c) Ternary local random variables
The purpose of this section is to characterize the set of TI nearest-neighbour marginals P h , P v for ternary local random variables (d = 3). For d = 2, all such P h and P v can be characterized simply by imposing (4.3). For d = 3, this condition is not sufficient. Take, for example, the distribution P h (0, 2) = P h (1, 0) = P h (2, 1) = P v (0, 2) = P v (1, 2) = P v (2, 1) = 1/3. It can be easily checked that this distribution satisfies equation (4.3). However, one can verify numerically that it does not belong to the first level of the hierarchy of approximations proposed in definition 3.1.
Luckily, this first approximation turns out to be enough to characterize the nearest-neighbour TI marginals.

Proposition 4.5. For d = 3, the nearest-neighbour marginals P h , P v admit a TI extension iff they are compatible with a distribution P 2×2 satisfying LTI.
Proof. The proof follows the same lines as the proof of Proposition 4.4. To find the extreme points/facets of the polytope of marginals admitting a P 2×2 LTI extension, we first use a linear programme to maximize random objective functions under such a set. The vectors thus obtained are either vertices of the polytope or lie on the facets. By computing the dual description of this set of vectors, we obtain a set of linear inequalities some of which may be facets of the polytope we seek. To check whether the inequalities are indeed facets, we use linear programming again but this time taking the vectors which define the inequalities (without the bounds) as objective functions. If an inequality is indeed a facet, then maximizing this objective function using the same constrains as before will give us the upper bound of the inequality as the value of the objective function. If it is not a facet, then this maximization will violate the upper bound of the inequality, at which time we can add the vector achieving this maximization to the list of potential vertices. By iterating this violation/addition procedure until no inequality can be violated when performing the second maximization, we obtain all the facets of the polytope, from which the extreme points can also be computed. In the d = 3 case, the polytope is defined by 98 extreme points, which fall into 10 equivalence classes after taking local permutations of outcomes and reflections into account. A representative from each class is given below, with P h (a, b) and P v (a, b) written as a vector whose ith coordinate in base 3 gives the values of a, b: P h,v (a, b) ≡ (0 a 0 b , 0 a 1 b , 0 a 2 b , . . . , 2 a 2 b ). It can be readily checked that each of these class representatives can be generated by symmetrizing the depicted deterministic distributions accompanying it.

The exact marginal problem in two dimensions: no-go theorems
In the view of the results in the previous section, one would imagine that the exact resolution of the marginal problem for scenarios with high local dimension d is merely a matter of computational power. Note that all the marginal sets characterized so far satisfy the following properties: (i) They allow us to solve ENERGY exactly.
(ii) When we represent them in probability space, they happen to be convex polytopes, i.e. sets determined by a finite number of linear inequalities. More broadly, they are semialgebraic sets. A subset S of R s is a basic closed semi-algebraic set if there exist a finite number of polynomials {F i } u i=1 onx ∈ R s and the slack vectorȳ ∈ R t such that All subsets of R t which one can characterize via linear programming [13] or the more general tool of semidefinite programming [14] fall within this category. A semi-algebraic set would be the union of a finite number of basic closed semi-algebraic sets. However, it is easy to show that any closed convex semi-algebraic set is also basic, so in the following we will use both terms interchangeably.
Extrapolating, one would expect that the marginal sets of 2D TI distributions with arbitrary d should retain at least one of the above features.
In the following pages, we show that this is not the case even when we aim at solving the simplest non-trivial marginal problem: the characterization of (P h , P v ).

(a) Undecidability of ENERGY
In this section, we will prove the following result:

Theorem 5.1. There exists no algorithm to solve ENERGY(F h , F v ) for arbitrary local dimension d and
The proof of this theorem, as well as the proof of the next no-go result, will rely heavily on certain mathematical results on tilings of the plane, so let us first introduce a basic vocabulary.
Let A be a finite alphabet. Any function f : Z 2 → A defines a tiling of the plane with the set of tiles A. Any pair of subsets T = (T h , T v ) of A × A defines a tiling rule. We say that a tiling f respects the tiling rule T if, for all (x, y) ∈ Z 2 , If the rule T is implicitly known, we call f a valid tiling. Certain rules T do not admit any valid tiling, e.g.: take A = {0, 1} and T h = {(0, 1)}, T v = A × A. The problem of deciding if a given rule T admits a valid tiling was first raised by Wang [21], and proven undecidable by Berger [22], who also showed the existence of tiling rules T admitting just aperiodic tilings. The construction used by Berger to derive the latter result uses 20 426 tiles. This number has been decreasing over the years as new aperiodic tiling rules requiring less tiles were found. The current record, held by Jeandel and Rao, uses only 11 tiles [23].
We will prove that ENERGY(F h , F v ) is undecidable by reducing it to the general tiling problem. Let T be a tiling rule, and define the input of ENERGY as We claim that ENERGY(F h , F v ) = 2 iff T admits a valid tiling. Consequently, the ENERGY problem is undecidable. First, suppose that T admits a valid tiling. We choose an n × n square of this tiling, and apply symmetrization to obtain the marginals P (1/n). That way, we obtain a sequence of TI marginals (P Conversely, suppose that ENERGY(F h , F v ) = 2. Then, there exists a 2D TI distribution P saturating the bound. Given P |n×n| , take any configuration a |n×n| such that P |n×n| = 0. Obviously, a |n×n| is a valid tiling for |n × n|. Now, consider P |(n+1)×(n+1)| and take any configuration a |(n+1)×(n+1)| such that a |n×n| is an inner square of a |(n+1)×(n+1)| and P |(n+1)×(n+1)| (a |(n+1)×(n+1)| ) = 0. That such a configuration must exist follows from the fact that a P |(n+1)×(n+1)| (a |n×n| , a ) = P |n×n| (a |n×n| ) = 0, where a denotes the local variables of the inner border of the |(n + 1) × (n + 1)| square. Iterating, we obtain a sequence (a |n×n| ) n of valid tiles for ever-growing squares with the particularity that, for all n, a |n×n| is contained in a |(n+1)×(n+1)| . The desired tiling of the plane is hence given by the function f assigning to each point (x, y), with |x| < n, |y| < n the corresponding tile in a |n×n| .

(b) The set of TI nearest-neighbours marginals is not a semi-algebraic set
The goal of this section is to prove the following theorem. Before proving this theorem, we want to remark that the value of d = 2947 comes from the number of tiles required by our construction below. It is not necessarily tight: for all we know, the critical value d which marks the transition from polytopes to curved sets could be as low as 4. It is possible that, similar to the history of the minimum set of Wang tiles discussed above, our upper bound d = 2947 be improved in the future. A consequence of theorem 5.2 is that, for d ≥ 2947, the set of nearest-neighbour TI marginals cannot be characterized by either linear or semidefinite programming [14].
Proof. The idea of the proof is as follows: given a finite alphabet A, we will define a tiling rule T , and we will consider the set P of all TI marginals (P h (a, b), P v (a, b) That is, we will consider the intersection between the set of TI marginals and two planes defined by integer coefficients. If the set of all TI marginals (P h , P v ) were semi-algebraic or a polytope, then so would P.
We will then define two linear 1-site parameters ω ≡ a∈Aω (a)P (0,0) (a) and η ≡ a∈Aη (a)P (0,0) (a), (5.5) and characterize the set S of feasible values (ω, λ) in P. We will find that the boundary of S contains both flat and a smoothly curved pieces, so we will conclude that the set of nearestneighbour TI marginals with variables of dimension d = |A| does not form a polytope. Similarly, by proving that S is not a basic closed semi-algebraic set, we will demonstrate that neither is the set of TI marginals. The proof relies heavily on the connection between aperiodic tilings and immortal points of dynamical systems first pointed out in [24] and later extended in [25]. Before giving our implementation which proves the statement above we will briefly review the general construction given in [24,25].
Let M be a 2 × 2 matrix with rational coefficients; and c, a rational vector in R 2 . We will consider a number u of unit squares {R i } u i=1 in R 2 , each of them defined by their integer corners U i = {(m i , n i ), (m i + 1, n i ), (m i , n i + 1), (m i + 1, n i + 1)}. Let R ≡ ∪ u i=1 R i . The tiles used in this construction are all derived from the prototile given in figure 3. A tile is labelled by four vectors in R 2 , called  t, b, l, r (top, bottom, left and right) and a region i ∈ {1, . . . , u}. The components of l and r are  The first rule can be seen to imply that, for any valid row of tiles k = 1, . . . , m of the form (i, t k , b k , l k , r k ), where t , b , are, respectively, the arithmetic means of the top and bottom vectors of the n tiles. Note that, by convexity, t ∈ R i , b ∈ conv(R) ( figure 4).
Let A be such that the set L of possible left vectors l, viewed as a subset of R 2 , is bounded and equal to the set of feasible right vectors. It follows that, for n sufficiently large, we can view each valid row as a vectorz ∈ R i undergoing the transformationz → Mz + c. Moreover, by the second tiling rule, given the validly tiled rectangle |m × n|, with m n 1, the sequence of top vector averages on each row j can be interpreted as an orbit inside R given by (z j = f −j (z 0 ) + O(n/m) : j = −n, . . . , n). In addition, the region of the jth row indicates which square R 1 , . . . , R u containsz j . Taking the limit n → ∞, n/m → 0, a valid tiling is only possible if f j (z 0 ) ∈ R for all j ∈ Z.z 0 is then called an immortal point of the dynamical system given by M, c and R.
The breakthrough in [24,25] was to realize that, conversely, there always exists a finite alphabet A such that any immortal point of the system (M, c, R) can be represented with a valid tiling of A. Letv ∈ R 2 , and, for any k ∈ Z, define A k (v) ≡ kv , where the floor is taken for each coordinate ofv. Now, denote B k (v) ≡ A k (v) − A k−1 (v). It can be shown that, ifv ∈ R j , then B k (v) ∈ U j . Similarly, N). Now, suppose that, for any i = 1, . . . , u, and any immortal pointv ∈ R i , A contains the set A of all tiles T k (v) with region i and top, bottom, left and right vectors of the form: for all k ∈ Z. It can be verified that these tiles satisfy the conditions f (t) T k+2 (v), . . . , T k+m (v) form a valid row for any m. We can extend this row below by placing another row of the form T k (f (v)), T k+1 (f (v)), T k+2 (f (v)), . . ., with Mv ∈ R i . Iterating, we can tile the whole plane in this fashion, provided thatv is an immortal point. From the properties of B k (v), the averages of the inputs of each row j tend to f j (v) in the limit n → ∞, n/m → 0.
It rests to see that A is a finite set. The top and bottom vectors of each tile belong to ∪ u i=1 U i , and so they can only take finitely many values. As for the left or right vectors of the tile, they must belong to the set L = {r k (v) :v ∈ R 2 , k ∈ Z}. Let l ∈ L and denote by l s its sth component. Using the where j + (j − ) ranges over all those j such that M s,j > 0 (M s,j < 0). Call m s the common denominator m s of the rational numbers {M s,j } j , c s (remember that M, c are assumed to be rational). From the definition of L , it follows that l ∈ L must satisfy ( m 1 0 0 m 2 )l ∈ Z 2 . Define then the set L = {l : m s l s ∈ Z, μ − s ≤ l s ≤ μ + s , s = 1, 2}. Clearly L is bounded and contains L . In order to generate A, we go through all regions i, j = 1, . . . , u and all combinations of top and bottom vectors t ∈ U i , b ∈ U j , and any left vector l ∈ L and verify that the right vector r = f (t) + l − b also belongs to L. If it does, then we add the tile (i, t, l, r) to the definition of A. That way, we end up with a finite alphabet that, via the above tiling rules, can describe all immortal orbits of the dynamical system (M, c, R). This is a simplification of the construction proposed in [25] to simulate a dynamical system where the transformation f may depend on the region i, i.e. z → M i z + c i .
To prove our result, we take M, c to be M = As we will see later, φ 0 is an irrational multiple of 2π . Hence, by applying M sequentially, we can induce a rotation arbitrarily close to any angle θ. It follows that a pointz is immortal in the dynamical system given by (M, c, R) iff z −ĉ 2 ≤ 2/5. We are interested in the following extensive quantities: Intuitively, the vector z defined at the top side of each row of the tiling is turning by an amount φ 0 from row to row. The witness (5.10) corresponds to the average ofz −ĉ's first coordinate in the region R 2 . That is, (ϕ + kφ 0 ), (5.13) where ϕ is the angle ofz at row j = 0 with respect to thex-axis; μ = z −ĉ 2 and χ O (p) is the characteristic function that equals 1 when p ∈ O and 0 otherwise. Since φ 0 is not congruent, we expect the above expression to converge to To justify equations (5.14) and (5.15), though, we need to prove a result that relates the integration of a function with the sampling over (ϕ + kφ 0 ) k . Lemma 5.3. Let φ 0 be such that mφ 0 = 0 (mod 2π ) for all m, and let g(φ) be any piecewise continuous bounded function in φ ∈ [−π , π ]. Then, for any ϕ ∈ [−π , π ], Since the process is, in fact, an average, |g R − g | < , and we conclude that equation (5.16) holds for all piece-wise derivable bounded functions. The general result can be obtained by noticing that, for any piece-wise continuous function g, there exist two sequences of piece-wise derivable functions (g u n ) n , (g d n ) n such that g d n (θ ) ≤ g(θ ) ≤ g u n (θ) for all n, θ and lim n→∞ π −π dθg d n (θ) = lim n→∞ π −π dθ g u n (θ ).
We have just shown that, in the limits n → ∞, n/m → 0, the vector of feasible values (Ω, H) in an n × m valid tiling is parametrized by equation (5.12). What does this have to do with TI marginals? Consider a random variable taking values in the tile set A, and a TI marginal (P h , P v ) satisfying the constraint: Consider also the linear functionals given by where t 1 (a) ∈ {0, ±1} (i(a) ∈ {1, 2}) denotes the second coordinate of the top side (the region) of tile a. By definition, for any n, m there exists an LTI distribution P n×m with nearest-neighbour marginals (P h , P v ). P n×m can be seen as a convex combination of tilings a k n×m with weight p k of the set n × m. Moreover, due to condition (5.19), all of them must be valid tilings. Hence, In the limit n → ∞, n/m → 0, the point (ω, η) belongs to the convex hull S of the curve (5.12). It can be verified that (d 2 η/dω 2 ) ≤ 0, i.e. the curve is concave. Hence the boundary of S is given by curve (5.12) and the segment joining its start and end points ( figure 6). It rests to show that any point of the shaded region in figure 6 is achievable by a TI marginal. Since the set of TI marginals is convex, it is enough to see that one can achieve the extreme points of the set S. Let 1 5 ≤ μ ≤ 2 5 and take any valid tiling a Z 2 of the plane describing a vector z ∈ R 2 with z −ĉ 2 = μ. Further take any increasing sequence of rectangles m (k) × n (k) , with lim k→∞ n (k) /m (k) = 0. Symmetrizing the deterministic distribution a m (k) ×n (k) , we obtain a sequence of TI marginals (P (k) h , P (k) h ) k which violate (5.19) by an amount O(1/n (k) ). By Weiestrass' theorem, there exists a subsequence of this sequence that converges to a pair of distributions, call them (P h , P v ). By closure of the set of TI marginals, (P h , P v ) admits a 2D TI extension. Moreover, it satisfies (5.19) and ω, η are given by equation (5.12).
As it is clear from figure 6, the upper piece of S's boundary is smoothly curved. In other words, S is not a polytope, and so neither is the set of nearest-neighbour TI marginals. In order to show that it is also not a semi-algebraic set, see equation (5.1), we will assume that it is and prove the result by contradiction.

Conclusion
In this work, we have studied the problem of deciding whether a number of distributions correspond to the marginals of a 2D TI system, what we called the MARGINAL problem. We found that this problem is exactly solvable in scenarios of low local dimension and nearest or nextto-nearest neighbour statistics. We also showed that a natural variant of the problem, where we also demand symmetry under reflection, is solvable for all local dimensions. For other scenarios, we proposed a general algorithm to approximately solve the MARGINAL problem, as well as its dual, the ENERGY problem, where the goal is to minimize a linear functional of TI marginals. We also proved several no-go theorems concerning these two problems. We showed that the ENERGY problem is undecidable in general, so we cannot expect to identify the sets of TI marginals exactly for arbitrary d. We find that for d high enough, those sets are neither real polytopes nor semi-algebraic sets. Our techniques to prove negative results relied on a correspondence, proposed by Kari [24,25], between aperiodic tilings and immortal points of dynamical systems. Perhaps because of this, our upper bounds on the minimal dimension over which the set of TI marginals ceases to admit a simple description seem very poor. It is an open question how to lower those bounds. Could it be that, already for dimensions of order 10, we can experience the transition from a rational polytope to a convex object where parts of the boundary are smoothly curved? And, could it be that, for dimensions small enough, the description of the sets ceases to be a polytope but nonetheless admits a practical description via semidefinite programming?
Finally, we hope our methods and results can be applied to the study of other thermodynamical quantities or out-of-equilibrium systems. An immediate follow-up to our work would be to relate TI marginals to the maximum entropy per site of the whole lattice configuration from which they originate. That way, we would be able to compute interesting thermodynamical quantities of TI systems, such as the free energy, at non-zero temperature. Such maximum entropy extensions of TI marginals have been studied on the 1D Euclidean lattice and the Bethe lattice [5], but nothing is known when the lattice is 2D Euclidean. In this regard, our results on the uniqueness of the extendibility of the extreme points of d = 2 TI marginals suggest that a full characterization of the set of achievable TI marginals plus maximum entropy in that scenario is on the horizon.
Data accessibility. This work does not have any experimental data. Authors' contributions. Both authors contributed equally to this work and gave final approval for publication. Competing interests. We have no competing interests. Funding. This work was supported by the FQXi grant 'The physics of events'.