Quantum Crystals and Spin Chains

In this note, we discuss the quantum version of the melting crystal corner in one, two, and three dimensions, generalizing the treatment for the quantum dimer model. Using a mapping to spin chains we find that the two--dimensional case (growth of random partitions) is integrable and leads directly to the Hamiltonian of the Heisenberg XXZ ferromagnet. The three--dimensional case of the melting crystal corner is described in terms of a system of coupled XXZ spin chains. We give a conjecture for its mass gap and analyze the system numerically.


Introduction
In this note, we study the quantization of the statistical system of a melting crystal in one, two, and three dimensions. In two and three dimensions, the melting crystal configurations correspond to the diagrams of random partitions, resp. plane partitions. Random partitions appear in many branches of mathematics and physics, such as Gromov-Witten theory, random matrix theory, Seiberg-Witten theory and others, see e.g. [1] for an overview. Also in string theory, random partitions are of importance. In [2,3], the partition function of the melting crystal corner given by the MacMahon function was shown to equal the partition function of the topological string A-model with target space X = C 3 .
In the spirit of stochastic quantization, we construct a quantum Hamiltonian whose ground state corresponds to the classical steady state distribution. The treatment is reminiscent of the work of Rokhsar and Kivelson on quantum dimer models [4] (the perfect matchings of the dimer model on the hexagonal lattice being in one-to-one correspondence with the configurations of the 3d crystal) and Henley's study of the classical and quantum dynamics for systems satisfying the detailed balance condition [5]. While the quantum dimer model only keeps track of the total number of configurations, the number of boxes of each crystal melting configuration is taken into account in our system. In the partition function, each configuration is weighted by q #boxes . The quantum dimer model at the Rokhsar-Kivelson point corresponds therefore to the limit q → 1 of the melting quantum crystal. The problem of the two-dimensional quantum crystal can be reformulated in terms of a one dimensional fermionic system. The system translates to an asymmetric, hard-core diffusion problem on a half-full infinite chain. In fact, the quantum Hamiltonian we derive is easily transformed into the Hamiltonian of the XXZ ferromagnetic Heisenberg spin chain.
The ground state of the partition problem corresponds to a particular sector (the half-filled) of the kink ground state of the infinite chain. Using the technology of basic hypergeometric series, we derive the limit shape [6], which corresponds to the integral of the magnetization in the limit q → 1. Also the general expression for n-point correlation functions is derived. Based on our findings in the two-dimensional case, we can express the 3d system as a system of interlaced spin chains. The quantum crystal corner is a rare example of a threedimensional system which is nonetheless fairly tractable. This is of course due to its very constrained nature. Expressed in terms of spin chains, it becomes manifest that its dynamics is not even fully two-dimensional.
Based on numerical evidence, we conjecture the mass gap of the 3d quantum crystal to be the same as in the 1d and 2d systems. More generally, all these systems can be interpreted as directed random walks on the graph of configurations, where going forward and going back has different weights. Depending on the dimensionality of the problem, the graph of configurations changes, but it always remains a directed acyclic graph. This more general framework unifies the three problems considered here and might also account for common properties, such as the value of the mass gap.
The plan of this note is as follows. In Section 2, we explain the general procedure of deriving the quantum Hamiltonian for a melting crystal. As a warm-up we specialize in Sec. 3  while Sec. 4.4 gives some numeric results for higher eigenstates. In Section 5, the threedimensional system is discussed and numeric results for the higher eigenstates of the 3d case are presented. The conjecture for the mass gap is stated in Sec. 5.3.
In Section 6, we end with concluding remarks and an outlook on the open problems.

The quantum Hamiltonian
In this section, we will explain the general framework of quantization we will be using. The quantum dimer model of Rokhsar and Kivelson [4] is an example thereof. In fact, it is a limit of the more general problem of quantum crystal melting.

