Quantum Revivals in Curved Graphene Nanoflakes

  Citation: de-la-Huerta-Sainz, S.; Ballesteros, A.; Cordero, N.A. Quantum Revivals in Curved Graphene Nanoflakes. Nanomaterials 2022,12, 1953. https://doi.org/ 10.3390/nano12121953 Academic Editor: Arthur P Baddorf Received: 23 April 2022 Accepted: 3 June 2022 Published: 7 June 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). nanomaterials Article Quantum Revivals in Curved Graphene Nanoflakes Sergio de-la-Huerta-Sainz 1, Angel Ballesteros 1and Nicolás A. Cordero 1,2,3,* 1Physics Department, Universidad de Burgos, E-09001 Burgos, Spain; [email protected] (S.d.-l.-H.-S.); [email protected] (A.B.) 2International Research Center in Critical Raw Materials for Advanced Industrial Technologies (ICCRAM), Universidad de Burgos, E-09001 Burgos, Spain 3Institute Carlos I for Theoretical and Computational Physics (IC1), E-18016 Granada, Spain *Correspondence: ncorder[email protected] Abstract: Graphene nanostructures have attracted a lot of attention in recent years due to their unconventional properties. We have employed Density Functional Theory to study the mechanical and electronic properties of curved graphene nanoflakes. We explore hexagonal flakes relaxed with different boundary conditions: (i) all atoms on a perfect spherical sector, (ii) only border atoms forced to be on the spherical sector, and (iii) only vertex atoms forced to be on the spherical sector. For each case, we have analysed the behaviour of curvature energy and of quantum regeneration times (classical and revival) as the spherical sector radius changes. Revival time presents in one case a divergence usually associated with a phase transition, probably caused by the pseudomagnetic field created by the curvature. This could be the first case of a phase transition in graphene nanostructures without the presence of external electric or magnetic fields. Keywords: graphene; curvature; quantum revivals; DFT; phase transition 1. Introduction The experimental isolation of a single graphitic layer (now known as graphene) by means of the so-called “Scotch tape method” by Geim and Novoselov [ 1 ] has undoubtedly opened a new field in science. The proof is that in recent years, over 1% of all the scientific publications included in the Web of Science ™ global citation database [ 2 ] are related to this one-atom-thick system. The outstanding properties of this 2D material have led to many application proposals. To cite just a few of them: nanocomposites for bone tissue engineering [ 3 ], composites for multifunctional applications [ 4 ], lubrication [ 5 , 6 ], solar cells [ 7 ], ultracapacitors [ 8 , 9 ], batteries [ 10 ], sensors [ 11 ], catalysis [ 12 ], nanomedicine [ 13 ], fuel cells [ 14 , 15 ] or even energy harvesting [16]. Quantum revivals consist in the temporal periodic reconstruction of a wave packet in systems with a commensurable discrete spectrum [ 17 ]. They have recently attracted attention due to their possible interest in quantum devices [ 18 – 21 ]. Quantum revivals have been studied in infinite graphene under magnetic fields [ 22 – 24 ] or circularly polarized radiation [25], as well as in graphene rings [26] and graphene quantum dots [27]. One of the common misconceptions about graphene is that it is flat. In fact, Peierls instability [ 28 , 29 ] creates ripples in free-standing graphene [ 30 – 35 ]. Atomistic Monte Carlo simulations based on a very accurate many-body interatomic potential for carbon [ 36 ] gave ripples with a size distribution that peaked around 80 Å, in agreement with experiments that yield results in the 50–100 Å range [ 31 ]. Besides, when graphene is grown or deposited on a substrate, nanobubbles with diameters between 80 Å and 1000 Å appear [37–41]. Rectangular graphene flakes either free [ 42 ] or subjected to planar bending [ 43 ] have been analysed in the literature. Conical graphene rings in a magnetic field have also been studied, but only in the continuum limit approximation in the vicinity of the Dirac points Nanomaterials 2022,12, 1953. https://doi.org/10.3390/nano12121953 https://www.mdpi.com/journal/nanomaterials Nanomaterials 2022,12, 1953 2 of 13 and very far from the vertex [ 44 , 45 ]. Both in-plane and out-of-plane bending result in the appearance of a pseudo-magnetic field. We present in this article the first (to the best of our knowledge) calculation of quantum revivals using Density Functional Theory (DFT) and use it to study the interplay between curvature and regeneration times in a non-planar graphene flake. 2. Materials and Methods Since Wallace did the first calculation of the electronic properties of graphene using the tight-binding method [ 46 ], many theoretical models have been used to study this system: molecular mechanics, molecular dynamics, semiempirical methods, Density Functional Theory, Hartree–Fock (HF), post-HF including correlation, Monte Carlo simulations, hybrid methods, continuum models, etc. We have used the Density Functional Theory (DFT) formalism [ 47 ] within the Local Density Approximation (LDA) [ 48 , 49 ] as implemented in the Gaussian 09 [ 50 ] and Gaussian 16 [ 51 ] suites of programs. We have selected DFT for its balance between accuracy and computational effort, and have chosen LDA because it gives better results than gradient corrected approximations (GGAs) for graphitic systems [ 52 – 54 ] and because it has been previously used to successfully study the interaction between carbon nanostructures and several small molecules and atoms [ 55 – 60 ] as well as—very recently—graphene nanoribbons [ 61 ]. We have selected the 6-31G** basis set [ 62 – 67 ] that adds to the 6-31G set d-type and p-type Cartesian–Gaussian polarization functions and is commonly used for carbon nanostructures calculations. In order to avoid non-trivial magnetic ground states induced by asymmetric edge extensions [ 68 ] we have used the hexagonal flake with zig-zag borders passivated with hydrogen atoms shown in Figure 1. The size of the flake (10 concentric hexagonal rings comprising 600 C atoms and 60 H atoms with a separation of approximately 47 Å between opposite vertices) has been selected as a reasonable compromise between the size of experimentally measured ripples and computational cost. Figure 1. Graphene flake used in the calculations (image generated using GausView 6 [69]). Nanomaterials 2022,12, 1953 3 of 13 We have studied (quasi-)spherical nanoflakes of different curvature radii with the three different bending possibilities depicted in Figure 2for a radius of 40 Å: forcing all 600 carbon atoms to lie on a spherical surface (upper panel), fixing only the 60 border atoms on the surface and allowing the rest to relax (middle panel) and forcing only the 12 vertex atoms to belong to the sphere and not imposing any condition on the rest (lower panel). Since carbon atoms tend to be in a sp 2 hybridization state, the nanoflake tends to flatten, as the restrictions are less demanding. Figure 2. Optimized geometries for a (quasi-)spherical graphene flake of radius 40 Å (image generated using GausView 6 [69]). 3. Results and Discussion We have studied both curvature energy and quantum revivals in this nanostructure. The first one for checking if the macroscopic continuum limit model that predicts that curvature energy is proportional to the Gaussian curvature (i.e., the inverse of the radius squared) [ 70 ] holds at this scale. The second one for trying to find trends in regeneration times. 3.1. Curvature Energy We present in Figure 3the dependence of the energy of the nanoflake E with the inverse squared radius of the sphere 1 R2 . Since we have taken the flat configuration as the energy origin, this plot represents the curvature energy. For the case of a perfect spherical surface that corresponds to our fixed surface calculations (shown in red in the figure), 1 R2 is precisely the Gaussian curvature. For the fixed borders (depicted in blue) and fixed vertices (painted in green) cases, curvature is no longer constant, but we can take the value of 1 R2 for the fixed atoms as an estimate of the curvature for comparison purposes. Logically, for a fixed value of R the geometry with only the vertices fixed is energetically more favourable than that with the borders fixed and this in turn is more stable than the one Nanomaterials 2022,12, 1953 4 of 13 with all the atoms forced to lie on a spherical surface. As expected, for the three geometries the curvature energy increases with curvature but the trend for the fixed surface case is somewhat unexpected. The macroscopic continuum result E∝1 R2 would lead to a straight red line, but that does not seem to be the case. Figure 3. Energy of a (quasi-)spherical graphene flake as a function of its curvature. Flat configuration is taken as energy origin. Lines are merely guides for the eye. To check if the macroscopic result holds at least for low curvatures, we present in Figure 4a log–log plot of the curvature energy versus radius. E=k1 R2⇒log E=logk−2log R, (1) and the red graph should present a constant slope equal to − 2. This is clearly not the case and the continuum approximation is not valid for so small a flake. Figure 4. Log–log plot of the energy of a spherical graphene flake as a function of its curvature. Flat configuration is taken as energy origin. Lines are merely guides for the eye. This deviation from the continuum behaviour is undoubtedly due to the difference between the classical physics laws governing the macroscopic world and the quantum ones Nanomaterials 2022,12, 1953 5 of 13 ruling at the nanoscale. This difference translates into two distinct aspects. The equilibrium geometry dictated by both sets of laws is different and, even for the same configuration, they lead to different energies. To check which one is more important in this case, we have made a series of non-self-consistent calculations with a hybrid Molecular Mechanics (MM)/DFT model. We have first used LDA calculations to determine a “semiclassical” Hooke potential for the C–C bond. To this end, we have used a graphene flake similar to the one in Figure 1, but with only 7 concentric hexagonal rings. We have changed the length of the central C–C bond and calculated the energy of the system for both stretched and compressed bond lengths using LDA. We have fitted the energies to a parabolic curve and determined a C–C pair potential for carbon nanoflakes. We have written an in-house code to determine the semiclassical equilibrium geometry for the spherical nanoflakes. This geometry is then used in a single-point calculation to determine the LDA energy of the nanoflake. We will label the results obtained by this procedure as “fixed surface (MM)”. We have included in Figure 3the energies resulting from this hybrid approach in gray. The energies are bigger that those corresponding to the self-consistent LDA calculation (in red) since, according to the variational principle, using a wave function different to that of the fundamental state leads to a higher energy. The greater the curvature, the bigger the energy increase. In order to check if the continuum model is valid for this approach, we have also included the corresponding results in Figure 4. The graph is now nearly a straight line, but there is still some non-linearity. Equation (1) seems to be approximately valid, but the factor in front of log Rdepends slightly on R, adopting the form E=k1 Rn⇒log E=logk−nlog R. (2) We present in Figure 5the value of n as R changes calculated using a 3-point finite differences method. Figure 5. Value of n in Equation (2) as a function of the radius of the spherical carbon nanoflake calculated with the hybrid MM/DFT approximation. It is possible to fit the results in this figure using the simplest Padé approximant n=a0+a1R −1+b1R. (3) The result of this fitting is presented in Table 1. Nanomaterials 2022,12, 1953 6 of 13 Table 1. Results of the fitting of the data in Figure 5to the Padé approximant in Equation (3). a0a1/Å−1b1/Å−1 −0.1818 0.0443 0.0219 We can use this Padé approximant to calculate the asymptotic behaviour of n. lim R→∞n=a1 b1=2.02. (4) This result is very close to 2, which is the value predicted in the continuum model. Therefore, the semiclassical MM/DFT approach tends to the classical continuum macroscopic limit, proving that the main quantum contribution to the deviation from this model is due to the small change in the equilibrium geometry. 3.2. Quantum Revivals We can write the initial state of a time-independent Hamiltonian ˆ H as a linear combination of its eigenfunctions: |Ψ(0)i= ∞ ∑ n=0 an|uni, (5) where anare constants and |uniare the eigenfunctions with energies En, ˆ H|uni=En|uni. (6) The temporal evolution of this state can be written as |Ψ(t)i= ∞ ∑ n=0 an|unie−i ¯hEnt. (7) Since we have calculated the energy of curved graphene nanoflakes using DFT (i.e., solving the Kohn–Sham equations for the system) we know their energy spectrum and can therefore study the time evolution of a wave packet in these systems. Let us consider a superposition of eigenstates of the Hamiltonian concentrated around a central energy level n0 characterised by an energy En0 . We can perform a Taylor expansion of the energy spectrum around En0: En=En0+E0 n0(n−n0) + 1 2!E00 n0(n−n0)2+1 3!E000 n0(n−n0)3+. . . . (8) Taking into account Equations (7) and (8), |Ψ(t)i= ∞ ∑ n=0 an|unie−i ¯hhEn0+E0 n0(n−n0)+ 1 2! E00 n0(n−n0)2+1 3! E000 n0(n−n0)3+...it. (9) Each term in the exponential (except the first one that is just a global phase) defines a characteristic time scale: TCl ≡2π¯h |E0 n0|is called classical time, (10) TRe ≡2π¯h |E00 n0|/2 is called revival time, and (11) TSup ≡2π¯h |E000 n0|/6 is called super-revival time (12) (see [17] for further details). Nanomaterials 2022,12, 1953 7 of 13 We are going to consider a Gaussian initial wave packet, an=1 σ√πe−(n−n0)2 2σ2, (13) where, we have selected the central level as the fourth level above the Highest Occupied Molecular Orbital (HOMO) n0=HOMO + 4 and, in order to get a sharply concentrated packet σ= 0.7 so that only five levels around n0 have a significative contribution (an>0.001). The easiest way of visualising wave packet regeneration is making use of the squared modulus of the so called autocorrelation function that measures the overlap of the wave packet at times 0 and t: |A(t)|2=|hΨ(0)|Ψ(t)i|2. (14) A typical case is presented in Figure 6. |A(t)|2 oscillates very fast, reaching a maximum every TCl inside an envelope with TRe periodicity. TSup is usually much larger that TRe and we will not consider it in this work. 0TCl TRe 0 1/4 1/2 3/4 1 Time A(t)2 Figure 6. A simple example of the time evolution of the squared modulus of the autocorrelation function in blue with its upper envelope in orange. Classical time is marked in red and revival time in green. In principle, calculating TCl is straightforward. One only has to search for the first maximum of |A(t)|2 . Determining TRe is not so easy. In this case, it is necessary to calculate the envelope of the function and calculate its first maximum. However, for some curvatures, the situation is not as clear as the one depicted in Figure 6. If TRe is not much bigger than TCl , both times interfere, the pattern is more complicated and it is difficult to determine both of them, especially TCl . It is then desirable to have another way of calculating these times. The solution is employing Equations (10) and (11). To do that, we have calculated an interpolating function by adjusting a parabolic curve to every three consecutive levels in the spectrum and used it to calculate the first two derivatives of the energy with respect to the level that appears in those equations. Therefore, we have two different ways of calculating these times. We will call the first one numerical (num.) and the second one analytical (an.). Experiments for measuring quantum regeneration times are based on the use of two consecutive laser pulses [ 17 ]. The first one, called pump, creates the initial wave packet, while the second one, called probe, measures its time evolution. By changing the delay Nanomaterials 2022,12, 1953 8 of 13 between the two pulses, different time scales can be explored. This pump–probe scheme was proposed by Alber, Ritsch, and Zoller [ 71 ], and was initially used to study atoms [ 72 , 73 ]. In this case, the second laser pulse was used to ionize the system and the photoionization signal was measured. The method was modified by Zewail (who was awarded the 1999 Nobel Prize in Chemistry for his studies in this field) and his group to study few-atom molecules, using the probe pulse to get the system to an upper fluorescent state and measure the fluorescence [ 74 , 75 ]. More recently, the scheme has been adapted so that the second pulse photoexcites the sample and the differential transmission spectra is analysed. This approach has made possible to study the wave packet evolution of CdSe quantum dots with a mean diameter of 6.4 nm (slightly larger than the diameter of the carbon nanoflake we have selected) [76]. 3.2.1. Classical Time We present in Figure 7classical times determined both numerically (points) and analytically (lines) for the four geometries considered: fixed surface (red), fixed borders (blue), fixed vertices (green) and fixed surface (MM) (gray). Analytical estimations are presented as lines for clarity purposes, but they can only be calculated at the same points as numerical ones. We have simply linked two consecutive points with a straight line. Figure 7. Classical times for a curved graphene nanoflake. Points correspond to numerical values while dotted lines represent analytical ones. The agreement between numerical and analytical estimations is good for all geometries for low curvatures (below approximately 10 −4Å−2 ), but in two of the three self-consistent calculations, this agreement is lost for high curvatures. The reason is that, as we will see in the following subsection, TRe decreases as the curvature increases and it becomes less that one order of magnitude bigger than TCl for the fixed surface and fixed border cases. Recalling Equation (10), TCl is proportional to the inverse of the first derivative of the energy with respect to the level index. This means classical time is related to the spreading of the energy spectrum. TCl decreases as the separation among energy levels increases. If we look at the analytical estimation for the classical times corresponding to the fixed surface and fixed border cases, they grow with curvature, while in the fixed vertices case this initial tendency breaks beyond 1.5 × 10 −4Å−2 and the graph becomes nearly flat. We Nanomaterials 2022,12, 1953 9 of 13 have to remember that this third case corresponds to the geometry with lower restrictions and the system can adapt itself better to the deformation of its fixed points. This means that both total energy (see Figure 3) and energy level separation are less sensitive to curvature. If we consider the numerical estimations for the three self-consistent geometries, TCl remains essentially flat since the slight decrease in energy spectrum spreading is compensated by the interference from TRe . In the fixed surface case this interference overcomes the tendency of TCl to grow and, in fact, it decreases for high curvatures. Finally, if we concentrate on the non-self consistent fixed surface (MM) case, the behaviour is simpler. Numerical and analytical estimations agree perfectly except for the high curvature regime because, as we will see in the next subsection, TRe is much higher than TCl . The spreading of the energy spectrum decreases monotonically with curvature and classical time increases in a linear way. 3.2.2. Revival Time We show in Figure 8revival times determined both numerically (points) and analytically (lines) for the four geometries considered: fixed surface (red), fixed borders (blue), fixed vertices (green) and fixed surface (MM) (gray). Analytical estimations are presented as lines for clarity purposes, but they can only be calculated at the same points as numerical ones. In this case, we have connected two consecutive points with a smoothed line. Figure 8. Revival times for a curved graphene nanoflake. Points correspond to numerical values while dotted lines represent analytical ones. In all cases, numerical and analytical estimates perfectly agree since superrevival times are much higher than revival ones and do not interfere with them. According to Equation (11), TRe is proportional to the inverse of the second derivative of the energy with respect to the level index. This means revival time is related to the non-linearity of the energy spectrum. TRe increases as the spectrum gets closer to being linear (i.e., energy levels tend to be equally spaced). Revival times for the fixed surface and fixed borders cases are very similar and decrease monotonically with curvature. TRe for the fixed vertices case coincides with them up to 1.5 × 10 −4Å−2 and from that point on it remains essentially constant. The situation is