New insights into the 1D carbon chain through the RPA

We investigated the electronic and structural properties of the infinite linear carbon chain (carbyne) using density functional theory (DFT) and the random phase approximation (RPA) to the correlation energy. The studies are performed in vacuo and for carbyne inside a carbon nano tube (CNT). In the vacuum, semi-local DFT and RPA predict bond length alternations of about 0.04 \r{A} and 0.13 \r{A}, respectively. The frequency of the highest optical mode at the $\Gamma$ point is 1219 cm$^{-1}$ and about 2000 cm$^{-1}$ for DFT and the RPA. Agreement of the RPA to previous high level quantum chemistry and diffusion Monte-Carlo results is excellent. For the RPA we calculate the phonon-dispersion in the full Brillouine zone and find marked quantitative differences to DFT calculations not only at the $\Gamma$ point but also throughout the entire Brillouine zone. To model carbyne inside a carbon nanotube we considered a (10,0) CNT. Here the DFT calculations are even qualitatively sensitive to the k-points sampling. At the limes of a very dense k-points sampling, semi-local DFT predicts no bond length alternation (BLA), whereas in the RPA a sizeable BLA of 0.09 \r{A} prevails. The reduced BLA leads to a significant red shift of the vibrational frequencies down to about about 1600 cm$^{-1}$, which is in good agreement with experimental estimates. Overall, the good agreement between the RPA and previously reported results from correlated wavefunction methods and experimental Raman data suggests that the RPA provides remarkably accurate results at moderate computational costs. It hence presents a useful addition to the repertoire of correlated wavefunction methods and its accuracy clearly prevails for low dimensional systems, where semi-local density functionals struggle to yield even qualitatively correct results.


I. INTRODUCTION
Carbyne -the infinite sp 1 hybridized carbon chain -is a model 1D material with fascinating properties that point to promising applications. Theoretical studies predict that carbyne's strength, elastic modulus and stiffness are greater than for other carbon allotropes like diamond, carbon nanotubes and graphene and, in fact, greater than for any known material. New composite materials may be able to exploit these properties in order to achieve superior functionality. 1 Another exciting application might be realized by using single chains of linear carbon as a channel for field-effect transistors, thus bringing the channel width to the extreme minimum of one carbon atom thickness. 2 Unfortunately, carbyne is extremely unstable and very reactive at ambient conditions, making its synthesis and detection experimentally very difficult. While for many years, the use of heavy end-capping groups for stabilization was very common to synthesize linear carbon chains, the lengths that could be achieved were very limited, e.g. 44 atoms in a 2010 study by Chalifoux and Tykwinski. 3 Another approach is to use carbon nanotubes (CNTs) as confining nanoreactors to synthesize 1D materials. In a 2016 study, Shi L. et al. reported a method for the bulk production of long acetylenic linear carbon chains protected by thin doublewalled CNTs. They showed that stable chains composed of more than 6,000 carbon atoms can be produced. 4 This important step towards the bulk production of carbyne is not only encouraging from an experimental perspective, but also shifts the focus of theoretical investigations towards linear carbon a) Electronic mail: georg.kresse@univie.ac.at chains confined by CNTs.
Theoretical studies on the ideal carbon chain in vacuum show that Peierls distortions 5 give rise to bond length alternations (BLA) that are associated with single and triple C-C bonds. The underlying reason is that the carbon chain is characterized by a half-filled band at the Fermi level, and a small distortion from the equidistant chain introduces a gap near the Fermi level, thus lowering the total energy of the electrons. Although conventional local and semi-local density functionals reproduce this effect qualitatively, there is a considerable quantitative discrepancy for the predicted BLAs between correlated wave function methods like coupled cluster singles and doubles with perturbative triples [CCSD(T)] or diffusion Monte-Carlo (DMC) and (semi) local functionals, such as LDA or PBE. Hybrid functionals, for instance HSE06, provide an improvement over the (semi) local functionals, but the carbyne in vacuo BLA from DMC is still around 50% larger than for HSE06. 6 Likewise, the phonon frequencies are not accurately predicted by (semi) local functionals. 6,7 These findings already indicate that, although carbyne is a prototypical seemingly simple 1D model system, one needs to go beyond standard semi-local DFT to quantitatively describe the phenomena exhibited by it.
For the theoretical study of carbyne inside CNTs, which is of particular interest due to the before mentioned pathway for the production of bulk carbyne, the situation is even more involved, since van-der-Waals (vdW) interactions between the CNTs and the linear carbon chain come into play. Wanko et al. reported that two distinct and separable effects lead to a red shift of 118-290 cm −1 in the Raman spectra of carbon chains inside CNTs compared to in vacuo frequencies. On the one hand, a charge transfer (CT) from the tube to the chain softens the chains Γ mode. On the other hand, vdW interactions lead to a similar softening that potentially even exceeds the one from the CT. 7 Since vdW interactions are not accounted for by (semi) local density functionals, the carbyne-CNT system clearly cannot be described accurately by semi-local functionals.
A promising method to tackle this materials problem is provided by the RPA. It has been shown that the RPA performs very well for different bonding types, including vdW interactions. [8][9][10] The availability of low-scaling O(N 3 ) RPA algorithms 11,12 as well as the possibility to calculate forces in the RPA 13 offers a potential route to accurately describe the carbyne-CNT system at moderate computational cost. Additionally, versatile tools for the RPA that would facilitate further investigations have been reported recently. On the one hand, RPA natural orbitals can be employed for the efficient calculation of accurate second-order Møller Plesset (MP2) or full configuration interaction (FCI) energies. 14 On the other hand, using the RPA to benchmark density functionals 15 might lead to the identification of a sufficiently accurate functional that can be used for large scale structure studies. Finally, the very recent introduction of a finite-temperature generalization of the low-scaling RPA algorithms 16 provides an elegant way to treat the partial occupancies that occur due to the charge transfer from the tube to the chain. These prospects lead us to believe that investigating the performance of the RPA for the carbyne-CNT system will yield a valuable contribution to further our knowledge on this intriguing material.