The general framework
Consider a classical system with a discrete set of configurations { C i }, such as a lattice model. LetH be the energy functional that associates to each configuration its energyH(C i ). If we consider the canonical ensemble, the probability of finding the configuration C i is given by its Boltzmann weight where β is the inverse temperature and Z is the partition function We want to quantize this system. To do so we introduce the separable Hilbert space H spanned by the orthonormal vectors |C i which are in one-to-one correspondence with the classical configurations. Each vector in |ψ ∈ H defines a probability measure on the C i as follows: where the scalar product is understood in l 2 .
Consider the state |Ψ 0 defined by It is immediate to verify that the measure associated to |Ψ 0 is precisely the Boltzmann distribution for the canonical ensemble, up to the normalization fixed by To be more explicit, consider a physical quantity Q that can be measured on the classical On the canonical ensemble, the statistical average of Q is given by To such a quantity we can associate an operatorQ in the Hilbert space H which is diagonal in the basis given by the |C i , For any such Q, the statistical average thus coincides with the expectation value ofQ on the state |Ψ 0 (which is equivalent to the above statement of the measure corresponding to the Boltzmann distribution): Thus the choice of normalization in Eq. (2.5) is meaningful.
Let now H be a Hamiltonian operator in the Hilbert space H that admits |Ψ 0 as its unique zero energy ground state, We call H the quantum Hamiltonian corresponding toH. Note in particular that β is now just a parameter in H. The quantum partition function which describes the quantum behaviour of the system contains on the other hand a new coupling constanth and is defined as For a given operatorQ the expectation value is defined as the following path integral: The classical limit is obtained forh → 0 and using a saddle point approximation one verifies that the main contribution to the expectation value comes from the ground state, hence reproducing the classical result: A geometric interpretation for the operator H can be found in terms of the classical configuration space. Consider the graph obtained by associating each of the C i to a node and linking two nodes with an edge weighted by l ij = exp −β H (C i ) −H(C j ) . Define the operator q on this graph as where q = e −β . The notation k i denotes the sum running over all configurations C k adjacent to C i . One can easily verify that Since this is the property we were looking for in the quantum Hamiltonian, we take q = H. There is some arbitrariness in the choice of q ; in fact one could have added longer range interactions and still preserved the annihilation of |Ψ 0 . The q operator can be interpreted either as a q-deformed Laplacian on the state graph or as the standard Laplacian plus a potential term: where J is the quantum energy scale (not to be confused with the classical coupling β of the classical HamiltonianH), and f (u) is a function defined on the graph. The choice of the minus sign in front of the Hamiltonian is needed to ensure that |Ψ 0 is the state of minimal energy. Note that in the β → 0 limit the potential term V q is proportional to the second variation of the classical Hamiltonian.
By following the above quantization description, the dynamics and a time direction are added to the initial classical problem, since the states in the Hilbert space satisfy the Schrödinger equation This dynamics corresponds to a random walk in the space of configurations where we introduce a (discrete) metric depending on the classical energy.
When this article was close to completion, we became aware of [7] which has some overlap with the treatment given in this section.

