# Quantum dislocations: the fate of multiple vacancies in two dimensional solid He

###### Abstract

Defects are believed to play a fundamental role in the supersolid state of He. We have studied solid He in two dimensions (2D) as function of the number of vacancies , up to 30, inserted in the initial configuration at Å, close to the melting density, with the exact zero temperature Shadow Path Integral Ground State method. The crystalline order is found to be stable also in presence of many vacancies and we observe two completely different regimes. For small , up to about 6, vacancies form a bound state and cause a decrease of the crystalline order. At larger , the formation energy of an extra vacancy at fixed density decreases by one order of magnitude to about 0.6 K. In the equilibrated state it is no more possible to recognize vacancies because they mainly transform into quantum dislocations and crystalline order is found almost independent on how many vacancies have been inserted in the initial configuration. The one–body density matrix in this latter regime shows a non decaying large distance tail: dislocations, that in 2D are point defects, turn out to be mobile, their number is fluctuating, and they are able to induce exchanges of particles across the system mainly triggered by the dislocation cores. These results indicate that the notion of incommensurate versus commensurate state loses meaning for solid He in 2D, because the number of lattice sites becomes ill defined when the system is not commensurate. Crystalline order is found to be stable also in 3D in presence of up to 100 vacancies.

###### pacs:

67.80.-s, 67.80.bd, 67.80.dj## 1 Introduction

Solid He is the subject of many experimental and theoretical studies[1, 2, 3] since the discovery of non classical rotational inertia[4] (NCRI), an expected manifestation of supersolidity.[5] Supersolidity is a new state of matter in which spatial order and off–diagonal long range order (ODLRO) are present at the same time, implying some form of superfluid properties, and was already predicted long ago.[6, 7] This novel state of matter attracts interest also in the prospect of finding it in cold atoms in optical lattices.[8] However the precise nature of solid He at low temperature is still elusive. Experimental evidence suggests that defects play an important role in NCRI, but which kind of disorder is responsible for the anomalous properties of solid He is still not clear. Many different defects have been considered, but none of the proposed models has shown to be able to capture all the phenomenology of supersolidity. For example the stiffening of He below the “transition” temperature[9] suggests dislocations as a candidate, but this possibility has difficulty in explaining the NCRI seen in solid He in Vycor[4] or in Aerogel.[10] Grain boundaries have been considered too,[1] but NCRI has been observed also in single crystals.[11]

Defected solid He systems have been studied by means of Quantum Monte Carlo (QMC) techniques.[12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] Such QMC methods are by construction equilibrium methods, so defects have to be stabilized by periodic boundary conditions (via a suitable choice of the simulation box) or by fixing the degrees of freedom of a number of atoms surrounding the interested defect. This constraint on the configurational space usually does not prevent the study of the physical properties of the defected system, and the results are generally in a good quantitative agreement with experimental observations; see, for example, the estimation of the vacancy activation energy[12] or the binding energy to the dislocation core of a He atom.[22] Exact QMC results, both at finite[17, 18, 19] and zero temperature[24], show that 2 and 3 vacancies form a bound state. A systematic study of many vacancies in solid He is still lacking; questions like whether the crystalline structure is stable in presence of a large number of vacancies, whether the vacancy–vacancy interaction is pairwise additive (or almost so), whether many vacancies lead to phase separation as suggested in Ref. [17, 19], or whether vacancies turn themselves into other kinds of defects are, according to us, still unanswered.

Since the earliest days of supersolidity quantum lattice models have been considered.[25] A terminology borrowed from lattice models is currently used for solid He in terms of a commensurate state (i.e. a perfect crystal in which there is an integer number of atoms per unit cell) or an incommensurate one (a crystal with non integer probability of occupation of the unit cell). QMC at K of solid He with one or few vacancies in the simulation cell beautifully[16] confirms the picture of an incommensurate state: the vacancy(-ies) is (are) mobile and the number of density maxima is equal to the number of particles plus the number of vacancies. But are we sure that we understand vacancies at a finite concentration in a macroscopic system? Is the option commensurate or incommensurate really appropriate for solid He? One should keep in mind that in a lattice model or in cold atoms in an optical lattice the lattice periodicity is externally imposed, whereas in solid He the same dynamical entities, the atoms, have to build up the periodicity as well as the delocalization required for supersolidity. Such “mermaid” aspect of the atoms is unique to a quantum solid. The purpose of this article is to address such issues with exact QMC methods at K in two dimensions (2D).