II. COMPUTATIONAL DETAILS
We calculated PBE 17 and RPA 11,12,16 total energies for different BLAs of carbyne in vacuo as well as inside a (10,0) (CNT). Throughout this work we used VASP 6.1 18,19 with its implementation of the projector augmented wave (PAW) formalism 20,21 , a 400 eV plane wave cut-off, and the potential C_s_GW. This potential uses three projectors for the 2s and 2p partial waves and a truncated all-electron potential as local potential. The radial cutoffs are 1.35 a.u., 1.75 a.u and 0.7 a.u. for the s, p, and local potential, respectively.

A. Geometries and k-points
As the "standard" cell, we employed an 11.200 Å × 11.200 Å × 12.875 Å tetragonal cell with periodic boundary conditions in which the carbyne's and CNT's axes were aligned along the z-axis. The choice of this cell resulted from preliminary calculations: The geometry of the carbyne chain was determined by varying the lattice constant in the periodic direction between 2.5 and 2.7 Å, and relaxing the internal coordinate (BLA) for each lattice constant. For the equilibrium lattice constant, a value of 2.575 Å and 2.590 Å was determined with PBE and RPA, respectively.
The experimental bond length of carbon in graphene is about 1.42 Å. For the (10,0) zig-zag tube, this yields a lattice constant of 4.26 Å. If the tube is three times replicated this yields a lattice constant of 12.780 Å. Thus, a 5-fold super cell of carbyne fits almost perfectly into the three times replicated tube. Since we wanted to keep the strain on the carbyne chain and the tube at a low level, we decided to fix the lattice constant to 12.875 Å. This corresponds to the five-fold supercell of carbyne in PBE, and hence the carbyne is strain free at the level of PBE. This lattice constant is also very close to the arithmetic mean of the experimental estimates for the three times replicated (10,0) tube and the five-fold replicated primitive carbyne chain in the RPA (12.95 Å). All the other coordinates were relaxed at the PBE level until the forces on all atoms were smaller than 10 −5 eV/Å.
In total there are 130 C-atoms in the filled standard cell-120 CNT C-atoms around 10 carbyne C-atoms, which are placed on the CNT's principal axis (Fig. 1). For the carbyne in vacuo calculations, we used a primitive 11.2 Å × 11.2 Å × 2.575 Å tetragonal cell with two atoms at (0, 0, BLA/2) Å and (0, 0, 1.2875) Å. With this setup, the number of k-points in z-direction needs to be five times as great as in the standard cell to allow a one-to-one comparison.
In the following, when we specify the k-point sampling, we will always refer to the number of k-points per primitive cell kpc. For example, 2 k-points in the standard cell correspond to 10 kpc. We used different Γ-centered 1 × 1 × k z k-meshes with k z ranging from 10 to 640 kpc for PBE and from 10 to 80 (40) kpc for the RPA in vacuo (in the CNT).