The quantum dimer model and quantum crystal melting
The quantum version of the dimer model was introduced in 1988 [4]. Given a dimer model on a bipartite graph, the perfect matchings form a basis for the Hilbert space of the quantum dimer model (QDM). Its (quantum) Hamiltonian takes the form where | and | represent the two possible dimer configurations in a plaquette. The first term is kinetic and runs over all pairs of perfect matchings |α β| which only differ by one elementary dimer flip move (plaquette move), in which a matching around a This construction can be generalized, as shown in [5]. The dynamics of any discrete stochastic classical model is described by the master equation where P α (τ) is the probability to be in configuration α at time τ, and W αβ is the transition rate to state α if the system is in state β. It was shown in [5] (see also [7]) that every discrete system fulfilling the detailed balance condition has a RK point with the same mapping of the eigenfunctions to the classical dynamics.
The system of a melting crystal corner can be mapped to the problem of stacking cubes in an empty corner of 3D space. The growth rules are the following. A cube can be added if three of its sides will touch either the wall or the faces of other cubes. This leads to a minimum energy configuration without any free-standing cubes. The partition function of (a) The empty room perfect matching (b) After one plaquette move (c) After five plaquette moves Figure 1: Empty room perfect matching and matchings after flip moves the melting crystal corner takes the following form: 20) where N(α) is the number of boxes in the configuration α and the rightmost expression is the so-called MacMahon function [8]. The convergence of the above expression is only guaranteed for 0 < q < 1, which we will assume henceforth. The configurations of the melting crystal corner are in one-to-one correspondence with the perfect matchings on an infinite hexagonal lattice. There is a unique (up to translation) empty room perfect matching which serves as the starting point. A plaquette move corresponds to adding or removing a box [2], see Figure 1. Nonetheless, the QDM explained above does not correspond directly to the quantum system for the melting crystal corner, since the number of boxes does not enter anywhere. In the following, we derive the quantum Hamiltonian for the crystal melting. Even though the following paragraph phrases everything in the language of the three-dimensional crystal melting problem, the procedure is completely general and applicable to any dimension.
To describe the (classical) dynamics of the crystal melting process, we express it as a random walk on the state graph (where each node corresponds to a configuration and two nodes are joined by a directed line if the endpoint can be reached from the initial one by adding a box -see Fig. 2). The dynamics is encoded in a master equation of the type (2.18).
We choose the following transition rates: In other words, where C αβ is the adjacency matrix for the state graph. One can verify that the (unique) stationary distribution for P α is where, as in Eq. (2.20), This is equivalent to saying that the system is in a Boltzmann distribution with classical In fact, one can easily verify not only that P (0) α is an equilibrium state, but that it also satisfies the detailed balance condition (2.19). It is customary to define the exit rate from state α as where deg ∓ (α) is the number of configurations connected to α containing one block more (less) than α. Note that in this way, the evolution can be described by the vector equation Instead of using the matrix W, one can define the symmetrized version one can verify that the ground state is now represented by the vectorφ Let us now switch to the | -notation used in the QDM. The ket | here represents a cube in the configuration that can be removed, while | represents a place where a cube can be added. Now, W αα can be written as where the sum runs over all places where a cube can be added or removed in the configuration. Similarly, the non-diagonal terms can be written as It is therefore natural to define the quantum Hamiltonian by acting on the Hilbert space generated by the orthonormal basis of the configurations |α . Notice that at the RK point V/J = 1, this coincides with W: Being proportional, H and W share the same eigenvectors with proportional eigenvalues: In particular, the ground state E 0 = 0 corresponds tõ which has to be interpreted now in the sense of a quantum superposition: where |α is the state corresponding to the perfect matching α. Note that we changed the normalization for the wave function as in Eq. (2.5), so that As as consequence of the detailed balance condition in Eq. (2.19), the ground state is frustration free, i.e. each local interaction is minimized separately.

The one dimensional problem
As a warm-up we briefly consider the one-dimensional analog to the crystal melting problem, i.e. a system where a set of blocks is added or removed along a line with integer lattice points starting from n = 0. The states are labelled by an integer number |n and the quantum Hamiltonian acts on them as follows (note that we placed ourselves at the V = J point): As shown in the previous section, the ground state is For n > 0, the Hamiltonian can be written in terms of the Laplacian on the line as Hence we can immediately read off the value for the mass gap. In fact the system is easily solvable and one finds the eigenvectors This set of eigenstates does not contain the first excited state |Ψ 1 which would correspond to k = 0. It satisfies This means that |Ψ 1 is harmonic: In one dimension, this is easily solved and using the equation for |0 as initial condition one finds

The two-dimensional problem: quantum partitions
After having discussed the simplest version of the problem in the last section, we now turn to the study of the growth of a two-dimensional crystal. When using the picture of squares being added to an empty corner of the plane, the growth rules are the following: a square can only be added when two of its sides will touch either the wall or the side of another square. It is easy to see that in this case the allowed crystal configurations are in one-to-one correspondence with random partitions. A partition of n is a non-increasing finite sequence of positive integers whose sum is n. We can visually represent a partition by its (Ferrers) diagram, for which we adhere to the Russian tradition, see Figure 4, where the corner of the plane is rotated by 135 degrees. The theory of partitions is at the same time an old subject of mathematics with many important results dating back to Euler, and one which has seen major progress in the last fifty years (see e.g. [9] for an overview).
Let { C i } be the set of all integer partitions (or equivalently Young diagrams). Since we are interested in the size of a partition, it is natural to define the "classical" Hamiltonian The classical partition function for 2d partitions is given by Using the procedure detailed in Sec. 2 we define the graph of configurations whose nodes are in one-to-one correspondence with the C i and where a line is drawn between two partitions if they differ by only one square. The Hamiltonian retains the form given in Eq. (2.33), acting on the Hilbert space generated by the orthonormal basis of the configurations |C i . The ket | now represents a square that can be removed from the partition, while | represents a place where a square can be added. The sum runs over all places where a square can be added or removed in the configuration.
The 2d case of random partitions is especially interesting for us because of the existence of a map between 2d partitions and the Neveu-Schwarz sector of the complex fermionic oscillator. It enables us to express the quantum Hamiltonian (4.3) in terms of fermionic operators.

