# Quasicrystals with open boundaries
###### tags: `masters projects`,`quasicrystal`

### Owners (the only one with the permission to edit the main text, others can comment)

Robert, Alptug, Laura, Frank

# Project description
Ordinarily, quasicrystal simulations are performed using periodic boundary conditions to avoid finite-size and boundary effects. However, this goes against a defining feature of quasicrystals, namely that there should be no long-range periodicity present. Instead, open boundary conditions could be implemented, where the crystal is confined by line tension. By keeping track of the effect this tension has on the crystal, we can compensate for it in the results of the simulation.

The aims of the project will be to (a) simulate a 2D quasicrystal with closed boundary conditions, (b) simulate a 2D quasicrystal with open boundary conditions, (c) compare the results of the different methods, and (d) exploit the open boundary conditions to measure the free energy of the quasicrystal.

## Quasicrystal vertices as points on a 4-dimensional lattice
Our system is comprised exclusively of squares and triangles, such that the only possible angles between our edges are $\frac{\pi}{2}$ and $\frac{\pi}{3}$, with a largest common denominator of $\frac{\pi}{6}$. This means that, after fixing a reference edge, all other edges have an orientation that is a multiple of $\frac{\pi}{6}$ from the reference edge, such that there are only $12$ possible orientations for all of our edges. All of these orientations can be described by only $4$ unit vectors:
$$
\mathbf{e}_i = \left( \text{cos}\left( \frac{\pi}{6}i \right), \text{sin}\left( \frac{\pi}{6}i \right) \right)
$$
with $i \in [0,3]$. Since our vertices will be unit distance apart, all vertices can be described by the sum of integer multiples of these unit vectors.

## Generating an initial quasicrystal configuration
### Stampfli inflation

We start with a $D \times D$ regular square lattice of vertices which we use as a seed for the quasicrystal through a method called Stampfli inflation. This method consists of "inflating" our lattice, meaning we scale up our coordinates by the inflation constant $\lambda = 2 +\sqrt{3}$, after which we replace every vertex with a dodecagonal wheel of vertices in either of the two possible orientations (It is the same configuration, but offset by $\frac{\pi}{6}$). We repeat this inflation process $I$ times to end up with a quasicrystal of $[I,D]$-type.

To scale up our 4-space coordinates, $\{n_i\}$, with $i \in [0,3]$ we need to figure out by what amount each of these components need to be increased in order to end up at the right position.
The parallel-space coordinates scale trivially, $\mathbf{x} = \lambda\mathbf{x}$, and are related to the 4-space coordinates according to

$$
\begin{gather}
\begin{pmatrix}
x\\
y
\end{pmatrix} =
\begin{pmatrix}
1 & \frac{\sqrt{3}}{2} & \frac{1}{2} & 0\\
0 & \frac{1}{2} & \frac{\sqrt{3}}{2} & 1
\end{pmatrix}
\begin{pmatrix}
n_0\\
n_1\\
n_2\\
n_3
\end{pmatrix}
\end{gather}
$$

Note the symmetry. Scaling up both sides, we find, for $x$

$$
\begin{align*}
\lambda x &= \lambda\left( n_0 + \frac{\sqrt{3}}{2}n_1 + \frac{1}{2}n_2 \right)\\
x' &= (2 +\sqrt{3})n_0 + (\sqrt{3} + \frac{3}{2})n_1 + (1 + \frac{\sqrt{3}}{2})n_2\\
n_0' + \frac{\sqrt{3}}{2}n_1' + \frac{1}{2}n_2' &= (2 +\sqrt{3})n_0 + (\sqrt{3} + \frac{3}{2})n_1 + (1 + \frac{\sqrt{3}}{2})n_2
\end{align*}
$$

where the primed coordinates are the coordinates after the transformation.

All these numbers live in the quadratic number field of $\mathbb{Q}(\sqrt{3})$, an extension of the rationals $\mathbb{Q}$. Someone better versed in number theory might be able to connect some dots here, since $\lambda$, which allows us to fit dodecagons in our lattice, is equal to the fundamental unit of the ring of integers of this field, which has a discriminant of $12$. Or maybe that is just a coincidence.

More importantly, this is a field of degree two, with a possible basis being $\{ 1, \sqrt{3} \}$. We can decompose our equation using this basis to find 2 equations directly relating the primed to the unprimed coordinates