B. RPA
In general, we performed the RPA calculations on top of PBE, i.e. we use occupied and virtual PBE orbitals to calculate the correlation energy perturbatively in the RPA 23 . For the RPA calculations of carbyne in vacuo the standard low scaling O(N 3 ) implementation in VASP 11,12 was used with 12 imaginary time/frequency grid points. However, for the filled standard cell, we used a FTRPA approach 16 that is particu-larly useful for systems with partial occupancies (in the oneelectron orbitals of the underlying mean-field ground state). In this method, the ground state is calculated with Fermi smearing, i.e. the occupancies follow a Fermi distribution that is controlled via the smearing parameter σ = k B T [eV ]. The RPA correlation is then calculated on an optimized imaginary time/frequency grid for the given σ . In this work, we calculated and compared BLA vs. FTRPA energy curves of carbyne inside a (10,0) CNT for σ [eV ] ∈ {0.03, 0.05, 0.10}, using 16 time/frequency grid points.

C. Phonon frequencies
The full phonon dispersion was calculated for carbyne in vacuo using DFT and the RPA, finite displacements and deriving from the forces the inter-atomic force constants. 13 For PBE, the primitive cell was replicated 20 times, whereas for RPA a more modest replication of 5 unit cells was chosen. To our knowledge, this is the first time that the phonon dispersion relation of carbyne was calculated using a correlated wavefunction method.
For carbyne inside the nanotube, calculating the phonon dispersion would be exceedingly demanding, since the translational symmetry is broken. Therefore, we only determined the highest frequency mode using a finite displacement corresponding to the longitudinal optical (LO) mode. To estimate ω Γ , we used the harmonic approximation for the energy E(b) as a function of the BLA b with the reduced mass of two carbon atoms µ C =6 u, and the BLA (energy) at equilibrium being b 0 (E 0 ). We used a polynomial fit for the total energy f E (b) and calculated the equilibrium BLA as and the phonon frequency ω Γ as For all PBE energy curves and for the in vacuo RPA energy curves, sixth-order polynomials were used, since the data were very smooth. For the less smooth FTRPA energy curves, we had to resort to third-order polynomial fits.

A. Preliminary tests
As mentioned above, the choice of our filled standard cell resulted from the ratio of the primitive vectors in z-direction, a z , for the relaxed chain in vacuum and the empty (10,0) CNT. Since the ratio is very close to 3:5, our choice of the filled standard cell provided an ideal starting point for the investigation of the carbyne-CNT system under periodic boundary conditions. To start with, we confirmed, on the PBE level, that our choice of the standard cell minimizes the total energy of the filled standard cell with respect to the cell parameters. Additionally we studied how the total energy of the carbyne-CNT system changes, when the carbon chain as a whole is shifted along the principal axis. The differences were below 1 meV and are thus small but not entirely negligible. Since we considered the geometry of the tube to be fixed in its equilibrium state, the chain's BLA is the only structure parameter left that significantly affects the energy of the system. Furthermore, we checked whether the primitive cell and the standard cell yield numerically identical results. In vacuo, the energy differences between a 5-fold super cell with N k-points and a primitive cell with 5N k-points and the same BLA and x-yvacuum distance was smaller than 0.2 meV. Also phonon frequencies at the Γ point are identical. Therefore, as expected, it is safe to compare results from the primitive cell of carbyne in vacuo with those of carbyne in the (10,0) CNT, as long as one uses the same kpc.