We find that size and box commensuration effects are very pronounced, so in order to be able to explore very large distances in a Monte Carlo calculation, up to 100 Å, we have mainly studied solid He in 2D. In this context, the 2D solid is also of interest as a simple model for out of registry solid He adsorbed on planar substrates, such as graphite or silica glasses. New experimental investigations[26] are under way to answer the question whether a supersolid state is present also in adsorbed He. We have studied 2D solid He systems at K with a number of atoms up and above thousand containing up to 30 vacancies. We find that crystalline order is stable also in presence of a large number of vacancies, at least in the range % of vacancy concentration studied by us. Multiple vacancies are highly correlated with a preference to form linear structures, but in presence of 10 or more vacancies the scars in the lattice are healed away transforming the vacancies mainly into dislocations. Such dislocations are not permanent structures, but are mobile and their number is fluctuating. This allows exchange of particles across the system, not only in the cores of dislocations, giving the possibility of establishing a well defined, at least local, phase[27] and indeed, in the simulated systems, we find a non decaying large distance tail in the one–body density matrix. This hints at the possibility that an extended system is supersolid but more work is needed to corroborate this conclusion.

We do not address here the origin of vacancies/dislocations, i.e. if the true ground state of solid He contains or not zero–point defects. Our finding is that solid He, at least in 2D, is much more complex than the picture derived from the notion of commensurate vs incommensurate state because the number of lattice sites is an ill defined quantity when dislocations are present. The number of lattice sites is a meaningful quantity, at least for the sizes we are able to simulate, only if the number of particles and the simulation box are such that the regular lattice exactly fits in or the misfit is limited to just one or very few vacancies. We have performed few simulations also in three dimensions (3D). Also in this case we find that, near melting, crystalline order is stable even in presence of 100 vacancies, but we have not yet performed a detailed characterization of the disorder in the system.

## 2 The Spigs Method

We employ the exact K Shadow Path Integral Ground State (SPIGS)[15] method, an “exact” QMC technique.[28, 29] SPIGS is an extension of the Path Integral Ground State (PIGS)[30] method. The aim of PIGS is to improve a variationally optimized trial wave function by constructing a path in the Hilbert space of the system which connects the given to the lowest energy state of the system, , constrained by the choice of the number of particle , the geometry and the boundary conditions of the simulation box and the density . During this “path”, the correct correlations among the particles arise through the “imaginary time evolution operator” , where is the Hamiltonian operator. For a large enough , an accurate representation for the lowest energy state wave function is given by , which can be written analytically by discretizing the path in imaginary time and exploiting the factorization property . In this way, turns out to be expressed in terms of convolution integrals which involve the “imaginary time propagator” for small , for which very accurate approximations can be found in literature.[31, 32] This maps the quantum system into a classical system of open polymers.[30] An appealing feature peculiar to the PIGS method is that, in , the variational ansatz acts only as a starting point, while the full path in imaginary time is governed by , which depends only on the Hamiltonian operator. We have recently shown that PIGS results for large enough are unaffected by the choice of both in the liquid and in the solid phase[28, 29] thus providing an unbiased exact K QMC method. Within SPIGS a shadow wave function (SWF)[33, 34] is taken as and this choice was shown to greatly accelerate convergence to . Another appealing feature of the SPIGS method is that it recovers the solid phase via a spontaneously broken translational symmetry[15] and there is no constriction on the atomic positions, so it is particularly indicated in studying crystals with defects.

## 3 Model System and Simulation Details

