# Elastic properties and thermodynamic anomalies of supersolids

Milan Rakic  Department of Physics, Blackett Laboratory, Imperial College
London, SW7 2AZ, UK  Andrew F. Ho  Department of Physics, Royal Holloway
University of London, Egham, Surrey TW20 0EX, UK  Derek K. K. Lee  Department
of Physics, Blackett Laboratory, Imperial College London, SW7 2AZ, UK

###### Abstract

We study a supersolid in the context of a Gross-Pitaevskii theory with a non-
local effective potential. We employ a homogenisation technique which allows
us to calculate the elastic moduli, supersolid fraction and other state
variables of the system. Our methodology is verified against numerical
simulations of elastic deformations. We can also verify that the long-
wavelength Goldstone modes that emerge from this technique agree with
Bogoliubov theory. We find a thermodynamic anomaly that the supersolid does
not obey the thermodynamic relation $\partial P/\partial
V|_{N}\,=-n\,\left(\partial P/\partial N|_{V}\right)$, which we claim is a
feature unique to supersolids.

##  I Introduction

A supersolid is a phase of matter that displays both crystalline order and
superfluidity in the form of non-classical rotational inertia (NCRI). A key
requirement is that both the continuous translational symmetry and the U(1)
global gauge symmetry of the system is spontaneously broken in the ground
state.

There have been several attempts to understand this quantum phase of matter
experimentally and to understand its properties theoretically. While it has
not been observed in bulk 4He [1, 2], there are hints that such a phase may
exist in the second monolayer of 4He on graphite [3] where NCRI has been
measured in a density regime near layer completion with an anomalous
temperature dependence of the specific heat capacity. More recent experiments
in ultracold dysprosium and rubidium atoms [4, 5, 6] have observed the
spontaneous breaking of continuous translational symmetry together with long-
range phase correlation.

The first considerations of a supersolid phase were made by Andreev and
Lifshitz [7] and Chester [8] who considered the possibility of a supersolid
phase in 4He. They argued that the superfluid fraction of a supersolid would
be reduced from 100% due to the coupling of the phonons of the crystalline
structure to the $U(1)$ phase of the condensate wavefunction. This was further
developed by Leggett [9].

Later attempts at theoretical work have taken a phenomenological symmetry-
based approach starting with the work of Nozières [10] and Dorsey et al. [11].
Son [12] has described how Galilean invariance puts constraints on the form of
the Lagrangian that describes the low-energy dynamics of supersolids.

There have also been approaches starting from microscopic Hamiltonians based
on Gross-Pitaevskii theory [13] and Bogoliubov theory [14]. It is natural to
ask whether these approaches give rise to the same predictions for the low-
energy properties of the supersolid system, such as elasticity and
superfluidity.

Moreover, there has been intensive study on the superfluid fraction and
excitation spectrum of the supersolid phase. However, there has been
comparatively little study on the elastic properties. In this paper we will
build on the homogenisation technique of Rica and co-workers [15, 16, 17] to
obtain an effective low-energy theory that agrees well with Bogoliubov theory
and numerical calculations. Importantly, the original formulation of the
homogenisation technique did not agree with the expected bulk modulus in the
superfluid phase, and lacked the internal symmetry constraints for a Cauchy
elastic tensor in the supersolid phase. The velocities of the long-wavelength
excitations in the previous work also did not agree with the results based on
Bogoliubov theory [14].

With some important corrections that will be derived in this paper, we will
show that the homogenisation approach can be reconciled with other techniques,
producing the predictions for all the elastic and superfluid properties of the
supersolid phase. More specifically, we provide a method with the Gross-
Pitaevskii approximation that calculates the elastic constants and superfluid
density (phase stiffness) which is valid for any particular ground-state
solution of a Gross-Pitaevskii Lagrangian. It should be noted that, in the
limit of zero temperature, where a homogeneous superfluid has a superfluid
fraction of 100%, a supersolid has a reduced superfluid fraction [9, 18]. It
is important to clarify that this is not an enhancement of the “normal
fraction” which is found in superfluids at non-zero temperatures and carries
entropy. Rather, it is indicative of a separate phase that spontaneously
breaks translational symmetry. In fact, it is a consequence of the fact that
the “phonons” of the crystalline structure couples to the $U(1)$ phase of the
condensate wavefunction. [19, 20, 21, 22, 23].

The outline of the paper is as follows. In section 2, we review the Gross-
Pitaevskii theory for Bose condensates and how a finite-range interaction can
give rise to a supersolid phase. In section 3, we lay the framework for the
homogenisation theory, specifying the deformation procedure and rigorously
defining all necessary steps to calculate a long-wavelength theory of the low-
lying excitations. Section 4 uses results from section 3 to derive elastic
constants analytically and finds that they agree with elastic constants that
we obtain numerically. In section 5, we will consider the additional $U(1)$
gauge and derive the supersolid fraction, as well as a coarse-grained
Lagrangian which is the analytic expansion around the ground-state of the
elastic strain and $U(1)$ fields. Section 6 uses the effective Lagrangian and
considers the additional energetic contributions of flow across the system,
and of work done by the strain on the surroundings. In doing so, we derive a
long-wavelength effective Lagrangian which describes the Goldstone modes of
the system. Section 7 solves for the excitation velocities and discusses some
of the key features of the theory. In section 8, we verify the excitation
velocities through use of Bogoliubov fluctuations in a Bloch theorem context,
and find excellent agreement between the two seemingly disparate theories. In
Section 9, we will discuss a thermodynamic anomaly in the compressibility that
we believe is unique to a supersolid.

We subsequently make the claim that since we have satisfied the symmetry
requirements for an effective theory imposed by Son [12], have verified our
elastic constants numerically, and have verified the velocities of the low-
lying excitations through two independent and seemingly unconnected
techniques; then we must have the correct effective Lagrangian for a
supersolid. Furthermore, the technique outlined in this paper can now be used
and applied directly to other more complicated and realistic systems.

##  II Meanfield Supersolids

In this paper, we study a Bose-Einstein condensate (BEC) of particles of mass
$m$ with a finite-range interaction $U(\mathbf{r})$ in a Gross-Pitaevskii
theory. The condensate wavefunction $\psi(\mathbf{r},t)$ is a complex-valued
function that can be written in number-phase representation (or Madelung form)
as $\psi(\mathbf{r},t)=\sqrt{\rho(\mathbf{r},t)}e^{i\phi(\mathbf{r},t)}$. The
Lagrangian of the system is

$\begin{split}\mathcal{L}=-\int_{\Omega}&\left[\hbar\rho\frac{\partial\phi}{\partial
t}+\frac{\hbar^{2}}{2m}\left(\rho(\nabla\phi)^{2}+\frac{1}{4\rho}(\nabla\rho)^{2}\right)\right.\\\
&\left.+\frac{1}{2}\,\rho(\mathbf{r})\int_{\Omega}U(\absolutevalue{\mathbf{r}-\mathbf{r}^{\prime}})\rho(\mathbf{r}^{\prime})\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime}\right]\mathop{}\\!\mathrm{d}\mathbf{r}\end{split}$
(1)

where $\Omega$ is the spatial domain of the system. In the absence of any
current or twisted boundary conditions, the phase $\phi$ is spatially uniform
in the ground state and can subsequently be set to zero. The ground state
density can be found by variational calculus and is given by

$\frac{\hbar^{2}}{4m}\left(\frac{(\nabla\rho)^{2}}{2\rho^{2}}-\frac{\nabla^{2}\rho}{\rho}\right)+\int_{\Omega}U(\mathbf{r}-\mathbf{r}^{\prime})\rho(\mathbf{r}^{\prime})\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime}=\mu\,.$
(2)

where $\mu$ is the chemical potential to enforce the constraint
$\int_{\Omega}\rho(\mathbf{r})\mathop{}\\!\mathrm{d}\mathbf{r}=N$, the total
number of particles in the system Equation 2 is a nonlinear equation which we
solve numerically via a Crank-Nicholson scheme. By evolving the Lagrangian in
imaginary time numerically, we obtain the (real) ground state wavefunction as
the steady-state solution, using a Runge-Kutta algorithm for the temporal
evolution. The numerics are performed on a 2D triangular grid in real space
with periodic boundary conditions in order to respect the hexagonal symmetry
of the spatially modulated condensate that we expect.

A necessary feature (as shown by Heinonen et. al [13]) for the spontaneous
development of a density wave in the condensate is the presence of an
interaction potential which at least has non-zero Fourier components in
$\mathbf{k}^{4}$ and upwards [13]. This rules out the use of the typical
contact potential for supersolid formation.

To illustrate this, consider a toy model of a soft-core repulsive interaction
with a finite-range of $a$ in two dimensions:

$U(\mathbf{r})=U_{0}\Theta(\absolutevalue{\mathbf{r}}-a)$ (3)

where $\Theta$ is the Heaviside step function. For this simple interaction,
the system is controlled by a single dimensionless parameter

$\Lambda=\frac{\pi U_{0}ma^{2}}{\hbar^{2}}na^{2}$ (4)

for a system of average number density $n$.

The phase diagrams for the soft-core system are described in detail in various
works [24, 25, 26]. We find that our results are in agreement, with small
$\Lambda$ corresponding to a superfluid phase and that for $\Lambda>37$ the
system spatially orders into the supersolid phase. This is in agreement with a
theoretical subcritical limit calculated by During et al [17]. This supersolid
phase is characterised by a triangular lattice in 2D, illustrated in Figure 1.

The bulk of this paper is concerned with an effective long-wavelength theory
for this supersolid phase where the Bose condensate is spatially modulated. We
will use these numerics to verify our formalism by calculating the elastic
constants which are dependent on the ground state.

Figure 1: A spatially modulated condensate. Left: real space density
distribution. Right: Fourier transform of density distribution implying non-
zero momenta in the condensate.

##  III Homogenisation for Elasticity

The spontaneously broken translational symmetry in a supersolid means that the
system must possess Goldstone modes which are the elastic modes. In this
section, we will build on the work of [15, 17] to derive the elastic response
of the system. The goal here is to consider an elastic deformation on the
ground state, e.g. stretching/shearing as described by some strain tensor, and
then to expand the Lagrangian to second order in said strain tensor and
recover the elastic moduli through the generalised Hooke’s law.

We follow the metholodology of Josserand et al [15], but our results have some
key differences. Namely, we believe some key terms were missed which skews the
elastic constants to ones that do not agree with simple calculations in the
superfluid limit and lead to physically inconsistent elastic moduli. More
specifically, the Cauchy elastic tensor does not obey internal symmetries
consistent with a triangular lattice, and possesses negative elastic moduli
which suggest an unstable state. Moreover, we conduct a different treatment of
the system with regard to the canonical ensemble and the degrees of expansion
in the strain tensor.

The analysis will take the following structure: we will first carefully define
the notion of an elastic deformation on the system and provide the basis for
the calculation. We then expand the Lagrangian to second-order in the elastic
strain tensor. We demand that the deformed system also be a ground state, and
therefore solve the Euler-Lagrange equations for the new system, which will
allow us to determine elastic constants.

We begin by describing a physical deformation by considering the displacement
of points in the material. The displacement can be written as

$\mathbf{r}^{\prime}=\mathbf{r}-\mathbf{u}(\mathbf{r})\,,$ (5)

moving points $\mathbf{r}$ (over some domain $\Omega$) in the undeformed
ground state to points $\mathbf{r}^{\prime}$ (with some domain
$\Omega^{\prime}$). This basis describes the deformed material in the lab
frame of an external observer. To clarify, this is considered an active
transformation which physically acts on the material and changes the domain of
the system in the lab frame. The basis $\mathbf{r}^{\prime}$ is a Cartesian
coordinate system in the lab frame of the external observer, but is a non-
Cartesian coordinate system in the frame of the material.

We then define an inverse passive transformation which maps the coordinate
system $\mathbf{r}^{\prime}$ to the coordinate system
$\mathbf{r}^{\prime\prime}$, which is a non-Cartesian basis in the material
frame that has the same domain as the inital $\mathbf{r}$ (domain $\Omega$).
This transformation is of the form

$\mathbf{r}^{\prime\prime}=\mathbf{r}^{\prime}+\mathbf{u}^{\prime}(\mathbf{r}^{\prime})\,,$
(6)

where the basis $\mathbf{r}^{\prime\prime}$ coincides identically with the
basis $\mathbf{r}$ in the material frame, i.e.
$\mathbf{r}^{\prime\prime}=\mathbf{r}$, but now the space has additional
components in its gradients, as well as a Jacobian associated with its volume
element. These extra curvatures are necessary to describe the full effect of
the displacement procedure on the system. We now relabel the basis
$\mathbf{r}^{\prime\prime}$ as $\mathbf{r}$, for convenience, where again we
must stress that even though the coordinate systems $\mathbf{r}$ and
$\mathbf{r}^{\prime\prime}$ coincide, the space is no longer Cartesian and has
an associated Jacobian and curvature. This procedure is illustrated
schematically in Figure 2.

Figure 2: A schematic of the deformation procedure. The leftmost image depicts
an undeformed system with a basis $\mathbf{r}$, which then undergoes an active
transformation given by some strain to become the middle system. The rightmost
image depicts the system after the passive transformation, such that the final
coordinate system $\mathbf{r}^{\prime\prime}$ coincides exactly with
$\mathbf{r}$, but now there are additional curvatures in the space.

As of now, the above transformation is completely general, but in order to
proceed we need to specify a particular class of transformation which allows
us to use linear response theory and elasticity theory. We make the
distinction that $\mathbf{u}({\mathbf{r}})$ is such that all gradients in
either space or time are small constants. More specifically, we set
$\partial_{i}u_{k}\sim\delta$ and $\partial_{t}u_{k}\sim\delta$, with $\delta$
an arbitrarily small constant that we choose to track up to
$\mathop{}\\!\mathcal{O}(\delta^{2})$. This is referred to as a uniform
deformation, in that the gradient of the deformation is constant everywhere
and is small, allowing us to treat it as a perturbative parameter.

We now refer to $u_{ik}\equiv\partial_{i}u_{k}$ as the strain tensor, and
endeavour to expand $\mathcal{L}$ in powers of $u_{ik}$. In order to employ
linear elastic theory we have to keep all terms up to
$\mathop{}\\!\mathcal{O}(\delta^{2})$ in the Lagrangian. We also need to
relate the gradient of the transformation $\mathbf{u}^{\prime}$ in the primed
deformed basis in the lab frame to the gradient of the transformation
$\mathbf{u}$ in the material basis. It is relatively easy to show that

$\partial_{i}^{\prime}u_{k}^{\prime}=\partial_{i}u_{k}+\partial_{i}u_{l}\partial_{l}u_{k}\equiv
u_{ik}+u_{il}u_{lk}$ (7)

up to $\mathop{}\\!\mathcal{O}(\delta^{2})$ in the strain tensor (n.b. we are
using the Einstein summation convention). This allows us now to express
derivatives in the primed (deformed lab frame) basis in terms of derivatives
in the unprimed (deformed material frame) basis to find

$\partial_{i}^{\prime}=\partial_{i}+\partial_{k}\left(u_{ik}+u_{il}u_{lk}\right)\quad.$
(8)

In order to simplify the expressions that we will derive in this paper, we
introduce notation for a set of key tensors that appear regularly.