Fermion operator formalism for 2d partitions
In this section, we will explain the map from partitions to fermions mentioned above. The empty partition is depicted in Figure 3.
all other anti-commutators zero, and the annihilation operators annihilate the vacuum, Note that the vacuum |0 corresponds to the bottom of the Fermi sea. In terms of fermions, the empty configuration (no squares) is represented by where each NW-SE diagonal line corresponds to a creator ψ * a . The state |half , corresponding to the filled Fermi sea (a picture also familiar from matrix models), obeys The partition with one square, see Figure 3.(b), corresponds to which corresponds to an annihilator for the SW-NE line at − 1 2 and a creator for the NW-SE line at 1 2 . On the projection to the horizontal line, the fermion at position − 1 2 was annihilated and a fermion was created at 1 2 , i.e. the fermion at − 1 2 hopped over to 1 2 . Adding squares to the partition corresponds to fermions hopping to the right, in this sense this can be seen as a diffusion process in which every position can be occupied only by one particle, (i.e. it is exclusive).
where n m = ψ * m ψ m is the fermion number operator. The first two terms are the kinetic hopping terms, while the last two terms are the potential terms, counting the number of possibilities to hop right and left per configuration. The different weights for hopping right and left ( √ q and 1/ √ q) make the system into an asymmetric diffusion process. The fermion number n f is conserved in this system, therefore states with differing fermion numbers belong to different superselection sectors S n f of the Hamiltonian.
Let us now go to the RK point V = J. The unique ground state, satisfying H |ground = 0, is written as in Eq. (2.37): where we notice that the number of squares in a given partition µ is given by

The Heisenberg XXZ spin chain
The Hamiltonian in Eq. (4.9) can be recast into a different form, which is more familiar in statistical physics. If we express the creation and annihilation operators in terms of the usual Pauli matrices σ k m , k = 1, 2, 3 at position m, Eq. (4.9) becomes 1 (4.12) Going to the RK point V/J = 1, this is precisely the Hamiltonian for the Heisenberg XXZ spin chain [10] with anisotropy parameter (4.13) 1 Note that our notation differs form the one generally used in condensed matter literature, where not √ q, but q is used which is related to q by q = √ q.
Note that the term proportional to σ 3 m − σ 3 m+1 is the only term in Eq. (4.12) not invariant under the exchange q → 1/q. It is a boundary contribution vanishing for the case of the infinite chain. In fact one could have noticed this symmetry before since for any given 2d partition one square more can be added than removed. It follows that the asymmetry between q and 1/q is only proportional to a trivial constant term. This implies also that we can restrict our attention in this case to V = J without loss of generality.
Since for the partitions, we are interested in the case 0 < q < 1, we have 1 < ∆ < ∞, which means that we are in the ferromagnetic regime. The case ∆ = 1 corresponds to the isotropic Heisenberg or XXX model. For ∆ = ∞, this is the one-dimensional Ising model. The XXZ model is one of the best studied systems in statistical mechanics and can be solved using the Bethe ansatz [11]. It was shown to correspond to the six-vertex or ice-type model.
That the finite XXZ chain admits the quantum group U √ q [SU(2)] as a symmetry was first pointed out by Pasquier and Saleur [12]. A large part of the literature is devoted to the anti-ferromagnetic case or the case −1 < ∆ < 1, but a fair amount of results exists on the ferromagnet, which is collected in [13].
We are specifically interested in the infinite volume chain. It admits four types of zero energy ground states. Two are translation invariant and correspond to all spins up and all spins down. The other two are the kink and the anti-kink [14,15]. This list of ground states was shown to be complete. These ground states are furthermore frustration free, i.e. they minimize each of the next-neighbour interactions separately. Over each of these ground states, a spectral gap exists with value [16] γ = −J q 1/2 + q −1/2 − 2 . (4.14) It was furthermore shown in [17] that droplet excitations exist over all four ground states. In our normalization, their energies are where n is the number of spins that deviate from the ground state and g = − log q. For n = 0, E(0) = 0, while for n = 1, the mass gap Eq. (4.14) is reproduced. These energy levels are symmetric under the exchange q → 1/q as shown in Fig. 5.