By studying the equation of state we find the freezing and the melting densities to be, respectively, Å and Å. As simulation box we choose then a rectangular box fitting a regular triangular lattice with lattice sites at Å. Periodic boundary conditions (pbc) are applied in both directions. Dealing with low temperature properties, He atoms are described as structureless zero–spin bosons, interacting through the HFDHE2 Aziz potential.[35] In order to avoid any approximation associated with the estimation of tail corrections due to the finite size of the simulation boxes, we have considered a truncated and shifted Aziz potential which goes to zero at Å. Since is well below the minimum size of each considered box, the results relative to different box sizes can be compared without any further correction. Our SPIGS computation is based on the pair–product approximation for the imaginary time projector[31] and the imaginary time step K has been chosen in order to ensure a good accuracy and a reasonable computational effort.[28]

An important parameter in the computation is the total projection time . Any trial wave function can be expanded in terms of the exact ground state and of the excited states of zero total linear momentum. If we call the minimum excitation energy for such excited states, has to be larger than in order to project out any excitation contribution in ; the rate of rejection of such contributions depends also on the quantity one is computing due to sensitivity to missing long range correlations in or to wrong short range correlations. In the perfect crystal, using K as energy unit, K was found[28] to be large enough to ensure convergence of both diagonal and off–diagonal observables when a SWF is used as . This is compatible with the criterion . In fact, in a perfect crystal, the lowest energy excitations are phonons; thus corresponds to the combination of two phonons of opposite with the smallest achievable value, i.e. in a box with pbc where is the largest side. Via a new analytic continuation method[36] we have estimated the longitudinal sound velocity for 2D solid He at Å to be KÅ. Then the minimum imaginary time requested for convergence is K in the largest simulation box considered here. This value is well below the used one K. Since the transverse sound velocity is typically similar to , the considered value is enough to ensure also the removal of transverse phonon contribution. The situation might be different in presence of vacancies, or other kinds of defects, because novel low energy excitation modes could be present. We have performed several tests in order to be confident that also for defected crystals the considered value of is large enough. The most straightforward one is to perform simulations with larger values and verify that the expectation values are compatible within the statistical uncertainty. Such a test turns out to be very time consuming and in practice it has been feasible only for a number of He atom , where we have increased the value of up to K finding good convergence for diagonal properties. An indirect test on the convergence is provided by the behavior of suitably chosen observables during the imaginary time path. If we define

(1) |

for , when is large enough we have and Eq. 1 gives the mixed expectation value as a function of the partial imaginary time . If has a plateau when , we are guaranteed that is large enough for convergence and the plateau represents the “exact” expectation value of . A couple of examples are shown in Fig.1 for the potential energy () and for the static structure factor at the lowest wave vectors (, where is the density fluctuation operator) along the simulation box axis.

As far as the potential energy is concerned, excitation modes are rapidly projected out and K is enough to achieve convergence both in the perfect and in the defected systems, as shown in Fig.1 by the presence pf a plateau of in both cases. Larger values of the imaginary time are needed to reach convergence for quantities that depend on long–range correlations, such as the static structure factor at small wave vectors () .[37] In the perfect crystal, convergence for is achieved at about K, as shown in Fig.1. In the defected crystal (especially when many vacancies are present as in the case considered in the figure) the chosen K value is just enough to get convergence. An interesting physical effect emerges from Fig.1, the static structure factor is strongly affected by the presence of defects especially along the nearest neighbor direction.

Two independent simulations are needed in order to calculate the formation energy of vacancies, one with particles (the regular ideal crystal, dubbed also perfect crystal), and one with less atoms (defected crystal). In this second case we say that we have vacancies, even if in general we can talk of vacancies only in the starting configuration of the simulation, where atoms are removed from the perfect lattice. In fact, in SPIGS, like in PIMC, there is no constraint on the atomic positions, so that vacancies are free to transform into different kinds of defects. In presence of vacancies, after removing particles from the ideal starting configuration, we rescale the dimensions of the simulation box to reset the system to the original density. This is performed to circumvent the need of correcting the energy due to density changes caused by the inclusion of vacancies. The formation energy of vacancies is proportional to the difference[12, 38] between the energy per particle of the defected crystal, , and of the perfect one, , i.e.

