# Interlayer correlation between two He monolayers adsorbed on both sides of -graphyne

###### Abstract

Path-integral Monte Carlo calculations have been performed to study the He adsorption on both sides of a single -graphyne sheet. For investigation of the interlayer correlation between the upper and the lower monolayer of He adatoms, the He-substrate interaction is described by the sum of the He-C interatomic pair potentials, for which we use both Lennard-Jones and Yukawa-6 anisotropic potentials. When the lower He layer is a C commensurate solid, the upper-layer He atoms are found to form a Kagomé lattice structure at a Mott insulating density of 0.0706 Å, and a commensurate solid at an areal density of 0.0941 Å for both substrate potentials. The correlation between upper- and lower-layer pseudospins, which were introduced in Ref. Kwon et al. (2013) for two degenerate configurations of three He atoms in a hexagonal cell, depends on the substrate potential used; With the substrate potential based on the anisotropic Yukawa-6 pair potentials, the Ising pseudo-spins of both He layers are found to be anti-parallel to each other while the parallel and anti-parallel pseudo-spin alignments between the two He layers are nearly degenerate with the Lennard-Jones potentials. This is attributed to the difference in the interlayer distance, which is Å with the Yukawa-6 substrate potential but as large as Å with the Lennard-Jones potential.

###### pacs:

67.25.bd, 67.25.bh, 67.80.B-## I Introduction

Among many substrates, graphite has long served as a test bed to investigate low-dimensional quantum fluids because of its strong binding of adsorbates. Up to seven distinct He layers were observed on graphite and each helium layer is considered to be a quasi-two-dimensional quantum system Zimmerli et al. (1992). The first He adlayer on graphite shows a commensurate-incommensurate solid transition as the helium coverage increases Greywall and Busch (1991); Greywall (1993); Crowell and Reppy (1996). Recently, a series of theoretical calculations have been performed to study the He adsorption on newly-synthesized (or -proposed) low-dimensional carbon substrates such as graphene Gordillo and Boronat (2009); Kwon and Ceperley (2012); Happacher et al. (2013), graphynes Kwon et al. (2013); Ahn et al. (2014), carbon nanotubes Cole et al. (2000); Gordillo (2008), and fullerene molecules Kwon and Shin (2010); Shin and Kwon (2012); Kim and Kwon (2013); Park and Kwon (2014). The phase diagrams of the He layers adsorbed on graphene were predicted to be very similar to those of the corresponding layers on graphite; the monolayer of He adatoms shows a C commensurate structure at the areal density of 0.0636 Å and goes through various domain-wall phases before crystallizing into an incommensurate triangular solid near its completion Kwon and Ceperley (2012); Happacher et al. (2013).

Graphyne is a two-dimensional (2D) network of - and -bonded C atoms Baughman et al. (1987); Coluci et al. (2004) which could be permeable to a He gas unlike graphene. Despite much experimental effort motivated by some promising theoretical predictions for graphyne as new Dirac materials Malko et al. (2012); Kim and Choi (2012); Chen et al. (2013) and high-capacity energy storage materials Zhang et al. (2011); Hwang et al. (2012, 2013), there has been no successful report yet for fabrication of extended 2D graphynes. However, some flakes or building blocks of finite-size graphynes have been synthesized Tobe et al. (2003); Haley (2008); Diederich and Kivala (2010), leading to a belief that graphynes will be fabricated in the near future. On the surface of -graphyne, which is the most stable graphyne structure according to quantum Monte Carlo calculations Shin et al. (2014), the He monolayer was predicted to exhibit a richer phase diagram than the corresponding layer on graphene or graphite, including various commensurate and incommensurate structures depending on the helium density Ahn et al. (2014). Recently one of us performed path-integral Monte Carlo (PIMC) calculations for He atoms adsorbed on a AB-stacked bilayer -graphyne Kwon et al. (2013), which is a hybridized honeycomb structure with each hexagon side consisting of one and two C atoms. It was found that the He monolayer was in a Mott insulating state at an areal density of 0.0706 Å while a commensurate solid was realized at 0.0941 Å. Introducing Ising pseudo-spin degrees of freedom for two degenerate configurations for three He atoms occupying a hexagonal cell (see Fig. 3 of Ref. Kwon et al. (2013)), this Mott-insulator to commensurate-solid transition was interpreted as a symmetry breaking process from a spin liquid of geometrically-frustrated antiferromagnets to a spin-aligned ferromagnet Kwon et al. (2013).