$\begin{split}&\epsilon_{ik}=-\frac{1}{2}\left(u_{ik}+u_{ki}\right)\quad;\quad\Delta_{ik}=\frac{1}{2}u_{il}u_{kl}\quad;\\\
&\omega_{ik}=\frac{1}{2}u_{li}u_{lk}\quad;\quad\chi_{ik}=\frac{1}{2}\left(u_{li}u_{kl}+u_{lk}u_{il}\right)\end{split}$
(9)

Since the deformation actively changes the domain of the system, there is an
associated change of volume. We refer to the deformations described in (5) and
(6) to define the differential volume elements as

$\int_{\Omega^{\prime}}\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime}=\int_{\Omega}\mathcal{J}_{\mathbf{r}^{\prime}\rightarrow\mathbf{r}}\mathop{}\\!\mathrm{d}\mathbf{r}$
(10)

with the Jacobian

$\mathcal{J}_{\mathbf{r}^{\prime}\rightarrow\mathbf{r}}=\det(\frac{\partial
r_{i}^{\prime}}{\partial r_{k}})=1+\epsilon_{ll}+M_{iklm}\,u_{ik}u_{lm}\quad.$
(11)

where $M_{iklm}\equiv(\delta_{ik}\delta_{lm}-\delta_{im}\delta_{kl})/2$. A
negative $\epsilon_{ll}$ corresponds to a reduction in volume. In two
dimensions, $M_{iklm}u_{ik}u_{lm}=u_{xx}u_{yy}-u_{xy}u_{yx}$. The quantities
$\epsilon_{ll}$ and $M_{iklm}u_{ik}u_{lm}$ are sometimes called the first and
second strain invariants as they are independent of the basis in which the
strain tensor is written. We stress now that the coordinate system
$\mathbf{r}$ describes the deformed system in the frame of the material: it
has the same domain as the undeformed material but is now a non-Cartesian
space.

To understand how this deformation affects the Lagrangian, we need to
understand how the density changes as a function of the deformation. After the
deformation occurs, the particles in the system will reorganise themselves in
such a way as to minimise the total energy of the new configuration, leading
to a new density described by $\rho(\mathbf{r}^{\prime})$. The reorganisation
of particles due to a deformation can be expressed as the pre-deformation
density plus a component, $\tilde{\rho}$, which completely accounts for all
changes. The component $\tilde{\rho}$ covers both local density changes at
length scales within a unit cell of the supersolid and also includes the
change in average density due to total number conservation and a change in
volume.

We consider that this new density is written in the $\mathbf{r}^{\prime}$
basis. We stress that this is not the density in the material frame, but
rather the density in the lab frame of the external observer. We can then use
the passive transformation to write the density in the
$\mathbf{r}^{\prime\prime}$ basis, which is the material frame that now has
some curvature. Without loss of generality we can express this new density at
some point $\mathbf{r}^{\prime}$ as

$\begin{split}\rho(\mathbf{r}^{\prime})&=\rho_{0}(\mathbf{r}^{\prime\prime})+\tilde{\rho}(\mathbf{r}^{\prime\prime})\,,\int_{\Omega^{\prime}}\\!\\!\\!\rho(\mathbf{r}^{\prime})\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime}=\int_{\Omega}\\!\\!\\!\rho_{0}(\mathbf{r}^{\prime\prime})\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime\prime}\end{split}$
(12)

where $\rho_{0}$ is the pre-deformation ground-state density profile. As of
now this is simply a mapping to a scalar function that has been subject to an
active transformation.

We stress that from this point we drop the $\mathbf{r}^{\prime\prime}$
notation and refer to the coordinate as $\mathbf{r}$ where now the coordinate
system is non-Cartesian and has a Jacobian associated with it.

The general form of $\tilde{\rho}$ is determined by minimising the energy of
the deformed system, subject to particle conservation. For linear elastic
theory, we only need to consider only a perturbative strain on the system and
expect that, since the strain is perturbative, the density change
$\tilde{\rho}$ must be analytic in the strain tensor. We truncate the Taylor
expansion to $\mathop{}\\!\mathcal{O}(\delta^{2})$ in the strain tensor and
use the ansatz

$\tilde{\rho}(\mathbf{r})=\rho_{1}^{ik}(\mathbf{r})\,\epsilon_{ik}+\rho_{2}^{iklm}(\mathbf{r})\,u_{ik}u_{lm}$
(13)

which we will justify later. This allows us to expand the new normalisation
condition for particle conservation.

$\int_{\Omega}(\rho_{0}+\tilde{\rho})\left(1+\epsilon_{ll}+M_{iklm}\,u_{ik}u_{lm}\right)\mathop{}\\!\mathrm{d}\mathbf{r}=\int_{\Omega}\rho_{0}\mathop{}\\!\mathrm{d}\mathbf{r}$
(14)

and after collecting powers of the strain tensor obtain normalisation
requirements for the different terms in the expansion

$\begin{split}\int_{\Omega}\left(\rho_{1}^{ik}+\rho_{0}\delta^{ik}\right)\mathop{}\\!\mathrm{d}\mathbf{r}&=0\\\
\int_{\Omega}\left[\rho_{2}^{iklm}+\rho_{0}\left(M_{iklm}-\delta_{ik}\delta_{lm}\right)\right]\mathop{}\\!\mathrm{d}\mathbf{r}&=0\,.\end{split}$
(15)

Due to the changing energy of the system (and indeed the change in average
density), we can expect to see a change in the chemical potential $\mu$ too.
In a similar spirit as that of equation (13), we can analytically expand the
new chemical potential in the strain tensor as

$\mu=\mu_{0}+\mu_{1}^{ik}\,\epsilon_{ik}+\mu_{2}^{iklm}\,u_{ik}u_{lm}\quad.$
(16)

where the first term is the chemical potential of the undeformed system, and
the constants $\mu_{1}^{ik}$ and $\mu_{2}^{iklm}$ encodes how the chemical
potential changes due to the deformation and are to be determined by the
least-action principle.

We additionally need to consider how the displacement field $u_{ik}$ changes
the interaction potential $U$. We note that the interaction is a function of
separation between particles, but importantly is originally written in terms
of the separation $\Delta r^{\prime}$ as measured in the real lab frame
coordinates: $U=U(\Delta r^{\prime})$ . For a uniform strain, we find
immediately that the lab-frame separation $\Delta\mathbf{r}^{\prime}$ can be
written in terms of the material-frame separation $\Delta\mathbf{r}$ such that

$\left(\Delta
r^{\prime}\right)^{2}=\left(\delta_{ik}+2\epsilon_{ik}+2\Delta_{ik}\right)\left(\Delta
r\right)_{i}\left(\Delta r\right)_{k}$ (17)

which is exactly equivalent to that of Landau & Lifshitz [27]. The tensor in
the above equation is sometimes referred to as the finite strain tensor. A
Taylor expansion of the interaction for small strain gives

$\begin{split}U(\Delta\mathbf{r}^{\prime})&=U(\Delta\mathbf{r})+\left(\epsilon_{ik}+\Delta_{ik}\right)f_{ik}(\Delta\mathbf{r})\\\
&\quad+\epsilon_{ik}\epsilon_{lm}W_{iklm}(\Delta\mathbf{r})\end{split}$ (18)

with

$\begin{split}f_{ik}(\mathbf{r})&\equiv\frac{r_{i}r_{k}}{|\mathbf{r}|}\frac{\partial
U(|\mathbf{r}|)}{\partial|\mathbf{r}|}\\\
W_{iklm}(\mathbf{r})&\equiv\frac{r_{i}r_{k}r_{l}r_{m}}{2|\mathbf{r}|^{2}}\left(\frac{\partial^{2}U}{\partial|\mathbf{r^{2}}|}-\frac{1}{|\mathbf{r}|}\frac{\partial
U}{\partial|\mathbf{r}|}\right)\;.\end{split}$ (19)

To summarise, we have now defined a perturbative deformation which changes the
geometry of our system, and have obtained the volume changes and curvature
changes due to said deformation. We have also derived a new density and
chemical potential, both of which we will solve for using the least-action
principle. We now have all the ingredients to formulate a theory of linear
response of the system to elastic deformations.

##  IV Elasticity Theory

We can now develop a linear elastic theory by expanding the Lagrangian to
second order in the strain tensor and finding the response, $\tilde{\rho}$, of
the density to a small applied strain using the principle of least action. The
full details are given in Appendix A where we calculated and collected all the
first- and second-order components of $u_{ik}$ in the Lagrangian. This should
be of the form:

$L(\Omega^{\prime})\simeq
L_{0}(\Omega)+|\Omega|\pi_{ik}u_{ik}-\frac{|\Omega|}{2}A_{iklm}u_{ik}u_{lm}$
(20)

where the first term is the ground-state Lagrangian, the second term contains
the coupling to the stress tensor $\pi_{ik}$, and the third term is related to
the Cauchy elastic tensor $B_{iklm}$, as carefully discussed in Bavaud et al.
[28]. ($|\Omega|$ is the total volume of the pre-deformed system.)

Let us consider the term linear in the strain tensor. The pressure tensor is
related to the stress tensor $\pi_{ik}$ by $P_{ik}=-\pi_{ik}$. We will be
studying a two-dimensional supersolid whose spatial modulation forms a
triangular lattice. The $C_{3}$ symmetry of the ground state dictates that
$P_{ik}=P\delta_{ik}$, i.e. the pressure is isotropic. Collecting the
$\mathop{}\\!\mathcal{O}(u_{ik})$ terms in the expansion of the Lagrangian
[see (75) and (76)], we can identify the pressure as

$\begin{split}P_{ik}=\int_{\Omega}&\left[\frac{\hbar^{2}}{4m}\left(\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}}-\nabla^{2}\rho_{0}\,\delta_{ik}\right)\right.\\\
&\left.\qquad\qquad\qquad-\frac{1}{2}\rho_{0}(f_{ik}*\rho_{0})\right]\frac{\mathop{}\\!\mathrm{d}\mathbf{r}}{|\Omega|}\end{split}$
(21)

which will be used in subsequent calculations. The above expression can be
directly calculated for any ground state density $\rho_{0}$. In the case of a
spatially homogeneous superfluid, it leads to what we expect:
$P_{ik}=\frac{1}{2}\mu_{0}n\delta_{ik}$. (See Appendix A.)

Let us now turn to the the second-order terms in the expansion (20). It is
important to note that the tensor $A_{iklm}$ is not the Cauchy elastic tensor
$B_{iklm}$ which provides the tensor in the Hooke’s law for the stress/strain
relationship between a deformed state stress and the strain that induced the
deformation, i.e.
$\pi_{ik}(\Omega^{\prime})=\pi_{ik}(\Omega)-B_{iklm}(\Omega)u_{lm}$. It can be
shown [28] that:

$B_{iklm}=A_{iklm}+P_{ik}\delta_{lm}-P_{im}\delta_{kl}$ (22)

which we will derive independently later in (VI) by considering work done by
the expansion on the surroundings.

For a supersolid with $C_{3}$ symmetry, we also expect the elastic tensor to
obey the indicial symmetries $B_{iklm}=B_{klim}=B_{lkmi}=B_{mlki}$ so that it
only contains two independent quantities — the bulk modulus $K=c_{44}$ and the
shear modulus $G=c_{66}$:

$B_{iklm}=K\delta_{ik}\delta_{lm}+G\left(\delta_{il}\delta_{km}+\delta_{im}\delta_{kl}\right)\,.$
(23)

It should be noted that $A_{iklm}$ is not indicially symmetric in the same way
as $B_{iklm}$. Elastic theory [27] in 2D posits that the bulk modulus $K$ is
given by $B_{xxyy}$, and the shear modulus $G$ is given by $B_{xyxy}$ (or all
equivalent indicially symmetric elements). For a triangular lattice or a
homogeneous system, the tensor also possesses the symmetry:
$B_{xxxx}=B_{xxyy}+2B_{xyxy}$. Note also that we expect the changes to the
chemical potential to have the same symmetry as the pressure, i.e.
$\mu_{1}^{ik}=\mu_{1}\delta_{ik}$. We will later make use of these properties
to simplify the expressions for the elastic constants.

We can deduce the elastic tensor by identifying the
$\mathop{}\\!\mathcal{O}(u_{ik}^{2})$ terms in the expansion of the Lagrangian
with $A_{iklm}$, and then using (22) to obtain $B_{iklm}$. (See Appendix (88)
and the derivative rules in Appendix (C)). We find that
$B_{iklm}=\absolutevalue{\Omega}^{-1}\int_{\Omega}b_{iklm}(\mathbf{r})\mathop{}\\!\mathrm{d}\mathbf{r}$
where

$\begin{split}b_{iklm}(\mathbf{r})=&+\frac{\hbar^{2}}{4m}\left(\frac{\partial_{i}\rho_{0}\partial_{l}\rho_{0}}{\rho_{0}}\delta_{km}+\frac{\partial_{k}\rho_{0}\partial_{l}\rho_{0}}{\rho_{0}}\delta_{im}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}}\delta_{lm}\right)\\\
&+\frac{1}{2}\left(f_{il}*\rho_{0}\right)\rho_{0}\delta_{km}-\frac{1}{2}\left(f_{ik}*\rho_{0}\right)\rho_{0}\delta_{lm}+\frac{1}{2}\left(f_{im}*\rho_{0}\right)\rho_{0}\delta_{kl}\\\
&+\left[\left(\left(U\delta_{ik}+2f_{ik}\right)\delta_{lm}+W_{iklm}\right)*\rho_{0}\right]\rho_{0}+\rho_{1}^{lm}\left[\frac{\hbar^{2}}{4m}\left(2\frac{\partial_{ik}\rho_{0}}{\rho_{0}}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}^{2}}\right)+\left(f_{ik}+\delta_{ik}U\right)*\rho_{0}+\mu_{1}^{ik}\right]\end{split}$
(24)

with
$(g*\rho)_{\mathbf{r}}\equiv\int\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime}g(\mathbf{r}-\mathbf{r}^{\prime})\rho(\mathbf{r}^{\prime})$
for any function $g$ and $\rho$ and where $f_{ik}$ and $W_{iklm}$ are defined
in (19).

We note that this result for the elastic constants depends on the ground state
density $\rho_{0}$ as well as the shift in the density $\rho_{1}$ and shift in
the chemical potential $\mu_{1}$ as defined in (13) and (16). In other words,
we need to deduce how the density profile has changed in response to the
applied strain. This is achieved by solving the Euler-Lagrange equations for
the variable $\tilde{\rho}$ as defined in (12). This involves solving the
following equation

$\begin{split}&\frac{\hbar^{2}}{4m}\gradient\cdot\left(\frac{\gradient\tilde{\rho}}{\rho_{0}}\right)-U*\tilde{\rho}\\\\[5.0pt]
+&\frac{\hbar^{2}}{4m}\left(\frac{(\nabla\rho_{0})^{2}}{\rho_{0}^{3}}-\frac{\nabla^{2}\rho_{0}}{\rho_{0}^{2}}\right)\tilde{\rho}\end{split}\text{\Large{\hbox{
=
}}}\begin{split}&-\left(\mu_{2}^{iklm}+\mu_{1}^{ik}\delta_{lm}+\mu_{0}M_{iklm}\right)u_{ik}u_{lm}+\epsilon_{ik}\left[-\mu_{1}^{ik}+\frac{\hbar^{2}}{4m}\left(2\frac{\partial_{ik}\rho_{0}}{\rho_{0}}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}^{2}}\right)\right.\\\\[5.0pt]
&\left.+\int_{\Omega}\left(f_{ik}(\mathbf{r}-\mathbf{r}^{\prime})+\delta_{ik}U(\mathbf{r}-\mathbf{r}^{\prime})\right)\rho_{0}(\mathbf{r}^{\prime})\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime}\right]\end{split}$
(25)

