$\alpha_{i}\in[0,1]$ for all $i$. Let us now define a vector $\bm{v}$ as

$\bm{v}:=\Big{(}\underbrace{1,\ldots,1}_{L},0\ldots,0\Big{)}.$ (222)

Using this $\bm{v}$, we can write the following

$\displaystyle\bm{v}\cdot
Q\bm{p}^{\downarrow}=\bm{v}\cdot\tilde{\bm{q}}^{\downarrow}=\bm{v}\cdot\bm{p}^{\downarrow}=(1-\epsilon)$
(223) $\displaystyle\Rightarrow$
$\displaystyle\sum_{i}\alpha_{i}\bm{v}\cdot\Pi_{i}\bm{p}^{\downarrow}=\bm{v}\cdot\bm{p}^{\downarrow}$

where we use Eq. (214) in the first line, and the convex decomposition of $Q$
from Eq. (221) in the second line. Note that Eq. (223) can be equivalently
written as

$\displaystyle\sum_{i}\alpha_{i}(\Pi_{i}^{T}\bm{v})\cdot\bm{p}^{\downarrow}=\bm{v}\cdot\bm{p}^{\downarrow},$
(224)

where $\Pi_{i}^{T}$ is the transpose of the permutation matrix $\Pi_{i}$.
Because the elements of $\bm{v}$ and $\bm{p}^{\downarrow}$ are arranged in a
decreasing order, the maximum value of
$(\Pi_{i}^{T}\bm{v})\cdot\bm{p}^{\downarrow}$ is
$\bm{v}\cdot\bm{p}^{\downarrow}$. Thus, we see that equality in Eq. (224)
holds if and only if $\bm{v}$ is invariant under $(\Pi_{i}^{T})$ for all $i$.
Therefore, $\Pi_{i}$ can be only of the form

$\Pi_{i}=\Pi_{i}^{L}\oplus\Pi_{i}^{L^{c}},$ (225)

where $\Pi_{i}^{L}$ and $\Pi_{i}^{L^{c}}$ are permutations that act trivially
for indices $i>L$ and $i\in\\{1,\ldots,L\\}$, respectively. Thus, we infer
using Eq. (221) that $Q$ is a block-diagonal matrix of the form

$Q=\tilde{B}\oplus B,$ (226)

where $\tilde{B}$ and $B$ are $L\times L$ and $(d-L)\times(d-L)$ bistochastic
matrices. Thus,

$\displaystyle\forall
i>L,\quad\tilde{q}_{i}^{\downarrow}=\sum_{j>L}B_{ij}p_{j}^{\downarrow}.$ (227)

Combining Eq. (220) and Eq. (227) we conclude that
$\tilde{\bm{q}}^{\downarrow}={\bm{q}}^{*\downarrow}$ with $\bm{q}^{*}$ defined
in Eq. (93). This finally implies $\tilde{\bm{q}}=\Pi\bm{q}^{*}$ with $\Pi$
being some permutation. Moreover, it is straightforward to show that the
minimal entropy of $\bm{q}^{*}$ is achieved for $B$ being the identity matrix,
since the entropy is increasing under application of any non-trivial
bistochastic matrix. ∎

##  Appendix C Eliminating the logarithmic term

We start with the following lemma that will be needed to prove our claim.

######  Lemma 8.

For a fixed value of $b>0$ and any $\bm{s}$, such that
$\|\bm{s}\|=\sqrt{\sum_{i=1}^{d}s_{i}^{2}}\leq b$, we have

$\log P^{N}_{\bm{k}(\bm{s})}=O(\log N),$ (228)

where $P^{N}_{\bm{k}}$ and $\bm{k}(\bm{s})$ are defined by Eqs. (126) and
(136), respectively.

###### Proof.

We start by using the definition,

$P^{N}_{\bm{k}(\bm{s})}=\binom{N}{k_{1}(s_{1}),...,k_{d}(s_{d})}p^{k_{1}(s_{1})}_{1}...p^{k_{d}(s_{d})}_{d},$
(229)

to write $\log P^{N}_{\bm{k}(\bm{s})}$ as

$\\!\\!\\!\log P^{N}_{\bm{k}(\bm{s})}\\!=\\!\log N!-\sum_{i=1}^{d}\log
k_{i}(s_{i})!+\sum_{i=1}^{d}k_{i}(s_{i})\log p_{i}.\\!$ (230)

Employing Stirling inequality,

$\displaystyle\log\sqrt{N}+N\log N-N$ $\displaystyle\qquad\leq\log N!\leq$
$\displaystyle\qquad\qquad 1+\log\sqrt{N}+N\log N-N,$ (231)

we first provide a lower bound for $\log P^{N}_{\bm{k}(\bm{s})}$,

$\displaystyle\log P^{N}_{\bm{k}(\bm{s})}\geq$
$\displaystyle(\log\sqrt{N}+N\log N-N)+\sum_{i=1}^{d}k_{i}(s_{i})\log p_{i}$
$\displaystyle-\sum_{i=1}^{d}(1\\!+\\!\log\sqrt{k_{i}(s_{i})}\\!+\\!k_{i}(s_{i})\log
k_{i}(s_{i})\\!-\\!k_{i}(s_{i}))$ $\displaystyle=$
$\displaystyle-\sum_{i=1}^{d}k_{i}(s_{i})\log\Big{(}\frac{k_{i}(s_{i})}{Np_{i}}\Big{)}-d$
$\displaystyle+\frac{1}{2}\Big{(}\log N-\sum_{i=1}^{d}\log
k_{i}(s_{i})\Big{)}.$ (232)

Recall that $k_{i}(s_{i})=Np_{i}+\sqrt{N}s_{i}$, which implies
$\sum_{i=1}^{d}s_{i}=0$. To simplify the above further, we lower bound the
first term by employing the inequality $\log(1+g)<g$ for $g>-1$ in the
following way:

$\displaystyle\sum_{i=1}^{d}k_{i}(s_{i})\log\Big{(}\frac{k_{i}(s_{i})}{Np_{i}}\Big{)}=\sum_{i=1}^{d}k_{i}(s_{i})\log\left(1+\frac{s_{i}}{\sqrt{N}p_{i}}\right)$
$\displaystyle\qquad\leq\sum_{i=1}^{d}k_{i}(s_{i})\frac{s_{i}}{\sqrt{N}p_{i}}=\sqrt{N}\sum_{i=1}^{d}(1+\frac{s_{i}}{\sqrt{N}p_{i}})s_{i}$
$\displaystyle\qquad\qquad=\sum_{i=1}^{d}\frac{s_{i}^{2}}{p_{i}}\leq\sum_{i=1}^{d}\frac{s_{i}^{2}}{p_{\min}}\leq\frac{b^{2}}{p_{\min}},$
(233)

where $p_{\min}=\text{min }\\{p_{1},\ldots,p_{d}\\}$. Moreover, observing that

$\log N\geq\left(\log N-\sum_{i=1}^{d}\log k_{i}(s_{i})\right)\geq-(d-1)\log
N,$ (234)

we can conclude that

$\log N-\sum_{i=1}^{d}\log k_{i}(s_{i})=O(\log N).$ (235)

Putting it together, we further simplify the bound given in Eq. (232) as

$\\!\\!\\!\log
P^{N}_{\bm{k}(\bm{s})}\geq-d-\frac{b^{2}}{p_{\min}}-\frac{(d-1)}{2}\log{N}=O(\log
N).$ (236)

Similarly, by employing the Stirling inequality from Eq. (231), we also prove
an upper bound for $\log P^{N}_{\bm{k}(\bm{s})}$ as follows

$\displaystyle\log P^{N}_{\bm{k}(\bm{s})}\leq$
$\displaystyle(1+\log{\sqrt{N}}+N\log N-N)+\sum_{i=1}^{d}k_{i}(s_{i})\log
p_{i}$
$\displaystyle-\sum_{i=1}^{d}\Big{(}\log\sqrt{k_{i}(s_{i})}+k_{i}(s_{i})\log
k_{i}(s_{i})-k_{i}(s_{i})\Big{)}$ $\displaystyle=$
$\displaystyle-\sum_{i=1}^{d}k_{i}(s_{i})\log\Big{(}\frac{k_{i}(s_{i})}{Np_{i}}\Big{)}+1$
$\displaystyle+\frac{1}{2}(\log N-\sum_{i=1}^{d}\log k_{i}(s_{i})).$ (237)

Using the inequality $\log(1+g)\geq\frac{g}{1+g}$ for $g>-1$, we have

$\displaystyle\sum_{i=1}^{d}k_{i}(s_{i})\log\Big{(}\frac{k_{i}(s_{i})}{Np_{i}}\Big{)}=\sum_{i=1}^{d}k_{i}(s_{i})\log\Big{(}1+\frac{s_{i}}{\sqrt{N}p_{i}}\Big{)}$
$\displaystyle\qquad\geq\sum_{i=1}^{d}k_{i}(s_{i})\frac{s_{i}}{s_{i}+\sqrt{N}p_{i}}=\sqrt{N}\sum_{i=1}^{d}s_{i}=0.$
(238)

The above inequality together with Eq. (235) imply that $\log
P^{N}_{\bm{k}(\bm{s})}\leq O(\log N)$ which completes the proof.

∎

Using Lemma 8, we will now be able to prove our claim that is captured by the
following result.

######  Lemma 9.

The following limits are equal:

$\displaystyle\lim_{N\rightarrow\infty}\sum_{\bm{s}}\Big{\\{}P^{N}_{\bm{k}(\bm{s})}\Big{|}\frac{1}{N}\log(P^{N}_{\bm{k}(\bm{s})})+\frac{\beta}{\sqrt{N}}\sum_{i}s_{i}E_{i}$
$\displaystyle\phantom{\lim_{N\rightarrow\infty}\sum_{\bm{s}}\Big{\\{}P^{N}_{\bm{k}(\bm{s})}\Big{|}}+F_{\rm
diss}^{N}\geq 0\Big{\\}}$
$\displaystyle~{}~{}=\lim_{N\rightarrow\infty}\sum_{\bm{s}}\Big{\\{}P^{N}_{\bm{k}(\bm{s})}\Big{|}\frac{\beta}{\sqrt{N}}\sum_{i}s_{i}E_{i}+F_{\rm
diss}^{N}\geq 0\Big{\\}}.$ (239)

###### Proof.

We start by introducing the following notation

$\displaystyle
A(N):=\sum_{\bm{s}}\Big{\\{}P^{N}_{\bm{k}(\bm{s})}\Big{|}\frac{1}{N}\log(P^{N}_{\bm{k}(\bm{s})})+\frac{\beta}{\sqrt{N}}\sum_{i}s_{i}E_{i}$
$\displaystyle\qquad\qquad\qquad\qquad\quad+F_{\text{diss}}^{N}\geq
0\Big{\\}},$ $\displaystyle
A(b,N):=\sum_{\bm{s}}\Big{\\{}P^{N}_{\bm{k}(\bm{s})}\Big{|}\frac{1}{N}\log(P^{N}_{\bm{k}(\bm{s})})+\frac{\beta}{\sqrt{N}}\sum_{i}s_{i}E_{i}$
$\displaystyle\qquad\qquad\qquad\qquad\quad+F_{\text{diss}}^{N}\geq 0\text{
such that }\|\bm{s}\|\leq b\Big{\\}},$ $\displaystyle
B(N):=\sum_{\bm{s}}\Big{\\{}P^{N}_{\bm{k}(\bm{s})}\Big{|}\frac{\beta}{\sqrt{N}}\sum_{i}s_{i}E_{i}+F_{\text{diss}}^{N}\geq
0\Big{\\}},$ $\displaystyle
B(b,N):=\sum_{\bm{s}}\Big{\\{}P^{N}_{\bm{k}(\bm{s})}\Big{|}\frac{\beta}{\sqrt{N}}\sum_{i}s_{i}E_{i}+F_{\text{diss}}^{N}\geq
0$ $\displaystyle\qquad\qquad\qquad\qquad\quad\text{ such that }\|\bm{s}\|\leq
b\Big{\\}},$
$\displaystyle\Omega(b,N):=\sum_{\bm{s}}\Big{\\{}P^{N}_{\bm{k}(\bm{s})}\Big{|}\text{
such that }\|\bm{s}\|\geq b\Big{\\}}.$ (240)

Our goal is to show that

$\lim_{N\rightarrow\infty}A(N)=\lim_{N\rightarrow\infty}B(N).$ (241)

From the definition it follows that

$\displaystyle A(N)-\Omega(b,N)$ $\displaystyle\leq$ $\displaystyle A(b,N)\leq
A(N),$ (242a) $\displaystyle B(N)-\Omega(b,N)$ $\displaystyle\leq$
$\displaystyle B(b,N)\leq B(N).$ (242b)

Taking the limit $N\rightarrow\infty$ of Eqs. (242a) and (242b), we have

$\displaystyle\lim_{N\rightarrow\infty}\Big{(}A(N)-\Omega(b,N)\Big{)}$
$\displaystyle\leq\lim_{N\rightarrow\infty}A(b,N)\leq\lim_{N\rightarrow\infty}A(N),$
(243a) $\displaystyle\lim_{N\rightarrow\infty}\Big{(}B(N)-\Omega(b,N)\Big{)}$
$\displaystyle\leq\lim_{N\rightarrow\infty}B(b,N)\leq\lim_{N\rightarrow\infty}B(N).$
(243b)

Now, let us define

$\displaystyle\lim_{N\rightarrow\infty}\Omega(b,N)=:\epsilon(b).$ (244)

As the multinomial distribution concentrates around mean for
$N\rightarrow\infty$, it follows that
$\lim_{b\rightarrow\infty}\epsilon(b)=0$. Therefore, taking the limit
$b\rightarrow\infty$ in Eq. (243a) we have

$\displaystyle\lim_{N\rightarrow\infty}A(N)\leq\lim_{b\rightarrow\infty}\lim_{N\rightarrow\infty}A(b,N)\leq\lim_{N\rightarrow\infty}A(N)$
$\displaystyle\qquad\Rightarrow\quad\lim_{b\rightarrow\infty}\lim_{N\rightarrow\infty}A(b,N)=\lim_{N\rightarrow\infty}A(N).$
(245)

Analogously, taking the limit $b\rightarrow\infty$ in Eq. (243b) we can show
that

$\lim_{b\rightarrow\infty}\lim_{N\rightarrow\infty}B(b,N)=\lim_{N\rightarrow\infty}B(N).$
(246)

Moreover, for any fixed $b$, by employing Lemma 8, we see that
$\frac{1}{N}\log(P^{N}_{\bm{k}(\bm{s})})=O(\frac{\log N}{N})$, which vanishes
as $N\rightarrow\infty$. Therefore, we have

$\lim_{N\rightarrow\infty}A(b,N)=\lim_{N\rightarrow\infty}B(b,N),$ (247)

and so by taking the limit $b\rightarrow\infty$, we arrive at

$\displaystyle\lim_{b\rightarrow\infty}\lim_{N\rightarrow\infty}A(b,N)=\lim_{b\rightarrow\infty}\lim_{N\rightarrow\infty}B(b,N).$
(248)

Combining the above with Eqs. (245)-(246) we have

$\lim_{N\rightarrow\infty}A(N)=\lim_{N\rightarrow\infty}B(N),$ (249)

which completes the proof. ∎

##  Appendix D Central limit theorem for multinomial distribution

######  Lemma 10.

The multinomial distribution with mean $\bm{\mu}=N\bm{p}$ and a covariance
matrix $\bm{\Sigma}$ can be approximated in the asymptotic limit by a
multivariate normal distribution $\mathcal{N}^{(\bm{\mu},\bm{\Sigma})}$ with
mean $\bm{\mu}$ and a covariance matrix $\bm{\Sigma}$.

###### Proof.

Assume $\bm{X}_{1}\ldots\bm{X}_{N}$ are independent and identically
distributed random vectors each of them with the following distribution

$\\!\\!\\!\text{Prob}(\bm{X}=\bm{x})=\begin{cases}\prod_{i=1}^{d}p_{i}^{x_{i}}&\quad\\!\text{if
$\bm{x}$ is unit vector},\\\ 0&\quad\\!\text{otherwise}.\end{cases}$ (250)

Then, the mean vector of $\bm{X}$ is $\bm{p}$ and the covariance matrix
$\frac{1}{N}\bm{\Sigma}=\text{ diag }(\bm{p})-\bm{p}\bm{p}^{T}$. Define
$\bm{S}_{N}:=\bm{X}_{1}+\ldots+\bm{X}_{N}$. Then

$\displaystyle\text{Prob}(\bm{S}_{N}=\bm{k})$ $\displaystyle=$
$\displaystyle\binom{N}{k_{1}\ldots k_{d}}p_{1}^{k_{1}}\ldots p_{d}^{k_{d}}.$
(251)

We thus see that a multinomial distribution arises from a sum of independent
and identically distributed random variables. Therefore, using the central
limit theorem, we obtain that the distribution of $\bm{k}$ approaches the
distribution $\mathcal{N}^{(\bm{\mu},\bm{\Sigma})}$ arbitrarily well for
$N\to\infty$, which completes the proof. ∎

##  Appendix E Entropy difference between uniformised and non-uniformised
embedding boxes

In this appendix, we show that there exists a permutation which transforms the
total final state $\hat{\bm{F}}^{N}$ from Eq. (94) into a state that is
uniform in almost all embedding boxes, leading to no dissipation up to higher-
order asymptotics. We start with the following observations about the total
initial and final states. First, since we consider the initial state composed
of $N$ independent systems in identical incoherent states
$\bm{P}^{N}=\bm{p}^{\otimes N}$, we note that
$\hat{\bm{P}}^{N}\otimes\hat{\tilde{\bm{G}}}^{N}$ has a polynomial number of
distinct entries with exponential degeneracy and exponentially many embedding
boxes. Second, note that according to Eq. (94) the entries of the embedded
total final state, $\hat{\bm{F}}^{N}$, are essentially given by the entries of
$\hat{\bm{P}}^{N}\otimes\hat{\tilde{\bm{G}}}^{N}$ (plus many equal entries
$(1-\epsilon)/K$). Thus, $\hat{\bm{F}}^{N}$ has polynomially many distinct
entries with exponential degeneracy and exponentially many embedding boxes.
Employing permutational freedom, we can then rearrange the entries of
$\hat{\bm{F}}^{N}$ so that they are uniform in almost every embedding box,
except poly($N$) of them. We will denote the probability distribution over the
embedding boxes by $\bm{q}$ and note that it is essentially equal to
$\bm{P}^{N}\otimes\,\tilde{\bm{G}}^{N}$. Moreover, we will denote by
$\bm{r}^{(i)}$ the normalised distribution within the $i$-th embedding box.

The next step is to write the entropy of the embedded total final state as the
entropy of probability distribution over different embedding boxes,
$H(\bm{q})$, plus the average entropy of normalised probabilities within each
box, $H(\bm{r}^{(i)})$,

$H(\hat{\bm{F}}^{N})=H(\bm{q})+\sum_{i}q_{i}H(\bm{r}^{(i)}).$ (252)

Note that whether we uniformise or not a given box, the entropy $H(\bm{q})$
does not change. Consequently, the entropy of the final state uniformised
within each embedding box takes the form of

$\displaystyle H(\hat{\bm{F}}_{\rm{uni}}^{N})$
$\displaystyle=H(\bm{q})+\sum_{i}q_{i}H(\bm{r}^{(i)}_{\rm{uni}}),$ (253)

with $\bm{r}^{(i)}_{\rm{uni}}$ representing the normalised and uniformised
probability within the $i$-th embedding box. Thus, the entropy difference
between the uniformised and non-uniformised distributions reads

$H(\hat{\bm{F}}_{\rm{uni}}^{N})-H(\hat{\bm{F}}^{N})=\sum_{i}q_{i}\left(H(\bm{r}^{(i)}_{\rm{uni}})-H(\bm{r}^{(i)})\right).$
(254)

Now, note that due to the previous argument, the above sum is performed only
over a polynomial number of boxes. Denoting the set with size poly($N$) as
$\Omega$, one can write

$\displaystyle H(\hat{\bm{F}}_{\rm{uni}}^{N})-H(\hat{\bm{F}}^{N})$
$\displaystyle=\sum_{i\in\Omega}q_{i}\left(H(\bm{r}^{(i)}_{\rm{uni}})-H(\bm{r}^{(i)})\right)$
$\displaystyle\leq\sum_{i\in\Omega}q_{i}H(\bm{r}^{(i)}_{\rm{uni}}).$ (255)

Let us analyse the right-hand side of the above equation. First, note that
$q_{i}$ is exponentially small, $q_{i}\propto\exp(-N)$; while
$H(\bm{r}^{(i)}_{\rm{uni}})$ is the entropy of a uniform state over the
dimension of the $i$-th embedding box (equal to $\prod_{i}D^{k_{i}}_{i}$), so
it will scale linearly in $N$:

$H(\bm{r}^{(i)}_{\rm{uni}})=\log\left(\prod_{i}D^{k_{i}}_{i}\right)\propto
O(N)\,.$ (256)

Therefore, we conclude that the entropy difference vanishes exponentially in
$N$

$H(\hat{\bm{F}}_{\rm{uni}}^{N})-H(\hat{\bm{F}}^{N})\propto\exp(-N).$ (257)

## References

  * Carathéodory [1909] C. Carathéodory, Untersuchungen über die grundlagen der thermodynamik, Math. Ann. 67, 355 (1909).
  * Giles [1964] R. Giles,  _Mathematical Foundations of Thermodynamics_ (Pergamon Press, 1964).
  * Seifert [2008] U. Seifert, Stochastic thermodynamics: principles and perspectives, Eur. Phys. J. B. 64, 423 (2008).
  * Seifert [2012] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Rep. Prog. Phys. 75, 126001 (2012).
  * Einstein [1906] A. Einstein, Zur theorie der brownschen bewegung, Ann. Phys. 324, 371 (1906).
  * von Smoluchowski [1906] M. von Smoluchowski, Zur kinetischen theorie der brownschen molekularbewegung und der suspensionen, Ann. Phys. 326, 756 (1906).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Rep. Prog. Phys. 29, 255 (1966).
  * Bonança [2008] M. V. S. Bonança, Fluctuation-dissipation theorem for the microcanonical ensemble, Phys. Rev. E 78, 031107 (2008).
  * Batalhão  _et al._ [2014] T. B. Batalhão, A. M. Souza, L. Mazzola, R. Auccaise, R. S. Sarthour, I. S. Oliveira, J. Goold, G. De Chiara, M. Paternostro, and R. M. Serra, Experimental reconstruction of work distribution and study of fluctuation relations in a closed quantum system, Phys. Rev. Lett. 113, 140601 (2014).
  * An  _et al._ [2015] S. An, J.-N. Zhang, M. Um, D. Lv, Y. Lu, J. Zhang, Z.-Q. Yin, H. T. Quan, and K. Kim, Experimental test of the quantum jarzynski equality with a trapped-ion system, Nat. Phys 11, 193 (2015).
  * Janzing  _et al._ [2000] D. Janzing, P. Wocjan, R. Zeier, R. Geiss, and T. Beth, Thermodynamic cost of reliability and low temperatures: Tightening landauer’s principle and the second law, Int. J. Theor. Phys. 39, 2717 (2000).
  * Horodecki and Oppenheim [2013] M. Horodecki and J. Oppenheim, (quantumness in the context of) resource theories, Int. J. Mod. Phys. B 27, 1345019 (2013).
  * Brandão  _et al._ [2015] F. G. S. L. Brandão, M. Horodecki, N. H. Y. Ng, J. Oppenheim, and S. Wehner, The second laws of quantum thermodynamics, Proc. Natl. Acad. Sci. U.S.A. 112, 3275 (2015).
  * Lostaglio [2019] M. Lostaglio, An introductory review of the resource theory approach to thermodynamics, Rep. Prog. Phys. 82, 114001 (2019).
  * Dahlsten  _et al._ [2011] O. C. O. Dahlsten, R. Renner, E. Rieper, and V. Vedral, Inadequacy of von neumann entropy for characterizing extractable work, New J. Phys. 13, 053015 (2011).
  * Yunger Halpern  _et al._ [2018] N. Yunger Halpern, A. J. P. Garner, O. C. O. Dahlsten, and V. Vedral, Maximum one-shot dissipated work from rényi divergences, Phys. Rev. E 97, 052135 (2018).
  * Alhambra  _et al._ [2016] A. M. Alhambra, L. Masanes, J. Oppenheim, and C. Perry, Fluctuating work: From quantum thermodynamical identities to a second law equality, Phys. Rev. X 6, 041017 (2016).
  * Halpern  _et al._ [2015] N. Y. Halpern, A. J. P. Garner, O. C. O. Dahlsten, and V. Vedral, Introducing one-shot work into fluctuation relations, New J. Phys. 17, 095003 (2015).
  * Tajima and Hayashi [2017] H. Tajima and M. Hayashi, Finite-size effect on optimal efficiency of heat engines, Phys. Rev. E 96, 012128 (2017).
  * Chubb  _et al._ [2018] C. T. Chubb, M. Tomamichel, and K. Korzekwa, Beyond the thermodynamic limit: finite-size corrections to state interconversion rates, Quantum 2, 108 (2018).
  * Chubb  _et al._ [2019] C. T. Chubb, M. Tomamichel, and K. Korzekwa, Moderate deviation analysis of majorization-based resource interconversion, Phys. Rev. A 99, 032332 (2019).
  * Korzekwa  _et al._ [2019a] K. Korzekwa, C. T. Chubb, and M. Tomamichel, Avoiding irreversibility: engineering resonant conversions of quantum resources, Phys. Rev. Lett. 122, 110403 (2019a).
  * Horodecki and Oppenheim [2013] M. Horodecki and J. Oppenheim, Fundamental limitations for quantum and nanoscale thermodynamics, Nat. Commun. 4, 2059 (2013).
  * Esposito  _et al._ [2009] M. Esposito, U. Harbola, and S. Mukamel, Nonequilibrium fluctuations, fluctuation theorems, and counting statistics in quantum systems, Rev. Mod. Phys. 81, 1665 (2009).
  * Parrondo  _et al._ [2015] J. M. R. Parrondo, J. M. Horowitz, and T. Sagawa, Thermodynamics of information, Nat. Phys 11, 131 (2015).
  * Alicki [1979] R. Alicki, The quantum open system as a model of the heat engine, J. Phys. A: Math. Gen. 12, L103 (1979).
  * Åberg [2013] J. Åberg, Truly work-like work extraction via a single-shot analysis, Nat. Commun. 4, 1925 (2013).
  * Skrzypczyk  _et al._ [2014] P. Skrzypczyk, A. J. Short, and S. Popescu, Work extraction and thermodynamics for individual quantum systems, Nat. Commun. 5, 4185 (2014).
  * Maxwell [1872] J. Maxwell,  _Theory of Heat_, Text-books of science (Longmans, Green, and Company, 1872).
  * Szilard [1929] L. Szilard, über die entropieverminderung in einem thermodynamischen system bei eingriffen intelligenter wesen, Zeitschrift für Physik 53, 840 (1929).
  * Bennett [1982] C. H. Bennett, The thermodynamics of computation—a review, Int. J. Theor. Phys. 21, 905 (1982).
  * Landauer [1961] R. Landauer, Irreversibility and heat generation in the computing process, IBM J. Res. Dev. 5, 183 (1961).
  * Alicki and Horodecki [2019] R. Alicki and M. Horodecki, Information-thermodynamics link revisited, J. Phys. A: Math. Theor. 52, 204001 (2019).
  * Wilde [2013] M. Wilde,  _Quantum Information Theory_, Quantum Information Theory (Cambridge University Press, 2013).
  * Narasimhachar  _et al._ [2019] V. Narasimhachar, J. Thompson, J. Ma, G. Gour, and M. Gu, Quantifying memory capacity as a quantum thermodynamic resource, Phys. Rev. Lett. 122, 060601 (2019).
  * Korzekwa  _et al._ [2019b] K. Korzekwa, Z. Puchała, M. Tomamichel, and K. Życzkowski, Encoding classical information into quantum resources, arXiv:1911.12373 (2019b).
  * Tomamichel and Hayashi [2013] M. Tomamichel and M. Hayashi, A hierarchy of information quantities for finite block length analysis of quantum tasks, IEEE Trans. Inf. Theory 59, 7693 (2013).
  * Boes  _et al._ [2020] P. Boes, N. H. Y. Ng, and H. Wilming, The variance of relative surprisal as single-shot quantifier (2020), arXiv:2009.08391 .
  * Alicki  _et al._ [2004] R. Alicki, M. Horodecki, P. Horodecki, and R. Horodecki, Thermodynamics of quantum information systems — hamiltonian description, Open Syst. Inf. Dyn. 11, 205 (2004).
  * Bhatia [1996] R. Bhatia,  _Matrix Analysis_, Graduate Texts in Mathematics (Springer New York, 1996).
  * Berry [1941] A. C. Berry, The accuracy of the gaussian approximation to the sum of independent variates, Trans. Am. Math. Soc. 49, 122 (1941).
  * Esseen [1942] C. G. Esseen, On the liapunov limit error in the theory of probability, Ark. Mat. Astr. Fys. 49, 1 (1942).
  * Esseen [1956] C. G. Esseen, A moment inequality with an application to the central limit theorem, Scand. Actuar J. 1956, 160 (1956).
  * Shevtsova [2011] I. Shevtsova, On the absolute constants in the berry-esseen type inequalities for identically distributed summands (2011), arXiv:1111.6554 .
  * Lostaglio  _et al._ [2015a] M. Lostaglio, D. Jennings, and T. Rudolph, Description of quantum coherence in thermodynamic processes requires constraints beyond free energy, Nat. Commun. 6, 6383 (2015a).
  * Lostaglio  _et al._ [2015b] M. Lostaglio, K. Korzekwa, D. Jennings, and T. Rudolph, Quantum coherence, time-translation symmetry, and thermodynamics, Phys. Rev. X 5, 021001 (2015b).
  * Ćwikliński  _et al._ [2015] P. Ćwikliński, M. Studziński, M. Horodecki, and J. Oppenheim, Limitations on the evolution of quantum coherences: Towards fully quantum second laws of thermodynamics, Phys. Rev. Lett. 115, 210403 (2015).
  * Woods  _et al._ [2019] M. P. Woods, N. H. Y. Ng, and S. Wehner, The maximum efficiency of nano heat engines depends on more than temperature, Quantum 3, 177 (2019).
  * Ng  _et al._ [2017] N. H. Y. Ng, M. P. Woods, and S. Wehner, Surpassing the carnot efficiency by extracting imperfect work, New J. Phys. 19, 113005 (2017).
  * Wang and Renner [2012] L. Wang and R. Renner, One-shot classical-quantum capacity and hypothesis testing, Phys. Rev. Lett. 108, 200501 (2012).
  * Buscemi and Datta [2010] F. Buscemi and N. Datta, The quantum capacity of channels with arbitrarily correlated noise, IEEE Trans. Inf. Theory. 56, 1447 (2010).
  * Datta  _et al._ [2013] N. Datta, M. Mosonyi, M. Hsieh, and F. G. S. L. Brandão, A smooth entropy approach to quantum hypothesis testing and the classical capacity of quantum channels, IEEE Trans. Inf. Theory. 59, 8014 (2013).
  * Li [2014] K. Li, Second-Order Asymptotics for Quantum Hypothesis Testing, Ann. Stat. 42, 171 (2014).


# The multinomial tiling model

Richard Kenyon111Department of Mathematics, Yale University, New Haven;
<EMAIL_ADDRESS>Cosmin Pohoata11footnotemark: 1

###### Abstract

Given a graph $\mathcal{G}$ and collection of subgraphs $T$ (called tiles), we
consider covering $\mathcal{G}$ with copies of tiles in $T$ so that each
vertex $v\in\mathcal{G}$ is covered with a predetermined multiplicity. The
_multinomial tiling model_ is a natural probability measure on such
configurations (it is the uniform measure on standard tilings of the
corresponding “blow-up” of $\mathcal{G}$).

In the limit of large multiplicities we compute asymptotic growth rate of the
number of multinomial tilings. We show that the individual tile densities tend
to a Gaussian field with respect to an associated discrete Laplacian. We also
find an exact discrete Coulomb gas limit when we vary the multiplicities.

For tilings of $\mathbb{Z}^{d}$ with translates of a single tile and a small
density of defects, we study a crystallization phenomena when the defect
density tends to zero, and give examples of naturally occurring quasicrystals
in this framework.

##  1 Introduction

The study of random tilings is a cornerstone area of combinatorics and
statistical mechanics. In its simplest form, the random tiling model consists
in the study of the set of tilings of a region (for example a subset of the
plane) with translated copies of a finite collection of shapes, called
prototiles. However even the simplest cases can lead to hard problems. The
mere existence of a tiling of a region in $\mathbb{R}^{2}$ with a prescribed
set of polyominos is an NP-complete problem [18], even if the prototiles
consist in just the $3\times 1$ and $1\times 3$ rectangles [1]. Enumerating
tilings is of course even harder.

However in the few cases where we _can_ analyze random tilings, like random
“domino” tilings (tilings with $2\times 1$ and $1\times 2$ rectangles) or
“lozenge” tilings (tilings with $60^{\circ}$ rhombi), we find very rich
behavior, with beautiful enumerative properties [11, 23, 8, 9, 20], phase
transitions [15], limit shapes [14], conformal invariance [12], and Gaussian
scaling limits [13]. Beyond these and other dimer models there are almost no
other cases we can analyze in detail. There are other cases where enumeration
is sometimes possible, like the $6$-vertex model [19], but for these models
very little is known about correlations, although they are sometimes predicted
in physics to be Gaussian in the scaling limit and “conformally
invariant”—such models were in fact the inspiration for conformal field
theory.

We study here a variant of the random tiling problem: the _multinomial tiling
problem_ , which is tractable in the sense that we can give exact generating
functions for enumerations, which in turn yield, in the limit of large
multiplicity, exact asymptotic expressions for growth rates and Gaussian
behavior for random tilings. This setting is quite general and works for
tilings in arbitrary graphs, not just plane regions. Furthermore we find all
of the phenomena discussed above: phase transitions, limit shapes,
crystallization phenomena, and conformal invariance (which we study in [16]).

It comes as an additional surprise that in certain situations our random
tilings form _quasicrystals_. Quasicrystals were first found in nature by
Schechtman et al [21]. Their physical and mathematical framework is still
debated, but examples of quasiperiodic tilings were first found by Berger [2]
and familiar examples like Penrose tilings are now well understood [6]; they
are sets of tiles which tile the plane but only in nonperiodic fashion. Our
quasicrystals arise from _random_ tilings; it is the statistical correlations
between tile densities that are quasiperiodic (even though there are periodic
points in the configuration space). For other examples of random
quasicrystals, see for example [10, 7].

Let $\mathcal{G}=(V,E)$ be a finite graph and let $T=\\{t_{1},\dots,t_{k}\\}$
be a collection of subgraphs of $\mathcal{G}$, called _tiles_. Let
$\mathbf{N}=\\{N_{v}\\}_{v\in\mathcal{G}}$ be nonnegative integers associated
to vertices of $\mathcal{G}$. Define a new graph $\mathcal{G}_{\mathbf{N}}$,
the “$\mathbf{N}$-fold blowup” of $\mathcal{G}$, to be the graph obtained by
replacing each vertex $v$ of $\mathcal{G}$ with $N_{v}$ vertices, and each
edge $uv$ with the complete bipartite graph $K_{N_{u},N_{v}}$. Now each tile
$t\in T$ can be lifted to a subset of $\mathcal{G}_{\mathbf{N}}$ in many ways:
if $t$ has vertices $v_{1},\dots,v_{k}$ then it has $N_{v_{1}}\cdots
N_{v_{k}}$-many lifts.

We consider _tilings_ of $\mathcal{G}_{\mathbf{N}}$ with lifts of copies of
tiles in $T$ (a tiling is a partition of the vertices of
$\mathcal{G}_{\mathbf{N}}$ into disjoint sets each of which is a lift of a
single tile of $T$). Let $\Omega(\mathbf{N})$ be the set of all tilings; we
call these _$\mathbf{N}$ -fold tilings_.

Let $w:T\to\mathbb{R}_{>0}$ be a positive real weight assigned to each tile.
An $\mathbf{N}$-fold tiling $m$ is assigned a weight $w(m):=\prod_{t\in
T}w_{t}^{m(t)}$ where $m(t)$ is the number of copies of $t$ used. The
_partition function_ for $\mathbf{N}$-fold tilings is defined to be

$Z(w,\mathbf{N})=\sum_{m\in\Omega(\mathbf{N})}w(m).$

We let $\mu=\mu(\mathbf{N})$ be the natural probability measure on
$\Omega(\mathbf{N})$, giving an $\mathbf{N}$-fold tiling a probability
proportional to its weight.

We note that $\mu$ is _not_ the same as the uniform measure on tilings of
$\mathcal{G}$ covering each vertex $N_{v}$ times; each such “multiple tiling”
of $\mathcal{G}$ can be typically lifted to a tiling of
$\mathcal{G}_{\mathbf{N}}$ in many ways.

###  1.1 Results

We compute a generating function for $Z(w,\mathbf{N})$ (Theorem 2.1), and the
asymptotic growth rate of $Z(w,\mathbf{N})$ as $\mathbf{N}\to\infty$, see (5).
This computation involves solving a nonlinear system of equations (4); however
the solution is realized as the unique minimizer of a convex functional
(Theorem 3.1).

In Theorem 5.1 we show that in the $\mathbf{N}\to\infty$ limit the tile
occupation fractions tend to a _Gaussian field_ governed by a discrete
Laplacian operator $\Delta$ on $\mathcal{G}$, the _tiling laplacian_.

For transitive graphs, when we vary the multiplicities, we obtain a _Coulomb
gas_ : defects in multiplicity interact via Coulombic potentials arising from
$\Delta$ (see Section 5.4).

Under certain conditions on transitive graphs, our random multinomial tilings
also undergo a crystallization phenomenon, where the correlations between
distant tiles no longer decay; the tiling freezes into a periodic or
quasiperiodic state. This occurs on $\mathbb{Z}^{2}$, for example, tiled with
translates of the $L$-triomino and a small density of singleton monomers
(Section 6.2.1). As the density of monomers tends to zero the correlation
length of the system tends to infinity, and the system _freezes_. There is a
spontaneous symmetry breaking, since there are three distinct crystalline
states (corresponding to the three distinct—up to translation—periodic tilings
of the plane with $L$ triominos).

For certain other polyominos we get similar freezing phenomena, and others we
don’t; the behavior depends on the presence and type of zeros of the
underlying _characteristic polynomial_ $p(z,u)$ on the unit torus
$\mathbb{T}^{2}\subset\mathbb{C}^{2}$. If the (isolated) zeros on
$\mathbb{T}^{2}$ are sufficiently “generic”, we show that the resulting tiling
will be a quasicrystal (Section 6.2.3). However there is a plethora of
nongeneric behavior for the roots of $p$ as the tile type varies, yielding a
similarly wide variety of behaviors for random tilings (Section 6.3).

Acknowledgments: We thank Wilhelm Schlag and Jim Propp for helpful
conversations. R.K. was supported by NSF DMS-1854272, DMS-1939926 and the
Simons Foundation grant 327929.

##  2 Combinatorics

It is convenient to generalize our definition of tile, to allow the vertices
of a tile to have multiplicity larger than one. The vertices of a tile $t$
then form a _multiset_ of vertices of $\mathcal{G}$, that is, a subset in
which each vertex $v$ has a nonnegative integer multiplicity $t_{v}$. We
identify a tile with its multiset. If we want to think about a tile as a
subgraph, we take all edges of $\mathcal{G}$ connecting vertices which have
positive multiplicity in $t$. A lift of a tile $t$ to
$\mathcal{G}_{\mathbf{N}}$ corresponds to a choice, for each
$v\in\mathcal{G}$, of $t_{v}$ distinct vertices of $\mathcal{G}_{\mathbf{N}}$
lying above $v$, along with the set of all edges of $\mathcal{G}_{\mathbf{N}}$
joining these vertices (that is, the induced subgraph of
$\mathcal{G}_{\mathbf{N}}$ on these vertices). Each lift of $t$ is a blow up
of (the subgraph underlying) $t$.

###  2.1 Generating function

We associate a variable $x_{v}$ to each vertex $v\in\mathcal{G}$. To each tile
$t\in T$ is associated the monomial
$x_{t}=\prod_{v\in\mathcal{G}}\frac{x_{v}^{t_{v}}}{t_{v}!}$. (Here the factor
$t_{v}!$ accounts for the indistinguishability of the vertices of the same
type in a lift of $t$.) Let $P=P(x_{1},\dots,x_{V})$ be the polynomial
$P=\sum_{t\in T}w_{t}x_{t}$. We call $P$ the _tiling polynomial_. The function
$F(X_{1},\dots,X_{V}):=\log P(e^{X_{1}},\dots,e^{X_{V}})$ is called the _free
energy_ (see Section 3.3 below).

######  Theorem 2.1.

Let $\mathrm{x}^{\mathbf{N}}=\prod_{v}x_{v}^{N_{v}}$. Then

$Z(w):=\exp(P)=\sum_{\mathbf{N}\geq
0}Z(w,\mathbf{N})\frac{\mathrm{x}^{\mathbf{N}}}{\mathbf{N}!}$

where the sum is over all vectors of nonnegative multiplicities. If we fix the
total number $K$ of tiles then the corresponding generating function is
$P^{K}/K!$.

###### Proof.

Suppose we use tile $t$ with multiplicity $K_{t}$. Label the abstract copies
of tile $t$ with labels $\ell\in\\{1,\dots,K_{t}\\}$. To place those tiles in
$\mathcal{G}_{\mathbf{N}}$, at each vertex $v$ of $t$, we must choose $K_{t}$
subsets of size $t_{v}$ (one of each label $\ell$) out of the $N_{v}$ vertices
of $\mathcal{G}_{\mathbf{N}}$ lying over $v$. Taking into account all tiles,
this is a multinomial coefficient at vertex $v$: it is
$\frac{N_{v}!}{\prod_{t}(t_{v}!)^{K_{t}}}$ total choices. We take the product
of these over all vertices, and then need to divide by $\prod_{t}K_{t}!$, the
set of choices of initial labellings. In total, the number of tilings with
tile multiplicities $K_{t}$ and vertex multiplicities $\mathbf{N}$ is

$\frac{\prod_{v}N_{v}!}{\prod_{t}K_{t}!\prod_{v,t}(t_{v}!)^{K_{t}}}.$

Multiplying by $\prod_{t}(w_{t}\prod x_{v}^{t_{v}})^{K_{t}},$ the tile weights
(and factors of $x$), dividing by $\mathbf{N}!$, this is

$\prod_{t}\frac{(w_{t}\prod_{v}\frac{x_{v}^{t_{v}}}{t_{v}!})^{K_{t}}}{K_{t}!}$

and summing over the $K_{t}$s gives the result. ∎

###  2.2 Feasible multiplicities

For a given graph $\mathcal{G}=(V,E)$ and tiling set $T$, not all
multiplicities $\mathbf{N}\in(\mathbb{Z}_{\geq 0})^{V}$ are feasible. The set
of feasible multiplicities
$\mathcal{M}_{\mathbb{Z}}=\mathcal{M}_{\mathbb{Z}}(T,\mathcal{G})$ is (just by
definition) the set of nonnegative integer linear combinations of vectors
$v_{t}=\sum_{v\in V}t_{v}e_{v}$, where $\\{e_{v}\\}_{v\in V}$ are the standard
basis vectors for $\mathbb{Z}^{V}$.

In other words we have $\mathcal{M}_{\mathbb{Z}}=D((\mathbb{Z}_{\geq 0})^{T})$
where $D:\mathbb{R}^{T}\to\mathbb{R}^{V}$ is the linear map defined by
$D(e_{t})=\sum_{v\in V}t_{v}e_{v}$. In the standard basis the matrix of $D$
(which we also denote $D$) is called the _incidence matrix_ of the tiling
problem: $D=(D_{v,t})$ where $D_{v,t}=t_{v}$, the multiplicity of $v$ in $t$.
Feasible multiplicities $\mathcal{M}_{\mathbb{Z}}$ are certain integer points
in a real polytopal cone $\mathcal{M}_{\mathbb{R}}\subset\mathbb{R}^{V}$;
$\mathcal{M}_{\mathbb{R}}=D((\mathbb{R}_{+})^{T})$.

Typically not all integer points in $\mathcal{M}_{\mathbb{R}}$ are in
$\mathcal{M}_{\mathbb{Z}}$. For example if all tiles have size $\delta$ then
necessarily $\mathbf{N}$ sums to a multiple of $\delta$. More generally if
$\phi$ is a homomorphism from $\mathbb{Z}^{V}$ to some abelian group, with the
property that $\phi(D(e_{t}))=0$ for all tiles $t$ then $\phi(\mathbf{N})=0$
as well. In the language of tilings these are called “coloring” conditions.

As a typical example of a coloring condition, suppose we wish to tile
$\mathbb{Z}^{2}$ or a subgraph of it with translates of bars of length $3$:
translates of $\\{(0,0),(1,0),(2,0)\\}$ and $\\{(0,0),(0,1),(0,2)\\}$. Let
$\phi:\mathbb{Z}^{\mathbb{Z}^{2}}\to\mathbb{Z}/3\mathbb{Z}$ be defined by
$\phi(e_{(x,y)})=x\bmod 3$. Note that $\phi$ applied to the translate of any
tile is zero:

$\phi(e_{(x,y)}+\phi(e_{(x+1,y)})+\phi(e_{(x+2,y)})=0=\phi(e_{(x,y)})+\phi(e_{(x,y+1)})+\phi(e_{(x,y+2)}).$

We conclude that $\phi(\mathbf{N})=0$ for any feasible multiplicity. This is
equivalent to saying that $\mathbf{N}$ must include an equal number of
vertices in each of the three translates of the sublattice
$\\{(x,y)\in\mathbb{Z}^{2}~{}|~{}x+y\equiv 0\bmod 3\\}$. The same argument
with $\phi_{-}(e_{(x,y)})=y\bmod 3$ gives another linear constraint on
$\mathbf{N}$.

###  2.3 Homology

The incidence map $D:\mathbb{R}^{T}\to\mathbb{R}^{V}$ is generally neither
surjective nor injective. Letting $D^{*}$ be its transpose with respect to the
standard bases, we write
$\mathbb{R}^{V}\cong\operatorname{Im}(D)\oplus\ker(D^{*})$ and
$\mathbb{R}^{T}\cong\operatorname{Im}(D^{*})\oplus\ker(D)$. These are
orthogonal decompositions with respect to the standard inner products. The map
$D$ is an isomorphism from $\operatorname{Im}(D^{*})$ to
$\operatorname{Im}(D)$, and likewise $D^{*}$ is an isomorphism from
$\operatorname{Im}(D)$ to $\operatorname{Im}(D^{*})$

We define
$H_{1}(T,\mathbb{R}):=\mathbb{R}^{V}/\operatorname{Im}(D)\cong\ker(D^{*})$.
Over the integers we define
$H_{1}(T,\mathbb{Z}):=\mathbb{Z}^{V}/\operatorname{Im}(D)$ to be the cokernel
of the map $D$. Colorings $\phi$ are then elements of $H_{1}$, that is, are
functions on vertices which sum to $0$ for each tile:
$D^{*}\phi(t)=\sum_{v}t_{v}\phi(v)=0$.

If all tiles have the same size $\delta$, then $H_{1}(T,\mathbb{Z})$ contains
a copy of $\mathbb{Z}/\delta\mathbb{Z}$; the corresponding coloring functions
are constant functions $f:V\to\mathbb{Z}/\delta\mathbb{Z}$.

For the above example with bars of length $3$, consider tilings of an $n\times
n$ grid, $n\geq 3$. Then $H_{1}(T,\mathbb{R})\equiv\mathbb{R}^{4}$: an element
of $H_{1}(T,\mathbb{R})$ is determined by its values on the lower left
$2\times 2$ square in the grid, which can be arbitrary reals.

The existence of nontrivial integer constraints has an effect on the long-
range behavior of random tilings, see Section 6 below.

###  2.4 Laplacian

The _tiling laplacian_ $\Delta:\mathbb{R}^{V}\to\mathbb{R}^{V}$ is the
operator $\Delta=DCD^{*}$, where $C$ is the diagonal matrix of tile weights
$w_{t}$. It has matrix $\Delta=(\Delta_{u,v})_{u,v\in V}$ with

$\Delta_{u,v}=\sum_{t\in T}w_{t}t_{u}t_{v}.$ (1)

Equivalently, for $f:V\to\mathbb{R}$ we have

$(\Delta f)(v)=\sum_{u}(\sum_{t}w_{t}t_{u}t_{v})f(u).$

The laplacian controls the covariances between tile densities, see Section 5.2
below.

###  2.5 Gauge equivalence

Tile weight functions $w,w^{\prime}$ on $T$ are said to be _gauge equivalent_
if there is a positive function $f:V\to\mathbb{R}_{+}$ such that for all $t\in
T$, $w^{\prime}_{t}=w_{t}\prod_{u\in t}f(u)$. We call $f$ a _gauge
transformation_.

######  Lemma 2.2.

For fixed multiplicities $\mathbf{N}$, gauge equivalent weight functions give
the same probability measure on multinomial tilings.

###### Proof.

Suppose $w^{\prime}$ is gauge equivalent to $w$, that is
$w^{\prime}_{t}=w_{t}\prod_{v}f(v)^{t_{v}}$. An $\mathbf{N}$-fold tiling $m$
for weights $w^{\prime}$ has weight

$\displaystyle\prod_{t}(w^{\prime}_{t})^{m(t)}$ $\displaystyle=$
$\displaystyle\prod_{t}\left(w_{t}^{m(t)}\prod_{v}f(v)^{t_{v}m(t)}\right)$
$\displaystyle=$
$\displaystyle\left(\prod_{t}w_{t}^{m(t)}\right)\prod_{v}f(v)^{\sum_{t}t_{v}m(t)}$
$\displaystyle=$
$\displaystyle\left(\prod_{t}w_{t}^{m(t)}\right)\prod_{v}f(v)^{N_{v}}.$

In particular its weight for $w^{\prime}$ is equal to its weight for $w$
multiplied by a constant independent of $m$. ∎

Note that if $w^{\prime}$ is gauge equivalent to $w$, then the gauge
transformation $f:\mathcal{G}\to\mathbb{R}_{+}$ from $w$ to $w^{\prime}$ may
not be unique: the set of functions $f$ satisfying $\prod_{v\in t}f(v)=1$ for
all $t$ is by definition the kernel of $D^{*}$, written multiplicatively (that
is, $\log f\in H_{1}(T,\mathbb{R})$).

##  3 Asymptotics

In this section compute the asymptotic growth of $Z(w,\mathbf{N})$ as
$\mathbf{N}\to\infty$.

###  3.1 Fixing the number of tiles

Given the multiplicities $\mathbf{N}$, is convenient to also fix the total
number of tiles $K$. If all tiles have the same size $\delta$, then the total
number of tiles $K$ is determined by the multiplicities $\mathbf{N}$: we have
$K=nN/\delta$ where $n=|V|$. More generally we proceed as follows.

We adjoin a new “dummy” vertex $v_{0}$ to $\mathcal{G}$, connected to all
other vertices. Let $\tilde{\mathcal{G}}=\mathcal{G}\cup\\{v_{0}\\}$ be this
new graph. We add to each tile a number of copies of the dummy vertex $v_{0}$
so that all tiles now have the same size $\delta$. Let $x_{0}$ be a variable
associated to the new vertex $v_{0}$, and let $P_{0}$ be the new tiling
polynomial; it is a homogenization of $P$, replacing a monomial $z$ by
$\frac{x_{0}^{m}}{m!}z$, where $m+\deg(z)=\delta$. The number $\delta$ is the
degree of $P_{0}$, and the size of every tile.

Let $N_{0}$ be an arbitrary multiplicity at $v_{0}$, and $M=\sum_{v}N_{v}$ be
the total multiplicity of the other vertices (not including $v_{0}$). For
tileability we need $M+N_{0}$ to be a multiple of $\delta$:

$M+N_{0}=K\delta.$ (2)

Note then that given the remaining multiplicities, the choice of $N_{0}$ is
linearly related to the number of tiles $K$.

We assume for the rest of the paper, unless explicitly stated, that all tiles
have the same size $\delta$. Notationally we can then use $\mathcal{G}$
instead of $\tilde{\mathcal{G}}$ and $P$ instead of $P_{0}$.

###  3.2 Saddle point

For each $v\in V$ let $\alpha_{v}\in\mathbb{R}_{+}$ be fixed. Let
$\vec{\alpha}=(\alpha_{v})_{v\in V}$. We suppose
$\vec{\alpha}\in\mathcal{M}_{\mathbb{R}}(\mathcal{G})$, that is,
$\vec{\alpha}$ is in the cone of feasible multiplicities. Take
$N_{v}\to\infty$ simultaneously for each $v$, in such a way that each
$\mathbf{N}\in\mathcal{M}_{\mathbb{Z}}(\mathcal{G})$, and
$\frac{N_{v}}{K}\to\alpha_{v}$. The quantity $\alpha_{v}$ is the (asymptotic)
fraction of tiles covering $v$, and

$\sum_{v}\alpha_{v}=\delta.$ (3)

From Theorem 2.1 we have

$\frac{K!}{\mathbf{N}!}Z(w,\mathbf{N})=[\mathrm{x}^{\mathbf{N}}]P^{K}.$

We extract the coefficient of $\mathrm{x}^{\mathbf{N}}$ of $P^{K}$ using a
contour integral:

$\frac{K!}{\mathbf{N}!}Z(w,\mathbf{N})=\frac{1}{(2\pi
i)^{V}}\int_{(S^{1})^{V}}\frac{P^{K}}{\prod_{v}x_{v}^{N_{v}}}\prod_{v}\frac{dx_{v}}{x_{v}}.$

For large $K$ we use the saddle-point method. The saddle point is located at
the critical point of the integrand, which is defined by the equations, one
for each $v\in\mathcal{G}$:

$\frac{\partial}{\partial x_{v}}(K\log P-\sum_{v}N_{v}\log x_{v})=0$

or in the large-$K$ limit

$\frac{x_{v}(P)_{x_{v}}}{P}=\alpha_{v}.$ (4)

Solutions to (4) are discussed in Theorem 3.1 below. Although positive
solutions always exist, they are not in general unique; however two positive
solutions differ only by a gauge equivalence in $H_{1}(T,\mathbb{R})$, and as
a consequence give rise to the same weight function and growth rate. (See
Section 3.5 below for an example with nonuniqueness.)

For a solution $\mathrm{x}=\\{x_{v}\\}_{v\in V}$ the growth rate of
$\frac{K!}{\mathbf{N}!}Z(w,\mathbf{N})$ is

$\displaystyle\sigma(w,\\{\alpha_{v}\\})$
$\displaystyle:=\lim_{K\to\infty}\frac{1}{K}\log\frac{K!}{\mathbf{N}!}Z(w,\mathbf{N})$
$\displaystyle=\log P(\mathbf{x})-\sum_{v}\alpha_{v}\log\mathrm{x}_{v}.$ (5)

We call $\sigma(w,\\{\alpha_{v}\\})$ the _exponential growth rate_ of the
multinomial tiling model.

Scaling so that $P=1$, the criticality equations (4) can be written: for all
$v\in V$,

$\sum_{t}w_{t}t_{v}x_{t}=\alpha_{v}.$ (6)

###  3.3 Critical gauge

######  Theorem 3.1.

For any $\vec{\alpha}\in\mathcal{M}_{\mathbb{R}}(\mathcal{G})$ and weight
function $w$ there is a unique gauge equivalent weight function $w^{\prime}$
with the property that for all $v$ the sum of weights of tiles containing
vertex $v$ (counted with multiplicity) is $\alpha_{v}$, that is
$\sum_{v}w^{\prime}_{t}t_{v}=\alpha_{v}$. A corresponding gauge transformation
$f:V\to\mathbb{R}_{>0}$ solves the criticality equations (4) with
$x_{v}=f(v)$.

We call $w^{\prime}$ of this theorem the (weight function in the) _critical
gauge_.

###### Proof.

Define $X_{v}=\log x_{v}$. The free energy $F(X_{1},\dots,X_{V}):=\log
P(e^{X_{1}},\dots,e^{X_{V}})$ is a smooth function of the $X_{i}$’s. It has
gradient lying in $\operatorname{Im}(D)$: its gradient is

$\nabla
F=(\frac{x_{1}P_{x_{1}}}{P},\dots,\frac{x_{T}P_{x_{T}}}{P})=\frac{1}{P}D(\sum_{t}w^{\prime}_{t}e_{t}),$

where $w^{\prime}_{t}=w_{t}\prod_{v\in t}x_{v}$.

Moreover we claim that $F$ is convex. If we interpret $P$ (after scaling so
that $P(1)=1$) as the probability generating function for a
$\mathbb{R}^{V}$-valued random variable $Y$, then the Hessian matrix $H$ of
$F$ is $H_{F}=(\frac{\partial^{2}\log P}{\partial X_{u}\partial
X_{v}})_{u,v\in V}$ is the covariance matrix of $Y$, hence positive
semidefinite. $F$ is strictly convex on directions in $\operatorname{Im}(D)$,
as these are directions where the variance is positive, and $F$ is constant on
directions in $\ker(D^{*})$, that is, those perpendicular to
$\operatorname{Im}(D)$.

Let $S(\vec{\alpha})$ be the Legendre dual of $F$: for
$\vec{\alpha}\in\mathcal{M}_{\mathbb{R}}\subset\operatorname{Im}(D)$, we have

$S(\alpha_{1},\dots,\alpha_{V})=\max_{X_{1},\dots,X_{V}}\left\\{-\log
P(e^{X_{1}},\dots,e^{X_{V}})+\alpha_{1}X_{1}+\dots+\alpha_{V}X_{V}\right\\}.$

Then $S$ is strictly convex and defined on all of
$\mathcal{M}_{\mathbb{R}}\cap\\{\sum_{v}\alpha_{v}=\delta\\}$. We have

$\alpha_{v}=\frac{\partial}{\partial X_{v}}\log P(e^{X_{1}},\dots,e^{X_{V}})$

so that the criticality equations are satisfied with
$\mathrm{x}_{v}=e^{X_{v}}$. Comparing with (5) we see that
$-S(\vec{\alpha})=\sigma(w,\vec{\alpha})$ is the growth rate function.

Here the maximizing $X_{v}$ are unique up to a global additive constant and up
to translations in $\ker D^{*}$; the latter correspond precisely to gauge
transformations not changing the tile weights. The former allow us to scale
all weights so that $P(1)=1$. After this scaling (4) or (6) say precisely that
the sum of weights of tiles containing $v$ (counted with multiplicity) is
$\alpha_{v}$. ∎

Since $\sigma$ is strictly concave, a solution to (4) can be found by
maximizing $\sigma$, written as a function of the $x$’s (even though $\sigma$
is not strictly concave when written as a function of the $x$’s since it is
constant on directions in $H_{1}$, it is strictly concave on orthogonal
directions.)

######  Corollary 3.2.

For the critical gauge $w^{\prime}$, tile probabilities are proportional to
tile weights, that is, the expected number of tiles of type $t$ is
$Kw^{\prime}_{t}=Kw_{t}\mathrm{x}_{t}$, where $K$ is the total number of
tiles.

One consequence of this corollary is that there is, for any choice of tile
probabilities (satisfying the necessary condition of summing to $\alpha_{v}$
at vertex $v$ for each $v$), a choice of tile weights $w_{t}$, unique up to
gauge, for which the multinomial tiling model has those tile probabilities.

###  3.4 Example

Consider tilings of $\mathcal{G}_{1}=\\{1,2,3,4,5\\}\subset\mathbb{Z}$ with
tiles consisting of single vertices and pairs of adjacent vertices: the tiles
are $T=\\{1,2,3,4,5,12,23,34,45\\}.$ We add dummy vertex $v_{0}$ and let
$\mathcal{G}=\mathcal{G}_{1}\cup\\{v_{0}\\}$ where $0$ is connected to all
vertices of $\mathcal{G}_{1}$. Suppose $N_{i}=N$ for $i\neq v_{0}$, and all
tile weights are $1$. Then

$P=x_{0}(x_{1}+x_{2}+x_{3}+x_{4}+x_{5})+x_{1}x_{2}+x_{2}x_{3}+x_{3}x_{4}+x_{4}x_{5}.$

We have $N_{0}+5N=2K$. Let $\alpha=N/K$ and $\alpha_{0}=N_{0}/K$ (note
$\alpha_{0}+5\alpha=2=\delta$).

The feasible range of $\alpha$ is $\alpha\in[\frac{1}{5},\frac{1}{3}]$: when
$\alpha=1/5$, $K=5N$ and we need to use only singleton tiles, and when
$\alpha=1/3$, $K=3N$ and we need to use the maximum proportion of long tiles
(which is two long tiles for every singleton tile); moreover the singleton
tiles must be $x_{1},x_{3}$ or $x_{5}$.

Solving the criticality equations (4) we find the tile probabilities

$\displaystyle x_{0}x_{1}$
$\displaystyle=\frac{1}{4}(1+3\alpha-\sqrt{1-10\alpha+41\alpha^{2}})$
$\displaystyle x_{0}x_{2}$ $\displaystyle=\frac{1}{2}(1-3\alpha)$
$\displaystyle x_{0}x_{3}$
$\displaystyle=\frac{1}{2}(1-7\alpha+\sqrt{1-10\alpha+41\alpha^{2}})$
$\displaystyle x_{1}x_{2}$
$\displaystyle=\frac{1}{4}(-1+\alpha+\sqrt{1-10\alpha+41\alpha^{2}})$
$\displaystyle x_{2}x_{3}$
$\displaystyle=\frac{1}{4}(-1+9\alpha-\sqrt{1-10\alpha+41\alpha^{2}})$

and the remaining probabilities are given by symmetry.

Tile probabilities are plotted in Figure 1.

Figure 1: Tile probabilities for $\alpha\in[1/5,1/3]$. $x_{0}x_{1}$ in blue,
$x_{0}x_{2}$ in orange, $x_{0}x_{3}$ green, $x_{1}x_{2}$ red, $x_{2}x_{3}$
purple.

###  3.5 Example

Here is an example with nontrivial homology. Consider tilings of a cycle of
length $4$: $\\{1,2,3,4\\}$ with dimers $\\{12,23,34,41\\}$. Then
$P=x_{1}x_{2}+x_{2}x_{3}+x_{3}x_{4}+x_{4}x_{1}=(x_{1}+x_{3})(x_{2}+x_{4})$.
The space $\operatorname{Im}(D)\subset\mathbb{R}^{V}$ is the orthocomplement
of the vector $(1,-1,1,-1)$, so $H_{1}(T,\mathbb{R})$ has rank $1$ and is
generated by this vector. The feasible $\vec{\alpha}$ are those which satisfy
$\sum_{v}\alpha_{v}=2$ and are in $\operatorname{Im}(D)$, that is, satisfy the
conditions $\alpha_{1}+\alpha_{3}=1=\alpha_{2}+\alpha_{4}$. The criticality
equations are

$\frac{x_{1}x_{2}+x_{1}x_{4}}{P}=\alpha_{1},\frac{x_{1}x_{2}+x_{2}x_{3}}{P}=\alpha_{2},\frac{x_{2}x_{3}+x_{3}x_{4}}{P}=\alpha_{3},\frac{x_{1}x_{4}+x_{3}x_{4}}{P}=\alpha_{4},$

which reduce to

$\frac{x_{1}}{x_{1}+x_{3}}=\alpha_{1},~{}~{}\frac{x_{2}}{x_{2}+x_{4}}=\alpha_{2}.$

Solutions are not unique: given any solution we can multiply $x_{1},x_{3}$ by
a constant $t$ and divide $x_{2},x_{4}$ by $t$ to get another solution. We
have

$F(X_{1},\dots,X_{4})=\log((e^{X_{1}}+e^{X_{3}})(e^{X_{2}}+e^{X_{4}})).$

This leads to the growth rate

$\displaystyle\sigma(\vec{\alpha})$
$\displaystyle=-\alpha_{1}\log\alpha_{1}-\alpha_{2}\log\alpha_{2}-\alpha_{3}\log\alpha_{3}-\alpha_{4}\log\alpha_{4}$
$\displaystyle=-\alpha_{1}\log\alpha_{1}-\alpha_{2}\log\alpha_{2}-(1-\alpha_{1})\log(1-\alpha_{1})-(1-\alpha_{2})\log(1-\alpha_{2})$
$\displaystyle=h(\alpha_{1})+h(\alpha_{2})$

where $h(p)$ is the Shannon entropy $h(p)=-p\log p-(1-p)\log(1-p).$

##  4 Dimers

A special case of the multinomial tiling model is the _multinomial dimer
model_ , where tiles are simply all pairs of adjacent vertices (also known as
“dimers”). A $1$-dimer tiling is then a perfect matching, also known as _dimer
cover_ of $\mathcal{G}$.

###  4.1 Bipartite graphs

For the dimer model, when $\mathcal{G}$ is bipartite, $V=B\cup W$, there is an
equivalent but perhaps more efficient method of computing $Z(w,\mathbf{N})$.
By Corollary 3.2 we need to look for a gauge function $x:V\to\mathbb{R}_{>0}$
such that

$\alpha_{\mathrm{b}}=\sum_{\mathrm{w}\sim\mathrm{b}}w_{\mathrm{w}\mathrm{b}}x_{\mathrm{w}}x_{\mathrm{b}}$
(7)

and likewise for white vertices. However we can use (7) to define
$x_{\mathrm{b}}$:

$x_{\mathrm{b}}=\frac{\alpha_{\mathrm{b}}}{\sum_{\mathrm{w}\sim\mathrm{b}}w_{\mathrm{w}\mathrm{b}}x_{\mathrm{w}}}.$

Then we only have equations involving the remaining half the variables: those
at the white vertices $\mathrm{x}_{\mathrm{w}}$. These equations are

$\alpha_{\mathrm{w}}=\sum_{\mathrm{b}\sim\mathrm{w}}\alpha_{\mathrm{b}}\frac{w_{\mathrm{w}\mathrm{b}}x_{\mathrm{w}}}{\sum_{\mathrm{w}\sim\mathrm{b}}w_{\mathrm{w}\mathrm{b}}x_{\mathrm{w}}}.$
(8)

In the standard case where $N_{v}\equiv N$, all the
$\alpha_{\mathrm{w}},\alpha_{\mathrm{b}}$ are equal and the saddle point
equations correspond to the property that the _sum of (critical) edge weights
at each vertex is $1$_. This is just a restatement of Cororollary 3.2 in this
setting, since the sum of edge probabilities at each vertex is $1$ for a
random dimer cover.

###  4.2 Aztec Diamond Example

The _Aztec diamond_ of order $n$ is a diamond-shaped subregion of
$\mathbb{Z}^{2}$ of horizontal diameter $2n-1$; see Figure 2, left panel for
the $n=4$ Aztec diamond. It is known to have $2^{n(n+1)/2}$ single-dimer
covers (see [8] and [9]). Consider $\mathbf{N}$-fold dimer covers with
$N_{v}\equiv N$ and $w_{t}\equiv 1$. The critical edge weights sum to $1$ at
each vertex, and have the property that around each square face the weights
$a,b,c,d$ satisfy $ac=bd$. The critical weights for $n=4$ are shown on the
left, and for general $n$ (scaled by $n(n-1)$) on the right in Figure 2.

Figure 2: Aztec diamond of order $4$ (left) in critical gauge. Critical gauge
edge weights for general $k$, multiplied by $k(k+1)$ (right). The weight on an
edge is a quadratic function of its $x,y$ coordinates and its parity. Out of a
white vertex at $x,y$ the edge weights E,N,W,S are $1/k(k+1)$ times
respectively $(x+y)(x-y),(x+y)(k+1-x+y),(k-x-y)(k+1-x+y),(k-x-y)(x-y)$.

We can work out the growth rate $\sigma$ in this case as follows. Since $P=1$,
and $\alpha_{v}=\frac{1}{k(k+1)}$ is a constant, we have

$\sigma=-\frac{1}{k(k+1)}\sum_{v}\log x_{v}=-\frac{1}{k(k+1)}\log\prod x_{v}.$

This product is the weight of any single dimer cover. The “all horizontal”
dimer cover has dimers of weight $\frac{1}{k(k+1)}$ times: $k^{2}$ for the top
row, $k(k-1)$ and $(k-1)k$ for the next row, and generally
$k(k-i+1),(k-1)(k-i+2),\dots,(k-i+1)k$ for the $i$ row, for $i$ from $1$ to
$k$, then repeating for the bottom half of the diamond. The total product of
edge weights of the dimer cover is

$\frac{(k^{k}\cdot(k-1)^{k-1}\cdots 2^{2}\cdot 1)^{4}}{(k(k+1))^{k(k+1)}}$

This yields for the exponential growth rate the remarkable value
$\sigma=1+O(\frac{\log k}{k^{2}})$.

Associated to a multinomial dimer cover with constant $N_{v}\equiv N$ of a
subgraph of $\mathbb{Z}^{2}$ is a _height function_ $h$ on the dual graph. The
height function is defined to be zero on a fixed face, and the change in
height across an edge $\mathrm{w}\mathrm{b}$ (when crossing the edge so that
the white vertex is on the left) is $-N/4$ plus the number of dimers on that
edge. In [3], see also [4], the authors prove a limit shape phenomenon for
single dimer covers: the (rescaled by $n$) height function for a random dimer
cover of $AD(n)$ converges with probability one as $n\to\infty$ to a nonrandom
piecewise analytic function on the rescaled diamond $|x|+|y|\leq 1$. For the
$N$-fold dimer cover discussed above we get a similar, but analytic, limit
shape. It is just the function $h(x,y)=x^{2}-y^{2}$. The general limit shape
phenomenon for multinomial tilings is discussed in [17].

###  4.3 Path example

For fixed $n>0$ take $\mathcal{G}$ to be the (bipartite) $n\times n$ honeycomb
graph of Figure 3, with edge weights $1$.

Figure 3: $3\times 3$ honeycomb graph.

There are $\binom{2n}{n}$ dimer covers: dimer covers correspond bijectively to
monotone lattice paths from $(0,0)$ to $(n,n)$. The bijection is obtained by
taking a dimer cover and shrinking all horizontal edges of $\mathcal{G}$ to
points.

Let us consider $\mathbf{N}$-dimer covers of $\mathcal{G}_{n}$, where
$N_{v}\equiv N$ and $w_{t}=1$. Index the white vertices $x_{i,j}$ as in the
figure. The criticality equations (8) are

$\frac{x_{i,j}}{x_{i,j}+x_{i-1,j}+x_{i,j-1}}+\frac{x_{i,j}}{x_{i,j}+x_{i+1,j}+x_{i+1,j-1}}+\frac{x_{i,j}}{x_{i,j}+x_{i,j+1}+x_{i-1,j+1}}=1$
(9)

and boundary conditions $x_{i,j}=0$ for $i<0$ or $j<0$ or $(i,j)=(n,n)$.

In the limit $n\to\infty$, there is a solution to (9) given by

$x_{i,j}=(i+j)!\binom{i+j}{i}.$

However we don’t know if this solution is the only one (since the graph is
infinite, unicity does not necessarily hold).

The edge probabilities for this solution are (for horizontal, NE, SE edges
respectively out of a white vertex $(i,j)$)

$\frac{i+j}{i+j+1},~{}~{}~{}\frac{j+1}{(i+j+1)(i+j+2)},~{}~{}~{}\frac{i+1}{(i+j+1)(i+j+2)}$

for the (horizontal, resp. NE, resp. SE) edge at $(i,j)$.

These edge probabilities give a unit flow on $\mathbb{N}\times\mathbb{N}$ from
$(0,0)$ to $\infty$, see Figure 4 left panel; the value on an edge represents
the flow from left to right along that edge.

Figure 4: Edge probabilities (left). Corresponding random walk probabilities
(right).

The value on an edge is the probability that a certain monotone random walk
uses that edge. The transition probabilities of this random walk are shown on
the right in Figure 4; this random walk is the _Polya urn_ 222An urn starts
out with one red and one green ball. A ball is selected at random and replaced
with another ball of the same color. This process is then repeated many times.
The resulting distribution of the number of red balls after $k$ steps is
uniform on $[1,k]$..

###  4.4 Example in higher dimensions

For a higher dimensional multinomial dimer example, consider the infinite
subgraph of $\mathbb{Z}^{3}$ in the slab $0\leq x+y+z\leq n$ where $n$ is
even, and $w\equiv 1$. To get a finite graph we can quotient by a cofinite
sublattice of $\\{x+y+z=0\\}$. Such a graph has a higher proportion of white
vertices than black vertices (assuming the origin is white); the density ratio
is $(n+1)/n$. Take $N_{\mathrm{w}}$ constant, and $N_{\mathrm{b}}$ a different
constant with $N_{\mathrm{w}}/N_{\mathrm{b}}=n/(n+1)$. Then a critical gauge
is given up to scale by: for $\mathrm{w}=(x,y,z),$
$x_{\mathrm{w}\mathrm{b}}=n-x-y-z$ if $b=\mathrm{w}+e_{i}$ and
$x_{\mathrm{w}\mathrm{b}}=x+y+z$ if $\mathrm{b}=\mathrm{w}-e_{i}$ (here
$e_{1},e_{2},e_{3}$ are the standard basis vectors).

There are analogous examples in all dimensions $d\geq 1$.

##  5 Multiplicities and fluctuations

###  5.1 Changing multiplicities

We compute the change in growth rate $\sigma$ (from equation (5)) under a
small change in the multiplicities $\alpha_{v}\to\alpha_{v}+d\alpha_{v}$. This
will be used below to compute tile covariances. Since
$\sum_{v}\alpha_{v}=\delta$, the sum of changes is necessarily zero:
$\sum_{v}d\alpha_{v}=0$.

Recall the incidence matrix $D=(D_{v,t})_{v\in V,t\in T}$, defined by
$D_{v,t}=t_{v}$. Differentiating (6) we find for each vertex $u$:

$\sum_{t}w_{t}t_{u}x_{t}(\sum_{v}t_{v}\frac{dx_{v}}{x_{v}})=d\alpha_{u},$

or

$\sum_{v,t}D_{u,t}w_{t}x_{t}D_{t,v}\frac{dx_{v}}{x_{v}}=d\alpha_{u}.$

Thus

$\sum_{v}\Delta_{u,v}\frac{dx_{v}}{x_{v}}=d\alpha_{u},$ (10)

where $\Delta=DCD^{*}$ is the tiling laplacian for the critical weights.
Equation (10) says that as a function of $v$, $\frac{dx_{v}}{x_{v}}$ is
harmonic with respect to the laplacian $\Delta$ at all vertices $u$ for which
$d\alpha_{u}=0$.

Equation (10) will have a solution if and only if $d\vec{\alpha}$ is in the
image of $\Delta$, which is the same as $\operatorname{Im}(D)$, since the
laplacian is invertible on $\operatorname{Im}(D)$, mapping it to itself (and
$\Delta$ is zero on $\ker D^{*}$). The solution is unique up to an element of
$\ker\Delta=\ker D^{*}=H_{1}(T,\mathbb{R})$; as discussed in Section 2.5 these
correspond to gauge transformations not changing the tile weights.

###  5.2 Covariance of tile densities

Let $X_{t}$ be the random variable counting the number of occurrences of tile
$t$ in an $\mathbf{N}$-fold tiling. We wish to compute the covariance
$\mathrm{Cov}(X_{t},X_{t^{\prime}})=\mathbb{E}[X_{t}X_{t^{\prime}}]-\mathbb{E}[X_{t}]\mathbb{E}[X_{t^{\prime}}]$
for two tiles $t,t^{\prime}$.

Note that $X_{t}$ is itself a sum of $\\{0,1\\}$-valued random variables,
$X_{t}=\sum_{i=1}^{M_{t}}X_{t}^{i}$, where the sum runs over all
$M_{t}=\prod_{v}\binom{N_{v}}{t_{v}}$ possible lifts of the tile $t$ to
$\mathcal{G}_{\mathbf{N}}$. It suffices to compute the covariance
$\mathrm{Cov}(X_{t}^{i},X_{t^{\prime}}^{j}).$

####  5.2.1 The case $t\neq t^{\prime}$.

Assume first that $t\neq t^{\prime}$. By the symmetry of
$\mathcal{G}_{\mathbf{N}}$, if $t$ and $t^{\prime}$ are disjoint this is
independent of $i$ and $j$:
$\mathrm{Cov}(X_{t}^{i},X_{t^{\prime}}^{j})=\mathrm{Cov}(X_{t}^{1},X_{t^{\prime}}^{1})$.
(If $t,t^{\prime}$ overlap, see below.) We have

$\displaystyle\mathbb{E}[X_{t}^{1}X_{t^{\prime}}^{1}]$
$\displaystyle=\Pr(X_{t}^{1}=1,X_{t^{\prime}}^{1}=1)$
$\displaystyle=\Pr(X_{t^{\prime}}^{1}=1|X_{t}^{1}=1)\Pr(X_{t}^{1}=1)$
$\displaystyle=\Pr\\!{}^{*}(X_{t^{\prime}}^{1}=1)\Pr(X_{t}^{1}=1),$

where the star denotes the probability measure on the graph
$\mathcal{G}_{\mathbf{N}^{*}}$ where we have reduced the multiplicities of
vertices $v$ by $t_{v}$. We thus have

$\mathbb{E}[X_{t}X_{t^{\prime}}]=M_{t}M_{t^{\prime}}\mathbb{E}[X_{t}^{1}X_{t^{\prime}}^{1}]=\mathbb{E}[X_{t}]\frac{M_{t^{\prime}}}{M^{*}_{t^{\prime}}}\mathbb{E}^{*}[X_{t^{\prime}}]=\mathbb{E}[X_{t}]\mathbb{E}^{*}[X_{t^{\prime}}],$

since $M_{t^{\prime}}=M^{*}_{t^{\prime}}$ when $t,t^{\prime}$ are disjoint.

When $t$ and $t^{\prime}$ overlap, the number of disjoint lifts of $t$ and
$t^{\prime}$ is

$M_{t,t^{\prime}}:=\prod_{v}\binom{N_{v}}{t_{v},t^{\prime}_{v}}=\prod_{v}\binom{N_{v}}{t_{v}}\binom{N_{v}-t_{v}}{t^{\prime}_{v}}=M_{t}\prod_{v}\binom{N_{v}-t_{v}}{t^{\prime}_{v}}.$

The calculation follows as before but we need to sum over the
$M_{t,t^{\prime}}$ lifts. We find

$\mathbb{E}[X_{t}X_{t^{\prime}}]=M_{t,t^{\prime}}\mathbb{E}[X_{t}^{1}X_{t^{\prime}}^{2}]=\mathbb{E}[X_{t}]\mathbb{E}^{*}[X_{t^{\prime}}],$

since the number of lifts of $t^{\prime}$ in $\mathcal{G}_{\mathbf{N}^{*}}$ is
exactly

$\prod_{v}\binom{N_{v}-t_{v}}{t^{\prime}_{v}}=\frac{M_{t,t^{\prime}}}{M_{t}}.$

Thus in either case, if $t\neq t^{\prime}$,

$\mathrm{Cov}(X_{t},X_{t^{\prime}})=\mathbb{E}[X_{t}](\mathbb{E}^{*}[X_{t^{\prime}}]-\mathbb{E}[X_{t^{\prime}}]).$

####  5.2.2 The case $t=t^{\prime}$.

Finally, if $t=t^{\prime}$ are the same tile, then we have a slightly
different computation. Let $X_{t,1}$ and $X_{t,2}$ correspond to disjoint
lifts of $t$. With $p=\mathbb{E}[X_{t,1}]$ and
$p^{*}=\mathbb{E}^{*}[X_{t,2}]$,

$\displaystyle\mathbb{E}[X_{t}^{2}]-\mathbb{E}[X_{t}]^{2}$
$\displaystyle=M_{t}\mathbb{E}[X_{t,1}]+M_{t}(M_{t}-1)\mathbb{E}[X_{t,1}X_{t,2}]-M_{t}^{2}\mathbb{E}[X_{t,1}]^{2}$
$\displaystyle=M_{t}p+M_{t}(M_{t}-1)pp^{*}-(M_{t}p)^{2}$
$\displaystyle=\mathbb{E}[X_{t}](1-p+\frac{M_{t}-1}{M_{t}}(\mathbb{E}^{*}[X_{t}]-\mathbb{E}[X_{t}])).$

####  5.2.3 Computation of $\mathbb{E}^{*}$.

We can now compute $\mathbb{E}^{*}[X_{t^{\prime}}]$ from the methods of
section 5.1. We need to change $\alpha_{v}=\frac{N_{v}}{K}$ to
$\frac{N_{v}-t_{v}}{K-1}$. Thus

$d\alpha_{v}=\frac{\alpha_{v}}{K-1}-\frac{t_{v}}{K-1}.$

Recalling $\mathbb{E}[X_{t^{\prime}}]=Kw_{t^{\prime}}x_{t^{\prime}},$ we have
(ignoring lower order terms)

$\displaystyle\mathbb{E}^{*}[X_{t^{\prime}}]-\mathbb{E}[X_{t^{\prime}}]$
$\displaystyle=(K-1)w_{t^{\prime}}x_{t^{\prime}}(1+\sum_{u}t^{\prime}_{u}\frac{dx_{u}}{x_{u}})-Kw_{t^{\prime}}x_{t^{\prime}}$
$\displaystyle=w_{t^{\prime}}x_{t^{\prime}}(-1+(K-1)\sum_{u}t^{\prime}_{u}\frac{dx_{u}}{x_{u}})$
$\displaystyle=w_{t^{\prime}}x_{t^{\prime}}(-1+(K-1)\sum_{u,v}D_{t^{\prime},u}\Delta_{u,v}^{-1}d\alpha_{v})$
$\displaystyle=w_{t^{\prime}}x_{t^{\prime}}(-1+\sum_{u,v}D_{t^{\prime},u}\Delta_{u,v}^{-1}(\alpha_{v}-D_{v,t}))$
$\displaystyle=w_{t^{\prime}}x_{t^{\prime}}(-1-(D^{*}\Delta^{-1}D)_{t^{\prime},t}+\sum_{u,v}D_{t^{\prime},u}\Delta^{-1}_{u,v}\alpha_{v}).$
(11)

Now note that

$\sum_{v}\Delta_{u,v}=\sum_{v,t}D_{u,t}w_{t}x_{t}D_{t,v}=\delta\sum_{t}D_{u,t}w_{t}x_{t}=\delta\alpha_{u},$

so $\Delta^{-1}_{u,v}\alpha_{v}=\frac{1}{\delta}1_{u}$ and thus
$\sum_{u,v}D_{t^{\prime},u}\Delta^{-1}_{u,v}\alpha_{v}=1$. The last sum in
(11) cancels the $-1$ and we have

$\mathbb{E}^{*}[X_{t^{\prime}}]-\mathbb{E}[X_{t^{\prime}}]=-w_{t^{\prime}}x_{t^{\prime}}\mathbb{K}_{t^{\prime},t},$

where $\mathbb{K}=D^{*}\Delta^{-1}D$, and so for $t\neq t^{\prime}$ at
critical gauge

$\mathrm{Cov}(X_{t},X_{t^{\prime}})=-Kw_{t}w_{t}^{\prime}\mathbb{K}_{t^{\prime},t}.$

Likewise when $t=t^{\prime}$,

$\mathbf{Var}(X_{t})=Kw_{t}(1-p-w_{t}\mathbb{K}_{t,t}).$

As long as $\delta>1$, we can ignore the $p$ for large $K$, and write

$\mathrm{Cov}(X_{t},X_{t^{\prime}})=Kw_{t}(I-C\mathbb{K})_{t^{\prime},t}$ (12)

where $I$ is the identity matrix. Note that $C\mathbb{K}=(C\mathbb{K})^{2}$ is
a projection matrix from $\mathbb{R}^{T}$ onto the subspace
$C\operatorname{Im}(D^{*})$, with kernel $\ker(D)$. Thus $I-C\mathbb{K}$ is
the complementary projection.

###  5.3 Gaussian fluctuations

For the multinomial tiling model with multiplicities $\mathbf{N}$, as before
let $X_{t}$ be the random variable representing the multiplicity of tile $t$.
We consider a limit $K\to\infty$ as in Section 3. Scaling $P$ to $1$, we can
consider $P$ to be the probability generating function (pgf) of a single tile.
Then $P^{K}$ is the pgf of placing $K$ i.i.d. tiles. The tile multiplicities
under this process are Poisson$(Kw_{t})$ random variables which tend in the
limit of large $K$ to Gaussian random variables which are independent except
for the constraint that their sum is $K$. When we impose the constraints on
the multiplicities $N_{v}$ this is an additional linear constraint (which
implies the first); the resulting random variable is thus also a joint
Gaussian. A (multidimensional) Gaussian is determined by its covariance
matrix. The covariance can in general be obtained from the original covariance
matrix and the constraint matrix. However in our case we have already computed
the covariance matrix above in (12).

######  Theorem 5.1.

In the limit of large $K$ the joint distribution of the $X_{t}$ tends to a
(multidimensional) Gaussian with mean $\mathbb{E}[\vec{X}]=K\vec{w}$ and
covariance matrix $\mathrm{Cov}(X_{s},X_{t})=Kw_{s}(I-C\mathbb{K})_{t,s}.$

Examples are explored in Section 6.

###  5.4 Coulomb gas

In this section we assume for simplicity that $\mathcal{G}$ is regular, $T$
and $w$ are symmetric under a transitive group of automorphisms of
$\mathcal{G}$, and $N_{v}\equiv N$. We also assume $t_{v}\in\\{0,1\\}$ for all
tiles and vertices. We have $K=nN/\delta$. Under these conditions the critical
weights are $x_{v}\equiv x$ where $x$ is a constant.

Suppose that we take a small perturbation of the multiplicities
$N_{v}=N+q_{v}$ where $\frac{q_{v}}{N}=o(1)$. We call $q_{v}$ the _$d$
-charge_ at $v$. It can be positive or negative; we’ll take the sum of
$d$-charges to be zero. Let us compute $\sigma$ as a function of the $q_{v}$.

We consider perturbations of $\sigma$ with respect to the multiplicities
$\alpha_{v}$. Each $\alpha_{v}$ is of the form
$\alpha_{v}=\delta/n+d\alpha_{v}$ where $d\alpha_{v}=\frac{\delta q_{v}}{nN}$.

To first order from (5) we have

$d\sigma(\vec{\alpha})=\sum_{v}\frac{P_{x_{v}}}{P}dx_{v}-\sum_{v}\alpha_{v}\frac{dx_{v}}{x_{v}}-\sum_{v}\log
x_{v}d\alpha_{v}=-\sum_{v}\log x_{v}d\alpha_{v}=0$

since $\log x_{v}$ is constant. To second order we have

$\sigma(\vec{\alpha}+d\vec{\alpha})-\sigma(\vec{\alpha})=$
$\sum_{v}(\frac{P_{x_{v}x_{v}}}{P}-\frac{P_{x_{v}}^{2}}{P^{2}})\frac{dx_{v}^{2}}{2}+\sum_{u\neq
v}(\frac{P_{x_{u}x_{v}}}{P}-\frac{P_{x_{u}}P_{x_{v}}}{P^{2}})dx_{u}dx_{v}-\sum_{v}\frac{d\alpha_{v}dx_{v}}{x_{v}}+\sum_{v}\frac{\alpha_{v}}{x_{v}^{2}}\frac{dx_{v}^{2}}{2}.$

Substituting $x_{v}=x,P_{x_{v}}=\alpha_{v}/x,P=1,P_{x_{v}x_{v}}=0$ and
$\alpha_{v}=\delta/n$ gives

$=\sum_{v}(\frac{\delta}{n}-\frac{\delta^{2}}{n^{2}})\frac{dx_{v}^{2}}{2x^{2}}+\sum_{u\neq
v}(\sum_{t}w_{t}t_{u}t_{v}-\frac{\delta^{2}}{n^{2}})\frac{dx_{u}}{x}\frac{dx_{v}}{x}-\sum_{v}d\alpha_{v}\frac{dx_{v}}{x}$

which we can write as

$=\frac{d\vec{x}^{t}}{x}\left(\frac{1}{2}\Delta-\frac{\delta^{2}}{2n^{2}}J\right)\frac{d\vec{x}}{x}-d\vec{\alpha}^{t}\cdot\frac{d\vec{x}}{x}$

where $J$ is the all-$1$’s matrix. But

$\frac{d\vec{x}^{t}}{x}J\frac{d\vec{x}}{x}=d\vec{\alpha}\Delta^{-1}J\Delta^{-1}d\vec{\alpha}=0$

(since $\Delta^{-1}$ has constant row and column sums,
$\Delta^{-1}J\Delta^{-1}$ is a multiple of $J$). We are left with

######  Theorem 5.2.

Under the above conditions on $\mathcal{G},\mathbf{N},T,w$ the second
derivative of $\sigma$ in direction $\vec{\alpha}$ (of sum zero) is

$d^{2}\sigma=-\frac{1}{2}d\vec{\alpha}\Delta^{-1}d\vec{\alpha}.$

If we fix the $d$-charges but allow them to move from vertex to vertex,
considering $\sigma$ as a function of position of these charges, we see that
$d$-charges interact with a potential defined by $\Delta^{-1}$. This is the
_Coulomb potential_ associated with $\Delta$.

####  5.4.1 Example: dimer case

Consider for example the case in which $\mathcal{G}$ is bipartite, and we have
the multinomial dimer model. Define the _charge_ $\tilde{q}_{v}$ of a vertex
of d-charge $q_{v}$ to be $\tilde{q}_{v}=d(v)q_{v}$ where $d(v)=1$ for a white
vertex and $d(v)=-1$ for a black vertex. Then note that
$\tilde{\Delta}_{u,v}=d(u)\Delta_{u,v}d(v)$ is the _standard_ graph laplacian.

Suppose we fix the charges $\tilde{q}_{v}$, but not their locations. We can
interpret $\sigma$ as a function of position as a potential energy which
causes like charges to repel and opposite charges to attract. That is,
$\sigma$ is larger when like charges are farther apart and opposite charges
are closer. The _force_ of repulsion/attraction is naturally given by the
gradient of the Green’s function $\tilde{\Delta}^{-1}$ (that is, the gradient
of the potential). For $\mathcal{G}=\mathbb{Z}^{2}$, for example this is a
“$1/r$” force for distant particles, where $r$ is the vector between them. For
$\mathcal{G}=\mathbb{Z}^{3}$ it is a “$1/r^{2}$” force for distant particles.
These conclusions are consistent with standard electrostatics in
$\mathbb{R}^{2}$ and $\mathbb{R}^{3}$, when the charges are far apart
(compared to the lattice spacing). In $2d$ these results also agree with
corresponding results obtained for the single dimer model by Ciucu [5].

##  6 Crystallization

In this section we restrict our graphs $\mathcal{G}$ to be subgraphs of
$\mathbb{Z}^{d}$, and we consider tilings with tiles $T$ which are translates
of one or more “prototiles”. In the simplest case we have only two prototiles
$\mathrm{t},\mathrm{t}_{0}$, where $\mathrm{t}_{0}$ is a _single vertex_. Then
a $\mathbf{N}$-fold tiling of $\mathcal{G}$ with $T$ is a tiling of
$\mathcal{G}_{\mathbf{N}}$ with lifts of translates of $\mathrm{t}$ which has
a number of _holes_ (which are locations which are covered by lifts of
translates of $\mathrm{t}_{0}$). We are interested in what happens when the
fractional density of holes goes to zero. Depending on the shape of
$\mathrm{t}$, the system will _crystallize_ (Theorem 6.1 below).

Before giving a general result, we work out some explicit examples, which are
of interest in their own right, and which illustrate the general situation.
First in dimension one we consider the case of triominos with no holes
(Section 6.1.1) and then with a positive fraction of holes (Section 6.1.2),
and then the $L$ triomino in $\mathbb{Z}^{2}$ (Section 6.2.1).

###  6.1 Examples in one dimension

####  6.1.1 Example: bars of length $3$

Consider tilings of the cycle $\mathcal{G}=\mathbb{Z}/n\mathbb{Z}$ with
translates of the triomino $\mathrm{t}=\\{0,1,2\\}$, that is
$T=\\{t_{0},\dots,t_{n-1}\\}$ with $t_{i}=\\{i,i+1,i+2\\}$ with cyclic
indices. We choose constant multiplicities $N_{v}\equiv N$. Then the total
number of tiles is $K=Nn/3$ and the fraction of tiles per vertex is
$\alpha=3/n$. Since the graph is regular, for the critical gauge we have
$x_{v}\equiv x,$ and $P=nx^{3}=1$, so the weight per tile is
$w_{t}=x^{3}=1/n$. The tile laplacian $\Delta$ satisfies, for a function
$f\in\mathbb{R}^{V}$,

$(\Delta f)(j)=\frac{1}{n}(3f(j)+2f(j-1)+2f(j+1)+f(j+2)+f(j-2))$

with cyclic indices.

It is convenient to use the Fourier transform to invert $\Delta$. For $z$ an
$n$th root of $1$ let $U_{z}$ be the subspace of $\mathbb{C}^{V}$ consisting
of $z$-periodic functions

$U_{z}=\\{f\in\mathbb{C}^{V}~{}|~{}f(x+1)=zf(x)\\}.$

Here $U_{z}$ is the eigenspace for translation by $1$ on $\mathcal{G}$ with
eigenvalue $z$. The operator $\Delta$ preserves each $U_{z}$ and its action on
$U_{z}$ is multiplication by

$\lambda_{z}=\frac{1}{n}(1+z+z^{2})(1+z^{-1}+z^{-2}).$

The matrix $D$ can likewise be written in a Fourier basis, if we identify a
tile $t_{i}$ with the location of its left endpoint $i$. The action of $D$ on
$U_{z}$ is then multiplication by $1+z+z^{2}$, and that of $D^{*}$ is
multiplication by $1+1/z+1/z^{2}$. We see that $\mathbb{K}=D^{*}\Delta^{-1}D$
is simply multiplication by $n$ on subspaces $U_{z}$ for which $z^{2}+z+1\neq
0$, and $0$ on subspaces $U_{z}$ for which $z^{2}+z+1=0.$

Suppose that $n$ is not a multiple of $3$. Then $\mathbb{K}$ is a scalar,
equal to multiplication by $n$. The covariance matrix is identically zero:
$\mathrm{Cov}(X_{t},X_{t^{\prime}})=0.$ This is not surprising since the tile
multiplicities $X_{i}$ are not random: we necessarily have $X_{i}=N/3$ for all
$i$.

Suppose that $n$ is a multiple of $3$. Then $\frac{1}{n}\mathbb{K}$ is a
projection matrix $\mathbb{P}_{3}$ onto the span of the subspaces $U_{z}$
where $z^{2}+z+1\neq 0$. The covariance matrix is
$\mathrm{Cov}(X_{t},X_{t}^{\prime})=\frac{K}{n}(I-\mathbb{P}_{3})$, where
$I-\mathbb{P}_{3}$ is the projection onto the (two-dimensional) span of the
$U_{z}$ where $z^{2}+z+1=0$. That is, in the standard basis on
$\mathbb{R}^{V}$,

$\mathrm{Cov}=\frac{K}{n^{2}}\begin{pmatrix}2&-1&-1&2\\\ -1&2&-1&-1\\\
-1&-1&2&-1\\\ 2&-1&-1&2\\\ &&&&\ddots\end{pmatrix}.$ (13)

Pairs of tiles $t,t^{\prime}$ at distance a multiple of $3$ are perfectly
correlated: we have a crystal. This is also not surprising since the
multiplicity at $j$ determines $X_{j}$ as a function of $X_{j-1},X_{j-2}$: we
have $X_{j-2}+X_{j-1}+X_{j}=N$, which implies $X_{j}=X_{j+3}$ for all $j$.

####  6.1.2 Example: bars of length $3$ and singletons

Let us consider a variant of the above model, where we also allow singleton
tiles, with weight $1$. As in Section 3.1 we include a dummy vertex $v_{0}$ in
$\mathcal{G}$; we can use the multiplicity $N_{0}$ of the dummy vertex to
control the number of singleton tiles. We have $N_{0}+nN=3K$, and
$\alpha_{0}=N_{0}/K$ and $\alpha=N/K$ satisfy $\alpha_{0}+n\alpha=3$. Also

$P_{0}=\frac{x_{0}^{2}}{2}\sum x_{i}+\sum
x_{i}x_{i+1}x_{i+2}=\frac{n}{2}x_{0}^{2}x+nx^{3}$

where we replaced $x_{i}$ with $x$ by circular symmetry. The critical weights
$w_{0},w$ are $w_{0}=\frac{\alpha_{0}}{2n}$ and $w=\frac{2-\alpha_{0}}{2n}.$

In this case the laplacian can be partially diagonalized. We write

$\mathbb{R}^{\mathcal{G}}=\left(\bigoplus_{z^{n}=1}U_{z}\right)\oplus U_{0}$

wjere $U_{z}$ is as above and $U_{0}\cong\mathbb{R}$ is the space of functions
on $v_{0}$. The laplacian $\Delta$ preserves each $U_{z}$ for $z\neq 1,0$, and
acts by multiplication by

$\lambda_{z}=w_{0}+w(1+z+z^{2})(1+z^{-1}+z^{-2})$

on these $U_{z}$. The laplacian does not preserve $U_{1}$ but does preserve
the sum $U_{1}\oplus U_{0}$. On the space $U_{1}\oplus U_{0}$, with basis
given by $e_{1}=(1,\dots,1,0)$ and $e_{0}=(0,\dots,0,1)$, the laplacian acts
as the matrix

$\begin{pmatrix}w_{0}+9w&2w_{0}\\\ 2nw_{0}&4nw_{0}\end{pmatrix}.$ (14)

More generally, for a tile of size $A$, the matrix would be

$\begin{pmatrix}w_{0}+A^{2}w&(A-1)w_{0}\\\
(A-1)nw_{0}&(A-1)^{2}nw_{0}\end{pmatrix};$

we’ll use this below. Likewise the kernel $\mathbb{K}=D^{*}\Delta^{-1}D$ has a
similar decomposition into the subspaces $U_{z}$ and $U_{1}\oplus U_{0}$. On
$U_{z}$ it acts as multiplication by

$\frac{(1+z+z^{2})(1+z^{-1}+z^{-2})}{w_{0}+w(1+z+z^{2})(1+z^{-1}+z^{-2})}.$

On $U_{1}\oplus U_{0}$, if we take a single triomino $X_{t}$ it corresponds to
the vector $\frac{3}{\sqrt{n}}e_{1}+0e_{0}$. Thus the contribution to the
covariance is $\frac{1}{wn},$ which is $\frac{9}{n}$ times the $1,1$ entry in
the inverse of (14).

Now using Theorem 5.1,

$\mathrm{Cov}(X_{s},X_{t})=Kw\delta_{s=t}-\frac{Kw}{n}-\frac{Kw^{2}}{n}\sum_{\stackrel{{\scriptstyle
z^{n}=1}}{{z\neq
1}}}\frac{z^{s-t}(1+z+z^{2})(1+z^{-1}+z^{-2})}{w_{0}+w(1+z+z^{2})(1+z^{-1}+z^{-2})}.$
(15)

Here the second term $-\frac{Kw}{n}$ is the component of
$-Kw^{2}(D^{*}\Delta^{-1}D)_{s,t}$ on the subspace $U_{1}$, that is
$-Kw^{2}(\frac{1}{nw})$.

The covariance (15) only depends on $s-t$; without loss of generality assume
$t=0$. Now (15) simplifies to

$\mathrm{Cov}(X_{0},X_{s})=\frac{Kw_{0}}{n}\sum_{\stackrel{{\scriptstyle
z^{n}=1}}{{z\neq 1}}}\frac{z^{s}}{w_{0}/w+(1+z+z^{2})(1+z^{-1}+z^{-2})}.$ (16)

As long as $w_{0}/w\ll 1$, to leading order this sum is localized on the
region where $z$ is close to one of the primitive cube roots of $1$. If
$\frac{1}{n^{2}}\ll\frac{w_{0}}{w}\ll 1$, we can approximate the sum with an
integral. Writing $z=e^{i(2\pi/3+2\pi j/n)},$ and $\varepsilon=w_{0}/w$, and
letting $u=2\pi j/n$, the contribution near $e^{2\pi i/3}$ is (up to a
$1+o(1)$ factor)

$\frac{Kw_{0}\exp(\frac{2\pi
is}{3})}{2\pi}\int_{-\infty}^{\infty}\frac{e^{isu}du}{\varepsilon+3u^{2}}=\frac{Kw_{0}\exp(\frac{2\pi
is}{3})}{2\sqrt{3\varepsilon}}e^{-|s|\sqrt{\varepsilon/3}}.$

The integral near the other primitive cube root of $1$ contributes the complex
conjugate, and we get

$\mathrm{Cov}(X_{0},X_{s})\sim\frac{Kw_{0}\cos(\frac{2\pi
s}{3})}{\sqrt{3\varepsilon}}e^{-|s|\sqrt{\varepsilon/3}}.$ (17)

We see the onset of the $3$-periodic correlations between the $X_{s}$ as
$\varepsilon\to 0$.

Note that the variances and covariances here are much larger than in the
$w_{0}=0$ case: here they are of order

$\frac{Kw_{0}}{\sqrt{\varepsilon}}=K\sqrt{ww_{0}}=\frac{K\sqrt{\alpha_{0}(2-\alpha_{0})}}{2n},$

which is (for constant $\alpha_{0}$) a factor of $n$ larger than that for the
pure triomino case above.

For smaller $\varepsilon$, write $\varepsilon=\beta/n^{2}$; in this case
expanding near the primitive cube roots of $1$ gives

$\displaystyle\mathrm{Cov}(X_{0},X_{s})$ $\displaystyle\sim$
$\displaystyle\frac{K\beta\cos(\frac{2\pi
s}{3})}{n^{2}}\sum_{j\in\mathbb{Z}}\frac{1}{\beta+12\pi^{2}j^{2}}$
$\displaystyle=$ $\displaystyle\frac{K\sqrt{\beta}\cos(\frac{2\pi
s}{3})}{2n^{2}\sqrt{3}}\coth(\frac{\sqrt{\beta}}{2\sqrt{3}})=\frac{K\cos(\frac{2\pi
s}{3})}{n^{2}}(1+\frac{\beta}{36}+O(\beta^{2})),$

where in the last line we used the Mittag-Leffler expansion for the hyperbolic
cotangent function

$\pi\coth(\pi z)=\frac{1}{z}+2\sum_{n=1}^{\infty}\frac{z}{z^{2}+n^{2}}.$

####  6.1.3 Quasiperiodic example in dimension $1$

On $\mathbb{Z}/n\mathbb{Z}$ let $\mathrm{t}$ be the tile with characteristic
polynomial $p(z)=2/z+1+2z$. It has two roots of modulus $1$ which are not
roots of unity. We tile with $\mathrm{t}$ and $\mathrm{t}_{0}$ the singleton,
and as before let $\varepsilon=w_{0}/w$. The covariance function is given by
the analogue of (16) where the denominator is replaced by
$w_{0}/w+|p(z)|^{2}$. For $1/n^{2}\ll w_{0}/w\ll 1$ we get

$\mathrm{Cov}(X_{0},X_{s})\sim\frac{Kw_{0}}{2\pi
i}\int_{S^{1}}\frac{z^{s}}{\varepsilon+p^{2}}\frac{dz}{z}.$

The integral can be localized near the two roots $e^{\pm i\theta_{0}}$ of $p$.
Letting $z=e^{i(\theta_{0}+x)}$ we get

$\mathrm{Cov}(X_{0},X_{s})\sim
2Kw_{0}\operatorname{Re}\left[\frac{e^{is\theta_{0}}}{2\pi}\int_{\mathbb{R}}\frac{e^{isx}dx}{\varepsilon+2p^{\prime}(\theta_{0})^{2}x^{2}}\right]=\frac{Kw_{0}\cos(s\theta_{0})e^{-|s|\sqrt{\varepsilon
a}}}{\sqrt{\varepsilon a}}$

where $a=2p^{\prime}(\theta_{0})^{2}=32\sin^{2}(\theta_{0}).$

For a similar example without multiplicity, we can take
$p(z)=1+z+z^{3}+z^{5}+z^{6}$.

###  6.2 Examples in higher dimensions

####  6.2.1 Example: $L$ triomino

For a $2d$ example, consider tilings of the torus
$\mathbb{Z}^{2}/n\mathbb{Z}^{2}$ with translates of the triomino
$\mathrm{t}=\\{(0,0),(1,0),(0,1)\\}$, and constant multiplicities
$\mathbf{N}\equiv N$. Let $\Gamma$ be the sublattice of $\mathbb{Z}^{2}$
generated by $(1,1)$ and $(3,0)$. Translates of $\mathrm{t}$ by $\Gamma$ tile
$\mathbb{Z}^{2}$. Translates of this tiling by the three cosets of $\Gamma$ in
$\mathbb{Z}^{2}$ give three distinct periodic tilings.

We have $K=Nn^{2}/3$, $\alpha=3/n^{2}$. The graph is regular, so for the
critical gauge we have $x_{v}\equiv x,$ and $w_{t}=1/n^{2}$.

For $z,u$ two $n$th roots of $1$, define $U_{z,u}$ to be the subspace of
$\mathbb{C}^{V}\cong\mathbb{C}^{n^{2}}$ consisting of $(z,u)$-periodic
functions, that is, functions $f:\mathbb{Z}^{2}/n\mathbb{Z}^{2}\to\mathbb{C}$
such that $f(i+1,j)=zf(i,j)$ and $f(i,j+1)=uf(i,j)$.

The action of $\Delta$ on $U_{z,u}$ is multiplication by

$\lambda_{z,w}=\frac{1}{n^{2}}(1+z+u)(1+z^{-1}+u^{-1}).$

If we index tiles according to the location of their lower left vertex, then
the matrix $D$ corresponds to multiplication by $1+z+u$, and $D^{*}$ by
$1+z^{-1}+u^{-1}$. Then $\mathbb{K}=D^{*}\Delta^{-1}D$ is multiplication by
$n^{2}$ of spaces $U_{z,u}$ for which $1+z+u\neq 0$, and zero on spaces
$U_{z,u}$ where $1+z+u=0$.

If $n$ is not a multiple of $3$, then there are no subspaces $U_{z,u}$ where
$1+z+u=0$, and so $\mathbb{K}=n^{2}I$ is a scalar multiple of the identity,
and $\mathrm{Cov}$ is the zero matrix: all tile multiplicities are determined
and constant. If $n$ is a multiple of $3$, then
$\mathrm{Cov}=\frac{K}{n^{2}}(I-\mathbb{P})$ where $I-\mathbb{P}$ is the
projection onto the span of $U_{\omega,\omega^{2}}\oplus
U_{\omega^{2},\omega}$ (here $\omega=e^{2\pi i/3}$). The covariance between
tiles is $2K/n^{4}$ if they lie in the same translate of $\Gamma$ and
$-K/n^{4}$ if they do not. Again this is a _perfect crystal_ ; the tiles
differing by translates in $\Gamma$ are perfectly correlated.

Let us now allow a small fraction of singleton tiles: $N_{0}+n^{2}N=3K$, and
$\alpha_{0}=N_{0}/K$. As before $w_{0}=\frac{\alpha_{0}}{2n^{2}}$ and
$w=\frac{2-\alpha_{0}}{2n^{2}}$. The laplacian $\Delta$ still preserves each
$U_{z,u}$ for $(z,u)\neq(1,1)$, and on $U_{z,u}$ acts by multiplication by

$\lambda_{z,w}=w_{0}+w(1+z+u)(1+z^{-1}+u^{-1}).$

On $U_{1,1}\oplus U_{0}$ it acts as the matrix

$\begin{pmatrix}w_{0}+9w&2w_{0}\\\ 2n^{2}w_{0}&4n^{2}w_{0}\end{pmatrix}.$

We have

$\displaystyle\mathrm{Cov}(X_{(0,0)},X_{(s,t)})$
$\displaystyle=Kw\delta_{s=t=0}-\frac{Kw}{n^{2}}-\frac{Kw^{2}}{n^{2}}\sum_{\stackrel{{\scriptstyle
z^{n}=1=u^{n}}}{{(z,u)\neq(1,1)}}}\frac{z^{s}u^{t}(1+z+u)(1+z^{-1}+u^{-1})}{w_{0}+w(1+z+u)(1+z^{-1}+u^{-1})}$
$\displaystyle=\frac{Kw_{0}}{n^{2}}\sum_{\stackrel{{\scriptstyle
z^{n}=1=u^{n}}}{{(z,u)\neq(1,1)}}}\frac{z^{s}u^{t}}{w_{0}/w+(1+z+u)(1+z^{-1}+u^{-1})}.$

Now assume $w_{0}/w\ll 1$ and is fixed as $n\to\infty$. The sum is then
localized to the region near $(z,u)=(\omega,\omega^{2})$ or
$(\omega^{2},\omega).$ Near the first point write $z=e^{i(2\pi/3+x)}$ and
$u=e^{i(4\pi/3+y)}$, and take $\varepsilon=w_{0}/w$. The contribution near
this point is

$\sim\frac{Kw_{0}\exp(\frac{2\pi
i(s+2t)}{3})}{(2\pi)^{2}}\iint_{\mathbb{R}^{2}}\frac{e^{i(sx+ty)}dx\,dy}{\varepsilon+x^{2}-xy+y^{2}}.$

The integral here for $(s,t)\neq(0,0)$ is a Bessel-K function (see the
appendix Section 8),

$=\frac{Kw_{0}\exp(\frac{2\pi
i(s+2t)}{3})}{\pi\sqrt{3}}B(\frac{2}{\sqrt{3}}\sqrt{\varepsilon(s^{2}+st+t^{2})}).$

Summing over both roots, and taking $\varepsilon$ small we have

$\mathrm{Cov}(X_{(0,0)},X_{(s,t)})=\frac{Kw_{0}\cos(\frac{2\pi}{3}(s+2t))}{\pi\sqrt{3}}\left(\log\frac{1}{\varepsilon}-\log\frac{2(s^{2}+st+t^{2})}{3}-2\gamma_{E}+O(\varepsilon)\right)$

where $\gamma_{E}$ is the Euler gamma.

For the variance, that is, when $(s,t)=(0,0)$,

$\mathbf{Var}(X_{0,0})=\frac{Kw_{0}}{(2\pi
i)^{2}}\iint\frac{dzdu}{(1+u)(z-r_{1})(z-r_{2})}$

where $r_{1},r_{2}$ are the roots of the denominator which is a quadratic
polynomial in $z$. Let $r_{1}$ be the root inside the unit circle; then using
residues this is

$=\frac{Kw_{0}}{2\pi i}\int_{S^{1}}\frac{du}{(1+u)(r_{1}-r_{2})}$
$=\frac{Kw_{0}}{2\pi}\int_{0}^{2\pi}\frac{d\theta}{\sqrt{3+6\varepsilon+\varepsilon^{2}+(4+4\varepsilon)\cos\theta+2\cos(2\theta)}}$

and splitting the integral into parts near $\theta=2\pi/3,4\pi/3$ and the
remainder, we arrive at

$=\frac{Kw_{0}}{\pi\sqrt{3}}(\log\frac{9}{\varepsilon})(1+o(1)).$

See Figure 5 for a plot of the covariances for small $\varepsilon$.

Figure 5: Covariances for the L polyomino for $\varepsilon=0.001$ (left) and
in the limit $\varepsilon=0$ (right). Here values are scaled to the range from
$1$ (orange) to $-1$ (dark blue) with $0$ being white.

####  6.2.2 Examples in $\mathbb{Z}^{2}$ with simple zeros.

The above example of the $L$ triomino can be generalized. Consider the case of
tiling with a polyomino $\mathrm{t}$ in $\mathbb{Z}^{2}$ and a small density
of singletons $\mathrm{t}_{0}$. The characteristic polynomial $p(z,u)$ of a
polyomino $\mathrm{t}$ is defined as

$p(z,u)=\sum_{(i,j)\in\mathrm{t}}z^{i}u^{j}.$

As in the previous section we have

$\displaystyle\mathrm{Cov}(X_{0,0},X_{s,t})$
$\displaystyle=\frac{Kw_{0}}{n^{2}}\sum_{\stackrel{{\scriptstyle
z^{n}=1=w^{n}}}{{(z,u)\neq(1,1)}}}\frac{z^{s}u^{t}}{w_{0}/w+p(z,u)p(1/z,1/u)}$
$\displaystyle\sim\frac{Kw_{0}}{(2\pi)^{2}}\iint_{(S^{1})^{2}}\frac{z^{s}u^{t}}{\varepsilon+p(z,u)p(1/z,1/u)}\frac{dz}{iz}\frac{du}{iu}.$
(18)

Suppose now that $p(z,u)$ has a finite number of roots
$\\{(z_{i},u_{i})\\}_{i=1,\ldots,k}$ on $\mathbb{T}^{2}=S^{1}\times S^{1}$. A
root $(z,u)$ is _simple_ if $\frac{zP_{z}}{uP_{u}}\not\in\mathbb{R}$; this
means the zero set of $p$ intersects the unit torus transversely at that
point. We suppose for the moment that all roots $(z_{i},u_{i})$ are simple.

For $w_{0}/w$ a small constant the integral (18) can be localized near each
root. Near a simple root $(z_{0},u_{0})$ the integral is approximated by

$\frac{Kw_{0}z_{0}^{s}u_{0}^{t}}{(2\pi)^{2}}\int_{\mathbb{R}^{2}}\frac{e^{i(sx+ty)}\,dx\,dy}{\varepsilon+ax^{2}+bxy+cy^{2}}$

where the quadratic form in the denominator is

$ax^{2}+bxy+cy^{2}=|z_{0}P_{z}(z_{0},u_{0})x+u_{0}P_{u}(z_{0},u_{0})y|^{2}.$

Since the root is simple this quadratic form is positive definite, and the
integral is defined and finite (for $(s,t)\neq(0,0)$). This integral is a
(linear image of a) Bessel-K function (see section 8). Summing over all roots,
the covariance is a superposition of these Bessel-type functions.

We also need the $s=t=0$ case, that is, the variance, which we cannot get
using Bessel functions. We have

$\mathbf{Var}(X_{0,0})=\frac{Kw_{0}}{(2\pi)^{2}}\iint\frac{1}{\varepsilon+p(z,u)p(1/z,1/u)}\frac{dz}{iz}\frac{du}{iu}.$

Localize to a ball $B$ of radius $M\sqrt{\varepsilon}$ around a simple root
$(z_{0},u_{0})$i (where $M$ is a large constant). Let
$z=z_{0}e^{is\sqrt{\varepsilon}}$ and $u=u_{0}e^{it\sqrt{\varepsilon}}$. The
contribution from this ball is

$\frac{Kw_{0}}{(2\pi)^{2}}\int_{B}\frac{ds\,dt}{1+as^{2}+bst+ct^{2}+O(\varepsilon^{1/2})}=\frac{Kw_{0}}{2\pi}\left(\frac{1}{\sqrt{4ac-b^{2}}}\log\frac{1}{\varepsilon}+O(1)\right).$

This is then summed over all roots.

######  Theorem 6.1.

Suppose $\mathrm{t}$ is a polyomino whose characteristic polynomial has only
simple roots $\\{(z_{j},u_{j})\\}_{j=1,\ldots,k}$ on $\mathbb{T}^{2}$, and
consider tilings of $\mathbb{Z}^{2}$ with translates of $\mathrm{t}$ and the
singleton $\mathrm{t}_{0}$. Fix a small $\varepsilon=w_{0}/w$, then the
covariance function is, up to a $1+o_{\varepsilon}(1)$ factor,

$\mathrm{Cov}(X_{0,0},X_{s,t})=Kw_{0}\sum_{j=1}^{k}\frac{z_{j}^{s}u_{j}^{t}}{(2\pi)^{2}}\int_{\mathbb{R}^{2}}\frac{e^{i(sx+ty)}\,dx\,dy}{\varepsilon+|z_{j}P_{z}(z_{j},u_{j})x+u_{j}P_{u}(z_{j},u_{j})y|^{2}}.$

####  6.2.3 Quasiperiodic example in dimension $2$

This is an example with simple roots which are not roots of unity. Consider
the “key” polyomino (Figure 6) with

$p(z,u)=1+z+z^{2}+z^{3}+z^{4}+z^{5}+z^{6}+u(1+z+z^{4}+z^{6}).$

It has $2$ simple roots on $\mathbb{T}^{2}$,
$(z_{0},u_{0}),(\bar{z}_{0},\bar{u}_{0})$, and $z_{0},u_{0}$ have arguments
$\theta,\phi$ which are not rationally related to each other, that is, there
are no integers $m_{1},m_{2}$ such that $z_{0}^{m_{1}}u_{0}^{m_{2}}=1$ except
$m_{1}=0=m_{2}$ (this requires a short Galois theory argument).

Figure 6: The key polyomino and its (quasiperiodic) covariance function in the
limit $\varepsilon\to 0$.

This implies that the resulting covariance function
$\mathrm{Cov}(X_{0,0},X_{s,t})$ is a quasiperiodic function of $(s,t)$ in the
$\varepsilon\to 0$ limit: specifically, to leading order

$\mathrm{Cov}(X_{0,0},X_{s,t})=Kw_{0}C(\log\frac{1}{\varepsilon})\cos(s\theta+t\phi)\Big{(}1+o_{\varepsilon}(1)\Big{)}$
(19)

for a constant $C$. See Figure 6 for a numerical plot.

A simpler quasiperiodic example, if we allow multiplicities, is the tile with
$p(z,w)=2+3z+4u.$ It also has $2$ simple roots (with arguments corresponding
to the angles of the $2,3,4$-triangle). The covariance formula is again of the
form (19) where $\theta,\phi$ are these angles.

Figure 7: Quasiperiodic covariance function (when $\varepsilon=0$) for the
$2,3,4$-weighted $L$ polyomino.

###  6.3 Other examples

Not all polyominos have roots on $\mathbb{T}^{2}$. For example the tile of
Figure 8 has the property that its characteristic polynomial has no roots on
$\mathbb{T}^{2}$.

Figure 8: A (nonconnected) polyomino with no crystal structure.

As a consequence its covariance function $\mathrm{Cov}(X_{0,0},X_{s,t})$
decays exponentially in $|s|+|t|$ even for $w_{0}=0$. This example is however
somewhat special since its characteristic polynomial
$p(z,u)=(1+z+z^{3})(1+u+u^{3})$ is a product of two $1d$ polynomials. A
genuinely 2d example which does not factor is not easy to find. Here is one:

$p(z,u)=1+z+z^{2}+z^{3}+z^{7}+z^{9}+z^{12}+z^{13}+z^{17}+u.$

For polyominos with multiplicity, an easy example is the one with
$p(z,u)=3+z+u$.

The third class of polyominos in $\mathbb{Z}^{2}$ has characteristic
polynomials with roots on $\mathbb{T}^{2}$ which are either _not simple_ or
_not isolated_. And generally for $d$-dimensional polyominos with $d>2$ the
roots on $\mathbb{T}^{d}$ are not typically isolated. It is harder to
formulate a general theory encompassing all these cases. We will simply
illustrate with a few examples.

####  6.3.1 The square polyomino

We consider the square polyomino with $p(z,u)=(1+z)(1+u).$ We evaluate the
integral (18). We first perform a contour integral over $u$. Assume $t\geq 0$.

$\displaystyle\frac{Kw_{0}}{(2\pi
i)^{2}}\int\frac{z^{s}u^{t}}{\varepsilon+(1+z)(1+u)(1+\frac{1}{z})(1+\frac{1}{u})}\frac{du}{u}\frac{dz}{z}$
$\displaystyle=\frac{Kw_{0}}{(2\pi i)^{2}}\int\frac{z^{s}u^{t}}{\varepsilon
zu+p^{2}}\,du\,dz$ $\displaystyle=\frac{Kw_{0}}{(2\pi
i)^{2}}\int\frac{z^{s}u^{t}}{(1+z)^{2}(u-r_{1})(u-r_{2})}du\,dz.$

Roots $r_{1},r_{2}$ of the denominator are real with product $1$; choose
$|r_{1}|<1<|r_{2}|$. We get

$=\frac{Kw_{0}}{2\pi}\int_{0}^{2\pi}\frac{z^{s}r_{1}^{t}}{(1+z)^{2}(r_{1}-r_{2})}d\theta=\frac{Kw_{0}}{2\pi\sqrt{\varepsilon}}\int\frac{\cos(\theta
s)r_{1}^{t}}{\sqrt{8+\varepsilon+8\cos\theta}}d\theta.$

This is an elliptic function. For small $\varepsilon$ the integral
concentrates near $\theta=\pi$, and is to leading order, for constant $s,t$,

$Kw_{0}\frac{(-1)^{s+t}}{2\pi\sqrt{\varepsilon}}\left(\log\frac{16}{\sqrt{\varepsilon}}-2h_{s}-2h_{t}+O(\varepsilon)\right)$

where $h_{m}=1+\frac{1}{3}+\frac{1}{5}+\dots+\frac{1}{2m-1}.$

We thus have

$\mathrm{Cov}(X_{0,0},X_{s,t})=Kw_{0}\frac{(-1)^{s+t}}{2\pi\sqrt{\varepsilon}}\left(\log\frac{16}{\sqrt{\varepsilon}}-2h_{|s|}-2h_{|t|}+O(\varepsilon)\right).$

For larger $(s,t)$, on the order $s,t=O(\varepsilon^{-1/2})$, the Fourier
coefficients decay exponentially at a rate determined by the component around
$(0,0)$ in the complement of the amoeba of $\varepsilon+|p|^{2}$, (the
_amoeba_ is the image of the zero set of $\varepsilon+|p|^{2}$ under the map
$(z,w)\mapsto(\log|z|,\log|w|)$). In this case the component around the origin
for small $\varepsilon$ tends to a square of width
$\frac{1}{2}\sqrt{\varepsilon}$, so the Fourier coefficients are of modulus of
order $\exp(-\frac{\sqrt{\varepsilon}}{2}\min\\{|s|,|t|\\}).$

Figure 9: Covariances for the square polyomino, when $\varepsilon=.001$
(left), and in the limit $\varepsilon\to 0$ (right).

####  6.3.2 The $+$ polyomino

We consider the polyomino with $p(z,u)=1+z+1/z+u+1/u.$ We need to compute
($Kw_{0}/n$ times) the Fourier series of $\frac{1}{\varepsilon+|p|^{2}}$. Let
$z=e^{i\theta}$ and $u=e^{i\phi}$. The polynomial $p=1+2\cos\theta+2\cos\phi$
vanishes on a whole curve on $\mathbb{T}^{2}$, where $\theta$ runs over the
range $\theta\in[\pi/3,5\pi/3]$. For small $\varepsilon$ the integral
concentrates near this curve. Let $(\theta_{0},\phi_{0})$ be a point on the
curve. For $\theta$ fixed, the denominator has the form
$\varepsilon+a(\theta)(\phi-\phi_{0})^{2}+O(\phi-\phi_{0})^{3}$ where
$a(\theta)=3-4\cos\theta-4\cos^{2}\theta$.

The contribution for this slice is (with $x=\phi-\phi_{0}$, and to leading
order)

$\frac{1}{2\pi}\int_{\mathbb{R}}\frac{dx}{\varepsilon+a(\theta)x^{2}}=\frac{1}{2\sqrt{\varepsilon
a(\theta)}}.$

We thus have to leading order (for fixed $(s,t)$ as $\varepsilon\to 0$)

$\mathrm{Cov}(X_{0,0},X_{s,t})=\frac{Kw_{0}}{2\pi\sqrt{\varepsilon}}\int_{\pi/3}^{5\pi/3}\frac{\cos(s\theta+t\phi)d\theta}{\sqrt{3-4\cos\theta-4\cos^{2}\theta}}$

where $\cos\phi+\cos\theta=1/2.$

Multiplied by $\sqrt{\varepsilon}$ these covariances still tend to zero as
$|s|+|t|\to\infty$, although the decay is slow, of order
$1/(s^{2}+t^{2})^{1/4}$. See Figure 10 for a numerical plot.

Figure 10: Covariances for the “plus” polyomino, in the limit $\varepsilon\to
0$.

##  7 Further directions

Because tilings are so diverse, there are many directions for further
research. Here are some ideas.

  1. 1.

How is the covariance function for a 3D polyomino different? Typically the
characteristic polynomial intersects the unit $3$-dimensional torus
$\mathbb{T}^{3}$ in a $1$-dimensional set. Is there an analogue of Theorem
6.1?

  2. 2.

Find interesting examples with multiple tiles in $\mathbb{Z}^{2}$.

  3. 3.

For tilings of a planar domain with copies of the $L$ polyomino and a small
density of monomers (or plus polyomino and no monomers), understand the
influence of the boundary on the tiling. Can one create situations where there
is coexistence of the multiple phases?

  4. 4.

What behavior do we expect for the Coulomb gas for a general tiling problem?
What about the $L$ triomino?

  5. 5.

Wang tiles (squares with colored edges, tiled so that adjacent tiles share the
same color) can be used to emulate any Turing machine. What is the
multinomial-tiling analog of such a computation?

  6. 6.

What is the natural multinomial analogue of the random partition model? What
is the limit shape of a random such partition in that model (in the
appropriate limit)?

##  8 Appendix: The Bessel-K function

The Bessel-K function, or modified Bessel function of the second kind, $B(s)$,
can be defined for $s\in\mathbb{R}\setminus\\{0\\}$ by the integral

$B(s)=\frac{1}{2\pi}\iint_{\mathbb{R}^{2}}\frac{e^{isx}\,dx\,dy}{1+x^{2}+y^{2}}.$

For $z\in\mathbb{C}$, $B(|z|)$ is the Green’s function for the massive
laplacian, that is, it satisfies

$(I-\Delta)B(|z|)=\delta_{0}$

where $\delta_{0}$ is the point measure.

The function $B(s)$ has a logarithmic singularity at the origin, with
expansion

$B(s)=\log\frac{1}{s}+\log 2-\gamma_{E}+O(s^{2}\log s).$

An integral of the form

$\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\frac{e^{i(sx+ty)}\,dx\,dy}{\varepsilon+ax^{2}+bxy+cy^{2}}$

where the quadratic form $ax^{2}+bxy+cy^{2}$ is positive definite, can be
converted into a Bessel-K integral with a linear change of coordinates,
yielding

$=\frac{1}{\sqrt{c-\frac{b^{2}}{4a}}}B\left(\sqrt{\frac{\varepsilon(cs^{2}-bst+at^{2})}{ac-b^{2}/4}}\right).$

## References

  * [1] D. Beauquier, M. Nivat, É. Remila, M. Robson, Tiling figures of the plane with two bars, Comput. Geom. 5 (1995), 1–25.
  * [2] R. Berger, The undecidability of the domino problem, Memoirs of the American Mathematical Society, 66 (1966), 293–357.
  * [3] H. Cohn, N. Elkies, J. Propp Local statistics for random domino tilings of the Aztec diamond. _Duke Math. J._ 85 (1996), no. 1, 117–166.
  * [4] H. Cohn, R. Kenyon, J. Propp, A variational principle for domino tilings, Journal of the American Mathematical Society 14 (2001), 297–346.
  * [5] M. Ciucu, Dimer packings with gaps and electrostatics _PNAS_ (2008) 105 (8) 2766–2772.
  * [6] N. G. de Bruijn, Algebraic theory of Penrose’s non-periodic tilings of the plane, I, II”, Indagationes Mathematicae, (1981), 43 (1): 39–66
  * [7] J. de Gier, B. Nienhuis, Integrability of the square-triangle random tiling model. _Phys. Rev._ E 55, 3926
  * [8] N. Elkies, G. Kuperberg, M. Larsen, and J. Propp, Alternating-sign matrices and domino tilings (Part I), J. Algebraic Combin. 1 (1992), 111–132.
  * [9] N. Elkies, G. Kuperberg, M. Larsen, and J. Propp, Alternating-sign matrices and domino tilings (Part II), J. Algebraic Combin. 1 (1992), 219–234.
  * [10] P. A. Kalugin, The square-triangle random-tiling model in the thermodynamic limit,  Journal of Physics A: Mathematical and General, Volume 27, Number 11.
  * [11] P. W. Kasteleyn, Graph theory and crystal physics, in Graph Theory and Theoretical Physics, Academic Press, London, 1967, pp. 43–110.
  * [12] R. Kenyon, Conformal invariance of domino tiling, Ann. Probab. 28 (2000), no. 2, 759–795
  * [13] R. Kenyon, Dominos and the Gaussian free field, Ann. Probab. 29 (2001), no. 3, 1128–1137.
  * [14] R. Kenyon, A. Okounkov, _Acta Math._ 199 (2007), no. 2, 263–302.
  * [15] R. Kenyon, A. Okounkov, S. Sheffield, Dimers and amoebae, Ann. Math., 163 (2006) 1019–1056.
  * [16] R. Kenyon, C. Pohoata, Conformal invariance in the multinomial tiling model, in preparation.
  * [17] R. Kenyon, C. Pohoata, Limit shapes in the multinomial tiling model, in preparation.
  * [18] L. A. Levin, Universal sequential search problems, Problems of Information Transmission, 9 (3): 265–266, 1973.
  * [19] E. H. Lieb. Exact solution of the problem of the entropy of two-dimensional ice, Physical Review Letters, 18(17):692, 1967.
  * [20] P. A. MacMahon, Combinatory Analysis Vol 2. (1916), Cambridge University Press.
  * [21] D. Shechtman, I. Blech, D. Gratias, J. Cahn, J. (1984). ”Metallic Phase with Long-Range Orientational Order and No Translational Symmetry”, Physical Review Letters 53 (20): 1951–1953.
  * [22] E. M. Stein, R. Shakarchi, Complex Analysis, Volume 2 of Princeton lectures in analysis, Princeton University Press, 2010.
  * [23] W. Temperley, M. Fisher, Dimer problem in statistical mechanics—an exact result, Philos. Mag. (8) 6 (1961) 1061–1063.