One interesting feature of a 2D carbon structure, such as graphene and graphyne, is that it can be suspended in the air Meyer et al. (2007) and He atoms can be coated on both sides. Noting that some new physics could emerge as a result of interlayer correlation between opposite-side He layers, some theoretical studies were recently done for the He adsorption on both sides of a single graphene sheet. Markić et al. found that the correlation between two He clusters adsorbed on opposite sides of graphene, Å apart from each other, was quite weak as evidenced by peakless pair distribution functions Markić et al. (2013). A weak correlation between two He systems on the opposite sides of graphene was also predicted by Gordillo’s diffusion Monte Carlo calculations, which showed that the phase diagram of the He monolayer on graphene would not be affected by the He adsorption on the other side Gordillo (2014). In this paper we report PIMC study of the He adsorption on both sides of a single -graphyne sheet. Because -graphyne is more porous than graphene, He atoms can penetrate through graphyne to allow physical exchanges among He atoms on opposite sides. This could result in stronger interlayer correlation than the corresponding systems on graphene. We find that He atoms in a Mott-insulating state form a 2D Kagomé lattice as a result of the interlayer correlation when the opposite-side He layer is a C commensurate solid, a ferromagnetic state in a pseudospin terminology. Effects of the interlayer correlation between two ferromagnetic C solids are found to depend on the substrate potential used; the parallel and the antiparallel pseudospin alignments between two He layers are nearly degenerate with the substrate potential based on the Lennard-Jones (LJ) He-C pair potentials while the antiparallel alignment is favored with the one described by the Yukawa-6 pair potentials. The vacancy formation in a He layer on -graphyne is also found to be affected by the presence of the opposite-side He layer.

## Ii Methodology

In this study, a single -graphyne sheet is set to be at . The He-graphyne interaction is assumed to be a pairwise sum of interatomic potentials between the carbon atoms and a He atom, which has been widely used to describe the interaction between a He atom and a carbon substrate Kwon and Ceperley (2012); Markić et al. (2013); Kwon et al. (2013); Gordillo (2014); Ahn et al. (2014). For the He-C interatomic pair potential, we employ two anisotropic potentials proposed by Carlos and Cole Carlos and Cole (1979, 1980), i.e., a 6-12 LJ potential and a Yukawa-6 potential. For the computational convenience our previous study for the He monolayer on bilayer -graphyne was done with only isotropic parts of the LJ pair potential. However, the original interatomic pair potentials of Carlos and Cole include anisotropic parts to fit helium scattering data from graphite surfaces. Even though the inclusion of the anisotropic parts of the interatomic potentials has little effect on quantum phases displayed by the He layer on one side of -graphyne, it allows some He atoms to be closer to the substrate, resulting in stronger correlation between two He layers on the opposite sides (the minima of the substrate potential made of the anisotropic pair potentials are deeper and closer to graphyne than the corresponding ones based on only isotropic parts of the pair potentials). This leads to our decision of using the substrate potentials based on fully-anisotropic interatomic pair potentials, which should give a better description of the interlayer correlation. Furthermore, since the LJ and the Yukawa-6 potentials used in this study were based on an interaction between helium and -bonded carbon atoms in graphite, we tested the sensitivity of our modelling of He-graphyne potentials to the well depth of the pair potentials. Although decrease in the well depth yields more fluctuations in He density distributions, the density modulations are found to change only little and our main results presented below are still, at least qualitatively, valid. For the He-He interaction, we use a well-known Aziz potential Aziz et al. (1992).

In the discrete path-integral representation, the thermal density matrix at a low temperature is expressed by a convolution of high-temperature density matrices with an imaginary time step of Ceperley (1995). While the isotropic parts of He-C pair potentials along with the He-He potential pair potentials are used to compute the exact two-body density matrices Ceperley (1995); Zillich et al. (2005) at the high temperature , their anisotropic parts are treated with the primitive approximation Ceperley (1995). This is found to give accurate description of both He-He and He-graphyne interactions with a time step of K. We employ the multilevel Metropolis algorithm to sample the imaginary time paths along with permutations among He atoms as described in Ref. Ceperley (1995). To minimize finite size effects, periodic boundary conditions are applied along the lateral directions.