which is an integro-differential equation that can be solved either by using
Bloch’s theorem (only applicable in periodic boundary conditions) or solving
directly via writing the linear operator as a matrix and using standard
iterative matrix solvers. More details on solving the equation and finding the
chemical potentials are provided in Appendix (79).

It must be noted that we actually have 8 equations to solve, as for each pair
of indices $(i,k)$ there is a different RHS source function, and for each of
these we need to solve the appropriate equation for both $\rho_{1}^{ik}$ and
$\mu_{1}^{ik}$. This generates 8 equations for a 2D system and 18 equations
for a 3D system. Once we have the solutions for both $\rho_{1}^{ik}$ and
$\mu_{1}^{ik}$, we can calculate the other functions in the expression for
(24) and obtain the elastic tensor $B_{iklm}$.

Let us check these results for the case of a homogeneous system, where the
modulation vanishes and we can check our results against the known results for
elastic moduli of a superfluid. We can see that the elastic tensor is greatly
simplified by the disappearance of the derivative terms
$\partial_{i}\rho_{0}$, and the chemical potential is analytically known:
$\mu=U*\rho_{0}$. We can find the convolutions $f_{ik}*\rho_{0}$ and
$W_{iklm}*\rho_{0}$, solve the equation (79) and use this to find that

$B_{iklm}=\mu n\delta_{ik}\delta_{lm}\quad.$ (26)

Using the definitions of $K$ and $G$ from the Cauchy elastic tensors (23), we
see here that we obtain the expected results of $K=\mu n$ and no shear
modulus: $G=0$.

We will now present the numerical calculation of these homogenisation results
for the elastic constants. The results are shown in Figure 3 (solid markers).
In our numerics, we see some small asymmetry, e.g. $B_{iklm}\approxeq
B_{lmik}$, in the results for the elastic tensor. The discrepancy with the
ideal form (23) decreases as we reduce the discretisation and increase the
system size in the numerics. With that in mind, we choose to define the bulk
modulus $K$ as the average of the values of $B_{xxyy}$ and $B_{yyxx}$, the
pressure $P$ as the average of $P_{xx}$ and $P_{yy}$, and the shear modulus
$G$ as the average of $B_{xyxy}$, $B_{xyyx}$, $B_{yxyx}$ and $B_{yxxy}$. We
then use the minimal and maximal values that the elastic moduli could have
taken in order to generate ‘error bars’ for the elastic constants. It should
be stressed that these error bars are not measures of some statistical error,
but rather a quantification of how much the elastic constants deviate due to
the crystal symmetry being only approximately triangular (as opposed to
exactly triangular) as a result of discretisation error in our numerics.

Figure 3: A comparison of the various elastic constants obtained via
homogenisation to those from numerics (homogenisation constants are shown in
solid markers, numerics in crosses). The top figure shows the bulk modulus and
shear modulus, the middle figure shows the pressure, and the bottom figure
shows the absolute value of $\mu_{1}$. The agreement between homogenisation
and numerics is very strong on the whole, but falls apart in the phase
transition region. This is likely due to numerical instabilities close to the
phase transition which affect the numerically obtained constants but not the
homogenisation constants.

We now turn to comparing these numerical results from homogenisation with a
direct calculation of the elastic constants by finding how the energy of the
ground state of the Gross-Pitaevskii equation changes when we deform the
simulation cell at fixed particle number (with periodic boundary conditions).
We begin with an undeformed ground-state obtained via numerically solving the
Gross-Pitaevskii with some given average density, area and interparticle
interaction. We then apply a strain to the system by changing the geometry of
the simulation cell by a strain tensor $u_{ik}$, whilst keeping the total
number of particles fixed, and solving the Gross-Pitaevskii equation for the
new geometry to find the deformed ground-state. For sufficiently small strain
$u_{ik}$, we expect that such a procedure will yield a ground state with a
Lagrangian of the form (20), where the tensor $u_{ik}$ is now a control
parameter in our computational experiment.

By numerically extracting the Lagrangian for a system which has been deformed
by some strain $u_{ik}$, we are able to find the pressure and elastic
constants by judicious selection of strains and subsequent fitting of the
coefficients. The full procedure is described in Appendix D, where we show the
extraction of elastic constants specifically for a 2D case in a triangular
geometry.

Our results are presented in Figure 3. These elastic constants should ideally
be the same as those obtained via homogenisation theory. In reality, they
suffer from numerical errors due to being unable to perform an arbitrarily
small strain, as well as convergence errors in the new ground-state when close
to the first-order phase transition between the homogeneous superfluid and the
supersolid. Scanning across the whole range of the interaction parameter
$\Lambda$, we find that, outside of a small region where the phase transition
occurs, the elastic constants obtained numerically are in excellent agreement
with those obtained analytically, with a typical discrepancy of 0.1%.

In Figure 3, we have also shown the shift in the chemical potential for unit
strain: $\mu_{1}^{ik}=\mu_{1}\delta_{ik}$.

In summary, we have developed an elastic theory for a supersolid using
homogeneration theory and calculated its elastic properties using the solution
of the Gross-Pitaevskii equation. Our homogenisation results agree well with
results from the numerical simulation of a strained supersolid.

##  V $U(1)$ Phase and Fractional Inertia

In the this section, we consider the application of a phase perturbation to
the meanfield supersolid, and analyse the coupling of the phase to the strain
field. We follow the metholodology of Josserand et al [15], but our results
have some key differences that we will point out.

We begin by specifying the nature of the phase perturbation. We consider
imposing a small superflow to the system in the form of a phase gradient. This
is applied to the supersolid in the deformed state, i.e the lab frame of the
deformed material. The system will respond to this long-wavelength
perturbation by generating phase fluctuations, $\tilde{\phi}$, at small
wavelengths within a unit cell of the supersolid. As with $\tilde{\rho}$ in
the previous section, $\tilde{\phi}$ is a correcting function to the
respective perturbation applied to the ground state. Namely, $\tilde{\rho}$ is
the function which allows the density to settle into the energetically optimal
configuration for a given strain, and $\tilde{\phi}$ is the function which
will allow the phase field to settle into the energetically optimal
configuration for a given applied phase gradient. We mathematically represent
this as

$\phi(\mathbf{r}^{\prime})=\phi_{0}(\mathbf{r}^{\prime})+\tilde{\phi}(\mathbf{r})$
(27)

where $\phi(\mathbf{r}^{\prime})$ is the total phase in the lab frame
$\mathbf{r}^{\prime}$, $\phi_{0}(\mathbf{r}^{\prime})$ is the applied
perturbative phase in the lab frame, and $\tilde{\phi}(\mathbf{r})$ is the
phase-correction provided by the material in the material frame of reference.
We express $\tilde{\phi}$ in terms of the material frame deformed coordinate
system $\mathbf{r}$ in keeping with linear response theory, similar to how we
expressed $\tilde{\rho}$ in terms of the material frame coordinate system
$\mathbf{r}$.

In order to apply perturbation theory to the phase field, we specify that
$\partial_{i}\phi_{0}(\mathbf{r}^{\prime})$ is small and constant within the
scope of the unit cell of the supersolid, in a similar way to the strain
tensor $u_{ik}$. We note that since the phase is applied in the deformed
frame, the total phase in the material frame can be expanded to obtain

$\phi(\mathbf{r}^{\prime})=\phi_{0}(\mathbf{r})-\mathbf{u}\cdot\gradient\phi_{0}(\mathbf{r})+\tilde{\phi}(\mathbf{r})$
(28)

which is the final expression for the phase field in the material frame of
reference. We can now proceed to expand the Lagrangian in terms of the phase
field, using the least action principle to find $\tilde{\phi}$ and the
associated energy.

Using the transformation rules for the Lagrangian (Jacobian, derivatives etc.)
and the $\phi$ expansion given by (28), the phase-dependent part of the
Lagrangian (1) can be expanded to $\mathop{}\\!\mathcal{O}(\tilde{\phi}^{2})$
as

$\begin{split}\mathcal{L}_{\phi}=&-\int_{\Omega}\Big{\\{}\hbar\rho_{0}\partial_{t}\phi_{0}\left(1+\gradient\cdot\mathbf{u}\right)\\\
&+\frac{\hbar^{2}\rho_{0}}{2m}\left[\left(\gradient\phi_{0}\right)^{2}+2\mathbf{A}\cdot\gradient\tilde{\phi}+(\gradient\tilde{\phi})^{2}\right]\Big{\\}}\mathop{}\\!\mathrm{d}\mathbf{r}\end{split}$
(29)

where

$\mathbf{A}=\gradient\phi_{0}+\left(\gradient\phi_{0}\cdot\gradient\right)\mathbf{u}+\frac{m}{\hbar}\partial_{t}\mathbf{u}\,.$
(30)

The vector field $\mathbf{A}$ is a kind of convective flow. First of all,
recall that the superfluid velocity is given by
$\mathbf{v}=(\hbar/m)\gradient{\phi_{0}}$. Also, we can define the material
derivative of a field as

$\frac{D}{Dt}=\partial_{t}+\mathbf{v}\cdot\gradient$ (31)

to describe the time derivative of a function which is co-moving with the
superfluid velocity field. We can rewrite $\mathbf{A}$ in term of these
variables:

$\mathbf{A}=\frac{m}{\hbar}\left(\mathbf{v}+\mathbf{v}\cdot\gradient\mathbf{u}+\partial_{t}\mathbf{u}\right)=\frac{m}{\hbar}\left(\mathbf{v}+\frac{D\mathbf{u}}{Dt}\right)\,.$
(32)

A point of subtlety is that the actual deformation we have applied is
$\mathbf{r}^{\prime}=\mathbf{r}-\mathbf{u}$, and so
$(\hbar/m)\mathbf{A}=\mathbf{v}-D(-\mathbf{u})/Dt$ is a relative velocity
between the superflow and the motion of the lattice deformations. This is
exactly the form predicted by Son [12] from symmetry based arguments. In fact,
all terms in $\mathcal{L}_{\phi}$ satisfy the Galilean symmetry requirements
imposed by Son. In this way we can show that we have microscopically satisfied
the necessary Galilean symmetry requirements imposed on any supersolid system
macroscopically.

Our phase Lagrangian $\mathcal{L}_{\phi}$ in (29) differs from the analogous
form in Josserand et al [15] by the existence of the term
$\hbar\rho_{0}\partial_{t}\phi_{0}\gradient\cdot\mathbf{u}$. This term is in
fact necessary to ensure that $\mathcal{L}_{\phi}$ satisfies Galilean symmetry
constraints as pointed out by Son [12]. We also note that this affects the
velocity of the excitations of the system, and we will see later in Section
VIII in that our results agree with the results obtained in Bogoliubov theory.

Figure 4: The reduction $\varrho/n$ in the superfluid fraction as a function
of interaction strength $\Lambda$, exhibiting a discontinuous jump at the
transition.

Let us now solve for the short-distance phase response $\tilde{\phi}$ within a
unit cell due to superflow and strain. We follow a similar procedure to Ref
[15]. By recognising that the externally imposed phase twist and strains
$\gradient\phi_{0}$ and $\partial_{i}u_{k}$ are defined as long-wavelength
relative to the unit cell, we can treat $\mathbf{A}$ as constant within a unit
cell. However, a uniform velocity with non-uniform density does not
necessarily obey local mass conservation in equilibrium. So we need a
correction to the phase field, $\tilde{\phi}$, to ensure that our system obeys
the continuity equation in the ground state:
$\gradient\cdot\left(\rho\gradient\phi\right)=0$. The Euler-Lagrange equation
for $\tilde{\phi}$ gives

$\tilde{\phi}=K_{i}A_{i}\quad\mbox{with}\quad\partial_{i}\rho+\mathbf{\gradient}\cdot\left(\rho\mathbf{\gradient}K_{i}\right)=0\,.$
(33)

This can be solved numerically in much the same way as (80). Physically, this
result expresses that the fact the phase twist (and therefore superflow) is
not uniform within a unit cell because it costs less kinetic energy to
introduce phase shifts in regions of low density compared to regions of high
density. Using this solution for $\tilde{\phi}$, it can be shown that the
phase-dependent Lagrangian can be written as:

$\begin{split}L_{\phi}=&\absolutevalue{\Omega}\frac{\hbar^{2}}{2m}\varrho_{ij}A_{i}A_{j}-\int_{\Omega}\Big{[}\hbar\rho_{0}\partial_{t}\phi_{0}\left(1+\gradient\cdot\mathbf{u}\right)\\\
&\qquad\qquad\qquad\qquad+\frac{\hbar^{2}\rho_{0}}{2m}\left(\gradient\phi_{0}\right)^{2}\Big{]}\mathop{}\\!\mathrm{d}\mathbf{r}\end{split}$
(34)

with

$\varrho_{ij}=-\frac{1}{\absolutevalue{\Omega}}\int_{\Omega}\rho_{0}\partial_{i}K_{j}\mathop{}\\!\mathrm{d}\mathbf{r}$
(35)

which can be interpreted as a correction to the superfluid phase stiffness
(see equation (37)). Figure 4 shows the evolution of this quantity as we
increase the interaction strength $\Lambda$ showing a discontinuous transition
at the superfluid-supersolid transition.

To restate the above, by carefully considering the dynamics of both an
applied/induced phase and a strain deformation (both of which are long-
wavelength relative to the unit cell), we calculate how the phase-dependent
dynamics of our system are dependent on the parameters $n$ and $\varrho_{ij}$
given by the form (34). $n$ is simply the density of particles in the system,
and corresponds to the phase-coherent BEC with components at both
$\mathbf{k}=0$ and $\mathbf{k}\neq 0$. We will later show that $\varrho_{ij}$
is the part of the system that corresponds to the pattern formation, i.e. it
only exists when $\psi_{\mathbf{k}\neq 0}\neq 0$. Consequently, we refer to
$\varrho_{ij}$ as being the “supersolid fraction” of the system. In this way
we have separated out the pattern response of the system to the coupled phase-
strain dynamical field from the standard superfluid response.

Finally, due to the $C_{3}$ symmetry, second-rank tensors representing bulk
properties should be isotropic. We therefore make the assumption that
$\varrho_{ij}=\varrho\delta_{ij}$, though in principle this is not necessary.
We also note that the dependency in (34) on $\phi_{0}$ and $\mathbf{u}$ can be
taken out of the integral (due to the same argument that $\mathbf{A}$ is a
constant in the unit cell). Combining these two simplifications, we find that
we can write

$\begin{split}\frac{L_{\phi}}{|\Omega|}=&-\hbar
n\dot{\phi}_{0}(1+\gradient\cdot\mathbf{u})\\\
&-\frac{\hbar^{2}}{2m}\left[n\left(\gradient\phi_{0}\right)^{2}-\varrho\left(\gradient\phi_{0}+\frac{m}{\hbar}\frac{D\mathbf{u}}{Dt}\right)^{2}\right]\end{split}$
(36)

where $n=\int_{\Omega}\rho_{0}\mathop{}\\!\mathrm{d}\mathbf{r}/|\Omega|$ is
the average particle density. We use this result, along with the elastic
expansion (24), to write out the full Lagrangian after the homogenisation
process as