$$
\left( n_0' + \frac{1}{2}n_2' \right) + \sqrt{3}\left[ \frac{1}{2}n_1' \right] = \left( 2n_0 + \frac{3}{2}n_1 + n_2 \right) + \sqrt{3}\left[ n_0 + n_1 + \frac{1}{2}n_2 \right]\\
\begin{cases}
\begin{align*}
n_0' + \frac{1}{2}n_2' &= 2n_0 + \frac{3}{2}n_1 + n_2 \\
\frac{1}{2}n_1' &= n_0 + n_1 + \frac{1}{2}n_2
\end{align*}
\end{cases}
$$

Similarly, we find 2 such equations for $y$, such that we can solve for the primed coordinates in terms of the original coordinates to find the following relationship.

$$
\mathbf{n}' = \Lambda \mathbf{n}
$$
where
$$
\Lambda =
\begin{pmatrix}
2 & 1 & 0 & -1\\
2 & 2 & 1 & 0\\
0 & 1 & 2 & 2\\
-1 & 0 & 1 & 2
\end{pmatrix}
$$

### Periodic boundary conditions
Before we start with open boundary conditions, we want to implement the more common periodic boundary conditions as a reference. This means that vertices near the left or bottom boundary (small $x$ or $y$ coordinates) need to be duplicated as mirror image near the right or top boundary. Where this edge is located, however, depends on what inflation step we are currently at. It is trivial for $I = 0$, but after that we find the following relationship between the primed coordinates of the image and the original coordinates
$$
x_i' = x_i + D\lambda ^I
$$
Where $x_i$ corresponds to either one or both of the parallel-space coordinates, depending on which of the boundaries the vertex is close to.

Translating this to 4-space coordinates is done by a procedure similar to what we used to inflate the coordinates before. The result for creating an image along the $x$-axis is
$$
\begin{gather}
\begin{pmatrix}
n_0'\\
n_1'\\
n_2'\\
n_3'
\end{pmatrix} =
\begin{pmatrix}
n_0 + \alpha D\\
n_1 + 2\beta D\\
n_2\\
n_3 - \beta D
\end{pmatrix}
\end{gather}
$$
where $\alpha$ and $\beta$ are the solutions to
$$
\lambda^I = \alpha + \beta\sqrt{3}
$$
Displacement along the $y$-axis can be found in a similar way, or through the symmetry considerations of $n_0 \leftrightarrow n_3$ and $n_1 \leftrightarrow n_2$.

:::spoiler Click for a (likely unimportant) tangent
These are the units of the (ring of integers $\mathbf{Z}[\sqrt{3}]$ of the) quadratic field $\mathbb{Q}(\sqrt{3})$, as $\lambda$ is the fundamental unit of this field.

I noticed that $\alpha$ and $\beta$, as functions of $I$, follow a sequence found on The On-Line Encyclopedia of Integer Sequences (A001075 and A001353), which also states that these solutions can be generated by the following recurrence relation
$$
A(I) = 4A(I-1) - A(I-2)
$$
for inflation step $I > 1$, where $\{ A(0), A(1) \}$ are given by $\{1, 2\}$ for $\alpha$ and by $\{0, 1\}$ for $\beta$. Furthermore, this sequence has the property that
$$
\lim_{I \rightarrow\infty}\frac{A(I)}{A(I-1)} = \lambda
$$
and the ratio between $\alpha$ and $\beta$ converges to the square root of three
$$
\lim_{I\rightarrow\infty}\frac{\alpha(I)}{\beta(I)} = \sqrt{3}
$$
:::

## Generating subsequent quasicrystal configurations (Monte Carlo)
To perform Monte Carlo simulations we need to generate new quasicrystal configurations using an update move. The move we will be using is called a *zipper* move, as described in Oxborrow, Henley (1983). First we need to select a random 'house', meaning a region that comprises an adjoining square and triangle. A new vertex in placed in the middle of the house and the bonds are changed to form two thin rhombi and a triangle. We will attempt to propagate these rhombi outward, by flipping them with the squares and triangles it touches in the direction of propagation. Once these two rhombi meet again, they are annihilated, and we have a new quasicrystal configuration. Thus, the first step is to find these 'houses'.
### Identifying houses
When a square and a triangle are adjacent, the shared vertices have an angle of $\frac{5\pi}{6}$ beteen them. In fact, because the only allowed angles in our system are $\frac{\pi}{2}$ and $\frac{\pi}{3}$, any time we come across two edges with an angle of $\frac{5\pi}{6}$ between them, it must be a house.

We can search through the list of neighbours of every vertex to see if there are bonds separated by $\frac{5\pi}{6}$, and when there are, that vertex defines the position of the house. There must also be a third bond in between them, which we use to define the orientation of the house. This third bond can be in one of two positions, but we only consider one of them, such that we do not double count the houses, and we can uniquely define each house by a single vertex ($i$) + orientation ($o$).
### Creating the middle vertex
The orientation of the house, labelled by an integer $o\in[0,11]$, determines where the new middle vertex must be placed relative to the existing vertex that defines the position of the house, labelled vertex $i$ with 4-space coordinates $\mathbf{n}_i$.

The middle vertex of the house, labelled $k$, is then given by the 4-space coordinates
$$
\mathbf{n}_k = \mathbf{n}_i + \mathbf{e}_m + \mathbf{e}_n
$$
with $m = o + 4 \pmod{12}$ and $n = o + 9 \pmod{12}$. Where the unit vectors with index greater than 3 can be constructed by linear combination of the first 4 unit vectors.


### Flip types
Once we have created this middle vertex, we have two rhombus shapes in our house. To use these rhombi to flip the tiles in the direction of propagation, we first need to determine the nature of these surrounding tiles. An 'A-type' flip is performed if one the edges of the rhombus is shared with a square (in the direction of propagation). If both are triangles, a 'B-type' flip is performed.

To check this, we consider the vertex $i$ that defines the position of the rhombus and check for neighbours in the $\mathbf{e}_o$ and $\mathbf{e}_{o+1}$ directions, with $o$ the integer defining the orientation of the rhombus. If a neighbour exists in either of these directions, then one of the adjacent shapes is a square and thus an A-type flip will have to be performed.

If we need to perform an A-type flip, there are two possibilities, depending on which of the rhombus' edges is shared with a square and which with a triangle. In the first case we place a new vertex that is $\mathbf{e}_{o+2} + \mathbf{e}_{o+9}$ seperated from the vertex $i$ and then remove the vertex $i$. The vertex defining the new rhombus after the flip is the neighbour of vertex $i$ in the $o$ direction. In the second case we place a new vertex that is $\mathbf{e}_{o+4} + \mathbf{e}_{o+11}$ seperated from the vertex $i$ and then remove the vertex $i$. The vertex defining the new rhombus after the flip is the neighbour of vertex $i$ in the $o+1$ direction.

For a B-type flip, we randomly choose between one of two possiblities. In the first we place the new vertex at $\mathbf{e}_{o+2} + \mathbf{e}_{o+9}$ from vertex $i$ and designate the neighbour in the $o+10$ direction as the new rhombus position. For the second possibility, we place the new vertex at $\mathbf{e}_{o+4} + \mathbf{e}_{o+11}$ and designate the neighbour in the $o+2$ direction as the new rhombus position. In both cases we remove vertex $i$ afterwards, just like for the A-type flip.

:::spoiler Example of path traversed by zipper:
:::

## Open boundary conditions
A potential problem with simulations like this is that simulations with periodic boundary conditions clash with one of the defining properties of real quasicrystals; namely that there should be long range periodicity. To get around this, we can discard the periodic boundary conditions and instead use open boundary conditions.

This, however, is not as straightforward as one would like. With open boundary conditions, there is no incentive for the quasicrystal to remain within a bounded region. Applying any type of Monte Carlo move can result in the the vertices near the edges to move away from the bulk. Thus, the quasicrystal falls apart.

Our attempt at getting around this will consist of imposing a "line tension". This tension imposes an energy cost (T.B.C.)

### Defects
To simulate real quasicrystals, we must allow for the presence of defects in our system. Zipper moves do place defects into our system, but only until they annihilate again. Neither creation nor annihilation of the rhombi preserves the number of vertices in our system. What we want is a way to introduce defects into the system without altering the number of vertices. We also want these defects to wander randomly throughout our system, as opposed to propagating in a specific direction.

To do this, we can make use of our newly acquired boundaries. The vertices along the edges of our system have some room to move around, allowing for defects (rhombi) to be introduced into the system, while preserving the total number of vertices.

An example of this is shown below:


Here, we consider the marked vertex at the edge of our system, which has an outside angle of $\frac{5\pi}{6}$, the same as a thin rhombus. This means we can pretend there actually is a rhombus attached to the outside of our quasicrystal, and propagate it into the system like we would in a zipper move. This brings a real defect into our system while preserving the total number of vertices.

Creating a defect like this will always increase the length of the contour of our system. This means that, under line tension, this move is energetically unfavourable compared to propagating an already existing defect through the system, as this incurs no energy penalty. Thus the creation of these defects can only happen at non-zero temperatures.

### Non-zero temperature
We introduce temperature into the system in the following way: we propose a random move, then check whether it has an associated energy cost (it does if it generates a new defect). If it does, we only accept it according to the probability $P = e^{-\beta\Delta U}$, where $\beta = \frac{1}{k_BT}$ and $\Delta U$ is the difference in energy between the configurations of our system before and after the proposed move. This is the Metropolis algorithm.

The energy follows that of linear surface-like tension: $U = \gamma L$, where $L$ is the total length of the contour of our system, and $\gamma$ some proportionality constant. The exact influence of this constant on the behaviour of our system is something still to be analysed. For $\gamma > 0$, this energy is minimised by minimising $L$. This means that this condition should favour configurations that keep the quasicrystal bounded within a small area.

~~Since any possible move can only change the length of the contour by 1 at most, the probability of allowing the creation of a new defect becomes $P = e^{-\beta k}$. Removing a defect and reducing the length of our contour (decreasing the total energy of our system) results in a $P > 1$, meaning these moves are always accepted in this Metropolis algorithm.~~

This is not true, since we can also allow angles greater than $\frac{5\pi}{6}$ to generate defects:


This increases the total contour by 2 unit lengths, because it destroys the existing rhombus when the new one is formed.

Also, consider the following situation:


Clearly, there are moves that introduce defects into the system while not actually increasing the length of the contour all! These "free" defects can enter and exit our system at no energy penalty. I am not sure whether this could be problematic for the stability of the quasicrystal.
As defects can freely exit the crystal in a similar matter, we might see a 'baseline' defect density that does not depend on $\gamma$ or $T

### Checks
0) Kies deeltje p