(2) |

When vacancies turn into dislocations this has the meaning of defect formation energy at fixed density. One can consider also the formation energy at a fixed lattice parameter; this is given by

(3) |

being the chemical potential, that for our system turns out to be K.

We check the presence of crystalline order by monitoring the static structure factor for the presence of Bragg peaks. In addition, we compute the particle coordination number, via a Delaunay triangulation[39] (DT) of the sampled configurations, in order to estimate the amount of disorder in the system.[40] In a perfectly ordered 2D triangular crystal each atom is linked to 6 other atoms in the DT. Atoms with number of coordination not equal to 6 are then a measure of the disorder in the crystal.[40] A conservation law exists for the coordination numbers:[40]

(4) |

where is the number of coordinated atoms, so that we can consider only with . Therefore we take as an estimate of disorder in the system (we never found 3-sided polygons, i.e. ). In a perfect () 2D quantum crystal the coordination is 6 only on the average, in fact, due to the large zero point motion, fluctuations are present such that atoms do not have always coordination 6. Thus, a more useful quantity measuring the net amount of disorder in a crystal with defects is the difference between the observed and that of the corresponding perfect crystal:

(5) |

When the number of vacancies is small, up to 6, we have also computed the vacancy–vacancy correlation function by recording the relative distances among vacancies during the Monte Carlo sampling. The determination of the vector positions of the vacancies in a crystalline configuration is far from being trivial due to the large zero point motion, to high vacancy mobility and because in our algorithm the center of mass is not fixed. Vacancy positions are obtained via the coarse–graining procedure illustrated in Ref. [24]. This method is efficient as long as we have few vacancies. For large the efficiency in recognizing the position of vacancies becomes extremely poor (in fact, as we shall see, vacancies turn into dislocations), then is no more possible to compute . Different definitions of vacancy positions, like the one in Ref. [18], give very similar results.[41]

A special care has been devoted in ruling out possible metastability effects. We have considered different starting configurations for vacancies, i.e. we have removed a number of atoms forming a compact cluster or a linear cluster, or the removed atoms come from random lattice positions. After long enough equilibration (our simulations are never shorter than Monte Carlo steps) we find agreeing results for systems with the same value. In order to rule out possible finite size effects and pbc bias we have considered systems of increasing sizes ( values). Again we find compatible results for systems of different sizes and with equal value, so no appreciable finite size effects have been detected. The only exception is the one–body density matrix that will be discussed below. Furthermore we have considered also crystals with no principal crystalline axis parallel to the box sides obtaining once more the same results. Moreover, we have considered systems where the starting configuration was obtained by removing particles from an equilibrated configuration of the perfect crystal and at the same time the positions of one (two) line(-s) of atoms, parallel to one (both) box side(-s) were kept fixed during the Monte Carlo sampling. Even in this cases, the results do not change, for instance, for large values, vacancies are found to turn into dislocations with the same features as in the fully mobile systems.

## 4 Results

We have studied a triangular crystal at Å with , 480, 960 and 1440 lattice positions of the perfect crystal. Our results for the formation energy at constant density are plotted in Fig. 2. The dependence of on is monotonic and systematically sublinear, confirming the existence of an attractive interaction among vacancies. Systems with the same and different have the same within the statistical uncertainty, i.e. has no significant dependence on , at least in the range % that we have studied. From the plot in Fig. 2 it is possible to recognize two different behaviors: for , deviates rapidly from the linear dependence and vacancies form a bound state as indicated by a (roughly) exponentially decreasing correlation function , plotted in Fig. 3. For the ratio remains practically constant with a value of about 0.61 K. This means that, when some vacancies are already present in the system, the creation of an additional vacancy has a very low cost, about one tenth of the cost of a single vacancy. Notice that in the linear regime , at fixed lattice parameter, is negative, there is a gain of 13.61 K for each additional vacancy.