B. BLA and phonons in vacuo
VASP uses an internal cut-off extrapolation to calculate the RPA energy. Here, we compare the BLAs and LO Γ frequencies obtained by using the RPA energies at the highest cut-off for the response function (here 266.7 eV) with those obtained using the cut-off extrapolated energies (see Fig. 2). We found that the cut-off extrapolation slightly increases the BLA and the energy gain related to the BLA. This also brings the RPA results in better agreement with results from other correlated wavefunction methods such as DMC 6 . The difference between the highest cut-off and cut-off extrapolation in the predicted BLA and the LO Γ frequency was, however, less than 5% and 1%, respectively, for carbyne in vacuo.
In Fig. 3 on the left hand side, we show the BLA as well as the phonon frequencies for the LO Γ frequency in vacuo. The results are very sensitive to the density of the k-point mesh, especially for the DFT (PBE) calculations. For example, for carbyne in vacuo, 32 and 64 kpc yield BLAs 27% respectively 5% larger than at 640 kpc, see Fig. 3 (a). Likewise, the LO frequency at the Γ point is strongly k-point dependent, see Fig.  3 (c). Typically about 120 kpc are required to obtain highly accurate results. Remarkably the RPA results are somewhat less sensitive to the k-point sampling. The differences between 80 kpc and 40 kpc are quite small, and although we have attempted to use 80 kpc wherever possible, for the carbyne in the nanotube we had to reduce the sampling to 40 kpc (see below).
In vacuo using PBE, our predicted BLA is just 0.037 Å (with 640 kpc), while the RPA yields 0.129 Å (with 80 kpc), see Fig. 3 (a). The RPA result is remarkably close to the 0.136(2) Å from DMC 6 . Considering that the k-point sampling was rather coarse in the DMC (equivalent to 16 kpc), the agreement is very satisfactory.
In Fig. 4 we show the full phonon dispersion relation for carbyne in vacuo, for the RPA and PBE. To be consistent with our other calculations, we have chosen a five-fold super cell with a lattice constant of 5× 2.575 Å for the RPA calculations, whereas the PBE calculations were performed for a 20 times repeated unit cell. To confirm the results, in particular at the Brillouin zone boundary, we performed an additional RPA calculation (although with less k-points) in a four-fold super cell at the RPA lattice constant 4× 2.590 Å. It appears that the slight variation in the lattice constant has only a very small effect on the phonon dispersion. Likewise this test calculation demonstrates that the RPA force constants decay sufficiently fast that a five times replicated unit cell yields reliable results (a five times replicated unit cell does not guarantee technically accurate vibrational frequencies at the Z point, since the corresponding phonon pattern is not commensurate with the supercell).
Both, Fig. 3 and Fig. 4, show that the LO Γ frequency of carbyne in vacuo is very small at the level of PBE (1219 cm −1 ). Compared to the DMC result of 2084(5) cm −1 the error is almost 40 %. As was the case for the BLA, the RPA . Again, it needs to be emphasized that the DMC results were obtained using supercells corresponding to just 16 kpc, and the wave function based results are extrapolated from finite molecules with 32 atoms. At such low k-point densities, we also obtain somewhat higher phonon frequencies using the RPA.
C. BLA and phonons in (10,0) nanotube Our PBE calculations for the filled standard cell confirm previous results 24,25 . The BLA vanishes for carbon chains inside single walled CNTs (SWCNT) when using (semi) local density functionals. However, we want to emphasize that this is again only true for sufficiently dense k-point grids (here ≥ 80 kpc, compare Fig. 3). Likewise, the LO frequency at Γ is strongly k-point dependent, see Fig. 3 (d).
The suppression of the BLA in PBE is fairly easy to understand by inspecting the band structures in Fig. 5. For the non-dimerized carbyne chain in the nanotube, the Fermi level of the combined system is well above the Fermi-level of the carbon chain in the vacuum (cross in Fig. 5 left panel). This causes a significant charge transfer from the nanotube to the carbyne chain. Since the energy gain by the Peierls distortion can only be realized, if the upwards moving band is unoccu- pied and the downwards moving band is occupied, a charge transfer to the chain effectively suppresses the Peierls distortion. Recall that the Peierls distortion was small to begin with in DFT. In contrast to the PBE calculations, we found that the carbyne's BLA does not vanish inside the tube using the RPA. This is a remarkable qualitative difference. Concerning the technical accuracy of this result, we have a few comments to make. As can be seen in Fig. 3 (a), the equilibrium BLA for the primitive cell at 40 kpc is only slightly lower (2%) than using 80 kpc in vacuum. We are therefore quite confident that the reported BLAs obtained in the (10,0) CNT are reasonably accurate for 40 kpc. Depending on σ , the BLA of carbyne inside a (10,0) CNT calculated using the FTRPA is 0.083, 0.089 or 0.101 [Å] for σ =0.10, 0.05 or 0.03 [eV]. Although, the data become smoother with increasing σ , these broadenings correspond to temperatures of about 1200 K, 600 K and 400 K. Hence a small width is certainly preferable in order to compare with experimental results measured around room temperature. However, the BLAs and LO Γ frequencies do not change qualitatively.
The band structure shown in Fig. 5 evaluated at the BLA of 0.09 Å-close to the RPA BLA -shows that the valence band edge of the nanotube and the conduction band edge of the carbyne are almost iso-energetic. This means that the charge transfer is small at the equilibrium RPA BLA.
However, this onset of charge transfer also makes it some-what difficult to make accurate predictions for the carbyne's LO Γ frequency. As the BLA becomes smaller, charge is progressively transferred from the tube to the carbyne chain. If the k-point density is not very large, the one electron occupancies can change somewhat discontinuously. Thus, the RPA calculations would certainly have benefited from denser k-point grids, however, more than 40 kpc for the FTRPA calculations of the filled standard cell would have been difficult with the available computational resources due to excessive memory requirements. Our best estimate of the tube induced redshift of the carbyne's LO Γ frequency is 345 cm −1 at 40 kpc. This value is quite close to the 290 cm −1 reported by Wanko et al. 7 , which was obtained by calculating the difference between the theoretical carbyne in vacuo frequency from DMC 6 and the experimentally observed frequencies of long linear carbon chains inside double-wall CNTs (DWCNTs) 4 . Since the PBE functional does not yield any BLA inside the tube, the energy surface becomes very soft, and the phonon frequencies are even further red-shifted to about 1100 cm −1 using 640 kpc. Wanko et al. 7 related this erroneous description to the lack of vdW interactions in the PBE functional, but we do not agree with this interpretation. Semi-local functionals are not accurate for 1D systems, in particular, in the prediction of Peierls instabilities. Because of self-interaction errors, they generally underestimate the magnitude of the BLA, both in vacuo as well as in the tube. In vacuum, the disagreement is only quantitative, but since the charge transfer from the tube to the chain decreases the tendency for a bond length alteration (see Fig. 5), PBE ultimately fails to yield any BLA inside the tube, whereas in the RPA the BLA is only decreased but does not vanish entirely.
A final comment on the technical challenges of the present calculations seems in place. We performed our calculations on the Vienna Scientific Cluster (VSC4) using Intel Skylake Platinum 8174 cores. A 80 kpc carbyne in vacuo RPA energy calculation took typically 45 min with 144 cores in parallel and used up to about 5.1 GB RAM per processor. For the filled standard cell and 20 kpc we used 192 cores in parallel that took typically 1 h 15 min for one FTRPA energy calculation and used up 4.6 GB RAM per processor. With 40 kpc the memory demand increased to about 8 GB RAM, and 80 kpc calculations were not possible using 192 cores. The RPA phonon calculation for the five-fold primitive cell of carbyne in vacuo with 40 kpc were most expensive and took 36 h on 384 cores and used 6.9 GB RAM per processor. All in all, these are very reasonable computational requirements, except for the fairly excessive memory demands of the low scaling RPA code.

