\section{\large{Supplementary Material for ``Competing Crystal Phases in the Lowest Landau Level''}}

\section{S1: Free energies of the crystal and liquid states}

The free energy of the type-1 $^{2p}$CFC at small temperatures is given by [\onlinecite{pr1}]:
\begin{equation}\label{eq:fws1}
 F^{\rm CFC}=E^{\rm CFC}-0.701\left(\frac{C^{\rm t}_{\rm classical}}{C^{\rm t}_{\rm CF}}\frac{2a_0}{\nu e^2/\epsilon}\right)^{4/3}\left(k_{\rm B}T\right)^{7/3}
\end{equation}
where $E^{\rm CFC}$ is the $T=0$ energy per particle of the CFC and $a_0$ is the lattice spacing. The free energy of the FQHE liquid is given by  [\onlinecite{pr1}]:
\begin{equation}
 F^{\rm FQHE}=E^{\rm FQHE}-(2\pi m_{\rm R}\ell^2)^{1/2}\frac{k_{\rm R}\ell}{\nu}(k_{\rm B}T)^{3/2}e^{-\Delta_{\rm R}/k_{\rm B}T}
\end{equation}
where $E^{\rm FQHE}$ is the $T=0$ energy per particle of the FQHE state with a type-2 $^4$CFC, where it is assumed that the temperature is sufficiently small that only the magnetoroton [\onlinecite{Girvin}] part of the dispersion is relevant; $\Delta_{\rm R}$ and $m_{\rm R}$ are the energy and mass of the magnetoroton, and $k_{\rm R}$ is the wave vector at the energy minimum. In our calculation of the temperature where the $^2$CF crystal melts into a FQHE state at $\nu=0.213$, we have assumed that the type-2 CFC is pinned and the lowest energy excitations are the magnetorotons of the 1/5 state. 
The parameters $\Delta_{\rm R}$, $m_{\rm R}$ and $k_{\rm R}$ are obtained by performing fits to the CF exciton dispersion in the vicinity of the lowest energy magnetoroton as shown in Fig.~\ref{fig:fit}. For the crystal states at $\nu=1/5$ and 1/7, we have $C^{\rm t}_{\rm CF}/C^{\rm t}_{\rm classical}\approx1.2$.


\section{S2: System size dependence of energies}
In Fig.~\ref{fig:totE} we show the energies of type-1 $^{2p}$CFCs as a function of the system size. For clarity, we mark the phase boundaries between the crystals with a thin vertical line and the filling factor of the transition. This figure shows that the type-1 $^{2p}$CFC phase boundaries are well converged for system sizes greater than 32 particles.

Fig.~\ref{fig:ediff} shows the energy difference between the type-1 $^2$CFC and the FQHE state with a type-2 $^4$CFC in the filling factor range $1/5\leq\nu\leq2/9$ for three different system sizes. (The typical error in the energy difference from Metropolis sampling is $\sim 3\times10^{-5} e^2/\epsilon \ell$.) This shows that the results for 96 particles accurately represent the thermodynamic limit. 

\section{S3: Calculating the Shear Modulus of the $^{2p}$CF crystals}

In this section we first show how the shear modulus can be obtained from the filling factor dependence of the energy per particle of the hexagonal lattice, under the assumptions that the dynamics of the crystal can be described in the harmonic approximation and that the total energy can be written as a sum of two-body interactions. The derivation proceeds along two steps. In the first step, we show how to calculate the shear modulus from the dynamical matrix. In the second step, we show that the second derivative of the energy per particle with respect to the filling factor is proportional to the shear modulus.  We then provide details of our procedure for obtaining the shear modulus from our finite system results for which the energy is known only at discrete values of $\nu$.

In what follows, we will set the magnetic length $\ell=1$. The variation in filling factor $\nu$ should be thought of in terms of a change in density while keeping the magnetic field $B$ constant. The hexagonal lattice sites are located at
\begin{equation}
\vec{R}_{n,m}=\sqrt{\frac{4\pi}{\sqrt{3}\nu}}\left(n+{m\over 2},\frac{\sqrt{3}}{2}m\right).
\end{equation}