The double regime behavior of is reflected also in other properties of the system, for instance in the static structure factor . For a finite crystal the height of the Bragg peaks has a strong dependence on the number of particle and shows a slow convergence (like in 3D[42] and in 2D) to the thermodynamic limit. Note that being at K, the crystalline order is stable also in 2D. Our computation of in the perfect crystal verifies this dependence as expected for a 2D crystal. This dependence of arises from missing phonon modes for in a finite box so that essentially the same contribution is expected to be present in the defected system. Therefore by taking the difference between in the defected system and the value in the prefect crystal such contribution cancels out to a large extent so we can compare results for crystal of different sizes and . Since the crystalline lattice relaxes around a vacancy, one expects that the heights of the Bragg peaks in decrease with increasing number of vacancies. This is observed, as shown in Fig. 4, only when is small. For large values we find a very different behavior. The main Bragg peak height has some broadening but its integrated intensity[43] does not show a significant dependence on , rather it oscillates around an almost constant value. In Fig. 5 we show the static structure factor in the plane obtained for , 5, 15 and 25 in the largest considered box with lattice sites. When the number of vacancies is small (Fig. 5a) we find Bragg peaks in the first reciprocal vector star as sharp as in the case of the perfect crystal, within the resolution due to the finite size of the simulation box, but with a decreased height. Otherwise, see Fig. 5c and d, for large , the Bragg peaks turn out to be slightly broadened, but their integrated intensity (see Fig. 4) shows no significant dependence on . The broadening of the main Bragg peak has a linear dependence on the concentration of defects , with Å.

A similar behavior is displayed also by the amount of disorder , plotted in Fig. 6. increases with for small values and saturates to an almost constant value for larger . By examining sampled configurations of the system in this regime, we find that it is no more possible to identify vacancy positions. Rather one finds dislocations, as shown by a DT of the position of the atoms, see Fig. 7a and b. Dislocation cores in 2D are point defects characterized by a couple 5–7 in the coordinations number of neighboring atoms in the DT. We find that dislocations do not prefer any of the principal lattice directions. Studying a large collection of independent configurations one finds that a large fraction of them contains a couple of dislocation cores (see Fig. 7). For example, for in the box, there is a probability of about 90% for two dislocation cores and 10% of finding more than two dislocation cores. In addition there is a probability around 3% of finding also some isolated vacancies besides dislocations (see Fig. 7c). Similar values are found for in the box.

We stress that, even if initially placed in a compact cluster configuration, vacancies do not separate into a vacancy rich region surrounded by a regular crystal. Rather they mostly reorganize themselves across the system giving rise to dislocations. In this regime of large the quantity that maintains direct physical meaning is , the number of atoms. When atoms rearrange themselves giving rise to dislocations, the number of lattice sites becomes ill defined because the lattice positions of the perfect crystal compatible with the simulation box and the pbc do no more represent equilibrium positions for the particles, and it is no more possible to recognize a regular lattice describing the equilibrium positions of the particles. Still we continue to use and as a convenient measure of deviation from the ideal crystal where crystalline order is perfect. In the regime of small , vacancies are strongly correlated and prefer to form linear configurations, as shown in Fig. 8.

These linear configurations can be considered as forerunners of the quantum dislocations found at larger . This behavior has some similarities with the behavior of colloidal crystals in 2D.[44, 45] We find that dislocations are not fixed structures, rather they are very mobile and their number changes along the simulation. As an example, in Fig. 9 we show the evolution of the DT of a configuration for vacancies in the box with sampled at different Monte Carlo times.