$\frac{L_{\mathrm{hom}}}{\absolutevalue{\Omega}}=L_{0}-\hbar
n\dot{\phi}_{0}(1+\gradient\cdot\mathbf{u})-P_{ik}u_{ik}-\frac{1}{2}A_{iklm}u_{ik}u_{lm}-\frac{\hbar^{2}}{2m}\left(n-\varrho\right)\left(\gradient\phi_{0}\right)^{2}+\hbar\varrho\gradient\phi_{0}\cdot\frac{D\mathbf{u}}{Dt}+\frac{1}{2}m\varrho\left(\frac{D\mathbf{u}}{Dt}\right)^{2}$
(37)

This effective Lagrangian contains all the linear response of the system due
to both strain and phase gradients. We can see that the phase stiffness is now
given by $n-\varrho$ and $\varrho$ is the inertia associated with the motion
of the modulated pattern of the density. In summary, we have followed Ref.
[15] and reproduced the same results for the superfluid phase stiffness of the
supersolid. Our Lagrangian differs from that reference in the dynamical term.
In the next section, we complete the development of an effective low-energy
theory for the supersolid by considering the free energy of the system in the
grand canonical ensemble. Then, we can finally proceed to discuss the
macroscopic dynamics of the system.

##  VI Free energy and Grand Canonical ensemble

In Sections III to V, we developed a homogenised theory of the supersolid in
the canonical ensemble at a fixed average particle density of $n$. We have
integrated out the intra-unit-cell dynamics of the supersolid and have now
obtained an effective Lagrangian $L_{\mathrm{hom}}(\phi_{0},\mathbf{u})$ which
describes the system under a given phase gradient and displacements
$\mathbf{u}$. We have assumed that all the unit cells responded in the same
way, i.e. the response has the same periodicity as the supersolid. This
effective Lagrangian obeys all the necessary requirements of Galilean
invariance, rotational invariance and $U(1)$ invariance as specified by Son
[12] and Andreev & Lifshitz [7], so we would expect the dynamics to be
phenomenologically correct.

In this section, we will consider the inter-cell dynamics and allow for
particle flow between unit cells leading to variations in the density $\delta
n$ across the system at length scales much larger than the size of the unit
cells of the supersolid.

It is important to note that the additional energy response due to flow
between cells is crucial to obtaining the correct dynamics of the system. If
we were to neglect this energy response, we would find that the dynamics of
the system would not match up with a Bogoliubov treatment. A similar phenomena
was observed in [29], where the authors found that the bulk modulus computed
via a homogeneous dilation technique was not in agreement with long-wavelength
phonon predictions unless charge redistribution was taken into account.

We assume that the inter-cell number fluctuation is small, and so we can
expand the Lagrangian analytically in $\delta n$ up to
$\mathop{}\\!\mathcal{O}(\delta n^{2})$. We now consider the ground state
Lagrangian density $\mathcal{E}(n)=-L_{0}/|\Omega|$ which is a function of the
average density $n$ of the undeformed system. Under a change of the local
density, this can be expanded as

$\mathcal{E}(n+\delta n)\simeq\mathcal{E}(n)+\mathcal{E}^{\prime}\delta
n+\frac{1}{2}\mathcal{E}^{\prime\prime}(\delta n)^{2}$ (38)

and similarly the local pressure can be expanded as

$P(n+\delta n)\simeq P(n)+P^{\prime}\delta n\,.$ (39)

These are the only terms we need to track in a theory that is an expansion up
to second order in the total power of all perturbative fields. The
coefficients, $P^{\prime}$ and $\mathcal{E}^{\prime\prime}$, can be calculated
numerically by solving Gross-Pitaevskii equation at several densities. We can
also use $\mathcal{E}^{\prime}=\mu_{0}$ and
$\mathcal{E}^{\prime\prime}=\partial\mu/\partial n$ if we have the density
dependence of the chemical potential. In fact, we have already calculated the
pressure derivative in terms of the chemical potential shift per unit strain,
$\mu_{1}$, as defined in (16) with $\mu_{1}^{ik}=\mu_{1}\delta_{ik}$. To see
this, we note that the pressure, $P=-(\partial E/\partial V)_{N}$, and the
chemical potential, $\mu=(\partial E/\partial N)_{V}$, are related by a
Maxwell relation so that

$P^{\prime}=\frac{\partial P}{\partial n}=V\left(\frac{\partial P}{\partial
N}\right)_{V}\\!\\!\\!=-V\left(\frac{\partial\mu}{\partial
V}\right)_{N}\\!\\!\\!=-\mu_{1}\,.$ (40)

since $\delta V/V=u_{kk}$ is the strain. Similarly, we can use
$\mathcal{E}^{\prime\prime}=\partial\mu/\partial n$ if we have already
obtained the density dependence of the chemical potential.

In general, the above analysis is valid for any tensorial quantity, but we
have elected to use the symmetry properties of tensors that belong to a
triangular lattice (i.e. $P_{ik}=P\delta_{ik}$ and
$\mu_{1}^{ik}=\mu_{1}\delta_{ik}$) to simplify our calculations. The following
analysis also uses this symmetry, but can be generalised to any tensorial
quantity.

Finally, we need to include the work done by the change in local density on
the local environment. Namely, if a part of the system expands due to an
increase in particle number, it does work on the surroundings which must
compress. The work done caused by a total volume change of $\delta V$ is:

$\delta W=-P\delta V=-P(\epsilon_{ll}+M_{iklm}u_{ik}u_{lm})|\Omega|\,.$ (41)

In order to describe the dynamics of these density functions, we need to add
this work done to the Lagrangian to obtain

$\frac{L_{\mathrm{dyn}}}{\absolutevalue{\Omega}}=\frac{L_{\mathrm{hom}}}{\absolutevalue{\Omega}}-P\left(\epsilon_{ll}+M_{iklm}u_{ik}u_{lm}\right)\,.$
(42)

Comparing with our form (37) for $L_{\mathrm{hom}}$ and the relation (22)
between the Cauchy tensor $B_{iklm}$ and $A_{iklm}$, we find that

$\displaystyle-$ $\displaystyle
Pu_{ll}-\frac{A_{iklm}}{2}u_{ik}u_{lm}\\!-\\!P(\epsilon_{ll}+M_{iklm}u_{ik}u_{lm})$
$\displaystyle\simeq-\frac{1}{2}B_{iklm}u_{ik}u_{lm}\,.$ (43)

Using this, we arrive at an effective quadratic Lagrangian that can describe
the dynamics of the low-energy excitations of the supersolid. This is the
principal result of our paper:

$\begin{split}\frac{L_{\mathrm{dyn}}}{\absolutevalue{\Omega}}=&-\mathcal{E}-\mu_{0}\delta
n-\frac{1}{2}\mathcal{E}^{\prime\prime}(\delta n)^{2}-\hbar
n\dot{\phi}_{0}\left(1+\gradient\cdot\mathbf{u}\right)-\hbar\delta
n\dot{\phi}_{0}+\mu_{1}\delta n\gradient\cdot\mathbf{u}\\\
&-\frac{1}{2}B_{iklm}u_{ik}u_{lm}-\frac{\hbar^{2}}{2m}\left(n-\varrho\right)\left(\gradient\phi_{0}\right)^{2}+\hbar\varrho\gradient\phi_{0}\cdot\dot{\mathbf{u}}+\frac{1}{2}m\varrho\dot{\mathbf{u}}^{2}\,.\end{split}$
(44)

In this process we note that we have directly identified the Cauchy elastic
tensor in two different ways, firstly by using the result from Bavaud [28],
and secondly by directly considering the work the excitations would need to do
on their environment.

To summarise, we have a coarse-grained Lagrangian to describe the low-energy
excitations of the supersolid. In the next section, we will investigate these
dynamics.

##  VII Effective supersolid Lagrangian

In the previous sections, we have derived an effective Lagrangian (44) by
integrating out the short-wavelength modes in the unit cell and coarse-
graining the system such that a unit cell is now considered a point in the
continuous field theory described by the fields $\\{\mathbf{u},n,\phi\\}$.
This Lagrangian formally describes a low-energy effective field theory wherein
the coupled Goldstone modes of the system are the low-lying excitations. We
now set out to calculate the normal modes by integrating out the field $\delta
N$ and solving the resulting set of coupled equations.

By minimising the Lagrangian for $\delta n$, we find

$\delta
n=\frac{\mu_{1}}{\mathcal{E}^{\prime\prime}}\gradient\cdot\mathbf{u}-\frac{\hbar}{\mathcal{E}^{\prime\prime}}\dot{\phi}$
(45)

which we substitute back into the Lagrangian to straightforwardly obtain

$\begin{split}L=&-\mathcal{E}-\hbar
n\dot{\phi}+\frac{\hbar^{2}}{2}\left[\frac{\dot{\phi}^{2}}{\mathcal{E}^{\prime\prime}}-\frac{n-\varrho}{m}\left(\gradient\phi\right)^{2}\right]\\\
&+\frac{1}{2}\left[m\varrho\dot{\mathbf{u}}^{2}-\left(B_{iklm}-\frac{\mu_{1}}{\mathcal{E}^{\prime\prime}}\mu_{1}\delta_{ik}\delta_{lm}\right)u_{ik}u_{lm}\right]\\\
&-\hbar\left(n-\varrho+\frac{\mu_{1}}{\mathcal{E}^{\prime\prime}}\right)\dot{\phi}\gradient\cdot\mathbf{u}\end{split}$
(46)

This Lagrangian should be able to describe the Goldstone modes for the system.
This effective Lagrangian is of the form predicted by Ye [30].

By counting the number of spontaneously broken continuous symmetries, we are
able to predict the number of Goldstone modes that emerge. Since we have
broken the $U(1)$ gauge symmetry, we expect to find one Goldstone mode. We
also have a triangular lattice, which means we have broken two continuous
translational and rotational symmetries, and so we expect to find two
Goldstone modes. Therefore we should find three Goldstone modes in total.

One can immediately find the equations of motion for both $\phi$ and
$\mathbf{u}$, but it serves as a useful check to ensure we recover the
Lagrangian for a Bogoliubov mode in the superfluid phase. In a superfluid
phase, one can show that $B_{iklm}=U_{k=0}n^{2}\delta_{ik}\delta_{lm}$,
$\mu_{1}=-U_{k=0}n$, $\mathcal{E}^{\prime\prime}=U_{k=0}$ and $\varrho=0$.
Consequently, the resultant Lagrangian is now

$L=-\mathcal{E}-\hbar
n\dot{\phi}+\frac{\hbar^{2}}{2}\left[\frac{\dot{\phi}^{2}}{U_{k=0}}-\frac{n}{m}\left(\gradient\phi\right)^{2}\right]$
(47)

which corresponds exactly with the superfluid Bogoliubov mode.

For the general supersolid case, we derive the equations of motion directly by
assuming the fields $\phi$ and $\mathbf{u}$ are of a plane wave form $\sim
e^{i\left(\mathbf{k}\cdot\mathbf{r}-\omega t\right)}$. The vector equation for
the variable $\mathbf{u}$ can be solved to provide both the shear and
longitudinal mode. The Euler-Lagrange equation for the variable $\mathbf{u}$
is immediately found to be

$\begin{split}m\varrho\ddot{\mathbf{u}}-&\left(K+G-\frac{\mu_{1}}{\mathcal{E}^{\prime\prime}}\mu_{1}\right)\gradient\left(\gradient\cdot\mathbf{u}\right)-G\nabla^{2}\mathbf{u}\\\
&-\hbar\left(n-\varrho+\frac{\mu_{1}}{\mathcal{E}^{\prime\prime}}\right)\gradient\dot{\phi}=0\end{split}$
(48)

where we have used (23) as we are analysing a triangular lattice.

We can now solve this equation by noting that the vector field $\mathbf{u}$
can be split into longitudinal and transverse elements
$\mathbf{u}=\mathbf{u}_{l}+\mathbf{u}_{t}$. Taking the curl of (48) isolates
the transverse motion ($\mathbf{u}_{t}$) giving:

$\gradient\times\left(m\varrho\ddot{\mathbf{u}}_{t}-G\nabla^{2}\mathbf{u}_{t}\right)=0$
(49)

This is a wave equation for $\gradient\times\mathbf{u}$ from which we obtain
the shear velocity

$c_{\textrm{shear}}=\sqrt{\frac{G}{m\varrho}}$ (50)

which is exactly equivalent to the shear velocity predicted in standard solid
elasticity.

We examine longitudinal dynamics ($\mathbf{u}_{l}$) by taking the divergence
of (48) we are examining the dynamics of $\mathbf{u}_{l}$. We find the two
coupled equations

$\hbar\left(\frac{\omega^{2}}{\mathcal{E}^{\prime\prime}}-\frac{n-\varrho}{m}\mathbf{k}^{2}\right)\Phi-i\omega\left(n-\varrho+\frac{\mu_{1}}{\mathcal{E}^{\prime\prime}}\right)\gradient\cdot\mathbf{u}_{l}=0\,,\\\
$ (51)
$\begin{split}\left[m\varrho\omega^{2}-\left(K-\frac{\mu_{1}}{\mathcal{E}^{\prime\prime}}\mu_{1}+2G\right)\mathbf{k}^{2}\right]\gradient\cdot\mathbf{u}_{l}&\\\
+i\omega\hbar\mathbf{k}^{2}\left(n-\varrho+\frac{\mu_{1}}{\mathcal{E}^{\prime\prime}}\right)&\Phi=0\end{split}$
(52)

in the variables $\mathbf{u}_{l}$ and $\phi$. It can be shown easily that if
we are in a superfluid regime the second equation is trivially satisfied, and
the first equation finds the Bogoliubov frequency. These equations can be
solved for $\omega$ as a function of $\mathbf{k}$ and leads to the following
dispersion relation

$\omega^{2}=\frac{A\pm D}{2m\varrho}\,\mathbf{k}^{2}$ (53)

with

$\begin{split}A&=K+2G+\left(n\mathcal{E}^{\prime\prime}+2\mu_{1}\right)\left(n-\varrho\right)\\\
D^{2}&=A^{2}-4\varrho(n-\varrho)\left[\mathcal{E}^{\prime\prime}(K+2G)-\mu_{1}^{2}\right]\,.\end{split}$
(54)

We note that, in the case of the homogeneous superfluid ($\varrho=0$), we
obtain the superfluid velocity directly from (47)

$\omega^{2}=\frac{n\,U_{k=0}}{m}\mathbf{k}^{2}$ (55)

which is exactly equal to the Bogoliubov velocity. In the opposite limit of
$\varrho=n$, the two modes have frequencies

$\omega^{2}=\frac{K+2G}{mn}\mathbf{k}^{2}\quad\text{and zero.}$ (56)

The zero mode corresponds to the collapse of the Bogoliubov mode and the
linear dispersing mode has a velocity as expected for a normal solid with
elastic moduli $K$ and $G$.

It is interesting to note that, although we expect the Gross-Pitaevskii
Lagrangian to struggle with describing a normal solid, the effective
phenomenological Lagrangian that emerges actually returns to that of a normal
solid in the limit $\varrho\rightarrow n$.

We can also now examine the nature of the normal modes of the system. By
writing (51) and (52) in terms of a matrix equation, we can find the
eigenvectors of the system for a given $\omega$. These eigenvectors are the
normal modes, and are some composite of the fields $\phi$ and
$\mathbf{u}_{l}$. We can then subsequently define an eigenvector angle
$\theta$ which describes the relative contribution of the fields $\phi$ and
$\mathbf{u}_{l}$ to the normal mode. For example, in the superfluid the
eigenvector angle is $\theta=0$ as the normal mode is purely $\phi$, i.e.
$\mathbf{e}_{n}\equiv\mathbf{e}_{\phi}$. A supersolid with normal mode
$\mathbf{e}_{n}=\left(0.5,0.5\right)$ would have $\theta=\pi/4$ etc. The
normal mode composition is illustrated in Figure 5.