Correlation functions for the kink states
Since the partitions correspond to the case of the kink centered in 0, we concentrate on the kink states. The kinks form an infinite family of ground states interpolating between spin up at −∞ and spin down at +∞. These kinks are the lowest energy eigenstates for the superselection sectors with fixed particle number into which the Hilbert space separates (note that [H, ∑ m σ 3 m ] = 0, which means that the action of the Hamiltonian does not change the particle number).

Grand canonical ground state
A useful representation for the kink state can be obtained in terms of the grand canonical ground state which represents a kink centered in x 0 = − log |z| / log q. Note that we shifted the spin chain such that the leftmost spin is at n = 1/2. The kink is exponentially localized [18]. In fact if one considers the magnetization profile m Ψ(z) (x) = Ψ(z)|σ 3 x |Ψ(z) , one finds that For q → 0, the system becomes an Ising model, where the kink is localized in a single point, i.e. the magnetization profile becomes a step function centered in x = x 0 . For q = 1, where the symmetry is the usual SU(2), the system is the XXX spin chain which only admits translation-invariant ground states (one could say that the kink is completely delocalized).
The grand-canonical kink (4.16) can be seen as a generating function for the kinks in the N-particle sectors: (4.18) In particular, our ground state corresponds to the half-filled infinite spin chain, i.e. the kink centered in the middle of the chain and is therefore precisely |Ψ L/2 .
Our strategy to compute correlation functions is the following: given an operator O that does not depend on z and leaves the fermion number invariant, one can expand the correlator on the grand-canonical kink as where we used the fact that the |Ψ N are orthogonal since they describe states with different particle numbers. It follows that the correlation function on the N-particle kink is given by This is particularly useful when O does not mix contributions from different points in which case the correlator on the grand-canonical kink can be found easily.
As a first example consider O = 1, i.e. the norm where we used the notation of the q-shifted factorials and w = z 2 . Using the q-binomial theorem, we immediately find that the norm of an N-particle kink is where the factor q N 2 /2 takes the shift of the spin chain into account.

One-point function
In general, other quantities of interest such as n-point correlation functions can be evaluated for Ψ(z) and then expanded in a series of w to get the contribution of the N-particle kink.
As an example, let µ(x) denote the profile of a partition µ, (i.e. a piecewise C ∞ function). The map to fermionic states is such that if µ (x) = −1, there is a spin up at positionx in the chain, while if µ (x) = 1, there is a spin down atx. Consider the operator σ 3 x acting on the pointx, If we define the magnetization for a state |Φ by where µ Φ (x) is the shape of the partition associated to |Φ . The evaluation of the one-point function for the generating function of the kinks yields By developing in a series (see App. A) one can show that the contribution of the N-particle kink is given by Shifting the center and sending N → ∞ one finds the exact expression for the shape of the kink: To study the q → 1 − limit one can introduce a q-difference equation satisfied by m ∞ (x), see App. A. It follows that the magnetization profile is The shape of the partition is hence given by or, in a more symmetric form (take µ φ (x) =: y), e (x−y)/2 + e −(x+y)/2 = 1 , (4.32) see Fig. 6. This reproduces the limit shape for 2d partitions given in [6]. (b) Limit shape Figure 6: Magnetization (a) and limit shape (b) for the ground state of the two-dimensional partition.

Two-point function
Given the factorized form of Ψ(z) in Eq. (4.16), one can easily evaluate two-point functions.
In particular, Let q x i = ζ i . We have the expansion We can now repeat the same construction as for the one-point function and find the contribution of the N-particle kink: In the N → ∞ limit, this becomes To take the q → 1 − limit, we define f (x) as in Eq. (A. 15) and, after the rescaling x i → −x i / log(q), (4.37) which shows that the two-point correlation function factorizes. The same holds true for the n-point functions as we show in the following section.