## Iii Results

The PIMC calculations were done with a fixed simulation cell with dimensions of Å, the same as in our previous study for He on bilayer -graphyne Kwon et al. (2013). We focus on the interlayer correlation between two He layers on the opposite sides of graphyne, which are either in a Mott-insulating state or a pseudospin-aligned commensurate solid state. The results obtained with two different substrate potentials, the LJ potential and the Yukawa-6 one, are presented separately below.

### iii.1 Lennard-Jones substrate potential

For PIMC calculations with the LJ substrate potential, we first prepare the -graphyne surface whose in-plane hexagon center is occupied by a single He atom and whose lower side is coated with a monolayer of He atoms constituting a C commensurate solid while each of the in-plane centers is occupied by a single He atom. The simulations for the He adsorptions on the upper side of the prepared graphyne surface begin from an initial configuration of He atoms being randomly distributed at the distances far away from graphyne. Figure 1 presents one-dimensional (1D) density distributions of He atoms, as a function of the vertical coordinate along the direction perpendicular to the graphyne surface, for two different combinations of particle numbers per simulation cell. Two distinct density peaks, which correspond to the upper and the lower He layers, are observed on the opposite sides of graphyne. Note that 36 and 48 He atoms per simulation cell correspond to the Mott-insulating density of 0.0706 Å and the C commensurate density of 0.0941 Å, respectively. The additional density peak at corresponds to the zeroth layer consisting of He atoms embedded onto the in-plane hexagon centers, which was also observed in our previous study for He adatoms on a bilayer -graphyne Kwon et al. (2013). Since the peak-to-peak distance between the upper and the lower layers is estimated to be about 4.8 Å, one can expect that the van der Waals interaction between He atoms on the opposite sides is weakly attractive (note that the Aziz potential we used for the He-He interaction has a minimum value at Å). In addition, the clear separation between the adjacent density peaks in Figure 1 suggests that exchange couplings among He atoms in different layers are nearly absent and any correlation between the upper and the lower layers, if it exists, should stem mostly from the weakly-attractive He-He interaction rather than particle exchanges.

Figure 2 shows two-dimensional (2D) density distributions of the upper-layer He atoms, while the lower-layer density peaks represented by the white stars constitute a C commensurate structure with all pseudospins being in the spin-up state (see Fig. 3 of Ref. Kwon et al. (2013)). Here a distinct density peak in each plot represents an occupancy of a single He atom. At an areal density of 0.0706 Å, every hexagonal cell of graphyne is seen in Fig. 2(a) to accommodate three upper-layer He atoms, which is a manifestation of a Mott-insulating state. Without the lower He layer in a pseudospin-aligned commensurate solid state, this Mott-insulating state is a nonmagnetic spin liquid of frustrated antiferromagnets in terms of pseudospin degrees of freedom (see Fig. 2(d) in Ref. Kwon et al. (2013)). However, in the presence of the pseudospin-aligned lower He layer, the upper-layer pseudospins are shown in Fig. 2(a) to be aligned in the same direction as the lower-layer ones. Our PIMC simulations at K have also produced the antiparallel pseudospin alignment between the two He layers. This is understood by the fact that the parallel alignment is energetically favored only by K per an upper-layer helium atom over the antiparallel alignment, i.e., two pseudospin alignments are nearly degenerate. We here note that the upper-layer He atoms in a pseudospin-aligned Mott insulating state of Fig. 2(a) constitute a 2D Kagomé lattice. This is also true when the upper-layer pseudospins are aligned in the opposite direction to the lower-layer ones. Therefore one can conclude that as a result of the interlayer correlation, the upper-layer He atoms form a Kagomé lattice structure at the Mott insulating density of Å when the lower He atoms constitute a pseudospin-aligned C commensurate solid.