An interesting point of note is that the slow longitudinal mode (which is
expected to die off in the limit of a solid) is a composite of both $\phi$ and
$\mathbf{u}_{l}$, but is mostly $\phi$ at all points in the phase diagram
(i.e. $\theta$ is small). More importantly, it appears there is a critical
point at which the slow longitudinal mode is ’saturated’, that is to say that
there is a global maximum in the mixing of the slow longitudinal mode (but
note that the mode is still mostly $\phi$). This feature appears to be absent
in the fast longitudinal mode. However, the fast longitudinal mode also
exhibits interesting behaviour as the mixing angle does not tend to $\pi/2$ or
$3\pi/2$ as one may initially expect (as we would expect the mode to tend to a
purely elastic phonon mode in the limit of a solid), but instead tends to
approximately $7\pi/8$. This is a curious feature of the fast longitudinal
mode, and we are currently investigating the origin of this behaviour.

The above holds for any system described by a Gross-Pitaevskii equation; we
can now use these results in the specific case of a soft-core interaction. We
find a ground state for the system numerically and by applying the procedure
detailed above we can determine the velocities of all modes that exist.
Importantly we note that our Lagrangian is tuned across the transition, i.e. a
single theory describes both the superfluid and supersolid regime and
continuously varies across a first order phase transition. Importantly,
symmetry arguments by Son [12] dictate the functional form of any Lagrangian
belonging to a supersolid and since our Lagrangian is a subset of those
allowed by symmetry, we argue that we completely capture the phenomenological
behaviour of the supersolid as described by a Hartree-Fock approximation.

Figure 5: Top: excitation velocities for all available sound modes of a system
governed by a Heaviside interaction of strength $\Lambda$ with a phase
transition at $\Lambda\sim 38$, as calculated via homogenisation. Bottom: the
eigenvector angle $\theta$ for the available longitudinal modes.

##  VIII Bogoliubov Dispersion

An alternative technique to finding excitations of the system is to directly
apply linearised oscillations in the number density and the phase. Namely,
since $\psi=\sqrt{\rho}e^{i\phi}$, we can apply the Bogoliubov fluctuations
$\rho\rightarrow\rho_{0}+\delta\rho$ and $\phi$. We expect to recover all the
Goldstone modes accessible to the system, and expect the quantitative
properties to be the same as those predicted in the elastic theory. More
precisely, due to the broken $U(1)$ phase and the broken continuous
translational symmetry, we expect to find three Goldstone modes, two of which
are longitudinal and one of which is a shear mode. Since this is a completely
independent technique to recovering the Goldstone modes, it serves as an
important check of our results.

Whilst both techniques access the speeds of sound, the elastic technique also
provides elastic moduli and fractional inertia, whereas the bogoliubov
technique provides a full band structure. In this way, the two techniques
complement each other, as they verify common ground but also provide
information that is inaccesible to the other.

By following a similar procedure to that outlined in [14], we are able to
obtain eigenvalue equations which give the excitation frequencies $\omega$ of
the Bloch waves in the system. These Bloch waves have a band structure, and
the lowest band structure at the $\Gamma$ point corresponds to the Goldstone
modes of the system whose propagation velocities were independently derived in
[14]. The band structure is fully described by the equations

$\begin{split}&\left[-\omega^{2}+\left(\hat{T}^{\mathbf{k}}+\hat{U}^{\mathbf{k}}\right)\hat{T}^{\mathbf{k}}\right](\Phi)^{\mathbf{k}}=0\\\
&\left[-\omega^{2}+\hat{T}^{\mathbf{k}}\left(\hat{T}^{\mathbf{k}}+\hat{U}^{\mathbf{k}}\right)\right](\delta\psi)^{\mathbf{k}}=0\end{split}$
(57)

with $\hat{T}$ and $\hat{U}$ defined as

$\begin{split}\hat{T}^{\mathbf{k}}_{\mathbf{G},\mathbf{G}^{\prime}}=&\left[\frac{1}{2}\left(\mathbf{G}+\mathbf{k}\right)^{2}-\mu\right]\delta_{\mathbf{G},\mathbf{G}^{\prime}}+U_{\mathbf{G}-\mathbf{G}^{\prime}}\rho_{\mathbf{G}-\mathbf{G}^{\prime}}\\\
\hat{U}^{\mathbf{k}}_{\mathbf{G},\mathbf{G}^{\prime}}=&2\sum_{\mathbf{G}^{\prime\prime}}\psi_{\mathbf{G}-\mathbf{G}^{\prime\prime}}U_{\mathbf{k}+\mathbf{G}^{\prime\prime}}\psi_{\mathbf{G}^{\prime\prime}-\mathbf{G}^{\prime}}\end{split}$
(58)

with $U_{\mathbf{q}}$, $\psi_{\mathbf{q}}$ and $\rho_{\mathbf{q}}$ the fourier
transforms of the interaction potential, the superfluid order parameter and
the density respectively. We expect that the band structure obtained from
solving these equations should find all the requisite Goldstone modes and that
they should have the same velocities as predicted via homogenisation.

Figure 6: Band structure for a supersolid at $\Lambda=50$. The three Goldstone
modes are highlighted as a fast longitudinal, slow longitudinal and shear
mode.

Upon recovery of the band structure, one finds exactly $(D+1)$ Goldstone modes
as expected. There is a band crossing at the $K$ point of the fast
longitudinal and shear modes, characteristic of a triangular lattice (n.b.
this is not actually a crossing but a band touching, as shown by the
$M\rightarrow\Gamma$ band structure).

For any particular given $\Lambda$, both the Bogoliubov and elastic techniques
predict all three speeds of sound (if supersolid, only one speed if
superfluid), and agree on the speeds with excellent precision. Subsequently,
by deriving the excitations of the system in two independent ways and
obtaining the same results to excellent numerical precision, we can reliably
say we have obtained the correct dynamics.

##  IX Results and Discussion

We can see that we have now verified the technique of homogenisation as a
means of recovering the low-lying excitations of the system. This was achieved
via means of two completely independent techniques, the results of which both
agree to excellent accuracy. An advantage of homogenisation over the
Bogoliubov technique is that it also provides the elastic constants of the
system, which cannot be obtained via Bogoliubov. Moreover, it provides the
fractional inertia of the system, which is a quantity that is of great
physical interest. A clear disadvantage however, is that homogenisation does
not provide the full spectrum of excitations, only the long-wavelength limit.
We stress that this is not an issue, as the Bogoliubov technique is completely
separate from homogenisation, and so by utilising both techniques one can
recover the full band structure from Bogoliubov and supplement it with
additional information in the long-wavelength limit from homogenisation.

The effective Lagrangian that is the end result of homogenisation elucidates
interesting dynamics. By considering that the flow excitations of the system
can be coupled to the strain excitations, we can see that the system is
capable of strain-induced flow. Furthermore, by considering the application of
an instantaneous strain to the system, we can see that flow will be induced
and contribute to the elastic response of the system. Alternatively, if we
apply a strain and keep it fixed for a long time so that the flow dissipates,
there will be a different elastic response. In this way we can see that the
elastic modulus has a time dependence, and that the system is capable of
dissipating energy via the flow of the system. This energy dissipation is
likely to be via the two longitudinal modes which are coupled crystal and
Bogoliubov phonons.

There is an interesting question of what the compressibility of the system is.
When one considers the elastic strain, it is natural to define the
compressibility of the system as $\kappa=-1/V\,(\partial V/\partial P)$.
Typically in a solid, one would simultaneously be able to say
$\kappa=1/n\,(\partial n/\partial P)$, where $n$ is the number density and
$\mu$ the chemical potential. However, in our analysis, we can show that these
two definitions do not agree, and are off by a considerable margin. These two
values are defined as $\kappa=1/B_{xxyy}$ and $\kappa=1/(n\mu_{1})$ in our
notation respectively.

The resolution of this seeming dilemma is to revisit the definition of the two
compressibilities. The compressibility $\kappa=1/B$ where $B=-V\,(\partial
P/\partial V\bigr{|}_{N})$ is the elastic bulk modulus of the system, taken at
constant particle number. This is the compressibility of the system as defined
via elastic theory, and is the one that we will take as the ’true
compressibility’, and consequently $B$ as the ’true bulk modulus’.
Importantly, the derivative is taken under constant particle number. As such,
to elucidate this a little further, we write

$B_{e}=-V\frac{\partial P}{\partial V}\biggr{|}_{N}$ (59)

where the subscript $e$ denotes that this is the elastic bulk modulus. We can
then write the compressibility as $\kappa_{e}=1/B_{e}$. We can subsequently
define the thermodynamic bulk modulus (using the typical definition with $n$
but now removing the factor of $V$ which is kept constant) as

$B_{t}=N\frac{\partial P}{\partial N}\biggr{|}_{V}$ (60)

and the thermodynamic compressibility accordingly.

The key consideration is that the $B_{e}\neq B_{t}$ due to the fact that we
are keeping different variables constant upon taking derivatives. If we want
to connect these two definitions, we need to use the cyclic relations of
thermodynamics

$\frac{\partial x}{\partial y}\biggr{|}_{z}=-\frac{\partial x}{\partial
z}\biggr{|}_{y}\frac{\partial z}{\partial y}\biggr{|}_{x}\quad.$ (61)

We can then write

$B_{e}=-V\frac{\partial P}{\partial V}\biggr{|}_{N}=V\frac{\partial
P}{\partial N}\biggr{|}_{V}\frac{\partial N}{\partial V}\biggr{|}_{P}\equiv
B_{t}\ \frac{1}{n}\frac{\partial N}{\partial V}\biggr{|}_{P}$ (62)

from which we can consider the following. In a regular material (solid, gas,
liquid) we would expect that the derivative of the particle number with
respect to the volume under constant pressure is simply the density of the
system, and we would recover the thermodynamic bulk modulus
$B_{t}=n\,(\partial P/\partial n\bigr{|}_{V})$ from (62). We can then examine
whether this is the case for our system, i.e. is $1/n\,(\partial N/\partial
V\bigr{|}_{P})=1$?

We can examine this quantity, henceforth denoted as $\mathcal{R}$ and called
the ’compressibility ratio’, by choosing an initial volume and particle number
$V_{0}$ and $N_{0}$ which will generate a $P_{0}$, adjusting the volume around
$V_{0}$ with some $\delta V$, and then varying the particle number
$N_{0}+\delta N$ until the system reaches $P_{0}$. We can then numerically
calculate $\mathcal{R}$ by finite differencing. If we expect the difference
between the two compressibilities to be due to the fact that
$\frac{1}{n}\frac{\partial N}{\partial V}\bigr{|}_{P}\neq 1$ then we would
expect $\mathcal{R}=\frac{\kappa_{t}}{\kappa_{e}}$, or equivalently

$\mathcal{R}=\frac{B_{e}}{B_{t}}\,.$ (63)

Upon performing the numerical determination of $\mathcal{R}$, we find that it
is within good agreement of the ratio of the two compressibilities as
determined via homogenisation. The percentage error between the homogenisation
and numerical $\mathcal{R}$ is typically between $3\%$ and $5\%$. We believe
this is an artifact of the numerical determination of $\mathcal{R}$, and that
as the algorithm and precision of simulations for determining $\mathcal{R}$ is
improved, the agreement will improve accordingly.

Figure 7: Compressibility ratio $\mathcal{R}$ for a system governed by a soft-
core interaction of strength $\Lambda$. The expected value of $\mathcal{R}=1$
is obtained in the superfluid phase, but is starkly differenti n the
supersolid regime.

While we have shown that the compressibility ratio obtained via homogenisation
is in good agreement with the numerical determination, we do not yet fully
understand the nature of this anomalous result. In our search across the
literature we have not found another example of this phenomenon, and we are
currently investigating the origin of this effect. A possible explanation for
the result is considering a pressure which can be written like $P(n,G(n,V))$,
where $G$ is the ordering wavevector. One can then subsequently show that the
difference between bulk moduli can be written as

$B_{e}=B_{t}-V\frac{\partial P}{\partial G}\biggr{|}_{n}\frac{\partial
G}{\partial V}\biggr{|}_{n}$ (64)

where it is clear that the last term is the anomalous term.

This term vanishes in ordinary materials in various ways. For gases and
liquids, there is no ordering vector, so this term never existed to begin
with. A regular solid has well-defined atomic positions. So, if the density is
specified, the ordering wavevector is also specified. Upon changing the
volume, there is no additional change in the ordering wavevector beyond that
caused by the change in density. Supersolids are unique in that the ordering
vector $G$ experiences a commensuration effect. For instance, in a one-
dimensional system, it needs to be an integer multiple of $2\pi/L$ with $L$
being the length of the system. If we stretch the system by $\delta L$ much
less than the lattice spacing, then the fractional change in the ordering
wavevector is $\delta G/G=-\delta L/L$ so that $L\partial G/\partial L=-G$. We
believe this is a feature unique to supersolids, due to their ability to exist
at a wide range of densities, each of which has their own ordering vector,
whereas regular solids only have one density at which they can exist (at some
given pressure) and therefore only have one ordering vector. This is (to the
best of our knowledge) a new phenomena which could be of great interest to the
characterisation of supersolids.

Some currently open questions are as follows. Is this effect unique to the
supersolid phase, or is it a more general phenomena? What kind of thermal
transport is occurring in the phonon modes of a supersolid? Does the notion of
the Landau criterion for superfluidity still hold in the supersolid phase, and
if so, how does it change now that we have a band structure? These questions
may potentially have avenues to solution via the techniques developed in this
paper.

The successful development of homogenisation techniques for the study of
supersolids is a significant step forward in the study of these systems. We
have shown that the homogenisation technique is able to accurately capture the
behaviour of the supersolid phase through comparison with well established
Bogoliubov techniques. We were able to obtain previously unknown elastic
properties of these materials as well as an effective Lagrangian which can be
used to study the dynamics of the system. In contrast with typical phases of
matter which have $\mathcal{R}=1$, we found that the compressibility ratio of
the supersolid phase is $\mathcal{R}>1$. This is a new phenomena which we
predict will arise in supersolids, and of which the origin is currently
unknown.