The dynamical matrix, $\Phi_{\alpha\beta}(\vec{k})$, of the classical 2D hexagonal system exposed to a perpendicular magnetic field is given by\cite{bm,mz}:
\begin{equation}\label{eq:dynmat}
 \displaystyle\Phi_{\alpha\beta}(\vec{k})=\sum_{\vec{R}\neq0}\left(1-\cos(\vec{k}\cdot\vec{R})\right)\frac{\partial^2V(R)}{\partial R_{\alpha} \partial R_{\beta}}=\left(\frac{\nu}{k}+C^{\rm L}-C^{\rm t}\right)k_{\alpha}k_{\beta} + C^{\rm t}k^2\delta_{\alpha\beta} + O(k^3)
\end{equation}
where the far right hand side of Eq.~(\ref{eq:dynmat}) gives the small $k$ limit, $R=|\vec{R}|$, $C^{\rm L}$ and $C^{\rm t}$ are the elastic constants which characterize the hexagonal lattice and $C^{\rm t}$ is the shear modulus. One can extract the shear modulus directly from the lattice sum in Eq.~(\ref{eq:dynmat}):
\begin{equation}
\displaystyle \Phi_{xx}((0,k_y))=\sum_{\vec{R}\neq0}\left(1-\cos(k_yR_y)\right)\frac{\partial^2V(R)}{\partial R_x \partial R_x}=k_y^2C^{\rm t}
\end{equation}
The long range interaction between $^{2p}$CFCs is Coulombic, while their short-range behavior gives rise to interesting elastic properties. In order to capture their short range behavior and to eliminate convergence issues, we define $V'\equiv V-V_{\rm MZ}$:
\begin{equation}
 \displaystyle C^{\rm t}_{\rm CF}=\frac{1}{k_y^2}\sum_{\vec{R}\neq0}^{\vec{R}_{\rm upp}}\left(1-\cos(k_yR_y)\right)\frac{\partial^2V'(R)}{\partial R_x \partial R_x} + C^{\rm t}_{\rm MZ}
\end{equation}
where $\vec{R}_{\rm upp}$ is chosen to evaluate the sum to arbitrary accuracy, and $V_{\rm MZ}$ is the MZ interaction given by\cite{mz} 
\begin{equation}
V_{\rm MZ}(R)=\frac{\sqrt{\pi}}{4}\frac{I_0(R^2/8)}{\cosh(R^2/8)}
\end{equation}
which is the Coulomb energy of two properly antisymmetrized gaussian wave packets in the lowest Landau level at distance $R$. 
Next, we expand $\cos({k_yR_y})$ in a Taylor series around $k_y=0$ and then take the limit $k_y\rightarrow0$:
\begin{equation}\label{eq:sm3}
 C_{\rm CF}^{\rm t}=\frac{1}{2}\displaystyle\sum_{\vec{R}\neq0}^{\vec{R}_{\rm upp}}(R_y)^2\frac{\partial^2V'(R)}{\partial R_x \partial R_x} + C^{\rm t}_{\rm MZ}
\end{equation}
Using the following identities (which rely on the hexagonal symmetry of the lattice):
\begin{equation}
 \frac{\partial^2}{\partial R_x^2}=\frac{(R_y)^2}{R^3}\frac{\partial}{\partial R} + \frac{(R_x)^2}{R^2}\frac{\partial^2}{\partial R^2}
\end{equation}
\begin{equation}
 \displaystyle\sum_{\vec{R}\neq0}(R_y)^4f(R)=\frac{9}{24}\sum_{\vec{R}\neq0}R^4f(R)
\end{equation}
\begin{equation}
 \displaystyle\sum_{\vec{R}\neq0}(R_y)^2(R_x)^2f(R)=\frac{3}{24}\sum_{\vec{R}\neq0}R^4f(R)
\end{equation}
we can rewrite Eq.~(\ref{eq:sm3}) as:
\begin{equation}\label{eq:sm4}
 C^{\rm t}_{\rm CF} = \frac{1}{16}\displaystyle\sum_{\vec{R}\neq0}\left(3R\frac{\partial V'}{\partial R}+R^2\frac{\partial^2 V'}{\partial R^2}\right) + C^{\rm t}_{\rm MZ}
\end{equation}
This completes the first step.


We now show explicitly how to calculate the shear modulus from the energy per particle. We note the following identities (when applied on a function of $R$):
\begin{equation}
 \frac{\partial}{\partial \nu}=\frac{\partial R}{\partial \nu}\frac{\partial}{\partial R}=-\frac{1}{2\nu}R\frac{\partial}{\partial R}
\end{equation}
\begin{equation}
 \frac{\partial^2}{\partial \nu^2}
 =\frac{1}{\nu^2}\left(\frac{3}{4}R\frac{\partial}{\partial R} + \frac{1}{4}R^2\frac{\partial^2}{\partial R^2}\right)\label{eq:finr}
\end{equation}
where we have used that $R\sim \nu^{-1/2}$.
Applying Eq.~(\ref{eq:finr}) to the energy per particle, we arrive at the final result:
\begin{equation}\label{eq:cfsm}
 \displaystyle \frac{1}{2}\nu^2\frac{\partial^2}{\partial\nu^2}(E_{\rm CF}-E_{\rm MZ})+C^{\rm t}_{\rm MZ}=\frac{1}{4}\nu^2\frac{\partial^2}{\partial\nu^2}\sum_{\vec{R}\neq0}V(R)'+C^{\rm t}_{\rm MZ}=C^{\rm t}_{\rm CF}
\end{equation}
Here, $E_{\rm MZ}$ is defined as the energy of a hexagonal crystal of {\em classical} particles interacting with the MZ interaction $V_{\rm MZ}$.  

This relation can be explicitly verified for the shear modulus of a hexagonal lattice of classical particles. The classical energy per particle for the 2d planar hexagonal lattice\cite{bm,mz} is given by $E=-0.782133\nu^{1/2}$, the shear modulus is given by $C^{\rm t}_{\rm classical}=0.0978\nu^{1/2}$ and the two are related by:
\begin{equation}
 \frac{1}{2}\nu^2\frac{\partial^2}{\partial\nu^2}\left(-0.782133\nu^{1/2}\right)=0.0978\nu^{1/2}
\end{equation}
We have also tested that the energies obtained by MZ\cite{mz} in which they assume two-body interaction, and confirmed that Eq.~(\ref{eq:cfsm}) exactly reproduces the shear modulus.



We next come to the issue of extracting the shear modulus from the discrete data points available to us from our finite system studies. For this purpose, we fit a smooth analytical curve to the energy per particle and calculate the shear modulus using Eq.~(\ref{eq:cfsm}). 
We show our smooth analytical fits to the $^{2p}$CFC energy relative to the Maki-Zotos energy $E_{\rm MZ}$ and their corresponding shear moduli for $N=128$ particles in Fig.~\ref{fig:ct}. Here, the energy $E_{\rm MZ}$ is given by the total energy of 128 classical particles located at the Thomson minima on the sphere that interact through the Maki-Zotos two-body interaction $V_{\rm MZ}$. The black curve in the right column panels of Fig.~\ref{fig:ct} is the shear modulus calculated from the strictly two-body Maki-Zotos interaction and is included for reference. 

To fit the CF energy, we use the functional form for the curves:
\begin{equation}
 \displaystyle f^{2p}_{\rm i}(\nu)=\left(\sum_{j=0}^{\rm m_{2p}+i}a_{2p,i}^j\nu^{3/2+j}\right)e^{-\nu^2}
\end{equation}
where $a_{2p,i}^j$ are found using a Least Squares fitting routine in Mathematica and $i$ is the fit number: $i=1,2,3$ or 4. The lines in Fig.~\ref{fig:ct} are drawn using these fit functions. In each fit function, referred to as `Fit $i$' in Fig.~\ref{fig:ct}, the degree of the polynomial is given by $m_{2p}+i$. The value of $m_{2p}$ is chosen separately for each crystal to give the best fit result. We perform fits of varying degree in order to check the consistency of our results, and take as our final result the average of four fits (given by the solid purple line in the right column panels of Fig.~\ref{fig:ct}). The weak oscillations in the shear modulus calculations are a result of the high degree of the polynomial fit. We note that the oscillations in the shear modulus fits are much smaller than the large overall structure that we observe near the transitions between $^{2p}$CFCs, implying that the latter is not an artifact of the fitting procedure. We have performed fits over the entire filling factor range that each $^{2p}$CFC exists, but have only shown the fits and their corresponding shear modulus for $\nu^*<1/2$ (or $\nu<{1\over 2p+2}$).


As a benchmark, in Fig.~\ref{fig:ctp0} we compare our calculated shear modulus for the $^0$CFC, which is a Slater determinant of otherwise uncorrelated electron wave packets, with that obtained previously by MZ. The two curves are very similar, but not identical. The difference arises because while the $^0$CFC is a fully antisymmetrized wave function, the MZ calculation imposes only pairwise antisymmetry, which is a good approximation only for small $\nu$. We believe that the shear modulus obtained in our calculation is more reliable. Interestingly, while Maki-Zotos results shows an instability of the crystal for $\nu\approx 0.45$, where the shear modulus becomes negative, our calculation eliminates that instability.  Of course, both are subject to the approximation that neglects three and higher body terms in the effective model mapping it into a crystal of classical particles.
For other values of $2p$, both the modified interaction between composite fermions and the antisymmetrization contribute to the deviation from the MZ shear modulus. Given that the distance between the neighboring lattice sites is large compared to the magnetic length, we expect that the assumption that the interaction is dominated by 2-body terms is accurate for $\nu^*<1/2$ or $\nu<1/(2p+2)$. However, as $\nu^*$ begins to approach unity, this approximation clearly breaks down (at $\nu^*=1$ we do not even have a crystal, but a liquid). Fortunately, the physically relevant regions (i.e. where the shear modulus is shown by solid lines) correspond to  $\nu^*<1/2$, with the exception of the filling factor $1/6<\nu<1/5$.




\section{S4: Density distribution of CF crystals}

Recently, Muraki and collaborators \cite{Muraki} have found that the spectral line shape of the NMR spectra contains information about the electron density distribution, because the Knight shift is proportional to the local electron density. We plot in Fig.~\ref{fig:pdens} the normalized probability of the CFC density as a function of density 
evaluated for a system with 96 particles. 
These density distributions are obtained at $\nu=0.215$ for type-1 crystals with $2p=0$ and $2$, and the type-2 crystal in the background of the $1/5$ FQHE state. From these plots it is clear that while the density distributions of the two type-1 CFCs are quite similar, they both qualitatively differ from the type-2 CFC, and it may be possible to identify the phase boundaries separating the type-1 and type-2 CFCs in the region $1/5<\nu<2/9$ by NMR measurements using the analysis of Ref. [\onlinecite{Muraki}]; in particular, an abrupt change in the dominant Knight shift is predicted at the phase transition. 



\begin{figure}[h]
\begin{center}
\includegraphics[width=0.40\textwidth]{MR_fit_1-5.pdf}
\includegraphics[width=0.40\textwidth]{MR_fit_1-7.pdf}
\end{center}
\caption{Fits to the CF exciton dispersion at $\nu=1/5$ (top panel) and 1/7 (bottom panel) near its minimum energy.  The blue dots are the energies extracted from Ref. [\onlinecite{cfr1}] and the dashed red lines are a least squares fit whose equation is shown in each panel.}
\label{fig:fit}
\end{figure}

\begin{figure}[h]
\begin{center}
\subfloat[$N$=32.]
{
\includegraphics[width=0.37\textwidth]{totE_N=32_G.pdf}
\label{fig:totE32}
}
\subfloat[$N$=48.]
{
\includegraphics[width=0.37\textwidth]{totE_N=48_G.pdf}
\label{fig:totE48}
}

\subfloat[$N$=64.]
{
\includegraphics[width=0.37\textwidth]{totE_N=64_G.pdf}
\label{fig:totE64}
}
\subfloat[$N$=96.]
{
\includegraphics[width=0.37\textwidth]{totE_N=96_G.pdf}
\label{fig:totE96}
}
\end{center}
\caption{Energy per particle of the type-1 $^{2p}$CFCs for $2p=0,2,4,6,8$ and 10 for several different system sizes. Energies are quoted relative to $E_{\rm fit}=-0.782133\nu^{1/2}+0.2623\nu^{3/2}+0.18\nu^{5/2}-15.1e^{-2.07/\nu}$. The colored dots mark the actual data points and the lines, which are guides for the eyes, are spline interpolations. The thin vertical lines mark the phase boundaries between the type-1 CFCs. We have evaluated a higher density of data points in the filling factor regions with more rapid variations of energy.}
\label{fig:totE}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.37\textwidth]{ediff.pdf}
\end{center}
\caption{Plots of the energy difference of the type-1 $^2$CFC and the FQHE state supporting a type-2 $^4$CFC in the filling factor range $1/5\leq\nu\leq2/9$ for three different system sizes. The uncertainty in the data points is $3\times10^{-5}$.}
\label{fig:ediff}
\end{figure}


\begin{figure}[h]
\begin{center}
\subfloat[$E-E_{\rm MZ}$ for $2p=0$.]
{
\includegraphics[width=0.35\textwidth]{ediff_fits_2p=0.pdf}
\label{fig:efitp0}
}
\subfloat[$C^{t}$ for $2p=0$.]
{
\includegraphics[width=0.35\textwidth]{cfc_ct_fits_2p=0.pdf}
\label{fig:ctp0}
}

\subfloat[$E-E_{\rm MZ}$ for $2p=2$.]
{
\includegraphics[width=0.35\textwidth]{ediff_fits_2p=2.pdf}
\label{fig:efitp1}
}
\subfloat[$C^{t}$ for $2p=2$.]
{
\includegraphics[width=0.35\textwidth]{cfc_ct_fits_2p=2.pdf}
\label{fig:ctp1}
}

\subfloat[$E-E_{\rm MZ}$ for $2p=4$.]
{
\includegraphics[width=0.35\textwidth]{ediff_fits_2p=4.pdf}
\label{fig:efitp2}
}
\subfloat[$C^{t}$ for $2p=4$.]
{
\includegraphics[width=0.35\textwidth]{cfc_ct_fits_2p=4.pdf}
\label{fig:ctp2}
}

\subfloat[$E-E_{\rm MZ}$ for $2p=6$.]
{
\includegraphics[width=0.35\textwidth]{ediff_fits_2p=6.pdf}
\label{fig:efitp3}
}
\subfloat[$C^{t}$ for $2p=6$.]
{
\includegraphics[width=0.35\textwidth]{cfc_ct_fits_2p=6.pdf}
\label{fig:ctp3}
}
\end{center}
\caption{Panels in the left column show the energies per particle of the $^{2p}$CFCs relative to the Maki-Zotos energy and the lines show the smooth analytical fits made to the data (for details, see text). The panels in the right column show the shear moduli calculated from the smooth analytical fits (dashed curves); the solid purple curve is the average from the four fits shown. The shear modulus calculated using the two-body Maki-Zotos interaction, shown by solid black line, is included for reference. All energies are calculated for a system with $N=128$ particles.}
\label{fig:ct}
\end{figure}

\clearpage

\begin{figure}[h]
\includegraphics[width=0.55\textwidth]{N=96_2Q=445_v=0_215_density.pdf}
\caption{Plot of the normalized probability of the density vs. density within a Wigner-Seitz cell on the sphere for type-1 and type-2 CFCs at $\nu=0.215$. The black arrow is a delta function centered at a density of 1 and represents the probability of the uniform-density FQHE liquid.}
\label{fig:pdens}
\end{figure}