The interlayer correlation between two ferromagnetic C commensurate solids is also analyzed. Figure 2(b) presents the 2D density distribution of the upper-layer He atoms at the areal density of Å, where they constitute a 4/3 commensurate solid. The upper-layer pseudospins are seen to be aligned in the same direction as the lower-layer ones. Similarly to the case of the Mott-insulating state, the antiparallel pseudospin alignment was also observed in our simulations. These two pseudospin alignments are more degenerate (the parallel alignment was found to be preferred by K per upper-layer He atom), than the parallel and antiparallel pseudospin alignments between a Mott-insulator and a C solid in Fig. 2(a). This can be understood by the fact that unlike He atoms inside a hexagonal cell, upper-layer He atoms at the vertices of the graphyne hexagons in a C solid state prefer the other sublattice sites over the ones occupied by the corresponding lower-layer helium atoms. We note that once a pseudospin alignment between the two 4/3 commensurate He solids is established, either parallel or antiparallel to each other, energy barrier is too large to reverse one alignment to the other.

Now we analyze the effects of the interlayer correlation on the formation of vacancies especially in a Mott insulator with the Kagomé lattice structure. Figure 3 shows the 2D density distribution of 35 upper-layer He atoms, one less than the Mott-insulating case, per simulation cell. So the upper-layer Mott insulator contains one vacancy per simulation cell while the lower He layer is the same 4/3 commensurate solid as in Fig. 2(a). One can see that every hexagonal cell, except one, is seen to accommodate three upper-layer He atoms and its pseudospin is aligned in the same direction as those of the ferromagnetic lower layer. As shown in Fig. 3, one cell involving only two upper-layer atoms does not show the clear pseudospin alignment. This tells us that the Kagomé lattice structure is sustained even with the creation of vacancies but those isolated vacancies are restricted at one triangle of this trihexagonal tiling structure without hopping to the neighboring sites because of high potential barrier provided by the graphyne surface.

### iii.2 Yukawa-6 substrate potential

Our PIMC simulations with the Yukawa-6 substrate potential start from an initial configuration of and He atoms being distributed randomly on the upper and the lower side of -graphyne, respectively. Figure 4 presents the 1D He density distributions as a function of the vertical coordinate for two different values of while is fixed to 48 per simulation cell. Unlike Fig. 1 for the LJ substrate potential, only two density peaks on the opposite sides of graphyne are observed without the zeroth layer consisting of He atoms embedded onto the in-plane hexagon centers. The absence of the zeroth He layer is attributed to the fact that the Yukawa-6 substrate potential is more slowly varying near the potential minima, that is, the hexagon centers, than the LJ substrate potential. Note that the Yukawa-6 He-C pair potential is less repulsive at short distances than the LJ interatomic pair potential Carlos and Cole (1980). The distance between two density peaks is Å, shorter than the corresponding distance for the LJ substrate potential. Furthermore, there is significant overlap between the two density peaks that are broader than those in Fig. 1. This indicates large quantum fluctuations of He adatoms along the vertical direction, which could result in frequent particle exchanges between these two layers.

Figure 5 shows 2D density distributions of the upper-layer He atoms, while the lower-layer density peaks are represented by the white stars. Even with the Yukawa-6 substrate potential, the lower He layer is seen to constitute a pseudospin-aligned a C commensurate solid at the areal density of 0.0941 Å. This provides another confirmation to the conclusion of Ref. Kwon et al. (2013) that most of quantum phases manifested in the He monolayer on the -graphyne surface, such as a Mott insulator, commensurate solids, and pseudospin degrees of freedom, are not sensitive to the specifics of the substrate potential but are determined mostly by the surface geometry. It is also shown in Fig. 5(a) that the upper-layer He atoms under the Yukawa-6 substrate potential are in a Mott-insulating state at the areal density of Å with each hexagonal cell accommodating three He atoms. Unlike Fig. 2(a), however, all upper-layer pseudospins are aligned in the opposite direction to those of the ferromagnetic lower-layer commensurate solid. We understand that the increase in the effective hard-core radii of He adatoms due to larger quantum fluctuations, along with a shorter interlayer distance, causes the antiparallel pseudospin alignment to be favored under the Yukawa-6 substrate potential. As observed with the LJ substrate potential, the upper-layer He adatoms in the pseudospin-aligned Mott insulating state constitute a 2D Kagomé lattice structure. The interlayer correlation that favors the antiparallel pseudospin alignment is more evident between two ferromagnetic C solids. We observe in Fig. 5(b) that the pseudospins of a upper-layer 4/3 commensurate solid are in spin-down state while the lower-layer pseudospins are in spin-up state. This antiparallel pseudospin alignment between the two He adlayers corresponds to the AB stacking of two triangular solids. So we conclude that with the Yukawa-6 substrate potential, the AB stacking is preferred to the AA stacking between two C triangular He solids while these two stacking orders are nearly degenerate with the LJ potential as discussed in Sec. III.1.