## References

  * Kim and Chan [2014] D. Y. Kim and M. H. Chan, Upper limit of supersolidity in solid helium, Physical Review B 90, 064503 (2014).
  * Day and Beamish [2007] J. Day and J. Beamish, Low-temperature shear modulus changes in solid 4 he and connection to supersolidity, Nature 450, 853 (2007).
  * Nyéki  _et al._ [2017] J. Nyéki, A. Phillis, A. Ho, D. Lee, P. Coleman, J. Parpia, B. Cowan, and J. Saunders, Intertwined superfluid and density wave order in two-dimensional 4 he, Nature Physics 13, 455 (2017).
  * Chomaz  _et al._ [2019] L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. Van Bijnen, A. Patscheider, M. Sohmen, _et al._ , Long-lived and transient supersolid behaviors in dipolar quantum gases, Physical Review X 9, 021012 (2019).
  * Guo  _et al._ [2021] Y. Guo, R. M. Kroeze, B. P. Marsh, S. Gopalakrishnan, J. Keeling, and B. L. Lev, An optical lattice with sound, Nature 599, 211 (2021).
  * Norcia  _et al._ [2021] M. A. Norcia, C. Politi, L. Klaus, E. Poli, M. Sohmen, M. J. Mark, R. N. Bisset, L. Santos, and F. Ferlaino, Two-dimensional supersolidity in a dipolar quantum gas, Nature 596, 357 (2021).
  * Andreev and Lifshits [1969] A. Andreev and I. Lifshits, Quantum theory of defects in crystals, Zhur Eksper Teoret Fiziki 56, 2057 (1969).
  * Chester [1970] G. Chester, Speculations on bose-einstein condensation and quantum crystals, Physical Review A 2, 256 (1970).
  * Leggett [1970] A. J. Leggett, Can a solid be” superfluid”?, Physical Review Letters 25, 1543 (1970).
  * Nozieres [2009] P. Nozieres, Superfluidity, crystallization and localization, Journal of Low Temperature Physics 156, 9 (2009).
  * Dorsey  _et al._ [2006] A. T. Dorsey, P. M. Goldbart, and J. Toner, Squeezing superfluid from a stone: Coupling superfluidity and elasticity in a supersolid, Physical Review Letters 96, 055301 (2006).
  * Son [2005] D. Son, Effective lagrangian and topological interactions in supersolids, Physical Review Letters 94, 175301 (2005).
  * Heinonen  _et al._ [2019] V. Heinonen, K. J. Burns, and J. Dunkel, Quantum hydrodynamics for supersolid crystals and quasicrystals, Physical Review A 99, 063621 (2019).
  * Kunimi and Kato [2012] M. Kunimi and Y. Kato, Mean-field and stability analyses of two-dimensional flowing soft-core bosons modeling a supersolid, Phys. Rev. B 86, 060510 (2012).
  * Josserand  _et al._ [2007a] C. Josserand, Y. Pomeau, and S. Rica, Patterns and supersolids, The European Physical Journal Special Topics 146, 47 (2007a).
  * Josserand  _et al._ [2007b] C. Josserand, Y. Pomeau, and S. Rica, Coexistence of ordinary elasticity and superfluidity in a model of a defect-free supersolid, Phys. Rev. Lett. 98, 195301 (2007b).
  * During  _et al._ [2011] G. During, C. Josserand, Y. Pomeau, and S. Rica, Theory of real supersolids, arXiv preprint arXiv:1110.1323 (2011).
  * Leggett [1998] A. Leggett, On the superfluid fraction of an arbitrary many-body system at t= 0, Journal of Statistical physics 93, 927 (1998).
  * Sepúlveda  _et al._ [2008] N. Sepúlveda, C. Josserand, and S. Rica, Nonclassical rotational inertia fraction in a one-dimensional model of a supersolid, Physical Review B 77, 054513 (2008).
  * Roccuzzo  _et al._ [2020] S. Roccuzzo, A. Gallemí, A. Recati, and S. Stringari, Rotating a supersolid dipolar gas, Physical Review Letters 124, 045702 (2020).
  * Aftalion  _et al._ [2007] A. Aftalion, X. Blanc, and R. L. Jerrard, Nonclassical rotational inertia of a supersolid, Physical Review Letters 99, 135301 (2007).
  * Tanzi  _et al._ [2021] L. Tanzi, J. Maloberti, G. Biagioni, A. Fioretti, C. Gabbanini, and G. Modugno, Evidence of superfluidity in a dipolar supersolid from nonclassical rotational inertia, Science 371, 1162 (2021).
  * Blakie [2023] P. B. Blakie, Superfluid fraction tensor of a two-dimensional supersolid, arXiv preprint arXiv:2308.14001 (2023).
  * Wessel and Troyer [2005] S. Wessel and M. Troyer, Supersolid hard-core bosons on the triangular lattice, Physical Review Letters 95, 127205 (2005).
  * Boninsegni and Prokof’ev [2005] M. Boninsegni and N. Prokof’ev, Supersolid phase of hard-core bosons on a triangular lattice, Physical Review Letters 95, 237204 (2005).
  * Boninsegni [2012] M. Boninsegni, Supersolid phases of cold atom assemblies, Journal of Low Temperature Physics 168, 137 (2012).
  * Landau  _et al._ [1986] L. D. Landau, E. M. Lifshitz, A. M. Kosevich, and L. P. Pitaevskii,  _Theory of elasticity: volume 7_ , Vol. 7 (Elsevier, 1986).
  * Bavaud  _et al._ [1986] F. Bavaud, P. Choquard, and J.-R. Fontaine, Statistical mechanics of elastic moduli, Journal of Statistical Physics 42, 621 (1986).
  * Sutton  _et al._ [1988] A. P. Sutton, M. W. Finnis, D. G. Pettifor, and Y. Ohta, The tight-binding bond model, Journal of Physics C: Solid State Physics 21, 35 (1988).
  * Ye [2008] J. Ye, Elementary excitation in a supersolid, Europhysics Letters 82, 16001 (2008).

##  Appendix A Expansion calculation

We start by expanding the non-interacting part of $\mathcal{L}$ using the
tools developed. We can express without expansion

$\left(\nabla^{\prime}\rho(\mathbf{r}^{\prime})\right)^{2}=\left[\left(\partial_{i}+\left(u_{ik}+u_{il}u_{lk}\right)\partial_{k}\right)\left(\rho_{0}(\mathbf{r})+\tilde{\rho}(\mathbf{r})\right)\right]^{2}$
(65)

where we subsequently drop the coordinates as it is clear everything is in the
unprimed basis $\mathbf{r}$. We can re-express the pre-factor

$\frac{1}{\rho(\mathbf{r}^{\prime})}=\frac{1}{\rho_{0}+\tilde{\rho}}=\frac{1}{\rho_{0}}-\frac{\tilde{\rho}}{\rho_{0}^{2}}+\frac{\tilde{\rho}^{2}}{\rho_{0}^{3}}+\mathop{}\\!\mathcal{O}(\tilde{\rho}^{3})$
(66)

to make it more algebraically pliable. This allows us to write

$\displaystyle\int_{\Omega^{\prime}}\frac{\hbar^{2}}{2m}$
$\displaystyle\frac{(\nabla^{\prime}\rho(\mathbf{r}^{\prime}))^{2}}{4\rho(\mathbf{r}^{\prime})}\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime}=\int_{\Omega}\Big{\\{}\frac{\hbar^{2}}{8m}\left(\frac{1}{\rho_{0}}-\frac{\tilde{\rho}}{\rho_{0}^{2}}+\frac{\tilde{\rho}^{2}}{\rho_{0}^{3}}\right)\left(1+\epsilon_{ll}+M_{iklm}u_{ik}u_{lm}\right)\times$
(67)
$\displaystyle\left[(\nabla\rho_{0})^{2}+2(\gradient\rho_{0}\cdot\gradient\tilde{\rho})-2\epsilon_{ik}\partial_{i}\rho_{0}\partial_{k}\rho_{0}+2\chi_{ik}\partial_{i}\rho_{0}\partial_{k}\rho_{0}-4\epsilon_{ik}\partial_{i}\rho_{0}\partial_{k}\tilde{\rho}+(\nabla\tilde{\rho})^{2}+2\Delta_{ik}\partial_{i}\rho_{0}\partial_{k}\rho_{0}\right]\Big{\\}}\mathop{}\\!\mathrm{d}\mathbf{r}$

which we can manipulate further. By isolating first and second order (in terms
of the strain tensor) parts of the non-interacting Lagrangian, and using the
following relation (obtained via integration by parts)

$\int_{\Omega}\frac{\gradient\tilde{\rho}\cdot\gradient\rho_{0}}{\rho_{0}}\mathop{}\\!\mathrm{d}\mathbf{r}=\int_{\Omega}\tilde{\rho}\frac{\left(\nabla\rho_{0}\right)^{2}}{\rho_{0}^{2}}-\tilde{\rho}\frac{\nabla^{2}\rho_{0}}{\rho_{0}}\mathop{}\\!\mathrm{d}\mathbf{r}$
(68)

we are able to find the first-order component

$\frac{\hbar^{2}}{4m}\left[\left(\frac{(\nabla\rho_{0})^{2}}{2\rho_{0}^{2}}-\frac{\nabla^{2}\rho_{0}}{\rho_{0}}\right)\tilde{\rho}+\epsilon_{ik}\left(\frac{(\nabla\rho_{0})^{2}}{2\rho_{0}}\delta_{ik}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}}\right)\right]\\\
$ (69)

and the second-order component

$\displaystyle\frac{\hbar^{2}}{4m}\left[\left(\Delta_{ik}+\chi_{ik}\right)\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}}\right]+$
(70)
$\displaystyle\frac{\hbar^{2}}{4m}\frac{1}{2}\left[-4\epsilon_{ik}\frac{\partial_{i}\rho_{0}\partial_{k}\tilde{\rho}}{\rho_{0}}\right.+\frac{(\nabla\tilde{\rho})^{2}}{\rho_{0}}+2\frac{\gradient\tilde{\rho}\cdot\gradient\rho_{0}}{\rho_{0}}\epsilon_{ll}-2\epsilon_{ik}\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}}\epsilon_{ll}$
$\displaystyle\qquad\quad+\frac{(\nabla\rho_{0})^{2}}{\rho_{0}}\,M_{iklm}u_{ik}u_{lm}-2\frac{\gradient\tilde{\rho}\gradient\rho_{0}}{\rho_{0}^{2}}\tilde{\rho}\left.+2\epsilon_{ik}\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}^{2}}\tilde{\rho}-\frac{(\nabla\rho_{0})^{2}}{\rho_{0}^{2}}\tilde{\rho}\epsilon_{ll}+\frac{(\nabla\rho_{0})^{2}}{\rho_{0}^{3}}\tilde{\rho}^{2}\right]$

respectively. We now turn to considering the interaction section of the
Lagrangian.

By first noting that the interaction potential is centrosymmetric (and indeed
we expect most effective potentials to be isotropic), we can write the
interaction as a function of radial distance. We can express the change in the
metric length of distance by simply applying the deformation definition to
$\absolutevalue{\Delta\mathbf{r}^{\prime}}^{2}$ and obtain the result

$\left(\Delta
r^{\prime}\right)^{2}=\left(\delta_{ik}+2\epsilon_{ik}+2\Delta_{ik}\right)\left(\Delta
r\right)_{i}\left(\Delta r\right)_{k}$ (71)

which is exactly equivalent to that of Landau & Lifshitz [27]. The tensor in
the above equation is sometimes referred to as the finite strain tensor.

Recalling the derivation of (19) we define the functions

$U(\Delta\mathbf{r}^{\prime})=U(\Delta\mathbf{r})+\left(\epsilon_{ik}+\Delta_{ik}\right)f_{ik}(\Delta\mathbf{r})+\epsilon_{ik}\epsilon_{lm}W_{iklm}(\Delta\mathbf{r})$
(72)

with

$f_{ik}(\mathbf{r})\equiv\frac{r_{i}r_{k}}{|\mathbf{r}|}\frac{\partial
U(|\mathbf{r}|)}{\partial|\mathbf{r}|}\quad;\quad
W_{iklm}(\mathbf{r})\equiv\frac{r_{i}r_{k}r_{l}r_{m}}{2|\mathbf{r}|^{2}}\left(\frac{\partial^{2}U}{\partial|\mathbf{r^{2}}|}-\frac{1}{|\mathbf{r}|}\frac{\partial
U}{\partial|\mathbf{r}|}\right)\;.$ (73)

This can be used to express the interaction component of $\mathcal{L}$ by
combining it with both the density mapping and the Jacobian to obtain

$\displaystyle
U=\frac{1}{2}\int_{\Omega}\left(\rho_{0}+\tilde{\rho}\right)_{r_{1}}\left[U(\mathbf{r}_{12})+\left(\epsilon_{ik}+\Delta_{ik}\right)f_{ik}(\mathbf{r}_{12})\right.$
$\displaystyle\left.+\epsilon_{ik}\epsilon_{lm}W_{iklm}(\mathbf{r}_{12})\right]\left(\rho_{0}+\tilde{\rho}\right)_{r_{2}}$
(74)
$\displaystyle\left(1+2\epsilon_{ll}+\epsilon_{ll}\epsilon_{kk}+2M_{iklm}u_{ik}u_{lm}\right)\;\mathop{}\\!\mathrm{d}\mathbf{r}_{1}\mathop{}\\!\mathrm{d}\mathbf{r}_{2}\quad$

where $\mathbf{r}_{12}=\mathbf{r}_{1}-\mathbf{r}_{2}$. We note that, since the
convolution is integrated, there is a term arising from
$\left(\mathcal{J}_{\mathbf{r}^{\prime}\rightarrow\mathbf{r}}\right)^{2}$.

Once again, this can be expanded and combined with the non-interacting
component of $\mathcal{L}$, and we can fully express the first and second
order components of the expansion. We make use of the Gross-Pitaevskii ground-
state equation (2), as well as the normalisation condition (14) to simplify
our expressions, and we can show that

$\mathcal{L}_{1}=\int_{\Omega}\left[\rho_{0}\epsilon_{ll}\left(U*\rho_{0}-\mu_{0}\right)+\frac{1}{2}\epsilon_{ik}(f_{ik}*\rho_{0})\rho_{0}+\frac{\hbar^{2}}{4m}\epsilon_{ik}\left(\frac{\left(\nabla\rho_{0}\right)^{2}}{2\rho_{0}}\delta_{ik}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}}\right)\right]\mathop{}\\!\mathrm{d}\mathbf{r}$
(75)

We can now make a direct connection between the linear component of the
expansion and the pressure tensor $P_{ik}$. Recalling that the first order
component of the expansion is the stress tensor, and that the stress tensor is
the negative of the pressure, we can write

$\begin{split}P_{ik}&=\frac{1}{|\Omega|}\int_{\Omega}\left[\left(\mu_{0}-U*\rho_{0}\right)\delta_{ik}-\frac{1}{2}f_{ik}*\rho_{0}-\frac{\hbar^{2}}{4m}\left(\frac{\left(\nabla\rho_{0}\right)^{2}}{2\rho_{0}^{2}}\delta_{ik}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}^{2}}\right)\right]\rho_{0}\mathop{}\\!\mathrm{d}\mathbf{r}\\\
&=\frac{1}{|\Omega|}\int_{\Omega}\left[\frac{\hbar^{2}}{4m}\left(\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}}-\nabla^{2}\rho_{0}\,\delta_{ik}\right)-\frac{1}{2}\rho_{0}(f_{ik}*\rho_{0})\right]\mathop{}\\!\mathrm{d}\mathbf{r}\\\
\end{split}$ (76)

The second line consists of the ground-state density profile only. This is our
first key result which is obtaining the pressure of any ground-state system
given a specified interaction. We can trivially calculate the expected
pressure for a homogeneous superfluid by noting that $\rho_{0}\equiv n$ and
that $\mu_{0}=U_{\mathbf{k}=0}n$ where
$U_{\mathbf{k}=0}=\int_{\Omega}U(\mathbf{r})d\mathbf{r}$. The pressure simply
reduces to
$-\frac{1}{2}n^{2}\int_{\Omega}f_{ik}(\mathbf{r})\mathop{}\\!\mathrm{d}\mathbf{r}$.
Integrating by parts gives the integral over $f_{ik}$ as
$-U_{\mathbf{k}=0}\delta_{ik}$ . So, we find that
$P_{ik}=\frac{1}{2}\mu_{0}n\delta_{ik}$, which is the expected result. We can
now move on to the more difficult task of solving for the second order
components.

