Extreme boundary conditions and random tilings

Standard statistical mechanical or condensed matter arguments tell us that bulk properties of a physical system do not depend too much on boundary conditions. Random tilings of large regions provide counterexamples to such intuition, as illustrated by the famous 'arctic circle theorem' for dimer coverings in two dimensions. In these notes, I discuss such examples in the context of critical phenomena, and their relation to 1+1d quantum particle models. All those turn out to share a common feature: they are inhomogeneous, in the sense that local densities now depend on position in the bulk. I explain how such problems may be understood using variational (or hydrodynamic) arguments, how to treat long range correlations, and how non trivial edge behavior can occur. While all this is done on the example of the dimer model, the results presented here have much greater generality. In that sense the dimer model serves as an opportunity to discuss broader methods and results. [These notes require only a basic knowledge of statistical mechanics.]


Contents
Perhaps one of the greatest strength of theoretical physicists is their ability to make simplifying assumptions which are obviously not true. These untrue assumptions are then used to make powerful predictions, some of which can be confirmed by experiments or even proved mathematically. For example, band theory assumes periodic boundary conditions to get a well defined momentum, and explains why certain materials are conductors, insulators or even topological insulators. To take an even older example, the hydrodynamic description of a glass of water relies on well defined conserved quantities such as mass density or momentum density. This is obviously not true, since the boundary breaks momentum conservation. However far from the boundary this assumption is much more reasonable to the leading order, which is why nobody is (and should) be worried about this. Similar assumptions are routinely made in classical statistical mechanics. For example, Onsager famously solved [1] the 2d Ising model, and found an exact formula for the free energy. What he did was assume periodic boundary conditions (or translational invariance), which simplify the calculations considerably. Then, his results for the free energy holds for any reasonable large chunk of the square lattice, since the free energy is well-known not to depend on boundary conditions. This simple but important observation is taught in any course on statistical mechanics. It also motivates introducing the dimer model, which we will discuss at great length in the lectures. One of the most striking property of the dimer model is that it totally violates the previous well-known fact of life. The free energy of the dimer model does depend, heavily, on boundary conditions.