We now try to create a single vacancy in the upper-layer Mott insulator, i.e., the Kagomé lattice structure, by putting only 35 He atoms on the upper side in our simulation. Figure 6 presents the 2D density distribution of the upper-layer He atoms along with the lower-layer density peaks of a pseudospin-up C commensurate solid. It is shown that there are still three upper-layer density peaks per every hexagonal cell, corresponding to a configuration of a Mott insulator with all pseudospins being in the down state. However, one lower-layer density peak or one white star is missing at a vertex of the graphyne honeycomb structure (see the bottom right corner), which prevents the lower-layer He atoms from forming a perfect C triangular lattice. This suggests that when a single vacancy is created in a upper-layer Kagomé lattice, one lower-layer He atom moves to the upper layer to form a perfect upper-layer lattice structure while a localized vacancy is created in the lower-layer 4/3 commensurate solid. This lower-layer vacancy is found at a vertex site on top of a carbon atom because it is a less favorable site for a He adatom than a site inside a hexagonal cell. The layer-to-layer hopping of a vacancy can be understood by a short interlayer distance and large quantum fluctuations along the vertical direction under the Yukawa-6 substrate potential.

## Iv Conclusion

According to our PIMC calculations of using two different He-substrate potentials, He atoms form distinct layers on both sides of a single -graphyne sheet. Regardless of the substrate potential used, the upper-layer He atoms form a 2D Kagomé structure at the Mott-insulating density of 0.0706 Å as a result of the interlayer correlation when the lower layer is a pseudospin-aligned C commensurate solid. Since the interaction of He atoms with a substrate or between themselves is similar to the corresponding interaction for He, the same Kagomé lattice structure is expected to be formed in the fermionic counterpart of the upper He layer , i.e., a He upper layer adsorbed on -graphyne, when its lower side is coated with the C commensurate helium solid. We speculate that some novel phenomena related with a geometrically-frustrated antiferromagnetism such as quantum spin liquids Balents (2010), could emerge in this He Kagomé lattice.

The interlayer correlation results in different stacking orders between two C commensurate triangular solids on the opposite sides of graphyne, depending on the substrate potential; with the Yukawa-6 potential, the AB stacking (an antiparallel pseudospin alignment between two He solids) is found to be favored but both AA (a parallel pseudospin alignment) and AB stacking configurations are nearly degenerate with the LJ substrate potential. This is attributed to the difference between two substrate potentials in the interlayer distance as well as in the magnitude of quantum fluctuations along the vertical direction. A more accurate He-graphyne potential would be required to draw a definite conclusion about the preferred stacking order of two commensurate triangular He solids on the opposite sides of -graphyne.

Recent theoretical studies done by Markić et al. Markić et al. (2013) and Gordillo Gordillo (2014) reported that the interlayer correlation between two He systems adsorbed on both sides of graphene was very weak and quantum phase diagram of one He layer would not be affected by the presence of the opposite-side He layer. On the other hand, our PIMC calculations have revealed some significant effects of the interlayer correlation on the structural properties of the He monolayers on -graphyne. This difference is understood to be due to much more porous nature of -graphyne than graphene.

###### Acknowledgements.

This work was supported by Konkuk University. We also acknowledge the support from the Supercomputing Center/Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2013-C3-033).## References