1) Kies richting, als neighbor in die richting, roteer om die neighbor n.

2) Roteer rond n in random richting, 1 stapje.

3) Als er al een deeltje zit, reject.

4) Update alle neighbors

5) Verwijder kruisingen

6) Check of alle voorgaande neighbors van p nog steeds minstens 1 neighbor hebben. Anders reject.

7) Check afstanden tussen deeltje p en alle deeltjes in de buurt. Moet ofwel 1/2 ofwel >= 1 zijn. Anders reject.

8) Bereken nieuwe perimeter length.

9) Check alle neighboring tiles van p. Check dat hun omtrek 3, 4, of gelijk aan perimeter length is.

10) Check dat alle oude neighbors bezocht worden tijdens check (9). Als er een neighbor is waarvoor dat niet zo is, check of die een neighbor deelt met p. Als dat ook faalt, reject.

11) Check acceptance rule: accept or reject.

### Area fractions:
The hyperslopes for the different tiles and their orientations are given by
$$
B_{S_1} =
\begin{pmatrix}
1 & 0\\
0 & -1
\end{pmatrix}
\qquad
B_{S_2} =
\begin{pmatrix}
-1/2 & -\sqrt{3}/2\\
-\sqrt{3}/2 & 1/2
\end{pmatrix}
\qquad
B_{S_3} =
\begin{pmatrix}
-1/2 & \sqrt{3}/2\\
\sqrt{3}/2 & 1/2
\end{pmatrix}
$$