The n-point function
The calculation of the n-point function can be simplified by introducing a recursive procedure to obtain the (n + 1)-point function from the n-point function.
The n-point function for the grand-canonical ensemble is written as where w = z 2 and ζ i = q x i +1/2 . It is convenient to develop it in a series in w: Passing to the (n + 1)-point function just amounts to multiplying by an extra factor of Developing in a series,  and yields where the indices i and j are understood in Z n . Now we can again use the same procedure outlined in App. A and find that on the Nparticle kink, (4.46) and in the N → ∞ limit this becomes: In the q → 1 − limit, this is We see that the system decorrelates completely.

Height function and numerics for higher eigenstates
We would like to form an idea about the excited eigenstates of the 2d quantum crystal. Since at this point, we do not have an exact expression, we resort to numerical calculations for finite systems containing partitions up to a certain number.
The height function is defined on the endpoints of the squares making up the partitions.
For a system allowing partitions fitting in a box of length N 0 , the height function is represented on a 1d integer lattice of length 2N 0 + 1. We can visually represent any eigenstate by plotting its average height. So we first have to express each individual partition in terms of a height function and for a given eigenstate sum all these heights with the probability coefficients of this eigenstate. More precisely, after the quantization of the system the height function becomes an observable and as such it has an expectation value for a given state |Ψ , where the sum extends over the partitions |α . We consider the contribution of each partition to the height function minus the profile of the empty corner, which has the form { . . . , 2, 1, 0, 1, 2, . . . }. A partition of N is given by N positive integers (a 0 , a 1 , . . . , a N−1 ) with ∑ N−1 i=0 a i = N. Each a i gives rise to a vector δ a i with all entries zero except for a sequence of a i twos from position i − a i + 1 to i. The sum of these N contributions gives the profile of the partition. The height function of a partition µ at position m is thus given by

Fermionic formalism
The general problem of the 3d melting crystal corner has already been explained in Section 2.2. We now want to express the 3d quantum Hamiltonian in terms of fermionic operators, analogously to the treatment of the 2d case in Section 4.1.
Not surprisingly, the quantum Hamiltonian for the three-dimensional partitions can be recast in terms of a system of Heisenberg XXZ chains. To do so, we first map the threedimensional partitions into a system of vicious walkers. Vicious walkers were introduced by Fisher in [19] as a system of N particles moving in non-intersecting random walks. The mapping to plane partitions is obtained as follows. Consider the finite problem in a cubic box of length N. The empty partition can be represented as an equilateral hexagon tessellated with rhombi. Choose the midpoints of the sides of these rhombi on the lower left side of the hexagon and join them with the corresponding midpoints of the rhombi on the upper right edge, following the sides of the cubes (Fig. 8(b)). One can verify that each partition corresponds to a realization of N non-intersecting directed paths, i.e. directed vicious walks on the lattice obtained by taking the midpoints of the rhombi in the tessellation of the empty partition. In particular, the empty partition is obtained as a set of parallel walks that first move to the right and then up.
It is convenient to represent the system on a square lattice by tilting the picture (Fig. 8(c)).
One can easily see that the shape of each vicious walk corresponds to a Young diagram obtained by parallel-slicing the three-dimensional partition. The set of N paths can be further mapped to a two-dimensional system of fermions living on the vertices of the lattice joined by a path (Fig. 8(d)). Let us now consider the quantum Hamiltonian of this system. Adding a cube corresponds to moving a fermion from (i, j) to (i − 1, j + 1). But this can only happen if the points (i − 1, j) and (i, j + 1) are already occupied. Similarly, removing a cube corresponds to jumping in the opposite way, respecting the same conditions. The Hamiltonian can therefore be written as is the Hamiltonian of the two-dimensional problem 2 . We find that the Hamiltonian describes an infinite system of parallel XXZ chains, each interacting with the two neighbouring ones. Equivalently, the system also represents the dynamics of a heap of dimers [20,21].
Interestingly enough, the same Hamiltonian is obtained when considering a diagonal slicing as shown in Fig. 9 instead of the vicious walker description which amounts to a parallel slicing 3 .
There is a direct mapping between the rhombus tiling making up the shape of the cubes (see Figure 8(b)) and the dimer model on the hexagonal graph (see also Figure 1). To the long diagonal of each rhombus we associate an edge covered by a dimer in the underlying lattice.
In this way, a perfect matching determines the orientations of the various rhombi in the tiling and therefore the shape of the 3d partition. From the dimer picture on the hexagonal lattice we can arrive directly at the fermion picture in Figure 9