IV. CONCLUSION
We have investigated the performance of the (FT)RPA for carbyne in vacuo and inside a SWCNT and found good agreement with experiments and other correlated wave function methods at moderate computational cost. Crucially, we found that a dense k-point sampling is absolutely required to obtain reliable and accurate results. Specifically, the PBE functional requires more than 128 k-points in the primitive cell to yield converged results. For example, for the carbyne chain inside the SWCNT, more than 80 kpc are necessary in order to observe that the BLA vanishes for PBE. Somewhat fortunately, the RPA does not show such a strong k-point dependence and reasonable results can already be expected around 40 kpc. This allowed us to calculate the full phonon dispersion relation of carbyne in vacuo. Remarkably, the results are quantitatively quite different from the PBE results, not only for the LO mode close to Γ point but throughout the entire Brillouine zone. The important point is that previous calculations for diamond and graphene showed only very small differences between semi-local functionals and the RPA 13 . The increased discrepancy in 1D is a strong indication that for low dimensional systems, beyond mean field methods are required to predict highly accurate results.
Our most interesting result is however the predicted BLA for carbyne inside a carbon nanotube. Using the RPA, our best estimate for the BLA is 0.09 Å in a (10,0) CNT. Previous assessments were mostly indirect and relied on establishing correlations between empty and filled tubes and the predicted and observed LO Γ frequency. The BLA predicted by the RPA seems to be in good agreement with the previous deductions, specifically the RPA predicts a LO Γ frequency of 1600 cm −1 in very reasonable agreement with the available experimental data (1800 cm −1 ). As to why our LO Γ frequency is too small, we can only speculate. First, the RPA might somewhat underestimate the phonon frequencies of the free carbyne chain (RPA 2000 cm −1 other correlated wave functions methods around 2080 cm −1 , albeit using only few k-points or rather short finite carbyne fragments). The origin for this small difference could be related to the approximations that the RPA makes, for instance, the neglect of second-order exchange or errors in the PBE density matrix. An alternative explanation is that the energy surface is very anharmonic at small BLAs, so that one might need to account for anharmonic effects in the theoretical calculations or beyond adiabatic coupling of phonons to electrons at the Fermi-level. Finally, the amount of charge transferred to the carbyne chain determines how strongly the BLA is suppressed and how much the LO phonon frequency changes compared to the free standing chain. The charge transfer might be different in double wall tubes, as the one used in experiments, than in our simplified theoretical model.
On the upside, the present theoretical study clearly demonstrates the capabilities of the RPA and the recently presented FTRPA method 16 . Calculations for carbyne inside a nanotube are simply not possible with any other correlated wave function method or diffusion Monte-Carlo. Using the RPA forces 13 , we were also able to predict the phonon dispersion relation of free standing carbyne in vacuo from beyond mean field methods, again, something hardly possible for other correlated wave function methods. Due to the lower symmetry inside the tube, calculations of the phonon dispersion relation for carbyne inside a SWCNTs are not yet possible, but might become possible in the not too distant future. We hope that our current results inspire more studies of 1D systems using the RPA and that the RPA can help to unravel the fascinating behavior of 1D materials.

CONFLICTS OF INTEREST
There are no conflicts to declare.

ACKNOWLEDGMENTS
Computations were performed on the Vienna Scientific Cluster VSC4.