Domino tilings and boundary conditions
Let us introduce the dimer (or domino tiling) model, on the square lattice. Edges connect nearest neighbors. Dimers are entities that cover edges of the lattice. The model then asks to cover the lattice with dimers, while following the following two constraints : a) A site cannot be touched by more than one dimer. This is called the hardcore constraint.
b) All sites have to be covered by exactly one dimer.
We then call dimer covering a valid configuration of dimers on the lattice. For example, all 5 possible coverings of the 4 × 2 square lattice are shown below.
We are of course interested in the thermodynamic limit, that is, dimer coverings of large M × L lattices. As we shall see shortly the dimer model is exactly solvable; for example, the partition function has been computed by Kasteleyn [2] and Temperley and Fisher [3]. They found that the corresponding free energy is in the limit L → ∞, M → ∞, L/M fixed. C is the Catalan constant, C/π 0.29156. Hence the number of possible coverings grows extremely fast when L, M are increased. Examples of dimer coverings on an L × L grid are shown in figure 1. In each picture, a dimer covering is picked uniformly at random, and drawn using a color code that we explain now. The lattice is bipartite, which means one can label lattice sites with two colors, black and white, in such a way that each black (white) lattice point has all its nearest neighbors of the opposite color. Then, we draw a vertical dimer in blue (resp. yellow) if its bottom part touches a black (resp. white) vertex. Similarly green (resp. red) dimers have their left part that touches a black (resp. white) vertex. At this stage the only purpose of this convention is to make the pictures look nicer. For the same reason we also make the dimers much thicker, so that they fill all space, and the underlying lattice cannot be seen anymore. Much later, in 1991, Ref. [4] studied dimer coverings of a peculiar region which they called Aztec diamond. They showed using combinatorial methods that the partition function, or number of dimer coverings, is given by the remarkably simple formula Therefore, there are much less available dimer coverings on the Aztec diamond that there would be on a regular square grid with the same area. We will see in the following that this is an effect of boundary conditions, so free energy does depend on boundary conditions for dimers. Even more remarkable is the following observation. Draw for reasonably large L a dimer covering out of the 2 L(L+1)/2 available ones, uniformly at random, as we did before. The corresponding pictures are shown in figure 2. The covering appears only random inside a region which looks roughly like disk for large L. Outside of the disk the orientations appear deterministic, with red dimers on the left, blue dimers at the bottom, etc. This is one of the Figure 2: The arctic circle theorem in action Dimer coverings of an Aztec diamond of order L, chosen uniformly at random. From left to right, L = 16, 128, 1024. As L increases, dimers appear totally frozen outside a region, which looks like a disk. There are, however, non trivial long range fluctuations inside the disk. most famous instance of what is called the limit shape phenomenon. We will give a more precise definition of limit shape later on, but for now, let us just introduce the terminology that we will use in these notes. In the deterministic region dimers are essentially frozen, so it is called frozen region. In contrast, dimer orientations fluctuate in the fluctuating region, sometimes also called liquid region. We will see later that correlations functions decay as power-laws, so the liquid region is critical. The interface between the frozen and liquid region is called an arctic curve. In the case discussed above, Ref. [5] proved that the arctic curve becomes an exact circle in the limit L → ∞, a result which now goes under the name arctic circle theorem.
There is no a priori reason why the arctic curve should be a circle or even smooth in this particular case, it just comes out of a long calculation. Lattice symmetries, however, do impose that the arctic curve be invariant under rotations by π/2. To illustrate this last point, one can generalise the model to include "interactions" between dimers. The simplest way to do that is to put weights favouring (or disfavouring) aligned dimers on a given plaquette [6,7]. In that case the partition function may be written as where N par counts the number of plaquettes with two dimers parallel to each other, that is, plaquettes of the form: or (6) Positive (negative) λ corresponds to attractive (respulsive) interactions between the dimers. For λ = 0 the arctic curve is not known analytically, but simulations shown in figure 3 clearly suggest that the arctic curve is very different from a circle. In fact, a closely related six vertex model has arctic curves which can be computed [8,9] and they are not even algebraic in general. See the bottom part of figure 3 for examples. Interacting dimers on the Aztec diamond with interaction set only on plaquettes of the even sublattice. This model can be mapped to the integrable six vertex model with domain wall boundary conditions (in which case a = b = 1 and ∆ = 1 − e λ ) [4]. Left: ∆ = 0.4, middle ∆ = 0, right ∆ = −6. The arctic curve is known but complicated away from the free point ∆ = 0, where we are back to a circle. Here we show in blue even flippable plaquettes, or odd empty plaquettes (this corresponds to c vertices in six vertex language). Note the appearance of a third region for negative ∆, where most vertices are of c type, and correlation functions decay exponentially. This region is usually called gas.
To finish this section let us mention that the pictures can be generated using a simple Markov chain Monte Carlo algorithm, which goes as follows in the free case. Start from any simple configuration, and then repeat the following lots of times: (i) Pick a plaquette uniformly at random.
(ii) If it has two horizontal (resp. vertical) dimers such as shown in (6), flip them to get vertical (resp. horizontal) dimers. Otherwise do nothing.
After lots of updates, the Markov chain will thermalize and show typical configurations, which were shown above. Of course, the number of necessary updates increases quickly with system size, so generating the best pictures might take a while for large L. There are more complicated (but more efficient) algorithms, see e. g. [7,10,11].
1.2 Detour: a one-dimensional classical Coulomb gas model The limit shape phenomenon may be illustrated by the following simple 1d model, which contains the basic ingredients needed to get a inhomogeneous density profile with frozen region. We consider N particles living in a box B = [−1, 1]. The probability of having the particles at positions x 1 , . . . , x N ∈ B is given by where β is a positive real number. This probability density function may be interpreted as a Boltzmann weight P (x 1 , . . . , is the electrostatic energy of a gas of N particles with 2d Coulomb interactions. For this reason it is often called Coulomb gas 1 . We will consider two versions of the model: • The first is the model as stated, with partition function This is well known from random matrix theory, where it is called (a limit of) the Jacobi ensemble.
• In the second we impose that the allowed positions for the particles be discrete. The particles now live in B L = {−1 + 2j L−1 j = 0, 1, . . . , L − 1}. The partition function reads This type of models usually go under the name discrete beta ensembles.
We are interested in the average distribution of the charges, in the limit N → ∞. In the second model, we further impose a fixed density d = N/L. Formally, the first may be recovered from the second by considering a low density limit d → 0. We study the density profile in both models separately. Before proceeding let us point out the following important fact: the interaction is long range so energy is of order N 2 , while entropy fluctuations are expected to be of order N . This is not the standard situation in thermodynamics, where there is typically a competition between energy and entropy. Here energy dominates, which means the average density profile (limit shape) can be obtained just by minimising energy. Hence, the precise value of β does not matter here. Of course, fluctuations on top of the limit shape are still there. They are more complicated, especially in the discrete case.

The continuous gas
The continuous gas may be treated using standard techniques [12]. Introducing the density , the energy may be rewritten, up to an unimportant additive constant as In the thermodynamic limit, ρ is expected to become a smooth function. Therefore, finding the equilibrium distribution of the charges boils down to minimizing the energy functional subject to the constraint Handling the constraint can be done by introducing a Lagrange multiplier λ. We consider the functional and write down the Euler-Lagrange (EL) equations The second equation gives back the constraint, while the first reads One can check that ρ(x) = λ 2π log 2 1 √ 1−x 2 is a solution to the above linear integral equation. Finally, normalization of ρ yields λ = 2N log 2. Hence we get the density profile This density, called the arcsine law, is integrable but unbounded near x = ±1. This means the particles tend to accumulate near the two boundaries, a non trivial effect of the longrange repulsion between them. Let us finally mention that this type of problem has been studied for a very long time, in relation to potential theory as well as orthogonal polynomials. For example, Stieltjes observed back in 1885 [13,14] that the positions of the particles that minimizes the electrostatic energy (8) on B = [−1, 1] coincide with the zeroes of known orthogonal polynomials 2 . In this context, the limiting distribution (18) was probably known even before. See also exercise 1.3.
2 More precisely, the equilibrium positions for the pdf N i=1 (1 + xi) a (1 − xi) b j>i (xi − xj) β on B N are the N roots of the polynomials P  ρ(x) x N = 16 limit shape Top: a typical configuration for N = 16 particles on B. Bottom: density profile for N = 16 particles (Monte Carlo simulation) and comparison with the thermodynamic limit (18).

The discrete gas
The discrete model can be treated using a similar method. Introducing the particle density ρ x = i δ x,x i , we expect that it becomes a continuous function ρ(x) normalized to B ρ(x) = 2d in the limit L → ∞, N → ∞, with fixed total density N/L = d. We still have to minimize the energy functional E[ρ] of (12) with the constraint (13). However, due to the discrete nature of the problem, density cannot exceed one, since two particles cannot sit on the same site. This means we get a second constraint This extra constraint is typical for discrete models. The previous solution (18) is unbounded, which means a new analysis is needed. The solution to minimizing (12) with the constraints (13) and (19) was found in Ref. [17]. Let us first give the result and comment on it later. The limit shape is given by As in the continuum the Coulomb interaction wants to make the density very large close to the edge. This is not possible, due to the constraint (19). Hence, the system compromises by having the maximum allowed density over a larger region In orthogonal on B with respect to the measure dµ(x) = (1 + x) p (1 − x) q dx (p, q > −1). The interested reader may consult [15] and references therein for a broader overview. The case α = β = 0 which we are dealing with here is degenerate [16], since the measure is not integrable; in this case the equilibrium positions of the particles are the zeroes of the polynomial (1 − x 2 )P this region the position of the particles become deterministic, we call this the frozen region in analogy to dimers. The other region is fluctuating. See figure 5 for an illustration. Let us now explain how to find the solution. The crucial point is that the second constraint (19) may be temporarily lifted by assuming that the density equals one in a given region x ∈ [−1, −r] ∪ [r, 1] close to the edges. Then, we have to minimise the new functional where the potential simply comes from the interaction of the charges in [−r, r] with the charges in the frozen region (in (21) there is also an additive constant E r which accounts for self-interaction in the frozen region). The constraint on the new equilibrium measure now read Now, depending on the value of r, the solution of the minimization problem without imposing the second constraint may well satisfy (24) anyway (e. g for r close to one it does not, while for r = 1 − d it trivially does, since there is no mass left in ρ r ). The solution to the minimization problem (21) with only constraint (23) can be found by writing down the Euler-Lagrange equation once again. Denoting by ρ r the solution, we choose the value In words, this is the largest r such that extra constraint is automatically satisfied. The solution to the minimization of (21) with constraint (23), (19) is then ρ(x) = ρ r 0 (x) if |x| < r, ρ(x) = 1 otherwise. This yields (20). For a rigorous proof, we refer to [17].
To finish this section, let us finally mention that some examples of limit shape problems in 2d statistical mechanics can also be mapped to a similar discrete coulomb gas problem in 1d, as was pointed out in [18]. In the following we will not pursue this direction, however, and aim for a more general hydrodynamic theory.
Plan of the rest of the lectures. It looks like a good idea to try and apply the general ideas we used to solve the discrete Coulomb gas model. However, it is not clear at this stage what we should minimize exactly to obtain the limit shapes in general. This question is addressed in section 2, where we introduce the height mapping and use that to understand which functional to minimize. Then, we actually compute this in section 3, for dimers on the hexagonal lattice. We use the transfer matrix formalism, and a mapping onto free fermions (see appendix A for a long description of free fermions techniques, if need be). The choice of the hexagonal lattice is motivated by technical simplicity, even though the square lattice is similar. Then we use that to find the limit shapes in section 4. Section 5 deals with exact lattice calculations, which allows to recover the arctic curves in a few selected cases. Finally, we discuss a few related and more complicated problems in section 6.1, and conclude.
1. By interpreting |Ψ(z 1 , . . . , z N )| 2 as a pdf for the N particles, what is the limit shape in the limit N → ∞? 2. The N particles are now constrained to live on the sites of a L × L square lattice with mesh a = 1 (the origin z = 0 is set at the barycenter). What is the limit shape in the limit N → ∞ with fixed density d = N/L 2 when d is not too large? What happens at higher densities?
Zeros of orthogonal polynomials [Stieltjes 1885] Consider the electrostatic energy considered before, with two extra charges at positions ±1, E(x 1 , . . . , for a, b ≥ 0. For large N , the limit shape does not depend on a, b.
1. Why? In the following, we set a = b = 1/4. 2. Write down a system of equations satisfied by the positions y 1 , . . . , y N that minimise the electrostatic energy.
Show that the previous system is equivalent to 2 . Show that p N (cos θ) = cos(N θ). 5. Find the zeroes of p N and recover the limit shape (18) in the limit N → ∞.

Variational principle
We show in this section that the right quantity to minimise is a variant of the free energy. To understand this, we first introduce an important ingredient, the height mapping.

The height mapping
Consider the dimer model on the square lattice. Remember the lattice is bipartite, which means one can label sites with two colors, black and white, in such a way that each black (white) lattice point has all its nearest neighbors of the opposite color. To each dimer configuration we associate a height configuration, as follows. Heights are discrete numbers (integers on figure 6) which live on plaquettes (or the dual lattice, which is also square). We pick a reference point, say the leftmost bottom plaquette, and set its height to zero. Then, turning counterclockwise around a black (resp. white) vertex, the height picks up +3 (resp. −3) when crossing a dimer, −1 (resp. +1) otherwise. A dimer configuration with the corresponding height configuration is shown in figure 6 on the left. To make the mapping clearer we use the following color code. We draw a vertical dimer in blue (resp. yellow) if its bottom part touches a black (resp. white) vertex. Similarly green (resp. red) dimers have their left part that touches a black (resp. white) vertex. Hence crossing a blue dimer from left to right, the height always picks −3, etc. Height mapping on the square lattice, in units of π/4. Left: heights corresponding to a dimer configuration picked uniformly at random. Observe the heights remain on average quite close to zero. Right: heights corresponding to a alternating boundary condition with zigzag blue dimers. Some other vertical dimer occupancies are automatically set as a result, those are shown in lighter blue. The corresponding mean horizontal slope is maximal, ∂ x h = π/2. The gradient in the vertical direction is also automatically set to ∂ y h = 0. In fact, in general the slopes satisfy |∂ x h| + |∂ y h| < π/2.
At the very end and for later convenience, all heights are remultiplied by π/4. The mapping is one to one, and has many nice properties that we will investigate in the following. For the moment let us just say that mapping to discrete heights has a long history in statistical mechanics, which dates back to Ref. [21,22]. In fact, the model studied in these references can be mapped onto dimers on the hexagonal lattice, which we will study later on.

Minimizing the free energy
From the height mapping, we know that dimers like to be in configurations where the height gradient is close to zero. A good example is provided by a rectangular domain, as illustrated in figure 6 on the left. Boundary conditions might spoil that, however. The reader can easily check that the Aztec diamond geometry of figure 2 can be cooked up by imposing the following boundary condition at the bottom, and a similar one at the top. In that case the slopes are maximal along the Aztec diamond, and alternate −π/2, π/2, −π/2, π/2 from one boundary to the other.
These extreme boundary condition play an important role, and are, as we shall see, responsible for the appearance of the arctic circle. Indeed, the number of available dimer coverings in a given region strongly depends on the average slopes imposed at the boundary, as can already be guessed by looking at right part of figure. 6.
A way to see this is to use the fact that all possible slopes are allowed when covering a torus with dimers, as can be seen in figure 7. So consider a × torus ( even), and associate a weight a = e µπ/ for yellow dimers, 1/a for blue dimers, b = e νπ/ for green dimers, 1/b for red dimers. The corresponding partition function may be written as Above, Z r,s precisely counts the number of dimer coverings in the (r, s) sector, where r is the mean horizontal slope ∂ x h and s the mean vertical slope ∂ y h . The total number of dimer coverings is Z(0, 0). Hence Z(µ, ν) encodes all the information about the Z r,s . For example, the generating function reads for a 8 × 8 torus (311853312 coverings in total, we will explain later how to calculate this): The corresponding free energy is, in the thermodynamic limit,  (27). Note that we refrain from identifying the bottom/top and leftmost/rightmost plaquettes here, for the sake of the argument.
The crucial point is that once r and s are fixed, the resulting free energy does not depend anymore on the details of the boundary conditions. So we are back the to the situation discussed in the introduction for the Ising model, provided r and s are fixed. This free energy (or 'surface tension') has many fascinating properties, and will play a central role in the following. For now let us just mention that it is minimum at (r, s) = (0, 0); in most cases it is also a convex function. With the free energy F (r, s) as a given, one way to understand the limit shape phenomenon is to consider a very large lattice, say L × L, and cut it in several much smaller (say square) cells of size × , where is still much larger than the lattice spacing (= 1 for us). Namely, our system is macroscopic, and we look at it at mesoscopic scales where (i) it can be described in the continuum (ii) it is uniform. The total system is then considered to be a collection of smoothly connected uniform cells, each with its own average boundary gradient ∇h, which is imposed by the surrounding cells.
The system will then try to minimize total free energy (maximize number of dimer coverings), by choosing the appropriate slopes in each mesoscopic cell. Hence we need to minimise the functional (or action) over some domain D. The Euler-Lagrange equations for this variational problem read where we have introduced the notation for the partial derivatives of F , and h = h(x, y). It also possible to add extra constraint in a similar fashion to what we did in section 1.2.
On physical grounds this variational (or hydrodynamic) principle is expected to hold under very general conditions, and has been used for a long time in the context of crystal surfaces in statistical mechanics (see e.g. [22,23]). For dimers the situation is even more favorable, since the validity of the variational principle is now a theorem [24].
Therefore, the limit shape is given by the solution to the PDE (30) with appropriate boundary conditions. Computing exactly the free energy can also be done [22,24], we will explain in the following section how to. It should be stressed, however, that even with an exact expression for the free energy, solving the PDE for arbitrary conditions is in general a difficult problem.

Fluctuations and free field theory
We have so far only discussed the limit shape, which gives the average height field at position x. There are of course fluctuations on top of that, which are, we argue here, described by a massless free field theory. Those fluctuations are in particular responsible for the power-law decay of correlation functions.
Homogeneous case. To discuss fluctuations, let us first consider the simpler case of the rectangular geometry (or really, any simply connected finite planar domain), for which ∇h is close to zero on average (see figure. 6), and the slopes in the neighborhood of (r, s) = (0, 0) dominate. Due to lattice symmetry considerations F (±r, ±s) = F (r, s) and F (r, s) = F (s, r), so F (10) and F (01) both vanish, and to lowest non trivial order we have So, neglecting higher order corrections, the free energy is determined up to a single unspecified parameter K > 0. K has many names 3 , in the following we call it the Luttinger parameter.
The Euler-Lagrange for the average height take a simple form in that case: so h is a harmonic function with appropriate boundary conditions on ∂D. The solution is sometimes called the classical (or zero mode) part of the field, and noted h cl . In fancier terms h cl (x, y) is the harmonic extension of h ∂ , its value at the boundary, to D. On top of that, fluctuations may be handled by writing where δh satisfies Dirichlet boundary conditions on ∂D (the factor 1/2 is set to match standard conventions later on). Plugging (32) and dropping F (0, 0) yields an action where the linear term drops out due to the Euler-Lagrange equations, combined with the fact that δh obeys Dirichlet boundary conditions (integrate by parts). The action for the fluctuations is The path integral formulation reads where δh satisfies Dirichlet boundary conditions, or, equivalently, where h = h cl at the boundary. The above is the euclidean action of a massless free (or gaussian) scalar field in two dimensions, in the following we will simply refer to it as the free field 4 . It is always desirable to visualize things, see figure 8 for two realizations of the discrete height field for dimers. The reader can then try to imagine what this becomes in the continuum limit.
Before heading to the inhomogeneous case several important remarks are in order.
• The argument we just provided is quite standard [25]; in fact, the very reason field theory techniques may be applied to statistical or condensed matter systems is that this type of reasoning often just works, and provides exact results for critical exponent and long range correlation functions. However, a proper derivation from a concrete lattice model is very often a difficult task.
• As can be seen from the figure, height field configurations look wilder and wilder as system size is increased. In particular, its variance at any point can be shown to diverge as ≈ log L. The free field is, in fact, a singular object. On a simply connected bounded planar domain D a possible mathematical definition is as follows. Consider (minus) the Laplacian −∇ 2 on D with Dirichlet boundary conditions. We write its normalized eigenfunctions as u k (x, y) and the corresponding eigenvalues as λ k (they are all strictly positive). Then introduce where the ξ k are independent centered gaussian random variables with variance E[(ξ k is the Green's function for the laplacian on D with Dirichlet boundary conditions, which means E[ϕ(x, y)ϕ(x , y )] = δh(x, y)δh(x , y ) , where the braket on the right denotes expectation value with respect to the path integral (37). At the level of correlations functions, δh and ϕ are essentially the same object.
The series (39) in fact, converges almost nowhere. Therefore, the free field has to be seen as a random distribution, not a random function 5 . The interested reader can have a look at References [26][27][28][29][30][31] for a mathematical treatment of the free field. For computations what matters is to be able to integrate against smooth test functions, so this is not a problem. It is, however, nice to keep in mind that the free field is fundamentally a singular object. transition [21,34] of BKT-type [35,36] to a phase where the height field becomes regular. In that case long range power law correlations are lost.
• The reader might be tempted to point out that the boundary heights are flat on average for the dimer model in a rectangular domain, so that it is reasonable to impose Dirichlet boundary conditions h ∂ = 0, and safely conclude that h cl (x, y) = 0 ∀x, y ∈ D since this is the only harmonic function satisfying the boundary conditions. This is not quite true, in fact the correct boundary average heights alternate (±π/8), and this slight change has an impact on some observables in the continuum limit. h cl is calculated in exercise 2.2.
• The action (36) is one of the simplest example of a conformal field theory. This becomes transparent when rewriting the action in terms of complex coordinates One can easily check that the action is invariant under any transformation z → g(z),z →ḡ(z), where g (ḡ) is any holomorphic (antiholomorphic) function. All such transformations preserve angles locally, they are conformal. The action (36) is then conformaly invariant. Conformal field theory is a vast subject, we will barely scratch the surface in these notes. We refer to [37][38][39] for reviews.
Inhomogeneous case. The inhomogeneous case is slightly more complicated, due to the fact that the classical solution solves a more complicated PDE (30). The free energy is still expected to be a strictly convex function, at least over a wide range of slopes. This means the determinant of the Hessian matrix is strictly positive. This allows to still define the Luttinger parameter, as 1 πK(r, s) where H is the Hessian matrix The (inverse) of the Luttinger parameter tells us how convex free energy is, so still measures stiffness. The smaller the K the more energy slope fluctuations on top of the classical solution will cost. By repeating the same arguments as before, we get and This can be written in covariant form as where g = H −1 is an emergent metric. It is a priori non flat, since H = H[r, s] depends on x, y through r = ∂ x h cl (x, y), s = ∂ y h cl (x, y). What we just wrote is an action for the fluctuation δh in a curved metric, which itself is determined from the limit shape h cl (x, y). A sample of the discrete height field is shown in figure 9, and illustrates what we just said. • In the covariant form of the action (45 , so the Luttinger parameter depends on position in general. This has important conceptual consequences, in particular it means that conformal invariance is broken in the general inhomogeneous case. For dimers we will see in fact that K = cst = 1, and explain this as a general property of models that map to free fermions. • The argument leading to (44) is appealing but oversimplified, for the following reason. As in the homogeneous case, we have used the EL equations to get rid of the linear term in the expansion. However, it is important to keep in mind how all those terms scale with L. The quadratic term is subleading (by 1/L) compared to the linear one, and EL only implies that the leading part of the linear term vanishes, but does not tell us anything about the lower order part, which would potentially break the quadratic nature of the action. As explained in Ref. [40], a more precise argument considers the dependence of the free energy on higher derivatives ∂ 2 h of the height field. Then, redoing the previous steps leads to a quadratic action identical to (44), albeit with an extra (1/L subleading compared to the limit shape) contribution to the average of the height field. This is not shocking if one compares to the homogeneous case, where the limit shape is zero to leading order, but the classical solution h cl is O(1) and still matters.
• This is as far as we can go without knowing the specific form of the free energy. We compute it in the next section in the simplest example, that is dimers on the honeycomb lattice.
• Of course, obtaining exact expressions for correlations on the lattice are also very much desirable, to check our hydrodynamic assumptions and their limitations. This approach is discussed in section 5.

Exercise 2.1 Sample of the free field
Using your favorite programing language, draw an (approximate) sample of the free field on a rectangle [0, π] 2 with Dirichlet boundary conditions.

Exercise 2.2
The even-odd effect [Inspired by [41][42][43][44][45]] We consider interacting dimers on the L x × L y square lattice, in the rectangular (figure 6) and cylinder geometries. L x and L y are assumed to be even. We wish to compute R(a) = Z(Lx,Ly)Z(Lx,Ly) Z(Lx,Ly+1)Z(Lx,Ly−1) in the limit L x , L y → ∞ with fixed aspect ratio a = L y /L x . To do that we assume that the only relevant contribution to the ratio of partition functions is that of the classical part e −S 0 [h cl ] , with the action (43). (meaning contributions from fluctuations are the same for even and odd). 1. Explain R(a) = O(1) from general arguments. Surprinsingly R(a) = 1, see below.
We deal with a cylinder of circumference L x and height L y first. 2. What are the possible height differences between top and bottom contributing to Z(L x , L y )? Same question for Z(L x , L y ± 1). 3. For each case, find the set of harmonic functions on the cylinder, that implement these height differences. [Hint: those are really simple functions.] Show then that We now treat the (more technical, use Maple/Mathematica) case of a rectangle. 5. Compute the average heights on all sides in the even and odd case. 6. Show that the real and imaginary parts of a holomorphic function are harmonic. 7. Show that the Dirichlet energy D dxdy(∇h e ) 2 −(∇h o ) 2 is invariant under conformal maps. The conformal map from the upper-half plane H = {z, Im z > 0} to D is given by the . For free dimers (K = 1) the result can also be extracted from the lattice [41].

Transfer matrix for dimers
We show here how the dimer model is equivalent to a system of free fermions, through the transfer matrix formalism. Before explaining the mapping a few remarks are in order. First, the method we present is not the only one. The most widely used for dimers relies on work of Kasteleyn [2], who was the first to solve the model, expressing e.g. the partition function for dimers on planar graphs as Pfaffians (see also [3]). In particular, many result on the mathematical side rely almost exclusively on Kasteleyn theory. The two methods are not that different, see e.g. Ref [46], and lead to essentially the same results. Our choice has the advantage of connecting to techniques better known to physicists, in particular free fermions, similar to the (simplified version [47,48] of) the Onsager transfer matrix [1] of the Ising model. This choice is also motivated by possible generalizations to interacting integrable systems such as six vertex model, where the transfer matrix method cannot really be avoided. The first mapping of dimers onto free fermions is due to Lieb [49]. The version we present here is slightly different, and leads in a transparent way to a Hermitian transfer matrix with conserved number of particles, in the spirit of Refs. [7,50,51]. For simplicity, we mainly focus on the dimer model on the hexagonal lattice, not square, even though the two are similar.

Reminder on the transfer matrix formalism
In this whole section we generalise the honeycomb dimer model slightly by putting alternating weights . . . , 1, u, 1, u, . . . for horizontal dimers along a given horizontal line (see figure 10, where the lattice is drawn as a brick wall). The first important observation is that a given dimer configuration on the lattice is uniquely determined by the occupancies along vertical edges, so we may completely ignore the occupancies of horizontal ones. Put now a 0 on vertical edges In thinner dark red are drawn the vertical edges not occupied by dimers, those become particles (1) in the following, while real dimers are holes (0). occupied by a dimer (shown in thick blue on the figure), and a 1 otherwise (thiner darkred). We see the ones as a collection of particles propagating upwards, from bottom to top. We then associate a basis vector to the vertical dimers occupancies along a given horizontal line.
For example for L = 6 in the picture, to the bottom line configuration 011001 we associate the vector |011001 . The scalar product is C|C = 1 if the two line configurations C and C are the same, zero otherwise. There are 2 L possible line configurations. Imagine now that it is possible to find a 2 L × 2 L matrix T , called the transfer matrix, such that C|T |C = u n if the configurations of two successive horizontal lines are compatible when stitched together -meaning they form a valid dimer covering-and zero otherwise. Here n is the number of horizontal dimers with a u weight when stitching the two horizontal lines. Then the partition function on the L × M lattice is simply Z = top|T M |bottom , where top| and |bottom are the bras and kets corresponding to the the top and bottom configurations respectively. Proof.
Each element in the sum is one for valid configurations, zero otherwise, so this just counts the number of dimer coverings compatible with the top and bottom boundary condition. It is also possible to require the bottom and top boundaries to coincide (periodic boundary conditions), in which case Z = Tr T M .
Here we actually need two transfer matrices T and T , since the rule changes depending on the parity of the row considered. Let us again consider the example in the figure.

Dimers as free fermions
It is important to realize that both transfer matrices (TMs) preserve the number of particles, so those are block diagonal, each block indexed by the number of particles (or ones).
Transfer matrix in the one particle sector. This one is not much more difficult. The corresponding matrix elements are Here i, j are the indices corresponding to the position of the (unique) one. A is similar.
Full Transfer matrix. For more particles the rules become more cumbersome, due to the dimer hardcore constraint, which prevents two 1 from occupying the same site. This is where the power of the free fermion formalism comes in. The reader not well acquainted with these techniques is invited to have a look at the self-contained introduction in appendix A, which has everything needed to understand the present lectures. If you already know (170), however, chances are you don't need to read it. First, let us represent vectors |C using the fermion formalism. We use free fermions operators to represent all states. For example for L = 4, we have |1101 = c † 1 c † 2 c † 4 |0 , where |0 = |0000 is the fermion vacuum. For any dimer configuration C, the associated ket reads in the configuration C. Note that since fermions anticommute, applying the operators in a different order might generate undesirable minus signs, e.g.
. Now the main claim is (see e.g. [25]) provided the logarithm of A makes sense. [For T just replace A by A in (47)]. Hence both TMs are the exponential of a Hamiltonian that is quadratic in the fermions operator, a free fermions Hamiltonian. Sketch of the proof. It goes in two steps. First, let us show that a sufficient condition for a good TM is that This obviously works for n = 0, 1 particles. One then needs to determine the action on states c † i 1 . . . c † in |0 , i 1 < . . . < i n in the n-particle sector, for n = 2, . . . , L. It is obtained by commuting T successively with all fermions operators using (49), and then using (48): The sum generates all possibilities for the particles to go upward-left or upward right. The jumps are not independent however, if two particles go to the same site i, the corresponding contribution is proportional to c † i c † i = 0, which means fermion operators effectively enforce the dimer hardcore constraint. [This observation also justifies the terminology free fermions: the particles interact only through the Pauli exclusion principle, which prevents two fermion from being in the same state (occupy the same site). In that sense free fermions are as free as fermions can ever be.] Importantly, nonzero contributions to the sum in (50) are still all ordered, since only nearest neighbor jumps are allowed in A. Hence all valid configurations are counted with the correct (+) sign.
The second step is to realize that (47) satisfies (48), (49). This is essentially formula (168) in appendix A, go there for the proof.

Various subtleties.
a) The reader might be worried that the matrices A, A contain Jordan blocks, so are not diagonalizable. While the derivation does not require A, A to be diagonalizable, this is still a nuisance. An easy fix is to consider T T and call it the transfer matrix. Using the identity (169) we have B (and log B) are now symmetric, so the Hamiltonian inside the exponential is a legitimate quantum mechanical Hermitian operator. In particular, one can use standard band theory techniques to examine its long range properties. All of them are solely determined from T T , which can be considered to be the true transfer matrix. In the following, we call T T the transfer matrix for the dimer model. b) Care must be taken when implementing periodic boundary conditions. Indeed, . Indeed, remember the correspondence with the stat mech model imposes that fermion operators be applied in order. PBC spoil that natural order when the number of particles is even. Hence the transfer matrix has to satisfy, T c † L = c † L + (−1)N −1 c † 1 T , whereN = L j=1 c † j c j is the total fermion number operator. In practice, this means it is necessary to consider separately even and odd fermion sectors, in order to keep the quadratic nature of the Hamiltonian. c) Diagonalization. As explained in appendix A.3, diagonalizing T T boils down to diagonalizing the L × L matrix B, which is a huge simplification. We obtain where the ε(k) are the eigenvalues of − 1 2 log B, Ω is the set of labels for those, and where v jk are the normalized eigenvectors of B, j v jk v jq = δ kq . This implies For dimers on the honeycomb lattice with open boundaries as described in this section, it is possible to get them explicitely. One possibility is to make the ansatz v jk ∝ sin(kj + γ), and check that those provide unnormalized eigenvectors of B provided γ = 0 and sin[k(L + 1)] + u sin[kL] = 0. Ω is then the set of the L solutions to the previous equation in the interval (0, π), and The minus signs in (52), (55) are here to mimic usual band theory, where one wants to minimize energy and look for ground states, the factor 2 to remember that our transfer matrix T T moves by two steps. The case of PBC turns out to be easier, and the reader can check that the set of allowed momenta is where α = 1 for evenN and α = 0 for oddN . This will be used in section 3.5.

Hamiltonian limit and quantum spin chains
We have just seen that dimers on the honeycomb lattice can be understood as a free fermion system, with dispersion (55). A similar result holds for the square lattice. For 0 < u < 1, the dispersion relation is analytic in [−π, π] so leads to a local Hamiltonian, that is, the hoppings decay exponentially fast with distance. Let us write T (u) the corresponding transfer matrix. It is easy to see that lim u→0 where H is a free fermion Hamiltonian with dispersion ε(k) = − cos k, that is a tight binding model with only next nearest neighbor hoppings. As is well known, this model also corresponds to the spin-1/2 XX chain [52]. Eq. (57) is called a Hamiltonian (or Trotter) limit, so the XX chain is a Hamiltonian limit of dimers on the honeycomb lattice (this works for the square lattice too, see figure 11). A way to think of this limit is to consider a finite L×M lattice, and notice that u essentially controls the hopping rates to the left and right. Then multiply the vertical size by an integer p and divide u by p. We have a L × (pM ) lattice with hopping rate u/p. The limit p → ∞ is nontrivial because while the hopping rate goes to zero, dimers are given more opportinities to hop since the number of application of the transfer matrix increases. It is also convenient to make the lattice spacing in the vertical direction proportional to 1/p, so that the distance between top and bottom stays the same in that limit.
The Hamiltonian limit is in a sense intermediate between discrete and continuous. The horizontal direction stays discrete, but the vertical one becomes continuous. Hamiltonian limits are key to the application of integrability techniques to quantum spin chains [53,54]. To finish this section let us mention the two most important examples of Hamiltonian limits: the XXZ spin chain is the Hamiltonian limit of the six-vertex model (interacting dimers), while the Ising chain in transverse field is the Hamiltonian limit of the 2d classical Ising model.

Connection to the height mapping
The height mapping for dimers on the honeycomb lattice is very similar to that of the square. Turning counterclockwise around black (resp. white) vertices, the height picks a factor +2 (resp. −2) when crossing a dimer, −1 (resp. +1) otherwise. As before, we use a color code to make the mapping more readable: vertical dimers always belong to the same sublattice, so we color them all in blue. Horizontal dimers are shown either in red or green. The rules are illustrated in figure 12. For future convenience, we also remultiply all heights by π/3 at the end.
The relation to the free fermion mapping is as follows. First, notice that the minimum slope in the horizontal direction correspond to all edges occupied by vertical dimers, so ∂ x h = − 2π 3 . The corresponding fermion density is n x = 0. The maximal slope is ∂ x h = π 3 and corresponds to fermion density n x = 1. In the vertical direction, the minimum slope is − π 2 while the maximal is π 2 . Notice again that the two slopes are not independent, ∂ y h takes values in ∂ x h is simply related to the fermion density, but how do we control ∂ y h? The simplest way to do that is to assign different weights to red and green dimers: starting from now we assign a weight b to red dimers, and a weight b −1 to green dimers. It is easy to see that the power of b in the expansion of the partition function allows to determine ∂ y h. Then, the corresponding transfer matrix has dispersion ε(k + iν), where b = e ν/Lx . Due to ε(−k) = ε(k) the transfer matrix is still normal and real, but not symmetric anymore. There is an alternative way to see this. Consider the fermion density n x (y) = e yH c † x c x e −yH in a Heisenberg picture. In the continuum limit 6 , it satisfies the continuity equation where j x (y) is the associated (imaginary 7 ) current. Since n x (y) = 1 π ∂ x h + 2 3 , we get ∂ y h = −πj x (y), so the current gives access to the vertical slope of the height field. Then, one can check that the total current operator J = x j x reads and has expectation value in the Fermi sea [−k F , k F ]. For the continuum Fourier operators c(k) and c † (k) we used the conventions (179).

Torus partition function and exact free energy for free fermions
We have now all the ingredients needed to compute the free energy. The generating function for all slopes on the L × M torus is given by This can be evaluated in closed form, using formula (167), and being careful about the point discussed in c). We find where and Ω α = (2m + α)π L , m = 1, . . . , L .
The reader can check that Ω 1 (resp. Ω 0 ) allows to generate all eigenvalues in the even fermion sector (resp. odd fermion sector), which has antiperiodic boundary conditions (resp. periodic boundary conditions). The leading asymptotic behavior may be determined by picking any 8 6 It is also possible to do the exercise in the discrete setting, and compute nx(y + 1) − nx(y). 7 It is important to realise that the imaginary current is not a hermitian operator, since it is i times the the real current. The consequences of this are not innocent. For example, an easy way to generate a real current is to consider a boosted Fermi sea [−k l , kr] with k l = kr. Our current operator takes imaginary values in a boosted Fermi sea, while we want it to take real values in order to set ∂yh. 8 Except e. g. when u = 1, ν = µ = 0 where Z − 0 vanishes identically, since one of the ε(k) is exactly zero. It can only happen to one term at a time, so does not matter. of the four terms (64). It is then easy to see that the only term that matter are those for which the argument of the exponential is strictly positive. Hence setting M = L for convenience, we obtain The Fermi sea FS is the set of k for which the real part of the integrant is negative, in our case it is a single symmetric interval [−k F , k F ]. Therefore, the generating free energy is simply the ground state energy corresponding to the dispersion ε(k + iν) − µ. One can check that it is real. The last step to get the asymptotic behavior of Z r,s is to choose µ, ν in such a way that the term e µr+νs Z r,s dominates in Z(µ, ν). In the thermodynamic limit the corresponding free energy F (r, s) is given by where µ and ν are determined from That is, F (r, s) is the Legendre transform of f (µ, ν), with r (resp. s) conjugate to µ (resp. ν). We finally obtain the following fundamental result where k F and ν are determined from (r, s) through Hence, the free energy is solely determined from the dispersion relation, extended to the strip [−π, π] + iR of the complex plane. From this formula it is also obvious that (71) is not specific to the honeycomb lattice, but applies to any one-band fermion problem, that satisfies ε(−k) = ε(k). The result can also be generalized to several bands problems, but we do not investigate this here. Let us now compute this in two important examples.
XX chain in imaginary time. Let us start with the simplest, which is the PNG droplet, or XX chain in imaginary time. The dispersion relation is ε(k) = − cos k, which makes explicit computations easy. We find F (r, s) = s π arcsinh s cos r − 1 π cos 2 r + s 2 .
It is shown in figure 13 on the left. Now r = k F − π/2 -notice the difference compared to dimers-so r ∈ [−π/2, π/2], while s ∈ R. The constraints between r, and s are relaxed in this limit, except for the fact that the free energy is infinite when r = ±π/2, s = 0.
Dimers on honeycomb. In this case integration is a little bit more tedious. We obtain after some algebra where recall ν solves (73) and dt is the dilogarithm. The free energy is shown in figure 13 on the right. To our knowledge, all explicitely known free energies can be expressed in terms of such dilogarithms [22,24,55].

K = 1 for free fermions
We are now in a position to compute the Hessian of the free energy, and check our previous claim that K = 1 in general for free fermions. The calculation presented below is essentially that of [56]. The first partial derivatives are where we have introduced z = k F + iν. Eq. (73), or 2s = i(ε(z) − ε(z)) was also explicitely used to derive the second equation. Now for the second derivatives: F (11) can be computed in two ways, either ∂ r F (01) (r, s) or ∂ s F (10) (r, s), so Hence where we have again used (73) to get the last line. So K = 1, independent on r and s.

Exercise 3.1 The log gas with free fermions
Let ψ(x 1 , . . . , i V (x i ) be a normalized wave function on the real line, x i ∈ R.

Consider the state
is a polynomial in x of degree at most k − 1. Under which condition do we have {d k , d † q } = δ kq for k, q ∈ {1, . . . , N }? We assume it is satisfied in the following.

Solving the minimization problem
From the previous sections, we have now all the necessary ingredients to solve the variational problem.

Complex Burgers equations
Recall the EL equations are ∂ x F (10) + ∂ y F (01) = 0, which read for free fermions, using (71,72,73), where k F and ν depend on position x, y. It is possible to express k F , ν in terms of r, s but this is not necessary. A convenient way to proceed is to introduce as we did before. In terms of those, the EL reads Since r, s are by definition partial derivatives of the height field, we also have the continuity equations ∂ y r − ∂ x s = 0 which may be also rewritten as ∂ y (z +z) − i∂ x (ε(z) − ε(z)) = 0. Combining with (88) yields the (conjugate of each other) equations [56] i∂ y z + ∂ x ε(z) = 0, −i∂ yz + ∂ x ε(z) = 0.
For ε(k) = k 2 /2 (free Fermi gas), this PDE is called complex Burgers equation, in the following we refer to (89) as a complex Burgers-type equation. An alternative but equivalent point of view is discussed in section 4.
The interpretation of these equations is very nice. Think of a free fermion system close to equilibrium, described by a (boosted) Fermi sea [−k l , k r ], where k l,r are the left (or right) Fermi momenta. The particle density is ρ = kr+k l 2π , while the current is J = ε(kr)−ε(k l ) 2π . Under unitary time evolution, they satisfy the continuity equations ∂ t k l,r + ∂ x ε(k l,r ) = 0. Since ∂ x ε(k) = (∂ x k)ε (k), this is essentially the statement that each quasiparticle with momentum k moves at a speed given by the group velocity v(k) = ε (k). Now z andz are the (respective) analytic continuations of k l , k r to the complex plane, in which case the continuity equations are continued by Wick rotation t = −iy, leading to (89,90). Hence (89,90) can be seen as Wick-rotated continuity equations.
In fact Ref. [56] proceeds in exactly the reverse order as we just did. Eqs. (89,90) are just assumed to hold in imaginary time, as the only sensible analytic continuation of the real time continuity equations. Then, it is possible to find the simplest lagrangian for which (89,90) are the EL equations. The reader can easily check that this is exactly the free energy (71). Here we derived all that starting from the lattice model.

Complex characteristics
Such equations may be "solved" by the method of complex characteristics. Without entering too much into specifics (see e.g. [57]), the solution z satisfies where G is an analytic function. This result is nothing more than a (nice) parametrisation of the whole set of solutions. As with all PDE's, it is very important to specify the boundary conditions; here those turn out to be enough to determine the desired analytic function G. Said differently, to each boundary condition (bc) there corresponds an analytic function. However, finding from a given bc is, in general, a difficult task. At this stage the reader might be worried that we did not achieve much from a practical perspective: we mapped the difficult problem of finding the limit shape onto the difficult problem of determining the analytic function. Fortunately the situation is not that bad, since there are a few interesting examples where G can be guessed. Those include the emptiness formation probability setup [56] (see also exercise 4.1), but the method can also be adapted to a domain wall geometry, as is discussed below. Note also that in this approach, if the density profile along any horizontal or vertical slice is known, then G is known, so everything is known.

Introducing the domain wall geometry
A simple way to generate limit shapes is to look at a 'domain wall geometry' similar to the one we used to generate the Aztec diamond for square lattice dimers. We consider lattice fermions on the infinite line Z, and impose the configurations |ψ = x<0 c † x |0 at the bottom and top boundaries. In height language these boundary conditions correspond to the maximum slopes ∂ x h = π 3 for x < 0, ∂ x h = − 2π 3 for x > 0. The boundary conditions are shown in figure 14 on the left, while a typical configuration is shown on the right. Specific examples with this geometry are worked out in the next subsection. We will see in particular that the corresponding arctic curve is an ellipse.

Examples
Here we discuss two examples with boundary conditions that correspond to a domain wall initial state. We focus on the XX chain first, and then proceed to dimers on the hexagonal lattice (assuming u < 1).
The XX chain in imaginary time (or PNG droplet [58,59]). Here ε(z) = − cos z, which gives ε (z) = sin z. One can check the right solution for the domain wall geometry takes the form R cos z = x + iy sin z.
[This corresponds to G(w) = arccos w R ]. It is sufficient to check that the boundary conditions are correct. Indeed, those are set at y = ±R, and read (for y = R) The rhs is real, so this imposes k F (x > 0) = 0 and k F (x < 0) = π, which is exactly the density corresponding to the domain wall boundary conditions. (93) is a quadratic equation, so can be easily solved. Provided x 2 + y 2 < R 2 it has two simple roots z and −z * , with From this we can compute essentially everything, as we demonstrate now. The density inside the arctic circle is simply given by ρ(x, y) = ∂xh π +1/2 = 1 π Re z = 1 π arccos . The height profile is given by and the corresponding (minimal) free energy is See figure 15 for pictures. The total free energy is S 0 [h cl ] = dxdyF (∂ x h cl (x, y), ∂ y h cl (x, y)) = −R 2 /2, which means the partition function scales as Z(R) ∼ e R 2 /2 . In fact, we will see in section 5 that Z(R) = e R 2 /2 exactly for all R.
Looking at the edge behavior is also worthwhile. Except near x = 0, y = ±R, the density profile vanishes with a square-root singularity. The height field in turn vanishes as x 3/2 , a result which is known more generally as the Pokrovsky-Talapov law [60]. Hilbert transform and dimers on the hexagonal lattice. For a general dispersion relation ε(z) = − p a p cos pz, v(z) = ε (z), it is easy to see that the right generalisation of (93) reads −Rṽ(z) = x + iyv(z) (98) whereṽ is obtained from v(z) = p pa p sin pz, by replacing all sin pz by − cos pz. This transformation is called Hilbert transform. Applying this to the dispersion (55) leads tõ We also get a quadratic equation, whose solution is inside the arctic ellipse X 2 + y 2 < R 2 , where X = 1−u 2 u x + Ru. It is also possible to compute everything as we did before, but we refrain from doing so.

Algebraic geometry
Recall the free energy may be put into the form where and µ, ν are determined from (72,73), namely, as the (double) Legendre transform of the generating function f (µ, ν). Let us take the honeycomb dimer model as an example. There is an alternative (equivalent) way of accessing the generating function, using the Kasteleyn approach [22,42], which leads to the more symmetric expression where P (z, w) = z + w + 1. The reader can check that the two expressions match, using dθ 2π log(a + be iθ ) = log max(|A|, |B|). In a sense, the transfer matrix approach does the integral over q for free in (103), but depending on context, the symmetric form can be nicer, and this is the way free energy is calculated in the mathematical literature, see e.g. [24]. The polynomial P essentially encodes the lattice, more complicated ones leads to higher order polynomials, for example P (z, w) = 1 + z + w − zw for the square lattice. In Kasteleyn's approach, P is related to the determinant of the Kasteleyn matrix.
In algebraic geometry context, f (µ, ν) as defined is (103) is called the Ronkin function of the polynomial P (z, w) = z + w + 1, and the free energy F (r, s) is then its Legendre dual. The dual takes values inside a polygon of allowed slopes, which is the Newton polygon of the polynomial P . One can also define the Amoeba of the algebraic curve P (z, w) = 0, as the set namely, the projection of the algebraic curve (in C 2 ) P (z, w) = 0 to R 2 by the map (z, w) → (log |z|, log |w|). Now, from general algebraic geometry machinery, the Ronkin function is convex, and linear in the complement of the Amoeba. The complement is an union of disjoint simply connected pieces, which correspond to frozen regions [22]. Under the Legendre duality each component of the complement is mapped to a single point of the Newton polygon. One can also show that the Hessian of the Ronkin function is constant det Hess(R) = π 2 for any point inside the Amoeba. This is interpreted as a Monge-Ampère equation. By Legendre duality this implies Hess(F ) = 1/π 2 , so K = 1 in the fluctuating region. This result happens to be a characterisation of algebraic curves known as Harnack curves, so the algebraic curves one gets from dimer models are always Harnack. The interested reader may have a look at references [57,61,62] for a much more precise discussion.
Hence, the statement "The free energy of the dimer model is the Ronkin function of a Harnack curve, so satisfies a Monge-Ampère relation for any point in the Amoeba" translates for us into the statement "Dimers map to free fermions, so bosonising the fermions we get K = 1 in the fluctuating region". A nice illustration of the two different types of jargons.

Edge behavior
The hydrodynamic solution gave z = z(x, y) = k F +iν from which the limit shape follows. Our aim is to provide a heuristic derivation of the universal edge behavior near the arctic curve. To do that, we look at the particular case where ν (or the current) is zero, which occurs at least in the middle at y = 0 in all the examples we discussed. Let us emphasize that the argument we present here may be generalized to ν = 0. However, the discussion would involve left/right ground states of non-Hermitian Hamiltonians, which would obscure the argument slightly. We now assume that z(x, y) = k F (x, y) is known from the hydrodynamic solution and happens to be real. The first claim is that the correlations around a given point (x, y) are those of the ground state of the Hamiltonian where µ is set such that ε(k F ) = µ, to ensure that k F be the Fermi momentum. Close to the arctic curve the density goes to 0 or 1, let us focus on the case of vanishing density. This means k F goes to zero, and it is possible to expand the dispersion around k = 0. Since ε(−k) = ε(k), we have ε(k) = ε(0) + 1 2 (0)k 2 + O(k 4 ), and we get the effective edge Hamiltonian Now how does k F depend on x, y close to the edge? The edge corresponds to the case where two (or more) solutions z s and −z * s to the hydrodynamic equation coalesce, so that we get a double root or higher. Generically this root will be only double, meaning ε (0) = 0. Writing x a (y, R) for the arctic curve, and setting x = x a −x, we get that, generically for some coefficient α that depends on x a /R. (for example for the domain wall XX chain α = √ 2). Therefore, k 2 F behaves linearly inx, close to the edge. Now what are free fermions with k 2 dispersion? This is just a free Fermi gas in the continuum, sometimes encountered via its relation to the Tonks-Girardeau gas in cold atom systems [63]. Making crudely the substitution k → −i d dx and undoing the Fourier transform yields the Hamiltonian that is free Dirac fermions {c(x), c † (x )} = δ(x −x ), in a linear potential. The change of variablesx = (R/α) 1/3 u leads to Therefore, a non trivial Hamiltonian emerges at distances ∼ R 1/3 from the edge. We expect correlations on such distance to be that of the ground state of (110). This is still a free problem, so diagonalizing H edge boils down to solving the single-particle eigenvalue equation u). The solutions are well-known to be Airy functions. It is an exercise to show that the (Dirac sea-like) ground state propagator for this Hamiltonian is which is known as the Airy kernel.
Let us now look at a region A = [s, ∞) of the infinite line. It is another nice exercise to compute the generating function Υ(α) = e α ∞ s c † (u)c(u)du , and show that it is given by the infinite series Now define the emptiness formation probability E(s) = lim α→−∞ Υ(α, s), which is, as its name indicates, the probability that the interval A = [s, ∞) be empty of particles. We have from the previous formula the exact series representation One can show that E(s) is smooth, positive, strictly increasing, and lim s→−∞ E(s) = 0, while lim s→∞ E(s) = 1. E(s) gives us information about the last (or rightmost) particle. Indeed E(s + ds) − E(s) is proportional to the probability that the rightmost particle lies in the interval [s, s + ds]. Hence p(s) = dE ds is actually the probability density function for the rightmost particle, so E(s) may be interpreted as the cumulative distribution for the rightmost particle.
The distribution that has the rhs of (113) as a cumulative distribution has a name. It is called the Tracy-Widom distribution [64]. We finish with a few remarks: • Exact series such as (113) or the one above are called Fredholm determinants. They are the continuum analog of the regular determinant, for operators. See e.g. (176,177) for a discrete analog.
• Proving that the distribution of the rightmost fermion does converge, after proper rescaling, to the T-W distribution requires of course more work than the heuristic argument we just gave. However, it still illustrates the physical mechanism through which T-W emerges, that is free fermions in a linear potential. We will see in the next section another mechanism through which T-W occurs.
• Convergence to T-W has been proved in all the models we considered so far. For the XX chain this was done by Praehoffer and Spohn [59], while dimers have been treated by Johannson (honeycomb [18] and square [65]).
• The attentive reader might have spotted a physical flaw with the argument we just made. We have implicitly assumed separation of scales throughout. This is fine in the bulk, but it is not clear that it still holds at the edge. In fact, one can show that the edge is exactly borderline with respect to separation of scale, since the density varies on distances that are comparable to inter particle distances. In that sense T-W is smoothly connected to the bulk, where separation of scale does hold. This also justifies exact calculations starting from the lattice, see section 5.
• The word generic near (108) is important. It is necessary for the density to vanish as square root to get a k 2 dispersion and a linear potential for the fermions, so R 1/3 behavior and T-W. If the arctic curve has cusps for example, edge correlations near (generic) cusps are described by a higher order kernel known as Pearcey kernel. We refer to [66] for a discussion of all these kernels. We discuss another example of higher order kernel in exercise 4.3.
• A plot of the T-W probability density function is shown in figure 17. As can be observed it resembles a gaussian, but it is slightly skewed.
• T-W scaling is part of a broader subject, which goes under the name Kardar-Parisi-Zhang (KPZ) equation [67] and KPZ universality class. Roughly speaking, the KPZ equation is a stochastic PDE that models interface growth [68,69], and it turns out the long time limit of this equation with certain initial conditions is exactly the T-W distribution. We refer to [25,[70][71][72][73] for reviews on this equation, the KPZ universality class, and related topics in random matrix theory. Emptiness boundary conditions [56] We consider fermions with dispersion ε(z) = − cos z. We have seen in the text, that solutions to complex Burgers may be put under the form cos z = G(x + iy sin z), where G is analytic. We wish to find the limit shape corresponding to emptiness boundary conditions, that is (i) vanishing density on an interval x ∈ [− , ], y = 0 of the full complex plane and (ii) density 1/2 at infinity.  Another higher order edge kernel [74,75] We consider the Hamiltonian where, semiclassically, µ depends position as µ = µ(x) = x/R for large R. The dispersion relation is ε(k) = − cos k + 1 4 cos(2k).
1. Where is the location of the right edge? 2. On which scales to you expect the distribution of the righmost fermion to be? 3. Compute the associated edge kernel and edge distribution.

Exact lattice calculations
The approach we have taken so far was variational or hydrodynamic: we showed how computing the limit shape boils down to solving PDEs, and found a few cases where this could be done explicitly. In turns out those precise cases can often be treated with other more direct methods. That is, computing explictely the two point function, and then recover our previous results by a careful asymptotic analysis. This approach is often more technical, but nicely complements hydrodynamics, by confirming its prediction, and also providing more information.
Our aim is to give a flavor how this can be done, on one of the simplest examples. We refer to [51,76,77] for similar calculations. We take the transfer matrices from before, and we try to compute the two point function for fermions in the domain wall geometry, where the top and bottom boundary are domain wall-like, with all particles packed to the left of the origin, set at x = 0. This is the precise geometry in which the hydrodynamic problem was solved in terms of Hilbert transform (section 4).
Of course as always in a free fermion calculation, if we know the propagator then Wick's theorem allows to reconstruct all higher order correlations. However, as we shall see, computing the two point function is more difficult than in standard condensed matter theory setups.
Let us now specify some notations. We take lattice sites to be half-integers, that is x ∈ Z+1/2. We focus on the special case where operators are measured at the same imaginary time y, but generalisation is straightforward. We wish to evaluate where expectation values are taken in the domain wall state |ψ = x<0 c † x |0 . We also have H = dk 2π ε(k)c † (k)c(k), and recall c † x = π −π dk 2π e −ikx c † (k), c † (k) = j∈Z+1/2 e ikj c † j . Using e (k)c † (k)c(k) c † (k)e − (k)c † (k)c(k) = e (k) c † (k), this may be rewritten as with This remaining term is unusual, and illustrates the extra layer of complexity associated to imaginary time problems -for a real time calculation, just set y = it and R = 0. It is very difficult to evaluate (116) in general; however, the special form of the domain wall state in which averages are taken allows for a small miracle.

A nice bosonization trick
Let us work out the case of nearest neighbor hoppings, for which ε(k) = − cos k. The main player in the calculation will be the operator Obviously, H = −(b + b † )/2. Think of a finite-size regularisation of the chain, e. g. with sites from −l to l. The commutator of b with b † is given by a telescopic sum, which simplifies to which means ψ|[b, b † ]|ψ = 1. Now expand e −2RH in power series. It is easy to check that ψ|[b, b † ]H p |ψ = ψ|H p |ψ provided l > p. Hence for any term in the power series, we can always choose l sufficiently large such that the commutator is scalar. For any finite R the series is expected to converge quite fast, which means we are allowed to assume [b, b † ] = 1 throughout. Hence the operator b is, effectively, a boson. It also happens to annihilate the domain wall state, b |ψ = 0. Now, recall the following formula which is a special case of the Baker-Campbell-Hausdorff identity 9 . Using this formula both in the numerator and denominator with α = R combined with e αb |ψ = |ψ , ψ| e αb † = ψ|, yields The last trick is to take derivative. Computing the commutator we obtain Integrating back we finally obtain

General dispersion relation
The case of general dispersion relation ε(k) = − n≥1 h n cos(nk) can be handled in a similar fashion. One introduces the set of operators which effectively satisfy the commutation relations Using this one can show in a similar fashion and The propagator finally reads whereε(k) is the Hilbert transform of ε(k). Let us insist again that this only holds for expectation values in the domain wall state. All what is left is to compute G 0 (k, q). We get from a direct calculation

Asymptotic analysis
The general method to evaluate the double integral in the limit R → ∞, x/R, x /R, y/R fixed is the stationary phase, or steepest descent method. The argument inside the exponential can have very large real and imaginary parts. Writing one expects the integral to be dominated, after proper contour deformation, by the region close to the points k c (resp. q c ) where the phase θ(k) (resp. −θ(q)) becomes stationary. The stationary points are the solution of the equation whose solution we denote by z and −z * . This equation is, in fact, identical to (98), which we obtained from the hydrodynamic approach. A full asymptotic analysis falls outside the scope of these lectures. However, let us just mention that what matters is the taylor expansion of the phase around the saddle points, that is the expansion Essentially computing the asymptotics boils down to computing a gaussian integral. The case of coinciding points x = x is more tricky, since the k and q saddle points might coincide. In that case one has to take into account the pole at k = q. See e.g. [78] for the details. Let us briefly comment on the edge behavior. The arctic curve corresponds to the points where the two solutions (assuming there are two) z and −z * become equal. This means the second derivative θ (z) vanishes, and it becomes necessary to expand to third order This naturally leads to the Airy kernel, since Airy functions may be alternatively defined as Ai(x) = R+i0 + dk 2π e ikx+ik 3 /3 . The subject would require longer exposition, but we have illustrated the two equivalent through which the Airy kernel emerges. Either with fermions in a potential that becomes linear in a Hamiltonian point of view, or through coalescence of two saddle points.

Exercise 5.1
Bosonizing Toeplitz determinants [76,79,80] A semi-infinite Toeplitz matrix is a matrix T = (T ij ) i,j∈N whose elements depend only on the difference i − j, T ij = g i−j . It is convenient to interpret the g l as Fourier coefficients of a periodic function, sometimes called symbol: We assume that g is sufficiently smooth, has a well-defined logarithm which we denote by ε(k) = log g(k), and also that π −π dk 2π ε(k) = 0. Consider the free fermions where the ε p are the Fourier coefficients of ε(k). We introduce a similar domain wall state |φ = ∞ j=0 c † j |0 as in the text, where notice the fermions are now located on nonnegative integers.

Show using bosonization that φ|e
3. Show using bosonization that where the g ± (k) = exp ±n>0 ε n e ikn are the Wiener-Hopf factors of g(k).
We now look at a 2N × 2N truncation of T , which we note T N . We want to evaluate det T N = det 0≤i,j≤2N −1 (g i−j ) when N → ∞. For this purpose, we introduce the state |ψ N = |x|<N c † x |0 where sites are now put on the half-integer line x ∈ Z + 1/2.
provided the series inside the exponential converges sufficiently fast. This result was first found by Onsager-Kaufman [79], and proved by Szegö [80] shortly thereafter (both used different techniques). 7. You are Onsager and Kaufman, and you just realized that the spin-spin correlations σ 0,0 σ n,n along the diagonal of the classical isotropic Ising model on Z 2 are given by a n × n Toeplitz determinant with symbol g(k) 2 = 1−αe −ik 1−αe ik , where 1/α = sinh 2 (βJ). This holds below the critical temperature, Jβ > Jβ c = arcsinh 1. What is the magnetisation exponent of the 2d Ising model? 46

Conclusion and related problems
We finish with a discussion of a few selected topics that go beyond the lectures, but still fit well with the spirit of the notes. We first examine the effects of interactions (section 6.1), how this (does not) affect the edge behavior in section 6.2, and finally explore the intricacies of the Wick rotation, section 6.3.

Interactions
For interacting systems, i.e. systems that cannot be mapped onto free fermions, the logic of the present notes still applies. The difference is that no exact formula for the free energy exists in general anymore. There are however deformations of the dimer model for which some analytical progress is possible. Those models are called integrable.
A discussion of all the intricacies of integrable models falls well outside the scope of the present notes, see [53,54,81] for reviews. We consider the case of the six vertex model, or even plaquettes interacting dimers, as explained in section 1. We parametrize the interaction term as e λ = 1 − ∆ or e λ = 1 − cos γ, depending on convenience. The only result that we need here is the following fact: in the same way that the free energy at the free fermion point could be determined from a simple ground state energy with some current, where ν and k F are determined from r, s, a similar expressions holds true for half-interacting dimers (or the six vertex model). Namely, the free energy is determined from the biggest eigenvalue of the transfer matrix (with appropriate particle number and current). The latter is given by [40,82] The λ j play a role similar to momenta for free fermions, and N/L controls r. The big difference is that the quantization condition is much more complicated. It is given by a set of equations called Bethe equations 10 . The free case corresponds to γ = π/2 for which the rhs simplifies to 1, and the λ k can be obtained explicitely. The reader can easily imagine that solving the Bethe equations is in general extremely difficult. The fact that (away from the free fermion point) the rhs is a complicated product over the positions of the particles has an important physical consequence, which usually goes under the name dressing: changing the number of particles (N ) affects all the rapidities, as illustrated in figure 18.
Dressing severely complicates asymptotic analysis, since any j f (λ j ) becomes in the thermodynamic limit f (λ)ρ(λ)dλ, where ρ is the density of Bethe roots. Now comes the big problem: the root density is only known in the case of zero current ν = 0, in which case it is a solution to a linear integral equation over a segment of the real line [53]. In the case of nonzero current these aspects have been investigated numerically in Ref. [40], and one finds that the Bethe roots condense on some non trivial curve in the complex plane.
A particular case that can be treated exactly is the so-called five vertex model, obtained from the six vertex model by setting one of the vertex weights to zero -it can also be seen as an (half-plaquettes) interacting version of the honeycomb dimer model [83]. This model is related to stochastic processes such as TASEP. The exact free energy can be computed, and limit shapes are also parametrized by analytic functions [55]. This is to date one of the most complicated model in which the variational/hydrodynamic program has been applied. This model is also the only one for which one can show that the Luttinger parameter K is not constant.
For the full six vertex model the hydrodynamic program has not been completed, the main bottleneck being determining the exact curve on which the roots densify. However, the arctic curve has been determined analytically, using the lattice approach. In a series of papers [84,85], Colomo and Pronko, and Colomo-Pronko-Zinn-Justin managed to compute exactly the lattice emptiness formation probability, which gives the distribution of the last particle. They managed to determine the precise location where this probability goes to zero in the thermodynamic limit. This location coincides with the arctic curve. Except at special values where γ/π is a rational number, this curve is non algebraic. Later, an attractive tangent method was also introduced to get the arctic curve in a slightly simpler way [86]. This method was also recently used to provide a proof [87] of the Colomo-Pronko formula for the curve in the special case ∆ = 1/2, the so-called combinatorial point.

Tracy-Widom at the edge
Tracy-Widom scaling is, in fact, also expected at the edge, for the following simple physical reason: near the edge the particle (or hole) density goes to zero, hence particles (holes) are diluted. For local interaction such as the plaquette terms we discussed previously, it is reasonable to assume that interactions become weaker and weaker. Hence particles become effectively free near the edge, and we expect the arguments presented in section 4.6 to hold, and T-W scaling at the edge. One can check numerically that this is the case: in 19 we show Monte Carlo simulations in interacting dimers and the six vertex model. We compute the lattice emptiness formation probability numerically, rescale appropriately, and compare to the T-W distribution. As can be seen the agreement is quite good. We also compute the skewness of the discrete distribution, and compare it to T-W, with excellent (and improving for larger L) agreement. Other similar checks in inhomogeneous quantum chains can also be found in [75], or with anharmonic chains in thermal equilibrium [88]. Showing this analytically in some generality in integrable models is not easy, despite the fact that the lattice emptiness formation probability is known exactly in the six vertex model. Note that the argument above does not assume integrability. However, the dilution argument can nicely be illustrated in Bethe-Ansatz integrable models, such as six vertex. Indeed looking back at the Bethe equations (137), the edge diluted limit corresponds to N L, for which the rhs can be considered a constant. Hence in this limit dressing disappears, and we are back to free fermions.
On the rigorous side, stochastics processes such are ASEP appear to be more tractable. For example, a proof of T-W scaling is known for ASEP with step initial conditions [89].

Wick rotation and inhomogeneous quantum quenches
We have seen that the expectation value of an observable O x in imaginary time is given by with either R and y discrete (dimers models, etc) or R and y continuous (XX chain in imaginary time). Our starting point will be the later case. The reader interested in quantum models might have spotted a similarity between the previous expression and regular time evolution steming from the Schrödinger equation. Indeed, for a quantum system prepared in a state |ψ , and let evolve in time with the Hamiltonian H, the wave function at time t is |ψ(t) = e −iHt , and the expectation value of O at time t becomes Formally the real time evolution may be recovered from setting y = −it in (138) and taking the limit R → 0 + . This procedure is the famous Wick rotation. This observation can be useful for two reasons.
a) First, exact calculation in the spirit of section 5 or using integrability techniques are quite algebraic in nature. For this reason they are often valid for any value of R, y ∈ C, which means the Wick rotation is perfectly justified for any finite time t. We can use this to derive exact highly nontriavial expressions for out of equilibrium quantities in a few selected case.
For example, the partition function Z(R) = ψ|e −RH |ψ of the six vertex model with domain wall boundary conditions is known exactly from the work of Korepin [90] and Izergin [91] (see also [92]). We can then take the Hamiltonian limit, perform the Wick rotation, and get an exact expression [93] for the amplitude ψ|e itH XXZ |ψ , where H XXZ is the Hamiltonian of the XXZ spin chain, an integrable generalisation of the XX chain. The modulus square of the amplitude is called return probability, or Loschmidt echo. This exact result is difficult to get from more direct approaches [94], which makes this method worthwhile.
b) At a more speculative level, it is tempting to try and Wick-rotate results already in the thermodynamic limit, following [95]. Of course, there is no mathematical justification for this. The statement that this can lead to wrong results is called Stokes phenomenon. Let us go back to our quench from a domain wall state. Hydrodynamics ideas can also be applied to out of equilibrium quantum integrable sytems, as was established recently [96,97]. This subject goes under the name generalized hydrodynamics. For a quench from a domain wall state a small miracle occurs, and it is possible to get the density profile exactly [98]. What is nice is that the Wick rotation of the arctic curve gives back precisely the location of the front, the simplest example being the free fermion point arctic circle x 2 + y 2 = R 2 which becomes x = ±t after Wick rotation.
This does not work for all observables, though. For example the return probability or the density profile have the crazy property of being nowhere continuous as a function of ∆ in the thermodynamic limit, which is clearly not the case in the original statistical problem.
A Self-contained reminder on free fermions techniques A.1 An explicit construction of lattice fermions Two level system. The two level system is obviously the most important genuine quantum system ever, so let us start from this one. We consider the Hilbert space H C 2 , and take the most down to earth approach, which is to work with two by two matrices (dim H = 2). We choose the two basis vectors |0 = 0 1 and |1 = 1 0 . |1 is interpreted as the presence of a particle, |0 as the absence of a particle (vacuum state). Pure states (or wave functions) are of the form Now let us introduce the two main heroes, and c † = 0 1 0 0 .
(143) † denotes hermitian conjugate. We use bra/ket notations, e. g. the bra 0| = (|0 ) † = (0 1) is a line vector. We have c † |0 = |1 , c |1 = |0 . Since c |1 = |0 , c destroys the particle, so we call it the annihilation operator. c † |0 = |1 creates a particle from the vacuum, so c † is the creation operator. We also have c † |1 = c |0 = 0 0 = 0, so that it is not possible to create two particles, or destroy a non existent particle. The zero vector is not to be confused with the vacuum. The rightmost hand side of the previous equation involves a slight abuse of notations; in the following we keep on writing 0 any vector/matrix that has all elements equal to zero.
Note also c † c = 1 0 0 0 , and c † c + cc † = I 2 = 1 0 0 1 , which will be useful in the following. Obviously, any matrix in M 2 (C) can be written as a linear combination of c, c † , c † c and cc † .
A collection of L two level systems. We now consider the Hilbert space H (C 2 ) ⊗L , where L is an integer ≥ 2. H has dimension dim H = 2 L . We want a set of (Dirac) fermionic operators, that is, a set of 2 L × 2 L matrices c i , c † i for i = 1, . . . , L that satisfy where I = I 2 ⊗ . . . ⊗ I 2 is the identity operator. ⊗ denotes tensor product, given by . . .
The chain of k − 1 tensor products of (−1) c † c in (149) is called a Jordan-Wigner string. One can readily check that the CAR (144), (145) are satisfied, using (147) lots of times. Of course, in practice, one often just uses the CAR without caring about an explicit representation. However, in the context of classical statistical mechanics, the above explicit construction turns out to be very useful.
• Exponentiation [follows from the fact that c † i c i is idempotent, (c † i c i ) 2 = c † i c i ]: exp τ c † i c i = 1 + (e τ − 1)c † i c i .
• Time evolution [take derivative]: A.3 How to diagonalize a free fermions Hamiltonian?
A free (lattice) fermions Hamiltonian is a 2 L × 2 L matrix that is quadratic in the fermions creation and annihilation operators (we assume hermiticity here, which is not necessary, strictly speaking): where A and B are L × L matrices (A is a Hermitian). This is of course a specific class of Hamiltonians, since we have at most 2L 2 free parameters while the Hilbert space size is 2 L . In this context, free really means quadratic in the fermion operators. In the following we explain how to diagonalise H in the special case B = 0, for simplicity. The procedure described below can be generalized to treat cases where B is a non zero matrix. Hamiltonians of the form conserve the number of particles, which means applying it on n-particle states c † i 1 . . . c † in |0 returns a sum over n particle states (any fermion destroyed by c j is immediately created back by c † i ). A is a hermitian L × L matrix, so can be diagonalized in an orthonormal basis. The corresponding eigenvalue equations read (assume no multiplicities for simplicity) L j=1 A ij u jk = k u ik , k = 1, . . . , L.
The eigenvalues are the k and the u jk are orthonormal, meaning L j=1 u * jk u jq = δ kq . Now introduce a new set of fermions as Then it is easy to show {f k , f † q } = δ kq and {f k , f q } = 0, so the new set of operators also obey the CAR. In terms of these the Hamiltonian reads Obtaining the spectrum becomes quite easy now. Obviously H |0 = 0. Using the anticommutation relations, Hf † k |0 = k f † k |0 . By induction, we obtain To get a nonzero eigenvector, the k i have to be pairwise distincts. Also, any permutation of the k i s gives back the same eigenvector up to a sign. Hence the spectrum is There are L n linearly independent eigenvectors in the sector with n particles, so the total number of eigenvalues is L n=0 L n = 2 L = dim H, as should be. An eigenstate with smallest eigenvalue is obtained by choosing all single particle energies k that are negative, |Ω = k, k <0 f † k |0 (irrespective of the order in which the product is taken).