Hexagonal lattice description
We can make use of the correspondence to the dimer model on the hexagonal lattice and describe the fermions in the plane more conveniently using an overcomplete basis. The plaquettes of a planar hexagonal lattice are labeled by a triplet of positive integers { n 1 , n 2 , n 3 } such that the center of each hexagon has coordinates n 1 e 1 + n 2 e 2 + n 3 e 3 . The e i are unit vectors that can be chosen to be e 1 = (− cos π 6 , − sin π 6 ) , e 2 = (cos π 6 , − sin π 6 ) , and e 3 = (0, 1) .
Note that the three unit vectors are not linearly independent and satisfy e 1 + e 2 + e 3 = 0 which implies that two triplets { n 1 , n 2 , n 3 } and { m 1 , m 2 , m 3 } represent the same point if We thus obtain a system where each fermion can jump from the point x to x + e 3 if both the points x − e 1 and x − e 2 are free (see Fig. 9(c)). Accordingly, we rewrite the Hamiltonian in Eq. (5.1) as (5.4) In this notation, the empty partition corresponds to the state |half = ∞ ∏ n 1 ,n 2 =0 ψ * n 1 e 1 +n 2 e 2 |0 . (5.5) A generic partition is expressed by where x i = α i e 1 + β i e 2 and α i , β i , and γ i are all positive. To be admissible, a set { x i , γ i } has to satisfy certain interlacing conditions to make sure it defines a valid plane partition.
These conditions are particularly simple in this notation once one realizes that (α i , β i , γ i ) are the coordinates of one of the corners of the topmost box in each column forming the three-dimensional partition. This means that for each column we have to verify that the two columns behind are taller, namely: The number of boxes in a given partition µ is The ground state of the quantum crystal expressed in fermionic language is given by where the sum runs over all the sets { x i , γ i } satisfying the condition given in Eq. (5.7).
This representation permits an explicit mapping between three-dimensional partitions and perfect matchings on the (bipartite) hexagonal lattice. Each perfect matching PM defines a unique height function on the plaquettes, up to a constant. This function h PM is defined as follows: ∓2 if e ∈ PM, where A and B are two plaquettes that share the edge e, the sign depends on the orientation 4 .
The empty room height function corresponds to the unique (up to translation) empty room perfect matching (See Figure 1.(a)) and is given by h 0 (n 1 e 1 + n 2 e 2 + n 3 e 3 ) = n 1 + n 2 + n 3 − 3 min { n 1 , n 2 , n 3 } + k , (5.11) where k is a constant which we set here to k = 0. The map from the partition µ expressed as in Eq. (5.6) to the height function h µ (y) and hence a perfect matching is obtained as follows:

Numerical results for higher eigenstates
As already explained in Sec. 4.4, we are interested in the average height function for the excited states. An analytic expression is at this point beyond our reach. We have nevertheless run some numeric simulations for a system containing plane partitions of integers of up to 12 (resulting in a total of 3462 three-dimensional configurations). Some of the results are shown in Figure 10. In particular one finds that the excited levels are more spread out than the ground state and we have some evidence that the first excited level has Z 3 symmetry and is non-degenerate, while the second excited level has only Z 2 symmetry and is thrice degenerate.

