Superfluid and supersolid phases of 4He on the second layer of graphite
aa r X i v : . [ c ond - m a t . o t h e r] M a y Superfluid and supersolid phases of He on the second layer of graphite
M.C. Gordillo
Departamento de Sistemas F´ısicos, Qu´ımicos y Naturales,Universidad Pablo de Olavide. E-41013 Seville, Spain
J. Boronat
Departament de F´ısica, Universitat Polit`ecnica de Catalunya, Campus Nord B4-B5, E-08034 Barcelona, Spain
We revisited the phase diagram of the second layer of He on top of graphite using quantumMonte Carlo methods. Our aim was to explore the existence of the novel phases suggested recentlyin experimental works, and determine their properties and stability limits. We found evidence of asuperfluid quantum phase with hexatic correlations, induced by the corrugation of the first Heliumlayer, and a quasi-two-dimensional supersolid corresponding to a 7/12 registered phase. The 4/7commensurate solid was found to be unstable, while the triangular incommensurate crystals, stableat large densities, were normal.
The light mass of Helium atoms and the strongCarbon-Helium interaction make He adsorbed ongraphite the most paradigmatic example of a two-dimensional (2D) quantum system. Its phase diagramwas extensively studied in the 90s, using a variety of ex-perimental techniques (see, for instance, Ref. 1). Theconsensus so far is that, at very low temperature, He indirect contact with the graphite surface is a √ × √ He layer, stable in the coverage range ∼ − . Those are total densities, including heliumatoms per surface unit both in the first and second lay-ers. Heat capacity [5, 6] and torsional oscillator [7, 8]experiments indicated that, after a promotion from thefirst to the second layer a quasi-two-dimensional liquidwas formed. Increasing the coverage, the liquid changesinto a commensurate phase (with respect to the first ad-sorbed Helium layer), and then to an incommensurateone before promotion to a third layer [5, 6].Remarkably, two recent experimental works have re-opened the doubts about the phase diagram of the sec-ond layer of He adsorbed on graphite. First, the calori-metric data of Ref. 9 suggests the existence of a liquidabove 0.175 ˚A − , followed upon an increasing of the he-lium coverage, of a stable phase in the 0.196-0.203 ˚A − range. That phase could be either a commensurate solid or a quantum hexatic phase. On the other hand, thetorsional oscillator data of Ref. 10 indicate a normal 2Dliquid between ∼ − and 0.1711 ˚A − (Ref. 10,supplementary information), followed by an arrangement showing a superfluid response from 0.1711 to 0.1996 ˚A − ,with a maximum at 0.1809 ˚A − . Since, according to pre-vious DMC calculations [11], those densities are abovethe stability limits of a liquid phase, that response wouldcorrespond to a quasi-two-dimensional supersolid. Thiswould be the first indication of a stable supersolid phasein He, after discarding that possibility in bulk [12].In this Letter, we revisit this problem from a theoret-ical microscopic point of view. Our aim is to clarify thenature of the stable phases of the second layer of He ongraphite, in the limit of zero temperature. Our resultsshow hexatic order [13–15] before crystallization into oneof the possible registered phases (7/12). In both cases,our measure of the superfluid fraction gives a finite value,larger for the hexatic but still very significant for the reg-istered solid. Therefore, on this layer we found two longpursued phases: a superhexatic [16] and a quasi-two di-mensional (registered) supersolid.Our zero-temperature first-principles study relies onthe diffusion Monte Carlo (DMC) method, used exten-sively in the past to analyze He phases in different ge-ometries [17]. The high numerical accuracy of DMC inthe estimation of the energy is crucial to disentangle thestability of different possible phases. This is specially rel-evant in the study of registered phases since the energydifferences between different commensurate solids is verytiny. In essence, the DMC method allows us to solve themany-body imaginary-time Schr¨odinger equation corre-sponding to the Hamiltonian describing the system [21].In the present case, H = N X i =1 (cid:20) − ¯ h m ∇ i + V ext ( x i , y i , z i ) (cid:21) + N X i 35 ˚A in the z direction andstacked in the A-B-A-B way typical of this compound, asin Ref. 3. We stopped at eight graphene sheets since toinclude a ninth one changes the energy per particle lessthan the typical error bars for that magnitude (around0.1 K) [3]. V ext ( x i , y i , z i ) sums, for each He atom i , allthe C-He pair interactions calculated using the accurateCarlos and Cole anisotropic potential [19]. The He-Heinteraction V He − He ( r ij ) is modeled using a standard Azizpotential [20]. The graphite substrate was consideredto be rigid, but the helium atoms in both the first andsecond layers were allowed to move from their crystallo-graphic positions.In order to reduce the variance and to fix the phaseunder study, DMC incorporates importance sampling byusing a guiding wave function. This wave function isdesigned as a first, simple approximation to the many-body system and is variationally optimized. In our case,we usedΦ( r , . . . , r N ) = Φ J ( r , . . . , r N )Φ ( r , . . . , r N ) × Φ ( r N +1 , . . . , r N ) , (2)with Φ J ( r , . . . , r N ) = N Y i Top panel: Full squares: Energy per He atom ina single layer of a triangular solid as a function of surfaceper helium atom. Open squares: Same for an arrangementcomprising two helium layers, the density of the first one being0.115 ˚A − . A least-squares fitting polynomial to the firstlayer data and the Maxwell construction line between the twoarrangements are also shown. Bottom panel: Data for theenergy of the second layer system after subtracting the valuesgiven by the Maxwell construction line. second layer and for the number of lattice points of thesolids. Therefore, no vacancies were considered in anysolid. a was chosen to minimize the total energy of thesystem. Notice that Eq. (5) can describe a translation-ally invariant Helium second layer by fixing a to zero.In most cases, N was fixed to 224 atoms, distributed ina simulation box including 14 × N was fixed to produce the desired density. Thismeant simulation cells including up to 356 He atoms.Calculating the energy per particle as a function ofdensity, one can establish the stability range of the dif-ferent phases. The DMC energies are shown in Figs. 1and 2 for different values of the surface area (the inverseof the surface density). The first issue to be addressedis the determination the first-layer solid density whichproduces the lowest total energy for a two-layer arrange-ment. To this end, we followed a procedure used previ-ously for graphene [11]. As in that work, the optimumdensity of the first layer turned out to be 0.115 ˚A − . Astandard Maxwell construction between a system with asingle Helium monolayer and a double-layered one withthat solid density (Fig. 1, top panel) and a translationalinvariant second layer, gives us a promotion density of0.113 ± − (corresponding to 8.86 ± ).This is similar to the experimental value 0.114 ˚A − re-ported in Ref. 10 and somewhat smaller than the otherexperimental measure, 0.118 ˚A − , that of Ref. 9. Tocalculate the lowest density limit for the bilayer arrange-ment, we subtracted the energy values for the Maxwellconstruction line from the direct simulation results. Theresults obtained are shown in the bottom panel of Fig.1. As one can see, this limiting value is 6.13 ± ,which corresponds to a density σ = 0 . ± . 002 ˚A − .A similar procedure allows us to establish the stabilitylimits for larger values of the total He density. The re-sults are depicted in Fig. 2. There, we can see the energyper atom for different phases. For total density valuessmaller than 0.185 ˚A − , the underlying Helium densitywas, as before, 0.115 ˚A − , while for that value up, theenergy was lowered by slightly compressing the first layerto a density 0.1175 ˚A − . In Fig. 2, the open circles corre-spond to the same translationally invariant double layerarrangement already discussed (Fig. 1). The full squarescorrespond to a bilayer comprising two incommensuratetriangular solid layers. The 4/7 and 7/12 commensuratephases with the triangular lattice of the first layer werealso considered by an appropriate choice of the set ofcrystallographic positions defining them. The fact that,in the figure, there are two points for each of those ar-rangements is due to the fact that we considered the twopossibilities for the first-layer densities discussed above.In any case, it is clear from Fig. 2, that the 4/7 arrange-ment is unstable. On the other hand, the 7/12 phaseis right on top of the Maxwell construction line betweenthe translationally-invariant phase and the double trian-gular solid one. This means we would have a first or-der phase transition between a translationally-invariantphase of density σ = 0 . ± . 002 ˚A − , and a 7/12registered solid of 0.182 ˚A − . The first layer of this ar-rangement would be then compressed up to σ = 0 . − , that is again in equilibrium with a triangular solidof density 0.188 ˚A − . This is similar to what happens ina He double layer system on graphite [23], where boththe 4/7 and 7/12 phases are stable.To better characterize the different stable phases, weneed to go further than to establish their stability limits.As indicated above, one of the issues raised in Ref. 10was the possible existence of a quasi-two-dimensional su-persolid and of a normal (non superfluid) liquid phase atsmaller densities. In order to study those claims, we es-timated the superfluid fraction σ s /σ on the second layerof the arrangements found to be stable. This is done byusing the usual winding-number estimator in the limit ofzero temperature [22, 24], σ s σ = lim τ →∞ α (cid:18) D s ( τ ) τ (cid:19) , (6) −96−94−92−90−88−86−84−82−80 5 5.2 5.4 5.6 5.8 6 E n e r gy p e r H e a t o m ( K ) Surface (Å )solidliquid4/77/12 FIG. 2. Same than in the previous figure for the two layersystem at higher densities. Data for different phases are dis-played. with τ the imaginary time in the Monte Carlo simu-lation. Here, α = N / (4 D ), D = ¯ h / (2 m ), and D s ( τ ) = h [ R CM ( τ ) − R CM (0)] i . R CM is the positionof the center of mass of the N He atoms on the secondlayer, taking into account only their x and y coordinates(where periodic boundary conditions apply). The resultsobtained for different cases are shown in Fig. 3. The errorbars in that figure correspond to one standard deviationcomputed from several independent Monte Carlo runs,and when not displayed, are similar to those shown.In Fig. 3, we see that the translationally invariantphase akin to a liquid found at low densities is a perfectsuperfluid, since the values for the superfluid estimatorare on top to the line corresponding to σ s /σ = 1. Onthe other hand, the superfluid estimator for a triangu-lar solid of σ = 0 . 196 ˚A − is zero within our numericalresolution . This situation is common to all other trian-gular lattice solids, irrespective of their total density, andmakes them normal solids. However, the superfluid frac-tion corresponding to a 7/12 structure with σ = 0 . − (open squares) has an intermediate value betweenzero and one, corresponding to a superfluid fraction of0 . ± . 1. This fraction is the same as the one for a σ = 0 . 182 ˚A − solid, whose data are not shown for sim-plicity. This means there is a supersolid stable phase inthe density range 0.182-0.186˚A − . Our results also sup-port the existence of superfluidity in the range between0.170 and 0.182 ˚A − , in which there is a mixture of a fullsuperfluid liquid-like phase in coexistence with the 7/12registered solid with the lowest density. The same canbe said of the 0.186-0.188˚A − interval, where the coexis-tence is with a normal triangular 2D-solid. S up e rf l u i d fr ac ti on t (K −1 ) s = 0.170Å −2 s = 0.165Å −2 FIG. 3. Superfluid density for different second layer arrange-ments. Full squares, a translationally-invariant phase with σ = 0 . 170 ˚A − . Open circles, same but for σ = 0 . 165 ˚A − .Open squares, results for a 7/12 phase with σ = 0 . 186 ˚A − .Full circles, same for a triangular solid with ρ = 0 . 196 ˚A − .The dotted line corresponds to a perfect superfluid with nonormal component. Previous studies on similar systems have always as-sumed the translationally-invariant phase to be an ho-mogeneous liquid. In this work, we have checked thatassumption by considering the possibility of the secondlayer being some kind of hexatic phase, as suggested inRef. 9. To study if this was so, and following Ref. 25,we write the pair distribution function g ( r ) as g ( r ) = ∞ X n =0 g n ( r ) cos( nθ ) (7)with n even and θ = r · e . Here, e is a unit vector ina reference direction. The hexatic order in a non-solidphase is associated to a periodic oscillation in the g ( r )component with algebraic decay at large distances [13–15, 25]. This hexatic order parameter is shown in Fig. 4for two translational invariant structures within the sta-bility range of that phase. It is important to stress thatthe guiding wave function used in DMC to describe thosearrangements, a product of Eqs. (3), (4) and (5), this lastone with a = 0, does not include any explicit hexaticcorrelation. What we see is a regular pattern of max-ima and minima extending to long distances with a slow −0.004−0.003−0.002−0.001 0 0.001 0.002 0.003 0.004 0.005 0 2 4 6 8 10 12 14 16 H e x a ti c o r d e r p a r a m e t e r Distance (Å) FIG. 4. Hexatic order parameter for translationally invariantphases of σ = 0 . 170 ˚A − (full squares) and σ = 0 . 163 ˚A − (open squares). Full circles correspond to a phase with σ =0 . 165 ˚A − , but in which the effect of the first solid layer havebeen smoothed out. When not displayed, the error bars aresimilar to the ones shown. decay that we cannot determine completely due to thefinite size of our simulation box. This ordering is simi-lar to the one found to be metastable in strictly 2D Heat larger densities ( σ > . 060 ˚A − ) [25] that the second-layer ones for the systems under consideration (0.047 and0.055 ˚A − ). This fact suggests it to be a consequence ofthe corrugation of the first layer solid substrate. This un-derlying structure supports a series of potential minimabetween every three Helium atoms with the right sym-metry to produce the observed order. To check if this isso, we calculated the same parameter for a structure inwhich the potential created by the second layer had beenaveraged over to make those potential minima disappear.The result is that the set of maxima and minima are notlonger present, indicating that the hexatic correlationsare indeed corrugation-induced.Our theoretical results are at least qualitatively com-patible with the available experimental data. For in-stance, both Refs. 9 and 10 suggest that the phase dia-gram of the second layer of He on graphite starts rightafter promotion with a gas-liquid coexistence zone, fol-lowed upon an increase of Helium density to a stableliquid-like region. This would be, at least for the low-est densities ( σ < . 170 ˚A − ), a normal fluid that wouldundergo a first-order phase transition to a commensuratephase. This last phase would change to a high-density tri-angular solid. Ref. 10 assigns a density range of 0.1711-0.1809 ˚A − to the liquid-commensurate transition, and of0.1809-0.1841 ˚A − for the stable registered phase. Bothof them are comparable with our suggestions: 0.170-0.182˚A − and 0.182-0.186 ˚A − , while the data of Ref. 9 isshifted further up in the density scale. In that entirerange, we see a superfluid response, first in the coexis-tence between the superhexatic [16] phase and the 7/12registered solid, and then in the stability range of thatphase itself, in accordance with the experimental data ofboth Refs. 8 and 10. On the other hand, previous DMCcalculations on graphene found the 7/12 commensuratesolid to be unstable [11]. This difference in the behaviorof those close related systems can be due to the delicateenergy balance needed for that structure to be seen: theintroduction of the exchanges in the description of thesupersolid decreases the energy enough for this commen-surate phase to emerge at T = 0 K. The effect of theadditional carbon layers might also play a role, as in thecase of the first-layer √ × √ − instead of the normal fluid found experimentally, can bedue to the lack of connectivity of the real substrate [8].Finally, the lack of that response beyond 0.188 ˚A − , isin agreement with the results of Ref. 8.Our work is carried out strictly at zero temperatureand compares very well with recent experiments per-formed in the mK regime. It would be very interesting tostudy the same system at finite temperature to determinethe thermal stability of the superhexatic and supersolidphases, even though those phases could be unstable attemperature values too low to be accessible to the pathintegral Monte Carlo method [2].We acknowledge financial support from MINECO(Spain) Grants FIS2017-84114-C2-2-P and FIS2017-84114-C2-1-P. We also acknowledge the use of theC3UPO computer facilities at the Universidad Pablo deOlavide. [1] L.W. Bruch,M.W. Cole and E. Zaremba. Physical Ad-sorption, Forces and Phenomena Dover, New York(1997). [2] P. Corboz, M. Boninsegni, L. Pollet, and M. Troyer,Phys. Rev. B 78, 245414 (2008).[3] M.C. Gordillo and J. Boronat, Phys. Rev. Lett. ,085303 (2009).[4] A. Noury, J. Vergara-Cruz, P. Morfin, B. Placais, M. C.Gordillo, J. Boronat, S. Balibar, and A. Bachtold, Phys.Rev. Lett , 165301 (2019).[5] D.S. Greywall and P.A. Busch. Phys. Rev. Lett 309 (1993).[7] P. A. Crowell and J. D. Reppy, Phys. Rev. Lett. , 3291(1993).[8] P. A. Crowell and J. D. Reppy, Phys. Rev. B , 2701(1996).[9] S. Nakamura, K. Matsui, T. Matsui and H. Fukuyama.Phys. Rev. B , 180501(R) (2016)[10] J. Nyki, A. Phillis, A. Ho, D. Lee, P. Coleman, J. Parpia,B. Cowan and J. Saunders. Nat. Phys. 455 (2017).[11] M.C. Gordillo and J. Boronat, Phys. Rev. B, ,155301 (2012).[13] B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. ,121 (1978).[14] D. R. Nelson and B. I. Halperin, Phys. Rev B , 2457(1979).[15] A. P. Young, Phys. Rev. B , 1855 (1979).[16] K. Mullen, H. T. C. Stoof, M. Wallin, and S. M. Girvin,Phys. Rev. Lett. , 4013 (1994).[17] J. Boronat, in Microscopic Approaches to Quantum Liq-uids in Confined Geometries , Vol. , ed. E. Krotscheckand J. Navarro (World Scientific, Singapore, 2002).[18] M.C. Gordillo and J. Boronat. Phys, Rev. Lett. , 339 (1980).[20] R.A. Aziz, F.R.W McCourt and C.C. K. Wong. Mol.Phys. , 8920(1994).[22] M.C. Gordillo, C. Cazorla, and J. Boronat. Phys. Rev. B84