It is important to note that since we have allowed $\tilde{\rho}$ to have both
first and second order components in $u_{ik}$, then even single powers of
$\tilde{\rho}$ must be taken into consideration for the
$\mathop{}\\!\mathcal{O}(u^{2})$ components. The above analysis only collected
the components of $\tilde{\rho}$ which were of $\mathop{}\\!\mathcal{O}(u)$,
but does not remove the singular powers of $\tilde{\rho}$ from the Lagrangian.
It is necessary to calculate the functional form of $\tilde{\rho}$ as there
will be terms like
$\int_{\Omega}\tilde{\rho}^{2}\mathop{}\\!\mathrm{d}\mathbf{r}$ that will
contribute to the perturbative Lagrangian, and as such we therefore need to
solve the equations of motion for $\tilde{\rho}$. It is not good enough to
simply know the normalisation condition as we did for the
$\mathop{}\\!\mathcal{O}(u)$ component.

The first step in solving for $\tilde{\rho}$ is to collect all the terms in
which it appears, regardless of first or second order splitting. We take all
the terms that contain $\tilde{\rho}$ and call them
$\mathcal{L}_{\tilde{\rho}}$, finding

$\displaystyle-\mathcal{L}_{\tilde{\rho}}=$
$\displaystyle\frac{1}{2}\frac{\hbar^{2}}{4m}\left[-4\epsilon_{ik}\frac{\partial_{i}\rho_{0}\partial_{k}\tilde{\rho}}{\rho_{0}}+\frac{(\nabla\tilde{\rho})^{2}}{\rho_{0}}+2\frac{\gradient\tilde{\rho}\cdot\gradient\rho_{0}}{\rho_{0}}\epsilon_{ll}-2\frac{\gradient\tilde{\rho}\cdot\gradient\rho_{0}}{\rho_{0}^{2}}\tilde{\rho}+2\epsilon_{ik}\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}^{2}}\tilde{\rho}-\frac{(\nabla\rho_{0})^{2}}{\rho_{0}^{2}}\tilde{\rho}\,\epsilon_{ll}+\frac{(\nabla\rho_{0})^{2}}{\rho_{0}^{3}}\tilde{\rho}^{2}\right]$
(77)
$\displaystyle+\mu_{0}\tilde{\rho}+\frac{1}{2}\left(\left(U*\tilde{\rho}\right)\tilde{\rho}+4\epsilon_{ll}\left(U*\rho_{0}\right)\tilde{\rho}+2\epsilon_{ik}\left(f_{ik}*\rho_{0}\right)\tilde{\rho}\right)$

which we can solve by the standard Euler-Lagrange technique. A key detail now
is that the new normalisation condition changes the constraints we apply.
Specifically, the equation we need to solve is

$\frac{\partial\mathcal{L}_{\tilde{\rho}}}{\partial\tilde{\rho}}-\partial_{i}\frac{\partial\mathcal{L}_{\tilde{\rho}}}{\partial\left(\partial_{i}\tilde{\rho}\right)}+\left(\mu_{0}+\mu_{1}\epsilon_{ll}+\mu_{2}^{iklm}u_{ik}u_{lm}\right)\left(1+\epsilon_{ll}+M_{iklm}u_{ik}u_{lm}\right)=0$
(78)

which can be shown via the strained normalisation condition (14).

A straightforward calculation follows of differentiating the $\tilde{\rho}$
dependent Lagrangian and collating all terms. We make use of the zeroth-order
equation of motion to simplify our expression and we find

$\begin{split}&\frac{\hbar^{2}}{4m}\gradient\cdot\left(\frac{\gradient\tilde{\rho}}{\rho_{0}}\right)-U*\tilde{\rho}\\\\[5.0pt]
+&\frac{\hbar^{2}}{4m}\left(\frac{(\nabla\rho_{0})^{2}}{\rho_{0}^{3}}-\frac{\nabla^{2}\rho_{0}}{\rho_{0}^{2}}\right)\tilde{\rho}\end{split}\text{\Large{\hbox{
=
}}}\begin{split}&-\left(\mu_{2}^{iklm}+\mu_{1}^{ik}\delta_{lm}+\mu_{0}M_{iklm}\right)u_{ik}u_{lm}+\epsilon_{ik}\left[-\mu_{1}^{ik}+\frac{\hbar^{2}}{4m}\left(2\frac{\partial_{ik}\rho_{0}}{\rho_{0}}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}^{2}}\right)\right.\\\\[5.0pt]
&\left.+\int_{\Omega}\left(f_{ik}(\mathbf{r}-\mathbf{r}^{\prime})+\delta_{ik}U(\mathbf{r}-\mathbf{r}^{\prime})\right)\rho_{0}(\mathbf{r}^{\prime})\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime}\right]\end{split}$
(79)

which contains both $\mathop{}\\!\mathcal{O}(u^{2})$ and
$\mathop{}\\!\mathcal{O}(u)$ terms. We can now begin to separate
$\tilde{\rho}$ into first and second order components via (13), and consider
solving the equation component by component. That is to say that we solve the
lowest order component first, and then use that solution to solve the next
order component, and so on.

The first order equation to solve is given by

$\begin{split}&\frac{\hbar^{2}}{4m}\gradient\cdot\left(\frac{\gradient\left(\rho_{1}^{ik}\epsilon_{ik}\right)}{\rho_{0}}\right)-U*\left(\rho_{1}^{ik}\epsilon_{ik}\right)\\\
&+\frac{\hbar^{2}}{4m}\left(\frac{(\nabla\rho_{0})^{2}}{\rho_{0}^{3}}-\frac{\nabla^{2}\rho_{0}}{\rho_{0}^{2}}\right)\left(\rho_{1}^{ik}\epsilon_{ik}\right)\end{split}\text{\Large{\hbox{
=
}}}\begin{split}\epsilon_{ik}&\left[-\mu_{1}^{ik}+\frac{\hbar^{2}}{4m}\left(2\frac{\partial_{ik}\rho_{0}}{\rho_{0}}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}^{2}}\right)\right.\\\\[5.0pt]
&\left.+\int_{\Omega}\left(f_{ik}(\mathbf{r}-\mathbf{r}^{\prime})+\delta_{ik}U(\mathbf{r}-\mathbf{r}^{\prime})\right)\rho_{0}(\mathbf{r}^{\prime})\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime}\right]\end{split}$
(80)

which we can solve numerically. We factor out $\epsilon_{ik}$ and solve the
set of equations where we have one for each index (in 2D we have 4 equations
to solve, in 3D there are 9). As this is an integrodifferential equation we
can solve the equation by writing the LHS as a matrix operator on an unraveled
vector and applying an iterative matrix solver algorithm which converges on
some pre-specified tolerance. We can introduce notation and write the first
and second order equations of motion as

$L\left(\rho_{1}^{ik}\right)=F_{1}\left(\rho_{0},\mu_{1}^{ik}\right)\,,\qquad
L\left(\rho_{2}^{iklm}\right)=F_{2}\left(\mu_{0},\mu_{1}^{ik},\mu_{2}^{iklm}\right)$
(81)

for brevity. It is important to note that we do not actually need to solve the
equation for $\rho_{2}^{iklm}$ as it only appears in the Lagrangian when
multiplied by a constant, thus we can directly use the normalisation condition
(14) and short circuit the task of solving it. We also do not need to find
$\mu_{2}^{iklm}$ for a slightly different reason, which is that it appears in
the Lagrangian to too high order and so is always discarded anyway. Thus our
only task now is to determine the first order correction to the chemical
potential $\mu_{1}^{ik}$ and the first order solution to $\tilde{\rho}$ which
is $\rho_{1}^{ik}$.

We first state $\rho_{1}^{ik}\equiv\rho_{1}^{ik}(\mathbf{r},\mu_{1}^{ik})$,
i.e. it is now a function of the chemical potential. We specifically choose
the functional form
$\rho_{1}^{ik}=\hat{\rho}_{1}^{ik}+\mu_{1}^{ik}\rho_{1}^{\prime}$, in a kind
of ’homogeneous + inhomogeneous’ fashion. Namely, we demand that

$L\left(\hat{\rho}_{1}^{ik}\right)=\frac{\hbar^{2}}{4m}\left(2\frac{\partial_{ik}\rho_{0}}{\rho_{0}}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}^{2}}\right)+\int_{\Omega}\left[f_{ik}(\mathbf{r}-\mathbf{r}^{\prime})+\delta_{ik}U(\mathbf{r}-\mathbf{r}^{\prime})\right]\rho_{0}(\mathbf{r}^{\prime})\mathop{}\\!\mathrm{d}\mathbf{r}^{\prime}$
(82)

i.e. $\hat{\rho}_{1}^{ik}$ satisfies the equation of motion without any
chemical potential considerations. Then we solve the equation

$L\left(\mu_{1}^{ik}\rho_{1}^{\prime}\right)=-\mu_{1}^{ik}\,\rightarrow\,L\left(\rho_{1}^{\prime}\right)=-1$
(83)

simultaneously, which we can do numerically as it is the same operator but
just a different RHS. We can state this in the notation of (81)

$L\left(\hat{\rho}_{1}^{ik}\right)=F_{1}\left(\rho_{0},\mu_{1}^{ik}=0\right)\quad,\quad
L\left(\rho_{1}^{\prime}\right)=F_{1}\left(\rho_{0}=0,\mu_{1}^{ik}=1\right)$
(84)

which we can readily substitute into the ’homogeneous’ and ’inhomogeneous’
solutions to see that

$L\left(\hat{\rho}_{1}^{ik}+\mu_{1}^{ik}\rho_{1}^{\prime}\right)=F_{1}\left(\rho_{0},\mu_{1}^{ik}\right)$
(85)

which is the ‘general’ solution we were looking for. The final task is to now
actually calculate $\mu_{1}^{ik}$ which we can do by using the normalisation
condition(14). Since we know what the first-order condition is, we can
rearrange to find that

$\mu_{1}^{ik}=-\int_{\Omega}(\rho_{0}\delta_{ik}+\hat{\rho}_{1}^{ik})\mathop{}\\!\mathrm{d}\mathbf{r}\Bigg{/}\int_{\Omega}\rho_{1}^{\prime}\mathop{}\\!\mathrm{d}\mathbf{r}$
(86)

which is the final ingredient.

We have now fully solved for $\tilde{\rho}$ (and discarded the parts that we
don’t need), so we can re-express $\mathcal{L}_{\tilde{\rho}}$ using
integration by parts. The calculation is relatively straightforward and makes
use of the base EoM (2), allowing one to find

$-\mathcal{L}_{\tilde{\rho}}=\frac{1}{2}\left(2\mu_{0}\rho_{2}^{iklm}u_{ik}u_{lm}+\rho_{1}^{ik}\left[\frac{\hbar^{2}}{4m}\left(2\frac{\partial_{ik}\rho_{0}}{\rho_{0}}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}^{2}}\right)+\left(f_{ik}+\delta_{ik}U\right)*\rho_{0}+2\mu_{0}\delta_{ik}+\mu_{1}^{ik}\right]\,\epsilon_{ik}\epsilon_{lm}\right)+\mathop{}\\!\mathcal{O}(\epsilon^{3})$
(87)

which contains all the second-order contributions from $\tilde{\rho}$.

We can combine this with all the second-order contributions that do not
contain $\tilde{\rho}$, and obtain a long expression

$\begin{split}-\mathcal{L}_{2}=\frac{1}{2}&\left[2\left(U*\rho_{0}\right)\rho_{0}M_{iklm}u_{ik}u_{lm}+\left(\left(U*\rho_{0}\right)\rho_{0}\delta_{ik}\delta_{lm}+2\left(f_{ik}*\rho_{0}\right)\rho_{0}\delta_{lm}+\left(W_{iklm}*\rho_{0}\right)\rho_{0}\right)\,\epsilon_{ik}\epsilon_{lm}\right.\\\
&\left.+\left(f_{ik}*\rho_{0}\right)\rho_{0}\Delta_{ik}+2\frac{\hbar^{2}}{4m}\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}}\left(\Delta_{ik}+\chi_{ik}\right)+\frac{\hbar^{2}}{4m}\left(-2\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}}\right)\delta_{lm}\epsilon_{ik}\epsilon_{lm}\right.\\\
&\left.+\frac{\hbar^{2}}{4m}\frac{\left(\nabla\rho_{0}\right)^{2}}{\rho_{0}}M_{iklm}u_{ik}u_{lm}+2\mu_{0}\rho_{0}\left(\delta_{ik}\delta_{lm}-M_{iklm}\right)u_{ik}u_{lm}\right.\\\
&\left.+\rho_{1}^{lm}\left[\frac{\hbar^{2}}{4m}\left(2\frac{\partial_{ik}\rho_{0}}{\rho_{0}}-\frac{\partial_{i}\rho_{0}\partial_{k}\rho_{0}}{\rho_{0}^{2}}\right)+\left(f_{ik}+\delta_{ik}U\right)*\rho_{0}+2\mu_{0}\delta_{ik}+\mu_{1}^{ik}\right]\,\epsilon_{ik}\epsilon_{lm}\right]\,.\end{split}$
(88)

which will be used in the calculation of the elastic constants.

To summarise, we have expanded our Lagrangian up to second order in the strain
tensor and collected all the terms. We’ve solved for $\tilde{\rho}$ such that
the new state is a ground state with the strained normalisation condition, and
subsequently re-expressed the second-order Lagrangian with that solution in
hand.

##  Appendix B Connecting Cauchy elastic tensor to expansion

This section closely follows the calculation done by Bavaud et al [28]. To
begin with, we define an analytic expansion in the free energy of our system
around a deformation
$x_{\alpha}^{\prime}=\left(\delta_{\alpha\beta}-u_{\alpha\beta}\right)x_{\beta}$
with a mapping $\Omega\mapsto\Omega^{\prime}$

$F(\Omega^{\prime})=F(\Omega)-\absolutevalue{\Omega}\pi_{\alpha\beta}(\Omega)u_{\alpha\beta}+\absolutevalue{\Omega}\frac{1}{2}A_{\alpha\beta\gamma\delta}(\Omega)u_{\alpha\beta}u_{\gamma\delta}$
(89)

where

$\pi_{\alpha\beta}(\Omega)=\left.\frac{1}{\absolutevalue{\Omega}}\frac{\partial
F(\Omega^{\prime})}{\partial\left(-u_{\alpha\beta}\right)}\right|_{u=0}=\left.-\frac{1}{\absolutevalue{\Omega}}\frac{\partial
F(\Omega^{\prime})}{\partial u_{\alpha\beta}}\right|_{u=0}$ (90)

and similarly

$A_{\alpha\beta\gamma\delta}(\Omega)=\left.\frac{1}{\absolutevalue{\Omega}}\frac{\partial^{2}F(\Omega^{\prime})}{\partial
u_{\alpha\beta}u_{\gamma\delta}}\right|_{u=0}\quad.$ (91)

To clarify, we map between spaces $\Omega\mapsto\Omega^{\prime}$ through the
transformation $\mathbf{x}\mapsto\mathbf{x^{\prime}}$.

It is important to note that the tensor $A_{\alpha\beta\gamma\delta}$ is not
the elastic tensor. It is the expansion of the free energy due to the
deformation, but the elastic tensor strictly connects the stress tensor to the
strain. We are searching for Hooke’s Law which maps strain in a deformed point
to strain in an undeformed point but now with the notation