The one–body density matrix has a central role since the presence of a plateau at large indicates the presence of ODLRO, i.e. the system has Bose–Einstein condensation (BEC). Even seems to show the signature of a double regime behavior depending on the number of vacancies. Some of our results for the one–body density matrix are reported in Fig. 10. As found in the 3D case[16], when only a single vacancy is present displays a non–zero plateau, which is the signature of ODLRO. As vacancies are added, the tail of one–body density matrix is depressed, even if it does not show a simple exponential decay. In order to be conclusive on the behavior of in the large distance range one should compute it in larger systems; we have faced these computationally intensive calculations in the range where vacancies give rise to quantum dislocations. Here (see Fig. 10b and Fig. 12c) displays a recognizable plateau at large distances (very evident in the largest system). In order to understand the origin of the observed off–diagonal contributions to in the dislocation regime it is useful to recall that in a SPIGS computation (which is the probability amplitude to destroy a particle in and to create one in ) is obtained by splitting one of the linear polymers in two half polymers,[15] one departing from and the other from . We find that the main contribution to the plateau comes from configurations where at least one of the half–polymers occupies a dislocation core, see Fig. 11(a,b).

This means that the system is able to transfer particles from a quantum dislocation to another. This is quite distinct from the observed superfluidity along the core of a screw dislocation[21] in 3D, in fact in 2D a dislocation core is a point defect and not a linear one as in 3D. Furthermore, the possibility of finding a half–polymer out of the cores like in Fig. 11(b,c) means that quantum dislocations are able also to induce vacancies in the surrounding crystal. Contrarily to diagonal properties, in the computed one–body density matrix we find presence of significant size effects, as can be inferred for example by the significant difference between the in the , system and in the , one, in the intermediate distance range (Å). This is shown more clearly, in Fig. 12 where we report in the full plane, one can notice the presence of ridges in the tail of in the smaller systems ( and , see panels (a) and (b)) that are no more present in the largest system (see Fig. 12c). Such ridges are due to a commensuration effect with the simulation box. We have evidence that in the and systems dislocation cores prefer to stay at a distance of about ( is the lattice parameter) that turns out to be comparable with . This is no more true in the box with , in which case has no ridge in the tail (see Fig. 12c).

In order to determine how the plateau of is affected by finite size effects one should perform computations for an even larger number of particles, bur is the maximum number that presently we can handle. Another reason for studying larger systems is to establish how the number of dislocations scales with the size of the system and with the number of vacancies. Presumably supersolidity and BEC are present in a macroscopic system only if the concentration of dislocations is finite. In addition there is the question if the projection time is large enough to have convergence to the exact result. On the basis of indirect tests performed on diagonal observables which depend on long–range correlations (see Fig. 1) we are rather confident that the used value of is large enough to get convergence, but we are not able to provide a direct test by using a larger value of again for computational limitations (a computation with K with K and atoms has to handle about degrees of freedom). Our results in the dislocation regime hint that ODLRO is present, but further tests are needed in order to firmly establishing this. In any case our results on off–diagonal properties in 2D systems enlighten a novel ability of dislocations in inducing exchanges even in their surroundings.

## 5 Conclusions

In summary, we find that crystalline order in 2D He is stable also in presence of a large number of vacancies and there is no tendency to phase separation. This is true also in 3D, we have studied up to 2548 particles and 98 vacancies and in all cases we find that crystalline order is present, as shown by the presence of Bragg peaks in the static structure factor. When, in the 2D system, the number of vacancies is of order of 10 and above, vacancies inserted in the initial configuration loose their identity and most of them become quantum dislocations. Such dislocations turn out to be very mobile and are able to induce exchange of particles across the system, which is necessary to supersolidity. The ability of dislocations to induce vacancies in the surrounding crystal could be relevant also for the 3D case. If this feature will be confirmed in 3D, the superfluidity would be not restricted to the dislocation cores,[21] but exchange processes necessary for supersolidity would be promoted by vacancies even in the bulk crystal far from dislocation cores. Simulations with numbers of particle orders of magnitude larger are needed to answer the question whether the number of dislocations cores will increase with the system size and with the number of vacancies to give a finite concentration of dislocations able to induce ODLRO even in a macroscopic 2D He crystal, but this is actually beyond our computational possibilities and we have to leave it for future investigations. Moreover, how much the phase correlation triggered by dislocations depends on the concentration of dislocations and its relevance for the 3D solid He are still open questions. Our results mainly indicate that the usual concepts of commensurate or incommensurate solid, borrowed from lattice models, is not appropriate for solid Helium, a system where the crystal lattice is not externally imposed, but is self-induced by the correlation among particles. In fact, the number of lattice sites is found to be an ill defined quantity when the crystal houses many vacancies/dislocations.