- Kwon et al. (2013) Y. Kwon, H. Shin, and H. Lee, Phys. Rev. B 88, 201403(R) (2013).
- Zimmerli et al. (1992) G. Zimmerli, G. Mistura, and M. H. W. Chan, Phys. Rev. Lett. 68, 60 (1992).
- Greywall and Busch (1991) D. S. Greywall and P. A. Busch, Phys. Rev. Lett. 67, 3535 (1991).
- Greywall (1993) D. S. Greywall, Phys. Rev. B 47, 309 (1993).
- Crowell and Reppy (1996) P. A. Crowell and J. D. Reppy, Phys. Rev. B 53, 2701 (1996).
- Gordillo and Boronat (2009) M. C. Gordillo and J. Boronat, Phys. Rev. Lett. 102, 085303 (2009).
- Kwon and Ceperley (2012) Y. Kwon and D. M. Ceperley, Phys. Rev. B 85, 224501 (2012).
- Happacher et al. (2013) J. Happacher, P. Corboz, M. Boninsegni, and L. Pollet, Phys. Rev. B 87, 094514 (2013).
- Ahn et al. (2014) J. Ahn, H. Lee, and Y. Kwon, Phys. Rev. B 90, 075433 (2014).
- Cole et al. (2000) M. W. Cole, V. H. Crespi, G. Stan, C. Ebner, J. M. Hartman, S. Moroni, and M. Boninsegni, Phys. Rev. Lett. 84, 3883 (2000).
- Gordillo (2008) M. C. Gordillo, Phys. Rev. Lett. 101, 046102 (2008).
- Kwon and Shin (2010) Y. Kwon and H. Shin, Phys. Rev. B 82, 172506 (2010).
- Shin and Kwon (2012) H. Shin and Y. Kwon, J. Chem. Phys. 136, 064514 (2012).
- Kim and Kwon (2013) B. Kim and Y. Kwon, J. Low Temp. Phys. 171, 599 (2013).
- Park and Kwon (2014) S. Park and Y. Kwon, Phys. Rev. E. 89, 042118 (2014).
- Baughman et al. (1987) R. H. Baughman, H. Eckhardt, and M. Kertesz, J. Chem. Phys. 87, 6687 (1987).
- Coluci et al. (2004) V. R. Coluci, S. F. Braga, S. B. Legoas, D. S. Galvão, and R. H. Baughman, Nanotechnology 15, S142 (2004).
- Malko et al. (2012) D. Malko, C. Neiss, F. Viñes, and A. Görling, Phys. Rev. Lett. 108, 086804 (2012).
- Kim and Choi (2012) B. G. Kim and H. J. Choi, Phys. Rev. B 86, 115435 (2012).
- Chen et al. (2013) J. Chen, J. Xi, D. Wang, and Z. Shuai, J. Phys. Chem. Lett. 4, 1443 (2013).
- Zhang et al. (2011) H. Zhang, M. Zhao, X. He, Z. Wang, X. Zhang, and X. Liu, J. Phys. Chem. C 115, 8845 (2011).
- Hwang et al. (2012) H. Hwang, Y. Kwon, and H. Lee, J. Phys. Chem. C 116, 20220 (2012).
- Hwang et al. (2013) H. Hwang, J. Koo, M. Park, N. Park, Y. Kwon, and H. Lee, J. Phys. Chem. C 117, 6919 (2013).
- Tobe et al. (2003) Y. Tobe, I. Ohki, M. Sonoda, H. Niino, T. Sato, and T. Wakabayashi, J. Am. Chem. Soc. 125, 5614 (2003).
- Haley (2008) M. M. Haley, Pure Appl. Chem. 80, 519 (2008).
- Diederich and Kivala (2010) F. Diederich and M. Kivala, Adv. Mater. 22, 803 (2010).
- Shin et al. (2014) H. Shin, S. Kang, J. Koo, H. Lee, J. Kim, and Y. Kwon, J. Chem. Phys. 140, 114702 (2014).
- Meyer et al. (2007) J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth, and S. Roth, Nature 446, 60 (2007).
- Markić et al. (2013) L. Vranješ Markić, P. Stipanović, I. Bešlić, and R. E. Zillich, Phys. Rev. B 88, 125416 (2013).
- Gordillo (2014) M. C. Gordillo, Phys. Rev. B 89, 155401 (2014).
- Carlos and Cole (1979) W. E. Carlos and M. W. Cole, Phys. Rev. Lett. 43, 697 (1979).
- Carlos and Cole (1980) W. E. Carlos and M. W. Cole, Surf. Sci. 91, 339 (1980).
- Aziz et al. (1992) R. A. Aziz, M. J. Slaman, A. Koide, A. R. Allnatt, and W. J. Meath, Mol. Phys. 77, 321 (1992).
- Ceperley (1995) D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
- Zillich et al. (2005) R. E. Zillich, F. Paesani, Y. Kwon, and K. B. Whaley, J. Chem. Phys. 123, 114301 (2005).
- Balents (2010) L. Balents, Nature 464, 199 (2010).