$$
B_{R_1} =
\begin{pmatrix}
1 & -2\sqrt{3} \\
0 & -1
\end{pmatrix}
\qquad
B_{R_2} =
\begin{pmatrix}
-2 & \sqrt{3} \\
-\sqrt{3} & 2
\end{pmatrix}
\qquad
B_{R_3} =
\begin{pmatrix}
1 & 0 \\
2\sqrt{3} & -1
\end{pmatrix}
$$

$$
B_{R_4} =
\begin{pmatrix}
1 & 0 \\
-2\sqrt{3} & -1
\end{pmatrix}
\qquad
B_{R_5} =
\begin{pmatrix}
-2 & -\sqrt{3} \\
\sqrt{3} & 2
\end{pmatrix}
\qquad
B_{R_6} =
\begin{pmatrix}
1 & 2\sqrt{3} \\
0 & -1
\end{pmatrix}
$$

$$
B_{T_1} =
\begin{pmatrix}
1 & 0\\
0 & 1
\end{pmatrix}
\qquad
B_{T_2} =
\begin{pmatrix}
-1 & 0\\
0 & -1
\end{pmatrix}
$$

$$
B_{T_3} =
\begin{pmatrix}
1 & 0\\
0 & 1
\end{pmatrix}
\qquad
B_{T_4} =
\begin{pmatrix}
-1 & 0\\
0 & -1
\end{pmatrix}
$$