We do not address the origin (intrinsic or extrinsic) of the vacancies/dislocations studied in the present paper. The issue if the ground state could contain zero point defects is still debated.[16, 3, 46, 47] In case the ground state of solid He contains defects, our findings suggest that the Andreev–Lifshitz–Chester scenario[6, 7] would in a certain sense revive, but in terms of ground state dislocations rather than vacancies. On the other hand, even if the ground state has no zero point defects, our findings suggest that in presence of extrinsic disorder this will manifest more in terms of fluctuating dislocations rather than vacancies, at least at low and in 2D.

This work was supported by the INFM Parallel Computing Initiative and by the Supercomputing facilities of CILEA. M. Rossi would thank W. Lechner for useful discussions.

## References

## References

- [1] N.V. Prokof’ev, Adv. Phys. 56, 381 (2007).
- [2] S. Balibar and F. Caupin, J. Phys.: Condensed Matter 20, 173201 (2008).
- [3] D.E. Galli and L. Reatto, J. Phys. Soc. Jpn., 77, 111010 (2008).
- [4] E. Kim and M.H.W. Chan, Nature 427, 225 (2004); Science 305, 1941 (2004); Phys. Rev. Lett. 97, 115302 (2006).
- [5] A.J. Leggett, Phys. Rev. Lett. 25, 1543 (1970).
- [6] A.F. Andreev and I.M. Lifshitz, Soviet Physics JETP 29, 1107 (1969).
- [7] G.V. Chester, Phys. Rev. A 2, 256 (1970).
- [8] V.W. Scarola, E. Delmer and S. Das Sarma, Phys. Rev. A 73, 051601(R) (2006).
- [9] J. Day and J. Beamish, Nature 450, 853 (2007).
- [10] N. Mulders, J.T. West, M.H.W. Chan, C.N. Kodituwakku, C.A. Burns and L.B. Lurio, Phys. Rev. Lett. 101, 165303 (2008).
- [11] A.C. Clark, J.T. West and M.H.W. Chan, Phys. Rev. Lett. 99, 135302 (2007).
- [12] F. Pederiva, G.V. Chester, S. Fantoni and L. Reatto, Phys. Rev. B 56, 5909 (1997).
- [13] D.E. Galli and L. Reatto, J. Low Temp. Phys. 124, 197 (2001).
- [14] D.E. Galli and L. Reatto, Phys. Rev. Lett. 90, 175301 (2003).
- [15] D.E. Galli and L. Reatto, Mol. Phys. 101, 1697 (2003); J. Low Temp. Phys. 134, 121 (2004).
- [16] D.E. Galli and L. Reatto, Phys. Rev. Lett. 96, 165301 (2006).
- [17] M. Boninsegni, A.B. Kuklov, L. Pollet, N.V. Prokof’ev, B.V. Svistunov and M. Troyer, Phys. Rev. Lett. 97, 080401 (2006).
- [18] B.K. Clark and D.M. Ceperley, Comp. Phys. Comm. 179, 82 (2008).
- [19] L. Pollet, M. Boninsegni, A.B. Kuklov, N.V. Prokof’ev, B.V. Svistunov, and M. Troyer, Phys. Rev. Lett. 101, 097202 (2008).
- [20] L. Pollet, M. Boninsegni, A.B. Kuklov, N.V. Prokof’ev, B.V. Svistunov, and M. Troyer, Phys. Rev. Lett. 98, 135301 (2007).
- [21] M. Boninsegni, A.B. Kuklov, L. Pollet, N.V. Prokof’ev, B.V. Svistunov and M. Troyer, Phys. Rev. Lett. 99, 035301 (2007).
- [22] P. Corboz, L. Pollet, N.V. Prokof’ev and M. Troyer, Phys. Rev. Lett. 101, 155302 (2008).
- [23] S.G. Soyler, A.B. Kuklov, L. Pollet, N.V. Prokof’ev and B.V. Svistunov, Phys. Rev. Lett. 103, 175301 (2009).
- [24] M. Rossi, E. Vitali, D.E. Galli and L. Reatto, J. Phys.: Conference Series 150, 032090 (2009).
- [25] K.S. Liu and M.E. Fisher, J. Low Temp. Phys. 10, 655 (1973).
- [26] Y. Shibayama, H. Fukuyama and K. Shirahama, J. Phys.: Conference Series 150, 032096 (2009) and preliminary results presented by J. Saunders during the Supersolid 2009 workshop in Banff.
- [27] P.W. Anderson, Phys. Rev. Lett. 100, 215301 (2008).
- [28] E. Vitali, M. Rossi, F. Tramonto, D.E. Galli and L. Reatto, Phys. Rev. B 77, 180505(R) (2008).
- [29] M. Rossi, M.Nava, L. Reatto and D.E. Galli, J. Chem. Phys. 131, 154108 (2009).
- [30] A. Sarsa, K.E. Schmidt and W.R. Magro, J. Chem. Phys. 113, 1366 (2000).
- [31] D.M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
- [32] M. Suzuki, Phys. Lett. A 201, 425 (1995).
- [33] S.A. Vitiello, K. Runge and M.H. Kalos, Phys. Rev. Lett. 60, 1970 (1988).
- [34] S. Moroni, D.E. Galli, S. Fantoni and L. Reatto, Phys. Rev. B 58, 909 (1998).
- [35] R.A. Aziz, V.P.S. Nain, J.S. Carley, W.L. Taylor and G.T. McConville, J. Chem. Phys. 70, 4330 (1979).
- [36] E. Vitali, M. Rossi, L. Reatto and D.E. Galli, arXiv:0905.4406.
- [37] L. Reatto and G.V. Chester, Phys. Rev. 155, 88 (1967).
- [38] Both and are ground state energies of two periodically repeated small systems, and their difference has been used to estimate the formation energy of extra vacancies in bulk, therefore is a derived quantity.[16]
- [39] The Delaunay triangulation (DT) for a set of points in the plane is a triangulation such that no point in is inside the circumcircle of any triangle in . DT corresponds to the dual graph of the Voronoi tessellation, i.e. the partitioning of the plane into convex polygons such that each polygon contains exactly one point in and every point of the plane closer to it than to any other point of . Then each point in is linked only to its nearest neighbors. The number of the sides of Voronoi polygons (of the links in DT) is the coordination number of the generating point.
- [40] B. Krishnamachari and G.V. Chester, Phys. Rev. B 61, 9677 (2000).
- [41] The algorithm proposed in Ref. [18] by construction provides always vacancy positions. For few vacancies the results with this algorithm agree with those obtained with the coarse graining method used here. When many vacancies are present (and turns into dislocations) the vacancy positions given by the algorithm of Ref. [18] are extremely erratic from one configuration to the following one and thus are not reliable.
- [42] E.W. Draeger and D.M. Ceperley, Phys. Rev. B 61, 12094 (2000).
- [43] The broadening of the Bragg peaks occurs on a small number of close points in the -grid (wave vectors in a finite system with pbc are quantized), thus we refer to integrated intensity as the sum of over the -points contributing to the Bragg peak.
- [44] A. Pertsinidis and X.S. Ling, Nature 413, 147 (2001); Phys. Rev. Lett. 87, 098303 (2001).
- [45] W. Lechner and C. Dellago, Soft Mat. 5, 2752 (2009).
- [46] P.W. Anderson, Science 324, 631 (2009).
- [47] M. Rossi, E. Vitali, D.E. Galli and L. Reatto, J. Low Temp. Phys. 153, 250 (2008).