A.4 More elaborate properties
Here is a collection of results that are very useful in practice. Almost all free fermions calculations make use of several of those at some point. We start with the most famous one. a) Wick's theorem [100]. Let the f j be linear combinations of the c i , c † i for j ∈ {1, . . . , 2n}, and H be a free fermion Hamiltonian. Then the thermal average T where Pf is the Pfaffian. For an antisymmetric matrix A = (A ij ) 1≤i,j≤2n , it is defined as Pf A = 1 2 n n! σ∈S 2n (−1) σ A σ(1),σ (2) . . . A σ(2n−1)σ(2n) , where the sum runs over all permutations σ of {1, . . . , 2n}, and (−1) σ is the signature of the permutation. The Pfaffian can be shown to be a square root of the determinant.
The factor 1/(2 n n!) may also be removed by imposing the ordering σ(2i) < σ(2i + 1), and σ(2i) < σ(2j) for j > i. As an example, the theorem yields We refer to [101] for a proof of (163). Note that as a particular case, expectation values in any ground state may be obtained by taking the limit β → ∞.
b) Wick's theorem (no pairings). We consider the special case where H is of the form (157), while for j = 1, . . . , n the f j are linear combinations of the c † i only, and the f j+n = g j are linear combinations of the c j only. Then We will need only this particular case in the limit β → ∞ in these notes.
c) Trace of exponential.
Proof : diagonalize the quadratic form in the exponential. d) General time evolution.
Proof : diagonalize the quadratic form in the exponential.
e) Product of exponentials.
Proof : use the Baker-Campbell-Hausdorff formula.
f ) Average of an exponential. Let P = (P ij ) 1≤i,j≤L be a L × L matrix. Then is any state where Wick's theorem applies we have e i,j P ij c † i c j = det(I + (e P − I)C).
C is the L × L matrix with elements C ij = c † i c j , I the L × L identity. Of course, it is also possible to combine (170) with (169) to compute the average of a product of exponentials.
The formula gives the full counting statistics as byproduct: for any subset A of {1, 2, . . . , L}: Proof : just pick P ij = λδ ij δ i∈A in (170) and check that only the block i, j ∈ A contributes to the determinant.
Proof of (170): Introduce a new set of fermions d k , d † k that diagonalise the quadratic form in the exponential on the lhs, as in section A.3. Denote by k the corresponding eigenvalues of P . The d k , d † k also satisfy the CAR. Then We have used in succession (153) Finally, multiplying by U on the left and U † on the right does not change the determinant, and gives (170).

A.5 Infinite lattice or continuum limit
So far we have introduced lattice fermions as 2 L × 2 L matrices, where L is the number of lattice sites. It is however very useful to consider generalisations where the matrices become infinite or operators in the continuum. The case of the infinite lattice, e. g. Z can easily be dealt with by considering a finite lattice j ∈ {−l, −l + 1, . . . , l} computing observables O l and then taking l → ∞.
We often make use of the following notation for the Fourier transform on the infinite lattice, where k ∈ [−π, π], in the main text. The obey the anticommutation relations {c(k), c † (k )} = 2πδ(k − k ). It is of course also possible to consider real space fermions, in which case we use the notation c † (x), which obeys the anticommutation relations {c(x), c † (x )} = δ(x − x ). In the main text, k and q always refers to momentum while x and y always refer to position, so it is not possible to get them confused.