We can then calculate the average hyperslope as a weighted average of the hyperslopes above
\begin{align*}
B &= \sum_{i=1}^3 \sigma_i B_{S_i} + \sum_{i=1}^4 \tau_i B_{T_i} + \sum_{i=1}^6 \rho_i B_{R_i}\\ \\
&=
\begin{pmatrix}
\sigma_1 - \frac{1}{2}(\sigma_2 + \sigma_3) & -\frac{\sqrt{3}}{2}(\sigma_2 - \sigma_3) \\
-\frac{\sqrt{3}}{2}(\sigma_2 - \sigma_3) & \frac{1}{2}(\sigma_2 + \sigma_3) -\sigma_1
\end{pmatrix} \\ \\
& \qquad +
\begin{pmatrix}
\tau_1 - \tau_2 + \tau_3 - \tau_4 & 0\\
0 & \tau_1 - \tau_2 + \tau_3 - \tau_4
\end{pmatrix} \\ \\
& \qquad +
\begin{pmatrix}
\rho_1 - 2\rho_2 + \rho_3 + \rho_4 -2\rho_5 + \rho_6 & \sqrt{3}(-2\rho_1 + \rho_2 - \rho_5 + 2\rho_6) \\
\sqrt{3}(-\rho_2 + 2\rho_3 - 2\rho_4 + \rho_5) & -\rho_1 + 2\rho_2 - \rho_3 - \rho_4 + 2\rho_5 - \rho_6
\end{pmatrix}
\end{align*}

Such that the determinant becomes
\begin{align*}
\det{B} = \quad&\rho_1(-\rho_1 - 2\rho_2 + 10\rho_3 - 14\rho_4 + 10\rho_5 - 2\rho_6 - 2\sigma_1 - 2\sigma_2 + 4\sigma_3) \\
+&\rho_2(-\rho_2 - 2\rho_3 + 10\rho_4 - 14\rho_5 + 10\rho_6 + 4\sigma_1 - 2\sigma_2 - 2\sigma_3) \\
+&\rho_3(-\rho_3 - 2\rho_4 + 10\rho_5 - 14\rho_6 - 2\sigma_1 + 4\sigma_2 - 2\sigma_3) \\
+&\rho_4(-\rho_4 - 2\rho_5 + 10\rho_6 - 2\sigma_1 - 2\sigma_2 + 4\sigma_3) \\
+&\rho_5(-\rho_5 - 2\rho_6 + 4\sigma_1 - 2\sigma_2 - 2\sigma_3) \\
+&\rho_6(-\rho_6 - 2\sigma_1 + 4\sigma_2 - 2\sigma_3) \\
-&\sigma_1^2 + \sigma_1\sigma_2 + \sigma_1\sigma_3 - \sigma_2^2 + \sigma_2\sigma_3 - \sigma_3^2 \\
+&\tau_1^2 - 2\tau_1\tau_2 + 2\tau_1\tau_3 - 2\tau_1\tau_4 + \tau_2^2 - 2\tau_2\tau_3 + 2\tau_2\tau_4 + \tau_3^2 - 2\tau_3\tau_4 + \tau_4^2.
\end{align*}
\begin{align*}
\det{B} &= \tau^2 - (\rho + \sigma)^2 + 12(\rho_1\rho_3 -\rho_1\rho_4 + \rho_1\rho_5 + \rho_2\rho_4 - \rho_2\rho_5 + \rho_2\rho_6 + \rho_4\rho_6 + \rho_3\rho_5 - \rho_3\rho_6 ) \\ & \quad + 6(\rho_1\sigma_3 + \rho_2\sigma_1 + \rho_3\sigma_2 + \rho_4\sigma_3 + \rho_5\sigma_1 + \rho_6\sigma_2) +3(\sigma_1\sigma_2 + \sigma_1\sigma_3 + \sigma_2\sigma_3) \\ &\quad - 4(\tau_1\tau_2 + \tau_1\tau_4 + \tau_2\tau_3 + \tau_3\tau_4)
\end{align*}