$\pi_{ik}(\Omega^{\prime})=\pi_{ik}(\Omega)-B_{iklm}(\Omega)u_{lm}\quad.$ (92)

An important additional feature of the elastic tensor $B_{iklm}$ is that, for
a sufficiently symmetric system, it has certain symmetry properties as
described by Landau & Lifshitz [27], whereas the expansion tensor $A_{iklm}$
does not necessarily have these symmetries.

We first define an additional two deformations, $u^{\prime}$ and
$u^{\prime\prime}$ such that

$x_{\alpha}^{\prime\prime}=\left(\delta_{\alpha\beta}-u_{\alpha\beta}^{\prime}\right)x_{\beta}^{\prime}\quad;\quad
x_{\alpha}^{\prime\prime}=\left(\delta_{\alpha\beta}-u_{\alpha\beta}^{\prime\prime}\right)x_{\beta}$
(93)

To clarify, we have an initial space $\Omega$ which we can deform with the
transformation $\mathbf{u}$ into the space $\Omega^{\prime}$. We can then
deform from $\Omega^{\prime}\mapsto\Omega^{\prime\prime}$ with the deformation
$\mathbf{u}^{\prime}$, and we can also deform to the $\Omega^{\prime\prime}$
space from the original space like $\Omega\mapsto\Omega^{\prime\prime}$ with
the deformation $\mathbf{u}^{\prime\prime}$. We can combine the deformations
to find the relation

$u_{\alpha\beta}^{\prime\prime}=u_{\alpha\beta}+u_{\alpha\beta}^{\prime}-u_{\alpha\sigma}^{\prime}u_{\sigma\beta}\quad.$
(94)

We then use the above relation to expand the stress tensor that connects the
$\Omega^{\prime}$ and $\Omega^{\prime\prime}$ space. Examining the expansion
of the free energy from the $\Omega^{\prime}$ space to the
$\Omega^{\prime\prime}$ space, we find

$F(\Omega^{\prime\prime})=F(\Omega^{\prime})-\absolutevalue{\Omega^{\prime}}\pi_{\alpha\beta}(\Omega^{\prime})u_{\alpha\beta}^{\prime}+\frac{\absolutevalue{\Omega^{\prime}}}{2}A_{\alpha\beta\gamma\delta}(\Omega^{\prime})u_{\alpha\beta}^{\prime}u_{\gamma\delta}^{\prime}\,.$
(95)

Expanding $\pi_{\alpha\beta}(\Omega^{\prime})$ and using (90) and (94), we
find

$\begin{split}\pi_{\alpha\beta}\left(\Omega^{\prime}\right)=&\left.-\frac{1}{\absolutevalue{\Omega^{\prime}}}\frac{\partial
F(\Omega^{\prime\prime})}{\partial
u_{\alpha\beta}^{\prime}}\right|_{u^{\prime}=0}\\\
=&\left.-\frac{1}{\absolutevalue{\Omega^{\prime}}}\frac{\partial
F(\Omega^{\prime\prime})}{\partial
u_{\gamma\epsilon}^{\prime\prime}}\frac{\partial
u_{\gamma\epsilon}^{\prime\prime}}{\partial
u_{\alpha\beta}^{\prime}}\right|_{u^{\prime}=0}\\\
=&\left.-\frac{1}{\absolutevalue{\Omega^{\prime}}}\frac{\partial
F(\Omega^{\prime\prime})}{\partial
u_{\gamma\epsilon}^{\prime\prime}}\left(\delta_{\gamma\alpha}\delta_{\beta\epsilon}-\delta_{\gamma\alpha}\delta_{\beta\sigma}u_{\sigma\epsilon}\right)\right|_{u^{\prime}=0}\\\
=&\left.-\frac{1}{\absolutevalue{\Omega^{\prime}}}\left(\delta_{\beta\epsilon}-u_{\beta\epsilon}\right)\frac{\partial
F(\Omega^{\prime\prime})}{\partial
u_{\alpha\epsilon}^{\prime\prime}}\right|_{u^{\prime}=0}\,.\end{split}$ (96)

We can now turn our attention to the derivative in the above expression.
Namely, since the derivative actually does not depend on $u^{\prime}$
explicitly (recall that we can map from $\Omega\mapsto\Omega^{\prime\prime}$
through $u^{\prime\prime}$ without needing $u^{\prime}$ at all), we can remove
the evalutation at $u^{\prime}=0$. As such, the derivative is now not a
number, but a function instead! This function describes the energy response of
the mapping $\Omega\mapsto\Omega^{\prime\prime}$ through $u^{\prime\prime}$
for any arbitrary $\Omega^{\prime\prime}$ and $u^{\prime\prime}$. That means
it must be equivalent to the function which maps from
$\Omega\mapsto\Omega^{\prime}$ through $u$. Thus we can write

$\left.\frac{\partial F(\Omega^{\prime\prime})}{\partial
u_{\alpha\epsilon}^{\prime\prime}}\right|_{u^{\prime}=0}=\frac{\partial
F(\Omega^{\prime})}{\partial u_{\alpha\epsilon}}$ (97)

and the final result from (96) becomes

$\pi_{\alpha\beta}\left(\Omega^{\prime}\right)=-\frac{1}{\absolutevalue{\Omega^{\prime}}}\left(\delta_{\beta\epsilon}-u_{\beta\epsilon}\right)\frac{\partial
F(\Omega^{\prime})}{\partial u_{\alpha\epsilon}}\quad.$ (98)

Since the derivative is no longer evaluated at $u=0$, we now substitute in
(89) to find

$\frac{\partial F(\Omega^{\prime})}{\partial
u_{\alpha\epsilon}}=-\absolutevalue{\Omega}\pi_{\alpha\epsilon}(\Omega)+\absolutevalue{\Omega}A_{\alpha\epsilon\gamma\delta}(\Omega)u_{\gamma\delta}\quad.$
(99)

We now also expand the $1/\absolutevalue{\Omega^{\prime}}$ term to find

$\frac{1}{\absolutevalue{\Omega^{\prime}}}=\frac{1}{\absolutevalue{\Omega}}\frac{1}{\det(\mathds{1}-\mathbf{U})}=\frac{1}{\absolutevalue{\Omega}}\left(1+\delta_{xy}u_{xy}+\mathop{}\\!\mathcal{O}(2)\right)$
(100)

where the indices $x$ and $y$ are tensor indices of the same nature as
$\alpha$, $\beta$ etc. It is important to note that we are discarding terms of
$\mathop{}\\!\mathcal{O}(u^{2})$ or higher as we are trying to derive Hooke’s
law, i.e. $F=kx$, and we do not want anharmonic terms.

We can now combine the results of (99) and (100) into (98) (discarding
$\mathop{}\\!\mathcal{O}(u_{ik}^{2})$ terms) to find

$\pi_{\alpha\beta}(\Omega^{\prime})=\pi_{\alpha\beta}-\left(A_{\alpha\beta\gamma\delta}+\delta_{\gamma\beta}\pi_{\alpha\delta}-\delta_{\gamma\delta}\pi_{\alpha\beta}\right)u_{\gamma\delta}$
(101)

where everything on the RHS is a function of $\Omega$. This is the expression
of Hooke’s Law as described in (92)! We can therefore find that

$B_{iklm}=A_{iklm}+\pi_{im}\delta_{kl}-\pi_{ik}\delta_{lm}$ (102)

with now the final task being to relate $\pi_{ik}$ to the pressure tensor
$P_{ik}$. Recalling that the stress tensor is just the negative of the
pressure, we can write the final expression as

$B_{iklm}=A_{iklm}+P_{ik}\delta_{lm}-P_{im}\delta_{kl}$ (103)

which is the Cauchy elastic tensor. For the case of $C_{3}$ symmetry, we have
the relation [27]:

$B_{xxxx}=B_{xxyy}+2B_{xyxy}\,.$ (104)

##  Appendix C Tensor differentiation

We have shown in previous sections how to expand a Lagrangian under a strain
tensor. Once we have that expansion, how do we recover the homogenised tensors
as first and second differentives: $\pi_{ik}=\partial F/\partial u_{ik}$ and
$A_{iklm}=\partial^{2}F/\partial u_{ik}\partial u_{lm}$?

We will now write out explicitly how expressions involving the tensors listed
in (9) differentiate. Consider first $F=F_{ik}u_{ik}$ with $F_{ik}=F_{ki}$,
its first derivative is: $\partial F/\partial u{ik}=F_{ik}$. Similarly, the
first derivative of $F=F_{ik}\epsilon_{ik}$ is $\partial F/\partial
u_{ik}=-F_{ik}$.

Let us now consider second derivatives. Suppose there is a term in the
Lagrangian of the form $F=F_{iklm}u_{ik}u_{lm}$. We have the result
$\partial^{2}F/\partial u_{ik}\partial u_{lm}=2F_{iklm}$. Consider now
$F=F_{ik}\chi_{ik}$ where
$\chi_{ik}=\left(u_{li}u_{kl}+u_{lk}u_{il}\right)/2$. Since $\chi_{ik}$ is
symmetric under $i\leftrightarrow k$, we can impose that $F_{ik}$ is also
symmetric under the same exchange.

$\displaystyle\frac{\partial^{2}\chi_{ik}}{\partial u_{cd}\partial u_{ab}}$
$\displaystyle=\frac{1}{2}\frac{\partial}{\partial
u_{cd}}\left(\delta_{la}\delta_{ib}u_{kl}+\delta_{la}\delta_{kb}u_{il}+\delta_{ka}\delta_{lb}u_{li}+\delta_{ia}\delta_{lb}u_{lk}\right)=\frac{1}{2}\left(\delta_{ib}\delta_{kc}\delta_{ad}+\delta_{kb}\delta_{ci}\delta_{ad}+\delta_{ka}\delta_{bc}\delta_{id}+\delta_{ia}\delta_{bc}\delta_{kd}\right)\,.$
(105)

This allows us to obtain

$\frac{\partial^{2}F}{\partial u_{ik}\partial
u_{lm}}=\frac{1}{2}\left(F_{bc}\delta_{ad}+F_{cb}\delta_{ad}+F_{da}\delta_{bc}+F_{ad}\delta_{bc}\right)=F_{kl}\delta_{im}+F_{im}\delta_{kl}\,.$
(106)

using the symmetry of $F_{ik}$ under exchange of indices. A similar
calculation can be done for other tensor contractions involving the tensors in
(9). This is listed in Table 1.

Table 1:  $F$ | $F_{iklm}\epsilon_{ik}\epsilon_{lm}$ | $F_{ik}\Delta_{ik}$ | $F_{ik}\omega_{ik}$ | $F_{ik}\chi_{ik}$
---|---|---|---|---
$\displaystyle{\frac{\partial^{2}F}{\partial u_{ik}\partial u_{lm}}}$ | $2F_{iklm}$ | $F_{il}\delta_{km}$ | $F_{km}\delta_{il}$ | $F_{kl}\delta_{im}+F_{im}\delta_{kl}$

It is important to briefly discuss the symmetry of each of these expressions.
For terms that look like $F_{iklm}\epsilon_{ik}\epsilon_{lm}$, these already
have all the symmetry properties we would want from our elastic tensors.
However, we can immediately see that the contribution from
$F_{ik}\Delta_{ik}$, $F_{ik}\omega_{ik}$ and $F_{ik}\chi_{ik}$ do not have the
requisite symmetries. We instead need to create these symmetries from
combinations of these tensors. This cannot be achieved a priori by adding in
whatever we like. However, if we expect elastic theory to hold, then we would
expect any expansion in any strain tensor (including ours) to lead us
eventually back to symmetric expressions. This is a good test to check whether
the expansion calculation has been done correctly.

##  Appendix D Numerical Deformations

In our numerics, we use a spatial discretisation consisting of a triangular
lattice as a grid so that the numerical scheme is consistent with the symmetry
of the supersolid phase under rotations. Therefore, each “pixel” is a
parallelogram and the simulation cell is also a parallelogram with basis
vectors $\mathbf{a}$ and $\mathbf{b}$ parallel to $(\sqrt{3}/2,1/2)$ and
$(0,1)$ in Cartesian coordinates. This appendix provides some algebraic
details on how deformations $\mathbf{V}$ of the simulation cell can be used to
extract the elastic tensors numerically using the energy change $F_{2}\equiv
Pu_{ll}+\frac{1}{2}A_{iklm}u_{ik}u_{lm}$ up to second order in the strain.

Let us write all positions and displacements in the basis of the triangular
grid: $\mathbf{r}=r_{a}\mathbf{a}+r_{b}\mathbf{b}$ and
$\mathbf{u}=u_{a}\mathbf{a}+u_{b}\mathbf{b}$. The strain tensor in this basis

$\mathbf{V}=\begin{pmatrix}\partial u_{a}/\partial r_{a}&\partial
u_{a}/\partial r_{b}\\\ \partial u_{b}/\partial r_{a}&\partial u_{b}/\partial
r_{b}\end{pmatrix}$

is related to the strain tensor in the Cartesian basis:
$(\mathbf{U})_{ij}\equiv u_{ij}$ by:

$\mathbf{U}=\mathbf{T}^{-1}\mathbf{V}\mathbf{T}\,,\qquad\mbox{with}\qquad\mathbf{T}=\begin{pmatrix}\frac{2}{\sqrt{3}}&0\\\
-\frac{1}{\sqrt{3}}&1\end{pmatrix}\,.$ (107)

In our numerics, we apply strains diagonal in the triangular-grid basis by
deforming the simulation cell in the directions of $\mathbf{a}$ and
$\mathbf{b}$. Using the above, we can calculate the energy changes we should
expect to see. There are three deformations that give linearly independent
results.

$\displaystyle\mathbf{V}=\begin{pmatrix}\epsilon&0\\\
0&0\end{pmatrix}\\!:F_{2}$
$\displaystyle=P\epsilon+\left(A_{xxxx}+\frac{A_{xyxy}}{3}\right)\frac{\epsilon^{2}}{2}$
$\displaystyle\mathbf{V}=\begin{pmatrix}0&0\\\
0&\epsilon\end{pmatrix}\\!:F_{2}$
$\displaystyle=P\epsilon+\left(A_{yyyy}+\frac{A_{xyxy}}{3}\right)\frac{\epsilon^{2}}{2}$
(108) $\displaystyle\mathbf{V}=\begin{pmatrix}\epsilon&0\\\
0&\epsilon\end{pmatrix}\\!:F_{2}$
$\displaystyle=2P\epsilon+\left(\frac{A_{xxxx}+A_{yyyy}}{2}+A_{xxyy}\right)\epsilon^{2}\,.$

For our system with $C_{3}$ symmetry, there are only two independent elastic
constants. Using (103) and (104), we see that $A_{xxxx}=A_{xxyy}+P+2A_{xyxy}$.
Therefore, we have three equations to solve for the elastic constants and the
pressure for a given $\epsilon\ll 1$.

We have calculated these quantities at several values of $\epsilon$ to check
for the self-consistency of this analysis. This gives a measure of the
numerical errors of this technique.

A good question to ask is whether or not we are allowed to a priori assume the
system is symmetric, and whether we are biasing towards a symmetric result.
This is a valid question, but really what we are doing is seeking to compare
the results we obtain with the analytical theory derived. Importantly, the
theory assumes nothing about symmetry, so if the numerics agree with the
theory, then this assumption is justified.