The mass gap
We have remarked already that the mass gaps of the 1-dimensional (see formula (3.4)) and two-dimensional quantum crystals (see Eq. (4.14)) are the same. For the 3d-quantum crystal, we have strong numerical evidence for the mass gap taking again at the same value, see Figure 11. We therefore conjecture that The mass gap of the quantum crystal is the same in all dimensions and equals to Why this should be true can heuristically be understood as follows. Consider the graph of configurations (Fig. 2), which is a directed, acyclic graph. This graph admits a height function, the height of a configuration being its number of boxes. The quantum crystal melting can be considered an asymmetric random walk on this graph, where a particle hops from one configuration to the next. This random walk corresponds directly to the 1-particle sector of an XXZ model living on the state graph. In one dimension, this corresponds directly to the XXZ model on the line with one particle moving around. In two dimensions, it corresponds to the half-filled sector of the 1d XXZ model and in three dimensions to the quantum crystal.
In all cases the Hamiltonian takes the form where deg − (α) and deg + (α) are the outdegree and indegree of the node α. We would like to argue that the mass gap is the same for the 1-particle sector of the XXZ model on any directed, acyclic graph satisfying the condition This is the case in the one-dimensional system in Sec. 3 (deg − (α) = deg + (α) = 1) and the two-dimensional system of Sec. 4 (deg − (α) = deg + (α) + 1). We conjecture it to be the same in three dimensions for the quantum crystal.
A similar result was already proven in [22] for the case of the graph of configurations being a tree.   Figure 11: The mass gap as a function of √ q. The red line corresponds to y = q 1/2 + q −1/2 − 2.
The blue points result from a numerical simulation in a cube of size 3 × 3 × 4 with a total of 4116 explored three-dimensional partitions.

Conclusions
In this note we have studied systems of crystal melting in one, two and three dimensions.
Generalizing the work of Rokhsar and Kivelson on the quantum dimer we obtain a quantum Hamiltonian whose unique ground state reproduces the classical Boltzmann distribution for the canonical ensemble.
The Hamiltonian becomes tractable once one realizes that it is possible to introduce a mapping to spin systems. More precisely, we found that the two-dimensional quantum crystal problem (which can also be phrased in terms of random partitions) is integrable, since it corresponds directly to the half-full sector of the XXZ ferromagnetic spin chain. In particular, the random height function can be interpreted as the integral of the magnetization of the underlying spin system. In a similar fashion, the three-dimensional quantum crystal can be described in terms of an infinite system of coupled XXZ spin chains. This, in turn, allowed us to carry out a numerical analysis, leading to the conjecture of a mass gap whose value does not depend on the dimensionality of the problem.
From this point on there is a variety of interesting questions to address, mainly concerning the three dimensional system. A rigorous proof for the value of the mass gap is required.
We would also like to obtain a better understanding of the full spectrum and explicit expressions for the higher-order eigenvectors in terms of three-dimensional partitions. Based on the number of conserved quantities, one is inclined to believe also the 3d quantum system to be integrable. If this is indeed the case is another a question worth studying. In this vein one could also try to generalize the Bethe ansatz for the XXZ spin chain to the constrained system of spin chains corresponding to the 3d case. Furthermore, we wonder whether the SU q (2) quantum symmetry of the XXZ spin chain carries over to the 3d melting quantum crystal.

A One-point function
In this appendix we describe in detail the computation for the one-point function for the two-dimensional kink introduced in Sec. 4.3.2.
The starting point is the one-point function for the generating function of the kinks.
Given the representation in terms of direct product over the spins (grand canonical ensemble) one finds: Ψ(z)|σ 3 x |Ψ(z) = 1 − w q x+1/2 2 (1 + w q x+1/2 ) (−w q 1/2 ; q) ∞ . (A.1) Using the development one can expand: The contribution of the N-particle kink is therefore where [ N k ] q is the q-binomial coefficient N k q = (q; q) N (q; q) N−k (q; q) k . (A.5) Using the identity (see [23]) N k q = (q −N ; q) k (q; q) k (−1) k q Nk q −( To study the q → 1 limit, introduce the function One can verify that f q (ζ) satisfies the q-difference equation and introducing the q-difference operator as in [23]: one obtains (1 − q) D q f q (ζ) = − 1 ζ 2 + 1 + ζ ζ 2 f q (ζ) . (A.14) In the q → 1 − limit, D q becomes the derivative with respect to ζ. Define Solving the differential equation for f (x), one obtains where E n (x) is the exponential integral E n (x) = f (x) = e x 1 + e x .
(A. 19) It follows that the magnetization profile is The shape of the partition is therefore given by 21) or, in a more symmetric form (take µ φ (x) =: y), e (x−y)/2 + e −(x+y)/2 = 1 .