Imposing 12-fold symmetry, $\sigma_i = \sigma/3$, $\tau_i = \tau/4$ and $\rho_i = \rho/6$ for all $i$, we find that this determinant vanishes.

Alternatively, the determinant is given by the weighted average of the determinants of the individual tiles' hyperslopes,
$$
\det{B} = \sum_{i=1}^3 \sigma_i \det{B_{S_i}} + \sum_{i=1}^4 \tau_i \det{B_{T_i}} + \sum_{i=1}^6 \rho_i \det{B_{R_i}}.
$$
Plugging these in, given by
$$
\det{B_{S_i}} = \det{B_{R_i}} = -1,\\
\det{B_{T_i}} = 1,
$$
we find that
$$
\det{B} = \tau - \sigma - \rho,
$$
again having used the 12-fold symmetry. Since these two expressions for the determinant must be the same, we find the following restriction for the area fractions of our tile types
$$
\tau = \sigma + \rho
$$

### Thermodynamic integration:
The free energy of our system is defined as
$$
F = -k_B T \log{Z}.
$$
Such that the derivative w.r.t. to the defect energy cost $\epsilon_R$ is given by
$$
\frac{d F}{d \epsilon_R} = -k_B T \frac{d}{d\epsilon_R}\log{\sum_m e^{-\beta U}} = \frac{\sum_m \frac{d U}{d \epsilon_R}e^{-\beta U}}{Z} = \biggl\langle \frac{dU}{d\epsilon_R} \biggr\rangle.
$$
Here we defined the partition function $Z$ as a sum of the Boltzmann factors over all the microstates labelled $m$, and $\langle \cdot \rangle$ denotes an ensemble average over these microstates. For our interaction potential,
$$
U = \gamma L + \epsilon_Rn_R,
$$
we get
$$
\frac{d F}{d \epsilon_R} = \langle n_R \rangle.
$$

Thermodynamic integration now lets us estimate the free energy at an arbitrary state $F(\epsilon_R)$ by starting at a state with a known free energy. In our case we can start from $F(\infty)$, as this correpsonds to the free energy of a pure square-triangle tiling, which is known. [source]
$$
F(A) - F(B) = \int_B^Ad\epsilon'_R \frac{d F}{d \epsilon'_R},\\
F(\epsilon_R) = F_\text{sq-tr} - \int_{\epsilon_R}^{\infty}d\epsilon'_R \langle n_R \rangle.
$$
We can estimate this integral numerically as a sum,
$$
\int_{\epsilon_R}^{\infty}d\epsilon'_R \langle n_R \rangle \approx \Delta\epsilon'_R\sum_{i=1}^N \langle n_R \rangle_i,
$$
where we do $N$ simulation runs at different values of $\epsilon'_R$, with spacing $\Delta \epsilon'_R$.

### Area fraction constraints for systems with "shields" and "eggs"
Analogous to before. First calculating the global average hyperslope

$$
B = \sum_{i=1}^3 \sigma_i B_{\text{Sq}_i} + \sum_{i=1}^4 \tau_i B_{\text{Tr}_i} + \sum_{i=1}^4 \varsigma_i B_{\text{Sh}_i} + \sum_{i=1}^6 \varepsilon_i B_{\text{Eg}_i}
$$
Calculating the determinant of this resulting matrix, we find that it vanishes after demanding 12-fold symmetry.
The determinants of the new tiles are given by
$$
\det{\text{Sh}_i} = -\det{\text{Eg}_i} = 7 - 4\sqrt{3} \qquad \forall i
$$
Which gives an average determinant of
$$
\det{B} = -\sigma + \tau + (4\sqrt{3}-7)(\varepsilon - \varsigma)
$$
Setting this equal to zero, and demanding that all area fractions add to one, we find
$$
\varsigma = \frac{1}{2} + (2\sqrt{3} + 3)\sigma - (2\sqrt{3} + 4)\tau \\
\varepsilon = \frac{1}{2} - (2\sqrt{3} + 4)\sigma + (2\sqrt{3} + 3)\tau
$$

## Alptug's trials to solve the quenching problem

All simulations start from $\gamma = 3 kT$ and $\epsilon_i = 0$

| Target $\epsilon_{T}$ | McSweeps | Hexagonal phase @ $\epsilon_{T}$ | Hexagonal phase @ MCSweep |
| --------------------- | -------- | -------------------------------- | ------------------------- |
| -2 | 1M | -1.36 | N/A |
| -1.36 | 1M | -1.254 | 920K |
| -1.25 |1M | | |