S Nemschokmichal

F Bernhardt

B Krames

J Meichsner

Laser induced fluorescence spectroscopy of N 2 (A 3 Σ + u ) and absolute density calibration by Rayleigh scattering in capacitively coupled rf discharges

Keywords: numbers: 52.25.Os, 52.25.Ya, 52.70.Kz, 52.80.Pi 52.25.Os Emission, absorption, and scattering of electromagnetic radiation 52.25.Ya Neutrals in plasmas 52.70.Kz Optical (ultraviolet, visible, infrared) measurements 52.80.Pi High-frequency and RF discharges

Laser-induced fluorescence spectroscopy of and absolute density calibration by Rayleigh scattering in capacitively coupled rf discharges

S Nemschokmichal, F Bernhardt, B Krames, J

Introduction

Neutral metastable excited atoms and molecules are important carriers of internal energy in a plasma. Due to their large lifetime, they can transfer a large amount of energy from the plasma to surfaces or keep this energy in the afterglow plasma. Especially for larger pressure, they gain importance because of the possibility for exoemission of seed electrons from surfaces, which are important for the ignition of atmospheric pressure discharges. Therefore, it is necessary to quantify metastable densities for the understanding of discharge development.

Taking into consideration the electronic states of molecular nitrogen, the first electronic excited state A 3 Σ + u is metastable with a radiative lifetime of seconds and an energy of 6.2 eV with respect to the ground state X 1 Σ + g . Just above (1.2 eV), two other triplet states (B 3 Π g , W 3 ∆ u ) exist [START_REF] Gilmore | Franck-Condon factors, r-centroids, electronic transition moments, and Einstein coefficients for many nitrogen and oxygen band systems[END_REF], but only the B 3 Π g state has visible optical transitions to the A 3 Σ + u state (1st positive system) [START_REF] Lofthus | The spectrum of molecular nitrogen[END_REF] which are accessible for laser induced fluorescence spectroscopy (LIF).

Examples for the LIF application to low pressure nitrogen discharges can be found in [START_REF] Benedictis | Time-resolved LIF spectroscopy on N 2 (A) metastable in a He/N 2 pulsed rf discharge[END_REF][START_REF] Benedictis | LIF measurement of N 2 (A 3 Σ + u , v = 4) population density in a pulsed rf discharge[END_REF][START_REF] Krames | In situ determination of absolute number densities of nitrogen molecule triplet states in an rf-plasma sheath[END_REF][START_REF] Krames | Spatial density profiles of the N 2 metastable state A 3 Σ + u , the B 3 Π g and C 3 Π u states in an asymmetric rf discharge[END_REF]. They differ in particular in the absolute density calibration, which can be done straight forward by calculation [START_REF] Benedictis | LIF measurement of N 2 (A 3 Σ + u , v = 4) population density in a pulsed rf discharge[END_REF] or by a comparison with scattered Rayleigh light [START_REF] Krames | In situ determination of absolute number densities of nitrogen molecule triplet states in an rf-plasma sheath[END_REF]. Other methods than LIF are recent measurements with threshold ionization mass spectrometry [START_REF] Matsuda | Detection of nitrogen metastable molecules by using the threshold ionization mass spectrometry[END_REF], cavity enhanced absorption spectroscopy [START_REF] Hancock | Absolute number densities of vibrationally excited N 2 (A 3 Σ + u ) produced in a low pressure rf plasma[END_REF] or cavity ring down spectroscopy [START_REF] Horikawa | Absolute densities of N 2 (A 3 Σ + u ), N( 4 S o ), and N( 2 D o ) in an inductively coupled nitrogen plasma source[END_REF][START_REF] Stancu | Time-resolved CRDS measurements of the N 2 (A 3 Σ + u ) density produced by nanosecond discharges in atmospheric pressure nitrogen and air[END_REF]. In particular, for atmospheric pressure discharges the LIF technique is still the most suitable one. Therefore, the simple one laser LIF spectroscopy has been improved to the so called optical-optical double resonance LIF spectroscopy, which makes use of both the first and second positive system [START_REF] Dilecce | Optical-optical double resonance LIF detection of N 2 (A 3 Σ + u ) in high pressure gas discharges[END_REF]. It was applied successfully to a dielectric barrier discharge [START_REF] Dilecce | N 2 (A 3 Σ + u ) density measurement in a dielectric barrier discharge in N 2 and N 2 with small O 2 admixtures[END_REF], where the absolute density calibration has been done by a comparison of the emission of the 2nd positive system, caused by the pooling reaction of two A 3 Σ + u molecules, with the NOγ system, caused by collisions of A 3 Σ + u with the NO ground state. This calibrated diagnostics has also been applied to a surface dielectric barrier discharge [START_REF] Ambrico | On the measurement of N 2 (A 3 Σ + u ) metastable in N 2 surface-dielectric barrier discharge at atmospheric pressure[END_REF][START_REF] Simek | N 2 (A 3 Σ + u ) behaviour in a N 2 -NO surface dielectric barrier discharge in the modulated ac regime at atmospheric pressure[END_REF]. Another possibility in atmospheric pressure discharges is the calibration directly by the pooling reaction, which has been done in a pulsed corona discharge [START_REF] Ono | Laser-induced fluorescence of N 2 (A 3 Σ + u ) metastable in N 2 pulsed positive corona discharge[END_REF][START_REF] Teramoto | Measurement of metastable in N 2 pulsed positive corona discharge with trace amounts of additives[END_REF].

The main goal of this work is the application of the LIF method together with the calibration by Rayleigh scattering in capacitively coupled radio-frequency discharges in nitrogen at pressures up to 1000 Pa which is above the typical low pressure range. Thereby, the most important challenge is the handling of the large quenching in comparison to the low optical emission. This is especially necessary for a further increase of the pressure up to atmospheric pressure discharges, which should be investigated afterwards.

The paper starts with the detailed description of the experiment which is designed to analyse the time response of the LIF signal. This analysis gives the optical emission to the ratio of the quenching by the comparison with the scattered laser light. Afterwards, the calibration by Rayleigh scattering is done. In the end, the measured axial density profiles in the capacitively coupled radio-frequency discharge are shown and discussed for different plasma processing parameters. Furthermore, the influence of small oxygen admixtures is discussed, which is known as a strong quencher of the metastable nitrogen state.

Experimental setup

A sketch of the experimental setup is shown in figure 1, which can be divided into the setup for the capacitively coupled radio frequency (cc-rf) discharge and the laser induced fluorescence (LIF) measurement, consisting of a Nd:YAG pumped dye laser and a combination of a monochromator and a photomultiplier as detector.

CC-RF discharge

The cylindrical stainless steel vacuum chamber (300 mm in diameter, 500 mm height) was evacuated by a turbo pump (Pfeiffer Vacuum, TMU-260, forepump: LABOVAC R PK 8 D) to a base pressure below 10 -4 Pa. Afterwards the chamber was filled with pure nitrogen (pureness of 99.999 %) to pressures from 40 to 1000 Pa by means of a mass flow controller (MKS PR4000F and MKS Mass-Flo R Controller). Furthermore, small admixtures of oxygen up to 2 % were performed by a second mass flow controller. To keep a constant nitrogen gas flow rate of 5 sccm for pure N 2 plasma and 10 sccm in the case of O 2 admixture, the total pressure was independently controlled by a diaphragm pressure gauge (MKS Baratron R Pressure Controller), a butterfly valve (MKS) and a process pump (TRIVAC D25BCSPFPE).

In the centre of the chamber, the discharge was ignited between two parallel plate stainless steel electrodes with diameter of 90 mm and distance of 10 mm. The lower electrode was powered via a matching network (ENI MW-10D) with a radio-frequency generator (ENI ACG-6B) at a fixed frequency of 13.56 MHz. The upper electrode was grounded as well as the walls of the vacuum chamber, wherefore the effective area of the grounded electrode was larger then the powered electrode and an asymmetric capacitively coupled discharge arised. Depending on the applied rf power of between 10 and 80 W as well as the total pressure the resulting self-bias voltage at the powered electrode was in the range from -15 to -350 V.

Laser induced fluorescence

The laser system for the excitation consists of a frequency doubled Nd:YAG laser (Quanta Ray GCR 130) with a maximal pulse energy of 170 mJ as pump laser and a tunable dye laser (LAS LDL 20505). The used dye was pyridine 1 which allowed the tuning of the laser wavelength of between 670 and 730 nm with a linewidth of about 0.6 cm -1 , beam diameter of 1.5 mm and pulse energy up to 30 mJ. In the experiment a maximum laser energy of 10 mJ was used to avoid stimulated emission and full saturation of the LIF signal.

The laser beam was guided parallel through the electrodes, crossing the discharge axis. To measure the density of the state, the laser wavelength was tuned to 687.44 nm, where the band head of the A 3 Σ + u , v = 0 → B 3 Π g , v = 3 (Q 12 and P 11 branches) is located. Perpendicular to the laser beam, the fluorescence was focused by means of a lens (focal length 150 mm) to the entrance slit of a 500 mm monochromator (ARC Acton Research Corporation, SpectraPro R , grating 1200 mm -1 , bandwidth 5 nm) and detected via a photomultiplier (PMT, Hamamatsu R928) afterwards.

The monochromator was tuned to a wavelength of 762 nm, which corresponds to the band head of the

B 3 Π g , v = 3 to A 3 Σ + u , v = 1 transition.
In addition, an RG9 filter was placed in front of the spectrometer to cut off the scattered laser light, which partially transmits the spectrometer even at the wavelength of 762 nm.

For the calibration procedure, the plasma was switched off and the spectrometer was set to 687.4 nm and the RG9 filter was removed to measure the intensity of the scattered laser light by pure nitrogen at room temperature. Because of the large intensity of the scattered light, a neutral density filter (NG3) was build in to avoid saturation effects in the PMT. Additionally, it should be mentioned here, that the linear laser polarization was perpendicular to the observation direction to get a sufficient large Rayleigh scattering intensity, considering the anisotropy of the differential scattering cross section.

The output of the PMT was connected via a tunable resistance R ≈ 1 kΩ with the ground, and the potential drop across the resistance was input to a boxcar integrator (SRS SR 250, input impedance: 1 MΩ). The gate delay of the boxcar integrator with respect to the laser pulse was varied via a pulse delay generator (SRS, DG535) to measure the time dependent LIF signal respectively the time dependent scattered signal for comparison. The voltage output of the boxcar integrator was measured with an A-D converter card (Keithley DAS-801) inside a PC.

U LIF [V] U LIF [V] U lsr [V] U lsr [V] Y LIF = 1.359 s = 0.368 V -1

Results and discussion

In this part the separated steps for the determination of absolute densities are discussed, starting with the consideration of the saturation effect, the analysis of the timedepending LIF signal, the comparison of the optics and the cross sections as well as the Rayleigh measurement itself. Afterwards, axial metastable density profiles for different plasma processing parameters are presented and interpreted in comparison with plasma induced optical emission. At last, the influence of oxygen admixtures on the absolute metastable density is discussed.

Saturation effect

In a first approach, the LIF signal is proportional to the laser pulse energy. This is true for small laser energies, but for larger ones the strong laser field causes also induced emission and therefore a smaller density of the excited state. The result is a saturation curve as shown in figure 2, where the measured LIF voltage U LIF is plotted versus the voltage of the laser energy detector U lsr , which is proportional to the laser pulse energy:

U LIF = Y LIF U lsr 1 + sU lsr . (1)
The dimensionless quantity Y LIF is the LIF signal per laser energy without the saturation effect. The parameter s describes the saturation level. It was calculated once from a fit, so that the LIF signal Y LIF could be calculated directly from U LIF for the following measurements. This was especially necessary because of the changing laser intensity and its temporal drift to smaller laser energies due to the heating of the laser after longer operation.

LIF of N 2 (A 3 Σ + u ) in ccrf discharges 6

Time response of LIF signal

The LIF signal depends strongly on the ratio of the fluorescence intensity at the observed wavelength of 762 nm and the remaining loss processes of the excited B 3 Π g , v = 3 state, such as the fluorescence at other wavelengths and the collisional quenching. In the case that the collisions of the considered B 3 Π g , v = 3 state with other species dominate and any additional production after the laser pulse at t 0 is negligibly, than the rate equation for its density is quite simple:

dn B 3 dt = - 1 τ eff n B 3 . (2)
The effective lifetime τ eff depends on the radiative lifetime τ r of B 3 Π g , v = 3 and the collision rates k q of the excited state with the quenchers q and their densities n q :

1 τ eff = 1 τ r + q k q n q . (3)
In this experiment, the short laser pulse duration (< 10 ns) can be neglected in comparison with the effective lifetime of the B 3 Π g , v = 3 state. At atmospheric pressure, this assumption is no longer valid and the time dependent laser pulse has to be taken into account in equation ( 2), too.

The solution of the rate equation is an exponential decay with the effective lifetime τ eff after the laser pulse which results also in an exponential decay for the LIF intensity

I LIF (t) = A B 3 →A 1 n B 3 (t), I LIF (t) = ÎLIF e - t-t 0 τ eff , (4)
where

A B 3 →A 1 = 6.446 • 10 4 s -1 [17] is the Einstein coefficient for spontaneous emission from B 3 Π g , v = 3 to A 3 Σ + u , v = 1.
In the experiment, the light passes the spectrometer and releases charges in the PMT, that drive an electric current via the external resistor to the ground. This is a simple RC circuit with a time constant of τ PMT = RC, which has been proven by varying the resistor R. Therefore, the voltage u LIF across the resistor has an additional exponential decay due to the electrical circuit of the PMT and the external resistor, which has to be calculated as a convolution of this exponential decay with the LIF intensity:

u LIF (t) ∝ t t 0 I LIF (t )e -t-t τ PMT dt , (5)
Calculation of equation ( 5) results in a difference of two exponential decays with the time constants τ eff due to the lifetime of the excited state and τ PMT due to the PMT circuit, which are in the microsecond respectively sub-microsecond time scale: Here, the ratio of the fluorescence rate, represented by the Einstein coefficient

u LIF (t) = ûLIF A B 3 →A 1 τ -1 eff τ eff τ eff -τ PMT × e - t-t 0 τ eff -e - t-t 0 τ PMT . (6)
A B 3 →A 1 to the overall loss rate τ -1 eff , is of importance to consider the losses of B 3 Π g , v = 3 by collisional quenching.
In the later calibration of the LIF signal by means of the Rayleigh scattering, the different temporal behaviour of the scattered light has to be considered. The scattered light appears immediately when the laser pulse takes place, so that there is no delay as for the fluorescence. In addition, the decay time of the PMT circuit is adjusted to be in the microsecond range, so that the laser pulse duration can be neglected again. Therefore, the voltage across the resistor shows a single exponential decay determined by the PMT circuit itself:

u scatt (t) = ûscatt e - t-t 0 τ PMT . (7)
As mentioned before, the voltages u LIF and u Ray across the resistance are not the measured ones, therefore, these voltages have additionally to be integrated, which yield the voltages at the output of the boxcar integrator, with a gate delay t GD and a gate width t GW :

U LIF/scatt (t GD ) = 1 t GW t GD +t GW t GD u LIF/scatt (t)dt. ( 8
)
The integrated voltage for the scattered signal remains an exponential decay with τ PMT as the time constant and another factor in front of the exponential function. The temporal signal for the LIF remains mainly the same in comparison with [START_REF] Krames | Spatial density profiles of the N 2 metastable state A 3 Σ + u , the B 3 Π g and C 3 Π u states in an asymmetric rf discharge[END_REF], but it starts with a finite value at t GD = 0 because of the finite gate width of the boxcar integrator. This is clearly recognizable in figure 3, where the scattered signal and the LIF signal are plotted versus the gate delay of the boxcar integrator.

In addition, the Levenberg-Marquardt algorithm is used to plot two fit functions in figure 3 for the LIF and scattering signal from equation [START_REF] Hancock | Absolute number densities of vibrationally excited N 2 (A 3 Σ + u ) produced in a low pressure rf plasma[END_REF]. The agreement of the exponential fit with the values from the measurement of the scattered light is very good and allows to find out the PMT decay time constant τ PMT . Following, the fit of the LIF signal is done and the effective lifetime τ eff and the signal strength ŷLIF were extracted. There, the signal strength ŷLIF represents a kind of original signal strength if there would not be any delay due to the lifetime of the excited state and the PMT circuit. For fitting the LIF signal only the measurements until 1 µs are taken into account to avoid the influence of backward processes, for example from the W 3 ∆ u state (see below).

τ -1 eff [MHz] τ -1 eff [MHz]
Considering equation ( 3) and assuming that the main quenching partner is the background gas (ground state X 1 Σ + g with the density n X = p/k B T gas ), then the effective decay rate τ -1 eff should depend linear on pressure. To proof this, the effective decay rates τ -1 eff are plotted versus pressure in figure 4.

First of all, one can see that the overall decay rate τ -1 eff is dominated by quenching and only the small intercept is caused by optical transitions. For example, the overall decay rate at 1000 Pa of τ -1 eff ≈ 6.5 MHz is about 100 times larger than the Einstein coefficient for the considered optical transition at 762 nm.

Obviously, the effective decay rates follow not a straight line. For comparison, some curves are drawn using quenching rates from literature, but these rates differ strongly from the measured values. One reason might be the collisional coupling of the B 3 Π g state to the W 3 ∆ u state, observed already in [START_REF] Rotem | Laser-induced fluorescence studies of collisional coupling of N 2 (B 3 Π g ) with N 2 (W 3 ∆ u ) and N 2 (A 3 Σ + u )[END_REF][START_REF] Sadeghi | Collisional coupling and relaxation of N 2 (B 3 Π g ) and N 2 (W 3 ∆ u ) vibrational levels in Ar and Ne[END_REF][START_REF] Piper | State-to-state N 2 (A 3 Σ + u ) energy pooling reactions. II. The formation and quenching of N[END_REF]:

B 3 Π g + X 1 Σ + g W 3 ∆ u + X 1 Σ + g . (9)
This could also explain the small increase in figure 3 of the LIF signal after 1 µs. In [START_REF] Sadeghi | Collisional coupling and relaxation of N 2 (B 3 Π g ) and N 2 (W 3 ∆ u ) vibrational levels in Ar and Ne[END_REF], the rate coefficient for the backward reaction in ( 9) is assumed to be much smaller than for the forward reaction, which might result in a larger B 3 Π g density late after the laser pulse.

The change of the effective lifetime has to be considered together with the variation of the plasma processing parameters. As it is obvious and shown, the effective lifetime

LIF of N 2 (A 3 Σ + u ) in ccrf discharges 9
depends on pressure. Furthermore, an effect due to the larger gas temperature at larger rf power is expected. Firstly, the mean thermal velocity of the background gas depends on temperature (v th ∝ T gas ), which yields in a first approximation to a temperature dependent quenching rate coefficient of k X proportional to T gas . In addition, the density of the most important quencher, the background gas in its ground state X 1 Σ + g depends on pressure by the ideal gas law for constant pressure: n = p/k B T gas . As a result, both effects give an effective decay rate τ -1 eff ≈ k X n X proportional to T -1/2 . The analysis of the rotational temperature from the

B 3 Π g , v = 2 → A 3 Σ + u , v = 0 and B 3 Π g , v = 3 → A 3 Σ + u , v = 1
bands provide a temperature growth of 200 K for a power increase from 10 to 70 W at 200 Pa. Unfortunately, the rotational temperature for low rf power (10 W) could not be determined and the gas temperature of 350 K slightly above the room temperature was assumed. For this, a factor of 550/350 ≈ 1.25 remains from the discussion above. To proof this, the rf power at 200 Pa was varied, but only a small increase of τ eff within the error range was observed. Therefore, the effective lifetime is taken to be constant for selected pressure. Nevertheless, this might be a reason for the deviations from the literature vaules in figure 4 for 700 and 1000 Pa, because at larger pressures also larger rf power was necessary for sustaining the discharge and stronger gas heating is expected.

Another influence on the effective lifetime might be due to the admixtures of oxygen. As it is known, the quenching of A 3 Σ + u by oxygen is large, but the quenching of the excited state B 3 Π g by oxygen is in the same order as the quenching by nitrogen. Consequently, the small admixtures of oxygen should not affect the effective lifetime. This was proofed experimentally, and there are no differences of the effective lifetimes with and without oxygen admixture within the error range.

Absolute density calibration

The absolute density calibration is done by the comparison of the LIF signal at 762 nm to the Rayleigh scattering signal at 687.4 nm to avoid the calculation of the interaction volume between the laser and the detection cone. The whole formula is given by

n A 0 = ŷLIF T Ray Q Ray T LIF Q LIF 3 2 σ Ray σ LIF n Ray ŷRay ( 10
)
for the density of the A 3 Σ + u , v = 0 state. First of all, it enters the signal strength ŷLIF , calculated from the integrated boxcar signal by use of the effective lifetime, the PMT decay time and the saturation factor. Furthermore, the ratios of the transmission optics at 687 nm (T Ray ) respectively 762 nm (T LIF ) and the quantum efficiencies of the PMT (Q Ray , Q LIF ) appear, as well as the ratio of the cross sections σ Ray and σ LIF . The additional factor 3/2 will be discussed below. At the end, the ratio of the known density of the Rayleigh scatterer n Ray to the Rayleigh signal strength ŷRay is included.

The ratio of the optics was determined by the comparison of the signal from a cold light lamp with known temperatures of between 2650 and 3200 K for the different wavelengths. This yielded a factor of 0.37 ± 0.02 for better quantum efficiency for the lower wavelength, partially compensated by the better transmission of the spectrometer at 762 nm (grating blazed for 750 nm). Additionally, a neutral density filter was used for the Rayleigh scattering experiment to avoid saturation inside the PMT.

T Ray Q Ray /T LIF Q LIF , because of the LIF of N 2 (A 3 Σ + u ) in ccrf discharges 10 U LIF [V] U LIF [V] σ LIF [10 -18 cm 2 ] σ LIF [10 -18 cm 2 ] λ lsr [nm] λ lsr [nm]
For the calculation of the LIF cross section σ LIF at the laser wavelength λ lsr = 2πc/ω lsr the convolution of the laser profile with the line profiles of the

A 3 Σ + u , v = 0 → B 3 Π g , v = 3 rotational band is essential: σ LIF (ω lsr ) = k F k ω σ k (ω)G(ω, ω lsr )dω. ( 11
)
The factor F k is the part of the rotation line at ω k of the whole

A 3 Σ + u , v = 0 → B 3 Π g , v = 3 band, which implies k F k = 1.
These factors F k are given by a numerical calculation which is based on data from high resolution Fourier spectrometry [START_REF] Roux | High-resolution Fourier spectrometry of 14 N 2 infrared emission spectrum: Extensive analysis of the B 3 Π g -A 3 Σ + u system[END_REF] and was already used for this calibration method [START_REF] Krames | In situ determination of absolute number densities of nitrogen molecule triplet states in an rf-plasma sheath[END_REF]. The first factor

σ k (ω) = πe 2 f ji 2ε 0 m e c 0 f (ω, ω k ) ( 12
)
in the convolution integral [START_REF] Dilecce | Optical-optical double resonance LIF detection of N 2 (A 3 Σ + u ) in high pressure gas discharges[END_REF] represents the frequency dependent LIF cross section for a line centered at ω k . It depends on the oscillator strength f ji = 1.404 • 10 -4 [START_REF] Laux | Arrays of radiative transition probabilities for the N 2 first and second positive, no beta and gamma, N + 2 first negative, and O 2 Schuhmann-Runge band systems[END_REF] and the Doppler broadened line profile f (ω, ω k ). As visible in equation [START_REF] Dilecce | Optical-optical double resonance LIF detection of N 2 (A 3 Σ + u ) in high pressure gas discharges[END_REF], the cross section for one line was convoluted with the laser profile G(ω, ω lsr ), which was assumed to be a Gaussian distribution with a FWHM of ∆ν ≈ 0.6 cm -1 , measured by a laser scan of a optogalvanic calibration lamp with neon.

To proof the calculation, the cross section for two rotational temperatures is plotted together with a measured laser scan in figure 5. The agreement with the simulation is very well, so it is possible to take the cross section from the simulation at the laser wavelength of 687.44 nm. Unfortunately, the cross section decreases with temperature due to the changing distribution of the rotational levels within the vibrational state.

Therefore, one has to take care about this effect if the power of the discharge increases significantly.

The calculation of the cross section for the Rayleigh scattering is less complicated. It can be done as by Miles et al. [START_REF] Miles | Laser Rayleigh scattering[END_REF] by

σ Ray = 8π 3 3ε 2 0 λ 4 α 2 6 -3ρ 0 6 -7ρ 0 ≈ 1.87 • 10 -27 cm 2 , (13)
where α = 1.76 7 • 10 -30 m 3 (cgs units) is the molecular polarizability of nitrogen and ρ 0 ≈ 0.01 the depolarization [START_REF] Bridge | The polarisation of laser light scattered gases[END_REF]. Another value for the molecular polarizability of 1.74 • 10 -30 m 3 (cgs) from Lide [START_REF] Lide | CRC Handbook of Chemistry and Physics[END_REF] yields a similar value of 1.81 • 10 -27 cm 2 for the cross section. In addition, more recent studies from Naus et al. [START_REF] Naus | Experimental verification of Rayleigh scattering cross sections[END_REF] yield 1.87 • 10 -27 cm 2 (theoretical) and 1.87 • 10 -27 cm 2 (experimental). Therefore, (1.84 ± 0.03) • 10 -27 cm 2 is taken as mean value for the Rayleigh cross section.

Although the calculation of the Rayleigh cross section is very simple, one has to take care because of the difference between the integral cross section and the differential cross section dσ/dΩ. For the LIF case the latter is just σ LIF /4π because of the isotropic fluorescence due to the fast rotation of the molecules. This is different for the Rayleigh scattered light, because it is scattered immediately and the polarization of the laser determines also the polarization of the scatterer. As a result, the differential cross section in the direction of the laser polarisation is very small and the cross section perpendicular to the laser polarization larger than just σ Ray /4π, namely 3/2σ Ray /4π [START_REF] Bridge | The polarisation of laser light scattered gases[END_REF][START_REF] Miles | Laser Rayleigh scattering[END_REF]. Therefore, the polarization of the laser light is important, which should be perpendicular to the observation direction as in this experiment. As a result, the anisotropic cross section of the Rayleigh scattering yields mainly a factor of 3/2 in comparison to the isotropic cross section for the LIF signal.

The last factor of equation ( 10) is the ratio of the density of the Rayleigh scatterer n Ray to the Rayleigh scattering signal represented by the original signal strength ŷRay . The scattering experiment is performed in pure nitrogen without plasma, so that the scatterer is the ground state of the nitrogen molecule. Unfortunately, one cannot measure the Rayleigh scattered light without the reflections of the laser at the walls of the vacuum chamber and the electrodes. Therefore, the pressure was varied between an empty chamber and atmospheric pressure, and the measured signal is plotted in figure 6. As expected, the measurement points are on a straight line with an intercept, which represents the pressure independent signal due to reflections at the wall. The slope gives the (inverted) ratio n Ray /ŷ Ray = (2.95 ± 0.07) • 10 22 cm -3 of the scatter density to the measured voltage.

Altogether, using equation [START_REF] Stancu | Time-resolved CRDS measurements of the N 2 (A 3 Σ + u ) density produced by nanosecond discharges in atmospheric pressure nitrogen and air[END_REF] for the example of 200 Pa and 10 W (figure 3) a density of (3.4 ± 0.6) • 10 12 cm -3 can be calculated for a distance of 6 mm from the powered electrode.

Absolute densities and axial profiles

In this section, the measured A 3 Σ + u , v = 0 densities, their parameter dependencies and axial density profiles are discussed. First of all, the pressure dependence of the maximum metastable density is shown in figure 7.

As it is expected, the metastable densities increase with pressure. This increase is steeper for low pressures and flatter for larger pressure. Furthermore, a jump is observed between 100 and 200 Pa, which might be just due to the different rf power of these measurements (see discussion below). The agreement with the literature is good, although the small discharge gap of 10 mm in this experiment differs strongly from the open discharge configuration with one planar electrode in [START_REF] Krames | In situ determination of absolute number densities of nitrogen molecule triplet states in an rf-plasma sheath[END_REF]. Therefore, the densities in [START_REF] Krames | In situ determination of absolute number densities of nitrogen molecule triplet states in an rf-plasma sheath[END_REF] should be larger due to reduced losses at the electrodes, but this might have been compensated by the missing factor of 3/2 from the comparison of the angular dependent Rayleigh cross section to the LIF cross section in equation [START_REF] Stancu | Time-resolved CRDS measurements of the N 2 (A 3 Σ + u ) density produced by nanosecond discharges in atmospheric pressure nitrogen and air[END_REF].

Axial density profiles of A 3 Σ u , v = 0 for various pressures are plotted in figure 8. In addition, the plasma induced optical emission of the transition from

B 3 Π g , v = 3 to A 3 Σ + u , v = 1 is shown for comparison.
The measurements were done for the lowest possible rf power to avoid heating of the background gas which makes the analysis more

LIF of N 2 (A 3 Σ + u ) in ccrf discharges 13 OES LIF ūOES [V] ūOES [V] n A0 [10 12 cm -3 ] n A0 [10 12
B 3 Π g , v = 3 → A 3 Σ + u , v = 1
and density profiles of A 3 Σ + u , v = 0 for various pressures. The powered electrode is at 0 mm.

difficult. The profiles are almost symmetric because of the small self-bias voltage at the powered electrode, which causes a nearly symmetric profile of the time averaged electric potential. Only for 40 Pa the asymmetric character is clearly visible for the plasma induced optical emission. The reason is the large sheath width at 40 Pa in front of the powered electrode, which is even larger than half of the discharge gap. For larger pressure, at 700 and 1000 Pa, the optical emission at 2.5 mm is slightly higher than at 7.5 mm which reflects the asymmetry of the discharge. For the density profiles of A 3 Σ + u , such an effect is not visible and the density is nearly constant in the middle of the discharge due to the large effective lifetime of this state. Only nearby the metallic electrodes, the de-excitation processes of A 3 Σ + u cause a steep decrease of the density. To compare the losses in the volume with the losses due to diffusion one can consider approximately a cylindrical volume in the discharge centre with the height of 5 mm corresponding to the difference between the axial positions of 2.5 and 7.5 mm, and the cross section area A. Within this region the metastable density is assumed to be constant. A comparison of the different quenching processes in the discharge centre is given in table 1. Despite the clearly larger ground state density, the quenching by the background gas is negligible and the pooling reaction dominates in pure nitrogen. Furthermore, the quenching by oxygen is also large, which will be discussed below. Besides the pooling reaction as the most important loss process in the volume, the diffusion and the loss at the electrodes is also important. Considering the assumed cylinder in the discharge middle, the diffusion in the radial direction is negligible because of the large electrode area in comparison to the small discharge gap, which implies a homogeneous discharge in the radial direction nearby the discharge axis. The axial

(A 3 Σ + u ) + N 2 (X 1 Σ + g ) → products n A 0 n X • 3.0 • 10 -16 cm 3 s -1 [28] 1.7 58 510 10 12 cm -3 s -1 N 2 (A 3 Σ + u ) + N 2 (A 3 Σ + u ) → products n 2 A 0 • 4 • 10 -10 cm 3 s -1 [29]
1.4 64 200 10 14 cm -3 s -1 n 2 A 0 • 1.5 • 10 -10 cm 3 s -1

[30] 0.5 24 76 10 [START_REF] Stewart | The reflection of metastable particles at a surface[END_REF]. The resulting diffusion current density j = D • dn A 0 /dz is nearly equal for 40 and 200 Pa due to the 1 : 1 increase of the metastable density with pressure (see figure 7). Between 200 and 1000 Pa, the metastable density only doubles. Hence, the density gradient is only slightly larger and the diffusion current decreases clearly due to the five times smaller diffusion coefficient. To estimate the diffusion losses in the cylinder, the diffusion current must be multiplied by the cylinder cross section area A and a factor 2 for the both sides. To get the losses per volume, this quantity is divided by the cylinder volume A • 5 mm, which cancels out the cylinder cross section area.

14 cm -3 s -1 N 2 (A 3 Σ + u ) + O 2 → products n A 0 n O 2 • 2.5 • 10 -12 cm 3 s -1 [29] 50 10 14 cm -3 s -1
(A 3 Σ + u ) molecule in N 2 is given from Dp = 2.1 • 10 4 Pa • cm 2 s -1
The comparison of the pooling reaction as the most important volume process with the diffusion shows, that for 40 Pa the diffusion is most important, whereas for 1000 Pa the pooling reaction dominates. A transition region exists at 200 Pa, where both processes have the same order of magnitude.

Due to the different loss processes one would expect different shapes of the density profiles. Therefore, the normalized densities are plotted in figure 9. Nevertheless, the shapes of the density profiles are the same. Only for 40 Pa some values differ due to the larger error for these measurements.

Another important discharge parameter is the rf power respectively the resulting self-bias voltage which is more characteristic for the rf plasma. Therefore, the optical emission from B 3 Π g , v = 3 to A 3 Σ + u , v = 1 and the densities of A 3 Σ + u , v = 0 are plotted again for different rf power in figure 10. Firstly, as it is expected, the asymmetric character of the optical emission profiles increases with the self-bias voltage and the maximum shifts toward the powered electrode, because electron impact excitation is the most important excitation process, located at the sheath edge [START_REF] Krames | Spatial density profiles of the N 2 metastable state A 3 Σ + u , the B 3 Π g and C 3 Π u states in an asymmetric rf discharge[END_REF]. In contrast to the optical emission, the A 3 Σ + u , v = 0 density profiles are still symmetric, whereas the absolute A 3 Σ + u , v = 0 density decreases with increasing self-bias voltage. This is a big difference to the measurements at lower pressures (and larger gaps) in [START_REF] Krames | Spatial density profiles of the N 2 metastable state A 3 Σ + u , the B 3 Π g and C 3 Π u states in an asymmetric rf discharge[END_REF] in a capacitively coupled rf discharge and [START_REF] Hancock | Absolute number densities of vibrationally excited N 2 (A 3 Σ + u ) produced in a low pressure rf plasma[END_REF] in an inductively coupled rf discharge, where the density increases with the self-bias voltage resp. rf power. One reason is the decrease of the LIF cross section due to the change in the rotational distribution caused by the temperature increase with rf power (cp. figure 5), which is not considered in the calibration due to the unknown absolute temperature. But this cannot explain the large decrease in figure 10. Reconsidering the results in [START_REF] Krames | Spatial density profiles of the N 2 metastable state A 3 Σ + u , the B 3 Π g and C 3 Π u states in an asymmetric rf discharge[END_REF], the increase is also related to an increase of the B 3 Π g density. This is not the case in our measurements, because the mean optical emission from B 3 Π g , v = 3 to A 3 Σ + u , v = 1 remains nearly constant. Hence, assuming that the cascade processes from the larger triplet states is most important one for the production of A 3 Σ + u , this means that the production rate should also be nearly constant with self-bias voltage. Therefore, the dependence of the A 3 Σ + u , v = 0 density on self-bias voltage is determined by the loss processes. These loss processes are also temperature dependent and increase with rf power. The rate coefficient of the pooling reaction is approximately proportional to √ T due to the larger thermal velocity. The other loss process, the diffusion, is determined by the diffusion coefficient D, which also depends on temperature. Firstly, D is proportional to 1/n, which results for constant pressure in a linear increase with temperature. Secondly, the factor √ T results again from the thermal velocity, so altogether D ∝ T 3/2 , which means a significant increase of the diffusion coefficient with temperature. Assuming again a temperature of 350 K at 10 W and 550 K at 50 W, the diffusion coefficient enlarges by a factor of 2. In summary, the temperature increase with power explains the decrease in metastable density.

LIF of N 2 (A 3 Σ + u ) in ccrf discharges 15 LIF n A0 / n A0 n A0 / n
LIF of N 2 (A 3 Σ + u ) in ccrf discharges 16 OES LIF ūOES [V] ūOES [V] n A0 [10 12 cm -3 ] n A0 [10 12 cm -3 ] |self-bias| [V] |self-bias| [V] 0.0 % O 2 0.5 % O 2 1.0 % O 2 2.0 % O 2 200

Small oxygen admixtures

As it is known from literature, the quenching of A 3 Σ + u by oxygen is much larger than by nitrogen itself, see table 1. For an admixture of 1 % oxygen the loss rate has the same order of magnitude as the pooling reaction in the discharge center. At the boundaries, it even exceeds the pooling reaction and becomes the dominant loss process. Furthermore, for large dissociation degrees of oxygen the losses by quenching due to atomic oxygen become dominante, because of the large quenching rate of 2.7 • 10 -11 cm 3 s -1 [START_REF] Herron | Evaluated chemical kinetics data for reactions of N( 2 D), N( 2 P ), and N 2 (A 3 Σ + u ) in the gas phase[END_REF]. However, the admixture of oxygen should decrease the A 3 Σ + u density significantly. Surprisingly, this is not always the case, as it is seen from figure 11, where the mean plasma induced optical emission and the A 3 Σ + u densities are plotted versus the self-bias voltage for various oxygen admixtures.

In the case of pure nitrogen plasma, the density of A 3 Σ + u decreases with rf power as it is discussed above. From 20 to 220 V, a factor of 2.5 is visible. In addition, the related optical emission intensity remains constant. In contrast to that, in the case of an oxygen admixture of 0.5 % the A 3 Σ + u density starts at a lower level and increases to a maximum with a larger A 3 Σ + u density compared with pure nitrogen. The 2.5 times lower density at small rf power is due to the larger quenching mentioned above. The increase for larger self-bias voltage is explainable by the increasing optical emission, which is a measure for the production rate due to the cascade process C 3 Π u → B 3 Π g → A 3 Σ + u . From 50 V selfbias to 210 V it increases by a factor of 5 with 0.5 % oxygen. Comparing the metastable densities for pure nitrogen and 0.5 % oxygen at 210/220 V, the twice larger value for the case with oxygen inside is explainable by the 2.5 times larger quenching and the 5 times larger production from B 3 Π g . For larger oxygen concentrations, the optical emission is similar to the case of 0.5 % oxygen admixture, so that the larger quenching rates result in a lower A 3 Σ + u density as expected.

Conclusions

The density of the lowest metastable A 3 Σ + u state of nitrogen in its vibrational ground state was successfully measured by laser induced fluorescence spectroscopy in combination with an absolute density calibration by Rayleigh scattering. The different time response of both signals allowed to measure the effective lifetime of the excited B 3 Π g state and to calculate densities without the large uncertainties of the quenching coefficients, even for pressures up to 1000 Pa.

The axial density profiles for the low pressure capacitively coupled radio-frequency discharges show similar shapes, although the dominating loss process for low pressures is diffusion to the electrodes and for larger pressure the collisional quenching within the volume. The increase of discharge power resp. self-bias voltage results in a lower metastable density, which can be explained by an increase of the rotational and kinetic (gas) temperature, respectively. Furthermore, the admixture of oxygen causes the well known decrease of the density due to the larger quenching, on the other hand it increases also the production of A 3 Σ + u due to the cascade processes from the higher triplet states so that the A 3 Σ + u density exceeds the value for the pure nitrogen case for some parameter sets.

Figure 1 .

 1 Figure 1. Experimental setup of the discharge and the LIF measurement.

Figure 2 .

 2 Figure 2. Saturation of the LIF signal depending on laser energy. The fit gives the saturation parameter s.

YFigure 3 .

 3 Figure 3. LIF voltage from the boxcar integrator depending on the gate delay of the boxcar for the LIF signal and the signal of the scattered laser light. The lines represent fits of the theoretical time responses.

Figure 4 .

 4 Figure 4. Effective decay rates τ -1eff of the excited B 3 Π g , v = 3 state depending on pressure. The straight lines represent quenching rate coefficients from literature (a[START_REF] Piper | State-to-state N 2 (A 3 Σ + u ) energy pooling reactions. II. The formation and quenching of N[END_REF], b[START_REF] Benedictis | Excitation and decay of N 2 (B 3 Π g , v) states in a pulsed discharge: Kinetics of electrons and long-lived species[END_REF]).

200Figure 5 .

 5 Figure 5. Comparison of the relative LIF intensity with a laser scan and the simulated LIF cross sections.

LIF of N 2 (n 3 200Figure 6 .

 236 Figure 6. Pressure dependent scattering signal. The slope represents the ratio of the signal due to Rayleigh scattering and the (ground state) density of nitrogen.

Figure 7 .

 7 Figure 7. Maximal A 3 Σ u , v = 0 densities for pressures from 40 to 1000 Pa and comparison with literature ( [4],•[START_REF] Krames | In situ determination of absolute number densities of nitrogen molecule triplet states in an rf-plasma sheath[END_REF]).

Figure 8 .

 8 Figure 8. Axial optical emission profiles ofB 3 Π g , v = 3 → A 3 Σ + u , v = 1and density profiles of A 3 Σ + u , v = 0 for various pressures. The powered electrode is at 0 mm.

Figure 9 .Figure 10 .

 910 Figure 9. axial density profiles of A 3 Σ + u , v = 0 for various pressures. The powered electrode is at 0 mm.

Figure 11 .

 11 Figure 11. Mean optical of B Π g , v = 3 → A Σ + u v = 1 and mean densities of A 3 Σ + u , v = 0 at 200 Pa depending on self-bias voltages for pure nitrogen and small oxygen admixtures.

Table 1 .

 1 Densities of important quenching partners and the corresponding quenching rates at 300 K for 40, 200 and 1000 Pa.

	ref.	40 Pa 200 Pa 1000 Pa

Table 2 .

 2 Quantification of the diffusion losses at discharge centre for 40, 200 and 1000 Pa.diffusion losses out of the cylinder can be estimated the density gradients at 2.5 mm and 7.5 mm. These gradients are determined by a linear fit in figure8(exemplary for 200 Pa) and are noted in table 2. In addition, the diffusion coefficient D for the excited N 2

	40 Pa 200 Pa 1000 Pa


R Chowdhury

S Adhikari
email: s.adhikari@swansea.ac.uk

F Scarpa

M I Friswell

Transverse vibration of single layer graphene sheets

Keywords: numbers: 62.25.-g Mechanical properties of nanoscale systems numbers: 68.55.-a Thin film structure and morphology numbers: 61.43.Bn Structural modeling: serial-addition models, computer simulation graphene, bending, molecular mechanics, frequencies

HAL is

Introduction

In 2003, Gan et. al. [START_REF] Gan | STM investigation on interaction between superstructure and grain boundary in graphite[END_REF] described the in-situ exfoliation of a single layer with an intersecting grain boundary in graphite using an STM operated in air. The exfoliation technique has led to the successful production of graphene, opening a new era in the field of nanoelectronics [START_REF] Novoselov | Electric Field Effect in Atomically Thin Carbon Films[END_REF][START_REF] Novoselov | Two-dimensional gas of massless Dirac fermions in graphene[END_REF][START_REF] Chowdhury | Graphene-based biosensor using transport properties[END_REF]. The very high in-plane stiffness of graphene sheets [START_REF] Huang | Thickness of graphene and single wall carbon nanotubes[END_REF] has suggested some possible use of graphite nanosheets for nanosensors and NEMS applications [START_REF] Jensen | An atomic-resolution nanomechanical mass sensor[END_REF][START_REF] Seoanez | Dissipation in graphene and nanotube resonators[END_REF], due to their extremely high surface to volume ratio, as well as large deflection capability under point loading [START_REF] Bao | Controlled ripple texturing of suspended graphene and ultrathin graphite membranes[END_REF][START_REF] Cranford | Folding multilayer graphene sheets[END_REF][START_REF] Lee | Measurement of the Elastic Properties and Intrinsic Strength of Monolayer Graphene[END_REF][START_REF] Scarpa | The bending of single layer graphene sheets: Lattice versus continuum approach[END_REF][START_REF] Schwartzbart | Molecular statical calculation of graphene sheet buckling[END_REF][START_REF] Duan | Nonlinear bending and stretching of a circular graphene sheet undera central point load[END_REF]. In this paper the out-of-plane or transverse vibration of single-layer graphene sheets (SLGSs) is considered. The vibration studies could be useful for graphene based mass and/or gas sensors [START_REF] Jensen | An atomic-resolution nanomechanical mass sensor[END_REF][START_REF] Arsat | Graphene-like nano-sheets for surface acoustic wave gas sensor applications[END_REF][START_REF] Arsat | Graphene-like nano-sheets for surface acoustic wave gas sensor applications[END_REF][START_REF] Rangel | Graphene Terahertz Generators for Molecular Circuits and Sensors[END_REF][START_REF] Mohanty | Graphene-Based Single-Bacterium Resolution Biodevice and DNA Transistor: Interfacing Graphene Derivatives with Nanoscale and Microscale Biocomponents[END_REF][START_REF] Moradian | Investigation of gas sensing properties of armchair graphene nanoribbons[END_REF][START_REF] Hadlington | Graphene sensor achieves ultimate sensitivity[END_REF].

The vibration of single and multiple layer graphene sheets has been investigated by several authors, using either continuum mechanics approaches [START_REF] Kitipornchai | Continuum model for the vibration of multilayered graphene sheets[END_REF], equivalent lattice structures made by atomistic-continuum models representing the C-C bonds [START_REF] Hashemnia | Vibrational analysis of carbon nanotubes and graphene sheets using molecular structural mechanics approach[END_REF], and molecular dynamics approaches combined with continuum mechanics for thickness identification [START_REF] Gupta | Elastic properties and frequencies of single-layer graphene sheets[END_REF]. The out-of-plane deformation of SLGS has been considered using the continuum mechanics models [START_REF] Duan | Nonlinear bending and stretching of a circular graphene sheet undera central point load[END_REF][START_REF] Hemmasizadeh | A method for developing the equivalent continuum model of a single layer graphene sheet[END_REF],together with continuum and truss-like structural assemblies [START_REF] Batra | Continuum models of multi-walled carbon nanotubes[END_REF][START_REF] Sakhaee-Pour | Elastic buckling of single-layered graphene sheet[END_REF][START_REF] Sakhaee-Pour | Elastic properties of single-layered graphene sheet[END_REF][START_REF] Odegard | Equivalent-continuum modeling of nano-structured materials[END_REF][START_REF] Chowdhury | Low frequency vibration of multiwall carbon nanotubes with heterogeneous boundaries[END_REF][START_REF] Tserpes | Finite Element modelling of single-walled carbon nanotubes[END_REF][START_REF] Scarpa | Effective elastic mechanical properties of single layer graphene sheets[END_REF][START_REF] Scarpa | Mechanical properties of non reconstructed defective single wall carbon nanotubes[END_REF]. In a recent paper a nonlinear mechanical model has been used [START_REF] Sadeghi | Nonlinear vibrational analysis of single-layer graphene sheets[END_REF] to take account for large deformations in SLGS. They observed higher resonance frequencies from the nonlinear model compared to the equivalent linear model. A molecular mechanics (MM) approach based on the computation of the Hessian matrix and its eigenvalues has been proposed by some authors to describe the structural dynamics of single wall carbon nanotubes [START_REF] Chowdhury | A molecular mechanics approach for the vibration of single walled carbon nanotubes[END_REF], and used to validate a lattice structural mechanics approach in nanoribbons [START_REF] Scarpa | Wave propagation and structural dynamics in graphene nanoribbons[END_REF]. In this work, we describe the behavior of the natural frequencies and modeshapes of single layer graphene sheets with various boundary conditions using the MM method, and compare the findings with continuum mechanics based on isotropic properties. Kirschoff-based plate formulations consider the graphene sheets as an isotropic continuum, while edge effects and the finite size of the sheets have been demonstrated to provide in-plane special orthotropic properties [START_REF] Huang | Thickness of graphene and single wall carbon nanotubes[END_REF][START_REF] Scarpa | Effective elastic mechanical properties of single layer graphene sheets[END_REF][START_REF] Reddy C D, Rajendran | Equilibrium configuration and elastic properties of finite graphene[END_REF][START_REF] Rajendran | Determination of Elastic Properties of Graphene and Carbon-Nanotubes Using Brenner Potential: The Maximum Attainable Numerical Precision[END_REF].

The main novelty of our paper is the molecular mechanics [START_REF] Adhikari | The calibration of carbon nanotube based bio-nano sensors[END_REF][START_REF] Chowdhury | Vibrational analysis of ZnO nanotubes: A molecular mechanics approach[END_REF], which is a higher fidelity model compared to the previously used models [START_REF] Hashemnia | Vibrational analysis of carbon nanotubes and graphene sheets using molecular structural mechanics approach[END_REF][START_REF] Sakhaee-Pour | Vibrational analysis of singlelayered graphene sheets[END_REF] where the structural mechanics (SM) approach is used. In the MM approach, the molecule first finds the lowest energy configuration from its pre-optimized state and then the subsequent calculations are performed. Finding the minimum potential energy surface is a key distinguishing feature in the MM approach. It can be noted that during the optimization the structures may change their configurations. Planar structure such as SLGS may deform to attain the lowest energy configuration. This is ignored in the SM approach [START_REF] Hashemnia | Vibrational analysis of carbon nanotubes and graphene sheets using molecular structural mechanics approach[END_REF][START_REF] Sakhaee-Pour | Vibrational analysis of singlelayered graphene sheets[END_REF], which only considers the initial geometry of the molecule. In the SM approach the Euler-Bernoulli beam model is used to represent the C-C bonds and subsequently the finite element method [START_REF] Bathe | Finite Element Procedures[END_REF] is used to discretise the equation of motion and obtain the natural frequencies. In the formulation presented in ref [START_REF] Hashemnia | Vibrational analysis of carbon nanotubes and graphene sheets using molecular structural mechanics approach[END_REF][START_REF] Sakhaee-Pour | Vibrational analysis of singlelayered graphene sheets[END_REF], the length and the diameter of the 'beams' representing the C-C bonds are almost similar (length = 0.142nm and diameter = 0.146nm). It is well known [START_REF] Timoshenko | Theory of Plates and Shells[END_REF] that for such 'short beams' the Euler-Bernoulli beam model used in [START_REF] Hashemnia | Vibrational analysis of carbon nanotubes and graphene sheets using molecular structural mechanics approach[END_REF][START_REF] Sakhaee-Pour | Vibrational analysis of singlelayered graphene sheets[END_REF] produce inaccurate results as the shear deformation is ignored. This can introduce error particularly for the calculation of higher frequencies. The results obtained using the MM approach do not have these drawbacks as neither the Euler-Bernoulli beam model, nor the finite element approach, is used. Contrary to SM approaches, the molecular mechanics method also allows the natural frequencies of the system in equilibrium to be computed without any requirement to use an equivalent thickness value for the nanostructure. In general, the introduction of the thickness concept in nanomaterials is highly contentious, leading to the well-known "Yakobson's paradox" [START_REF] Dibiasio | Difference between bending and stretching moduli of single-walled carbon nanotubes that are modeled as an elastic tube[END_REF][START_REF] Huang | Thickness of graphene and single-wall carbon nanotubes[END_REF], responsible for the high scattering of Young's moduli and Poisson's ratio results in open literature.

We will show that the MM approach is able to capture the equivalent anisotropic properties and their influence in the structural dynamics of the graphene sheets, and therefore provide a valid prediction tool to simulate the resonance behavior of graphenebased NEMS devices. The paper is organized in the following way. The continuum mechanics approach for the frequency analysis of graphene sheets is presented in section 2. Section 3 will be centered on the analysis and calculation of the frequencies using a molecular mechanics model. The numerical results and discussion will be presented in section 4. Finally, the major conclusions of this paper will be drawn in section 5 based on the results and analyses in section 4.

Vibration of single layer graphene

A single layer graphene sheet (shown in Fig. 1) may be approximated by a thin elastic plate [START_REF] Scarpa | The bending of single layer graphene sheets: Lattice versus continuum approach[END_REF]. The equation of motion of the transverse free vibration of a thin elastic plate [START_REF] Timoshenko | Theory of Plates and Shells[END_REF][START_REF] Soedel | Vibrations of Shells and Plates[END_REF] can be expressed as A rectangular single layer graphene sheet (SLGS) and its mathematical idealization using a thin continuum plate. The SLGS is assumed to be of dimension a × b.

D ( ∂ 4 w ∂x 4 + 2 ∂ 2 w ∂x 2 ∂ 2 w ∂y 2 + ∂ 4 w ∂y 4 ) + ρ ∂ 2 w ∂t 2 = 0 (1)
Here w ≡ w(x, y, t) is the transverse deflection, x, y are coordinates, t is the time, ρ is the mass density per area and the bending rigidity is defined by

D = Eh 3 12(1 -ν 2 ) (2)
E is the Young's modulus, h is the thickness and ν is the Poisson's ratio. We consider rectangular graphene sheets with cantilevered (clamped at one edge) and bridged (clamped at two opposite edges) boundary conditions. Following Blevins [START_REF] Blevins | Formulas for Natural Frequency and Mode Shape[END_REF], the natural frequency (in rad/s) of a rectangular plate of dimension a × b can be expressed as

ω ij = { π 4 D a 4 ρ } 1/2 { G 4 x + G 4 y ( a b ) 4 +2 ( a b ) 2 [νH x H y + (1 -ν)J x J y ] } 1/2 (3)
where i, j = 1, 2, 3, . . . are mode indices. The values of the coefficients G x , H x , J x and G y , H y , J y depend on the boundary conditions and the mode indices i, j. The first set of coefficients depends on the boundary conditions of the edges of width (side b) while the second set of coefficients depends on the boundary conditions of the edges of length (side a). The boundary conditions on the two edges of length (side a) are free. In this paper we consider the lower modes of vibration. For the first three modes the coefficients G y , H y and J y are given in Table 1. The coefficients G x , H x and J x for both boundary 2. General expressions of the coefficients for the higher values of i and j are given in [START_REF] Blevins | Formulas for Natural Frequency and Mode Shape[END_REF]. The values given in Tables 1 and 2 will be used to obtain the natural frequencies and compare with the molecular mechanics simulation described in the next section.

Molecular simulation approach

Since atomic configurations can have significant impact on the mechanical properties of single layer graphene sheets, zigzag and armchair models are adopted in this study. The zigzag and armchair models of graphene sheets under consideration are:

• Zigzag sheet clamped at one edge (Cantilevered condition)

• Zigzag sheet clamped at two opposite edges (Bridged condition)

• Armchair sheet clamped at one edge (Cantilevered condition)

• Armchair sheet clamped at two opposite edges (Bridged condition)

Different atomic configurations and boundary conditions can be considered in an unified manner within the scope of molecular mechanics. The general expression of total energy is a sum of energies due to valence or bonded interactions and non-bonded interactions [START_REF] Rappe | a full periodictable force-field for molecular mechanics and molecular-dynamics simulations[END_REF]

E = N B ∑ 0 1 2 k IJ (r -r IJ ) 2 + N A ∑ 0 k IJK (C 0 + C 1 cos θ + C 2 cos 2θ) + N T ∑ 0 1 2 V ϕ (1 -cos (nϕ 0 ) cos (nϕ)) + N I ∑ 0 V ω ( C I 0 + C I 1 cos ω + C I 2 cos 2ω ) + N nb ∑ 0 R IJ [ -2 ( x IJ x ) 6 + ( x IJ x ) 12 ] + N nb ∑ 0 q I .q J ε.x (4)
N B , N A , N T , N I and N nb are the numbers of the bond-, angle-, torsion-, inversion-and the non bonded-terms, respectively. k IJ and k IJK are the force constants of the bondand angle-terms, respectively. r and r IJ are the bond distance and natural bond distance of the two atoms I and J, respectively. θ and θ 0 are the angle and natural angle for three atoms I -J -K, respectively. ϕ and ϕ 0 are the torsion angle and torsion natural angle for three atoms I -J -K -L, respectively. V ϕ , n, V ω , ω are the height of the torsion barrier, periodicity of the torsion potential, height of the inversion barrier and inversion-or out-of-plane-angle at atom I, respectively. C I 0 , C I 1 and C I 2 are the Fourier coefficients of the inversions terms. x and x IJ are the distance and natural distance of two non bonded atoms I and J. R IJ is the depth of the Lennard-Jones potential. q I and ε are the partial charge of atoms I and dielectric constant. For the general nonlinear case, the bend function should have a minimum θ = θ 0 , with the second derivative at θ 0 equal to the force constant (k IJK ). The Fourier coefficients of the general angle terms C 0 , C 1 and C 2 are evaluated as a function of the natural angle θ 0 :

C 2 = 1 4sin 2 θ 0 C 1 = -4C 2 cos θ 0 C 0 = C 2 ( 2cos 2 θ 0 + 1 ) (5)
The bond stretching force constants (k IJ ) are atom based and are obtained from generalization of Badger's rules. The assumption is that the bonding is dominated by attractive ionic terms plus short-range Pauli repulsions [START_REF] Rappe | a full periodictable force-field for molecular mechanics and molecular-dynamics simulations[END_REF]. The force constant (in units of (kcal/mol)/ Å2 ) then becomes

k IJ = 644.12 Z * I Z * J r 3 IJ (6)
The Z * I is the effective atomic charges, in electron units. Similarly, the angle bend force constants (k IJK ) are generated using the angular generalization of Badger's rule. The force constant (in units of kcal/mol.rad 2 ) then becomes [START_REF] Frisch | Gaussian 09[END_REF]:

k IJK = 644.12 Z * I Z * J r 5 IJ [ 3r IJ r JK ( 1 -cos 2 θ 0 ) -r 2 IK cos θ 0 ] (7)
The torsional constant (kcal/mol) is defined as

V ϕ = 5 √ U I U J [1 + 4.18 ln (BO JK )] (8)
where, BO JK is the bond order for Atom-J and Atom-K. U I and U J are the atomic constants defined with UFF sp 2 . Regarding the inversion term, the coefficients are

C I 0 = 1, C I 1 = -1 and C I 2 = 0 for sp 2 atom type.
In this study, we used the UFF model [START_REF] Rappe | a full periodictable force-field for molecular mechanics and molecular-dynamics simulations[END_REF], wherein the force field parameters are estimated using general rules based only on the element, its hybridization and its connectivity. Hybridization determines the type of bonding of the carbon atoms with its neighbours. The sp 3 hybridization corresponds to the well-known tetrahedral configuration in which carbon binds to four neighbours giving rise to three-dimensional inter-connectivity of carbon atoms that is found in diamond. The sp 2 bonding in which carbon atoms bind to three neighbours also known as trigonal hybridization gives planar structures found in graphite and graphene. The sp 2 -hybridized carbon atoms, which are at their energy minimum in planar graphite (or graphene), must be bent to form the closed sphere (fullerenes) or tube (CNT), which produces angle strain. The characteristic reaction of fullerenes is electrophilic addition at 6,6-double bonds, which reduces angle strain by changing sp 2 -hybridized carbons into sp 3 -hybridized ones. The change in hybridized orbitals causes the bond angles to decrease from about 120 • in the sp 2 orbitals to about 109.5 • in the sp 3 orbitals. This decrease in bond angles allows for the bonds to bend less when closing the sphere or tube, and thus, the molecule becomes more stable. The force field functional forms and parameters used in this study are in accordance with [START_REF] Rappe | a full periodictable force-field for molecular mechanics and molecular-dynamics simulations[END_REF]. The calculation of frequency and their validation for CNTs were detailed in [START_REF] Chowdhury | A molecular mechanics approach for the vibration of single walled carbon nanotubes[END_REF]. In the following section we provide the methodology of the frequency calculation.

Calculation of the natural frequencies

We start with the Hessian matrix f CAR , which holds the second partial derivatives of the potential E with respect to the displacement of the atoms in cartesian coordinates (CAR) [START_REF] Frisch | Gaussian 09[END_REF][START_REF] Mcintosh | The determination of WilsonDecius F matrix elements from Cartesian force constants[END_REF]:

f CAR ij = ( ∂ 2 E ∂ξ i ∂ξ j ) Opt (9)
This is a 3N × 3N matrix (N is the number of atoms), where ξ 1 , ξ 2 , ξ 3 ,. . ., ξ 3N denote the displacements in cartesian coordinates, ∆x 1 , ∆y 1 , ∆z 1 , . . ., ∆z N . The () Opt refers to the fact that the derivatives are taken at the equilibrium positions of the atoms. These force constants are then converted to mass weighted cartesian coordinates (MWC) [START_REF] Mcintosh | The determination of WilsonDecius F matrix elements from Cartesian force constants[END_REF]

f M W C ij = f CAR ij √ m i m j = ( ∂ 2 E ∂c i ∂c j ) Opt ( 10
)
where,

c 1 = √ m 1 ξ 1 = √ m 1 ∆x 1 , c 2 = √ m 1 ξ 2 =
√ m 1 ∆y 1 and so on. f M W C is diagonalized, yielding a set of 3N eigenvectors and 3N eigenvalues. The next step is to translate the center of mass to the origin, and determine the moments and products of inertia, with the goal of finding the matrix that diagonalizes the moment of inertia tensor. Using this matrix we can find the vectors corresponding to the rotations and translations. Once these vectors are known, we know that the rest of the normal modes are vibrations. The center of mass R COM is found in the usual way:

R COM = ∑ α m α r α ∑ α m α ( 11
)
where the sums are over the atoms, α. The origin is then shifted to the center of mass r COMα = r α -R COM . Next we have to calculate the moments of inertia (the diagonal elements) and the products of inertia (off diagonal elements) of the moment of inertia tensor (I).

I =      I xx I xy I xz I yx I yy I yz I zx I zy I zz      =      ∑ m α (y 2 α + z 2 α ) - ∑ m α (x α y α ) - ∑ m α (x α z α ) - ∑ m α (y α x α ) ∑ m α (x 2 α + z 2 α ) - ∑ m α (y α z α ) - ∑ m α (z α x α ) - ∑ m α (z α y α ) ∑ m α (x 2 α + x 2 α )      (12)
The sums appearing in the above expression is over the index α. This symmetric matrix is diagonalized, yielding the principal moments (the eigenvalues I ′ ) and a 3 × 3 matrix (X), which is made up of the normalized eigenvectors of (I). The eigenvectors of the moment of inertia tensor are used to generate the vectors corresponding to translation and infinitesimal rotation of the molecular system. A Schmidt orthogonalization is used to generate N vib = 3N -6 remaining vectors, which are orthogonal to the six rotational and translational vectors. The result is a transformation matrix D which transforms from mass weighted cartesian coordinates q to internal coordinates S = Dq, where rotation and translation have been separated out. Now that we have coordinates in the rotating and translating frame, we need to transform the Hessian, f M W C , to these new internal coordinates (INT) [START_REF] Frisch | Gaussian 09[END_REF][START_REF] Mcintosh | The determination of WilsonDecius F matrix elements from Cartesian force constants[END_REF]. Only the N vib coordinates corresponding to internal coordinates will be diagonalized, although the full 3N coordinates are used to transform the Hessian. The transformation is straightforward as follows:

f IN T = D T f M W C D (13)
The N vib × N vib submatrix of f IN T , which represents the force constants internal coordinates, is diagonalized yielding N vib eigenvalues λ = 4π 2 ω 2 , and N vib eigenvectors.

If we call the transformation matrix composed of the eigenvectors L, then we have

L T f IN T L = Λ ( 14
)
where Λ is the diagonal matrix with eigenvalues λ i . At this point, the eigenvalues need to be converted to frequencies (in Hz) as

ω i = √ λ i 4π 2 (15)

Results and discussions

The resonant frequencies of single layer graphene (SLG)-based resonators depend on the geometric configurations. The atomic structures of SLGS could also exert significant influence on their vibration behaviours. Thus, in this work, we analyze two groups of SLG resonators, i.e., three zigzag SLGS (10,0), (14,0), (18,0) and three armchair SLGS [START_REF] Scarpa | The bending of single layer graphene sheets: Lattice versus continuum approach[END_REF][START_REF] Scarpa | The bending of single layer graphene sheets: Lattice versus continuum approach[END_REF], [START_REF] Rangel | Graphene Terahertz Generators for Molecular Circuits and Sensors[END_REF][START_REF] Rangel | Graphene Terahertz Generators for Molecular Circuits and Sensors[END_REF], [START_REF] Kitipornchai | Continuum model for the vibration of multilayered graphene sheets[END_REF][START_REF] Kitipornchai | Continuum model for the vibration of multilayered graphene sheets[END_REF], with increasing length. In this study, we computed our results using bridged (atoms at the two sides along width are restrained) and cantilevered (atoms at one side along width are restrained) boundary conditions. The computational results of the first four vibrational frequencies of these zigzag SLGS are calculated and presented in Fig. 2 and Fig. 3, respectively, for bridged and cantilevered boundary condition. Similarly, Fig. 4 and Fig. 5, respectively, presents the first four vibrational frequencies of these armchair SLGS, for bridged and cantilevered boundary condition. The widths of the SLGS are given in the figure captions.

Dependence of the length

As shown in Fig. 2 and Fig. 3, for SLGS with the length rising from around 20 Å to 120 Å, the fundamental frequencies are in the ranges of 100-3000 GHz (ref. 4) for the zigzag SLGS with bridged and cantilevered boundary conditions, respectively. While for the armchair SLGS, the variation of frequencies is between 40-1070 GHz (ref. Table 5) and 2-415 GHz (ref. Table 6) for bridged and cantilevered boundary conditions, respectively, with the length rising from around 40 Å to 210 Å. The trend of the frequency changes with length are generally in accordance with that given in the literature [START_REF] Hashemnia | Vibrational analysis of carbon nanotubes and graphene sheets using molecular structural mechanics approach[END_REF]. The discrepancy is primarily a result of the different end constraints, geometric configurations of SLGS and the differences in the simulation approaches. Recall that here the molecular mechanics approach is used as opposed to the finite element method employed in [START_REF] Hashemnia | Vibrational analysis of carbon nanotubes and graphene sheets using molecular structural mechanics approach[END_REF][START_REF] Sakhaee-Pour | Vibrational analysis of singlelayered graphene sheets[END_REF].

4-1300 GHz (ref. Table

Dependence of the geometric configuration

For both zigzag and armchair SLGS, the frequencies of all five modes generally decrease with increasing length. The curves of frequency becomes steeper for the SLGS of smaller sheet length (≤ 50 Å). This shows that the dependence on the aspect ratio is stronger for the frequencies of shorter SLGS. In the meantime, it is also seen from Fig. 2 and Fig. 4 that, for a given length the frequencies of SLGS always decline with rising sheet width. This issue can be further clarified from Tables 3 to 6. The frequencies of small-width SLGS are always higher than the corresponding frequencies of large-width SLGS. This is especially so when the length is relatively small. However, the effect of width on the frequencies diminishes for SLGS with larger widths. As an example, for zigzag SLGS (ref. 3. Thus we see that when the aspect ratio of SLGS grows, the difference in frequency due to the variation of length decreases significantly whereas the ratio between the frequencies remains almost unchanged. This observation suggests that the influence of width on the vibration frequency of SLGS does not significantly change with increasing length.

Here the decreasing frequencies with increasing length and width observed in Fig. 2, and Fig. 4 can be attributed to the fact that SLGS of larger length and width possess lower dynamic structural stiffness in both longitudinal and transverse directions. In particular, for SLGS of small width their transverse stiffness is high. The frequency of such SLGS thus becomes more sensitive to their longitudinal rigidity, which finally leads to stronger effect of the length for SLGS with smaller width.

Based on the molecular mechanics method, mode shapes of the SLGS are obtained. The first six mode shapes of zigzag and armchair sheets are demonstrated in the supplementary document (see supplementary materials). The first mode shape plays a significant role in the design of the nanomechanical resonators. It is perceived that the SLGS with different boundary conditions has a sinusoidal and/or hyperbolic sine and cosine configuration. These configurations guarantee the ease of detection of any small deflection in the SLGS. In addition, mode shapes of the SLGS, in contrast to the CNTs, are not changed by the length or aspect ratio.

Dependence of the atomic structure

Next we examine the effect of atomic structures on the frequency of SLGS. To this end we consider the calculated frequencies in Table 5 and Table 3 respectively, for armchair and zigzag SLGS with the bridged boundary condition. Similarly, we consider the calculated frequencies in Table 6, and Table 4 respectively, for armchair and zigzag SLGS with the cantilevered boundary condition. The results in Fig. 6 show that the chirality does not have a significant influence on the natural frequencies of vibration. It is shown that for almost same width and length the fundamental frequencies of zigzag SLGS are only higher than those of armchair SLGS for bridged case, whereas for cantilevered case, it is almost comparable. In the context of single wall carbon nanotubes (SWCNTs), the difference between the frequencies of the two types (e.g., zigzag and armchair) are also not very large [START_REF] Chowdhury | A molecular mechanics approach for the vibration of single walled carbon nanotubes[END_REF]. This may be expected as SWCNTs are effectively rolled up SLGS. The frequency of SLGS is primarily determined by their geometry, i.e., length and the aspect ratio, and cannot be substantially changed by varying their atomic structure. This finding demonstrates that the continuum models can produce a good approximation for the vibration of SLGS with different atomic structures.

Comparison with the continuum theory

In this section we investigate whether the vibrational frequencies obtained from the simple continuum plate model are comparable with the vibrational frequencies computed using the molecular mechanics approach. The density per unit area on the SLGS is computed from the total mass divided by the total area on the SLGS. We use the first natural frequency to obtain the bending rigidity D using the natural frequency equation ( 3) with the values of the coefficients corresponding to i = 1 and j = 1 in Tables 1 and2. We then use this value of D to compute the higher natural frequencies to understand if the simple plate model is applicable. The first natural frequency is not shown in the figures because the constant D is calculated using the first natural frequency.

In Fig. 7, the second and third bending mode frequencies are compared for the bridged boundary condition. We used (10,0) zigzag SLGS as an example. The results for the other types of SLGS used in this study show similar behavior. The analytical results shown in Fig. 7 are obtained from equation [START_REF] Novoselov | Two-dimensional gas of massless Dirac fermions in graphene[END_REF] with the values of the coefficients 1.

The trend in the variation of the frequencies with respect to the length is similar for both the methods. The difference between the two theories is more prominent for SLGS with smaller dimension. This is expected as the continuum theory may not be very suitable for SLGS in this case. The continuum model tends to overestimate the MM predictions (17 % and 37 % for the second and third mode respectively). The continuum mechanics formulation assumes an isotropic equivalent material for the graphene. However, edge effects have been demonstrated to play a significant role in the static [START_REF] Scarpa | Effective elastic mechanical properties of single layer graphene sheets[END_REF][START_REF] Reddy C D, Rajendran | Equilibrium configuration and elastic properties of finite graphene[END_REF][START_REF] Rajendran | Determination of Elastic Properties of Graphene and Carbon-Nanotubes Using Brenner Potential: The Maximum Attainable Numerical Precision[END_REF][START_REF] Reddy | Edge elastic properties of defect-free single-layer graphene sheets[END_REF] and dynamic [START_REF] Scarpa | Wave propagation and structural dynamics in graphene nanoribbons[END_REF] mechanical properties of SLGS, leading to an equivalent orthotropic, rather than isotropic material model for the graphene. Although the dimensions of the SLGS considered in this work are dissimilar, we notice a general agreement in terms of magnitude between the eigenvalues calculated with our MM approach, and the results using the MM3 potential in [START_REF] Gupta | Elastic properties and frequencies of single-layer graphene sheets[END_REF]. Figure 8 shows equivalent plots for the SLGS with cantilevered boundary condition. The analytical results shown in Fig. 8 are obtained from equation ( 3) with the values of the coefficients corresponding to i = 2, 3 for the clamped-free case in Table 2 and j = 1 in Table 1. The results predicted by the two approaches agree more closely for this boundary condition. The results obtained here shows that the boundary condition has an effect on the accuracy of the predictions from the continuum theory.

Conclusions

The vibrational properties zigzag and armchair single wall graphene sheets (SLGS) are studied. A molecular mechanics based approach is used to estimate the frequencies. We used the UFF model, wherein the force field parameters are estimated using general rules based only on the element, its hybridization, and its connectivity. Two types of boundary conditions are considered, namely, cantilevered and bridged. First five Natural frequencies are calculated for three zigzag, namely, (10,0), (14,0) and (18,0) and three armchair, namely, [START_REF] Scarpa | The bending of single layer graphene sheets: Lattice versus continuum approach[END_REF][START_REF] Scarpa | The bending of single layer graphene sheets: Lattice versus continuum approach[END_REF], [START_REF] Rangel | Graphene Terahertz Generators for Molecular Circuits and Sensors[END_REF][START_REF] Rangel | Graphene Terahertz Generators for Molecular Circuits and Sensors[END_REF] and [START_REF] Kitipornchai | Continuum model for the vibration of multilayered graphene sheets[END_REF][START_REF] Kitipornchai | Continuum model for the vibration of multilayered graphene sheets[END_REF] SLGS. The natural frequencies of SLGS decrease with length but they are generally insensitive to the atomic structure. Results obtained from the molecular mechanics approach are compared with the same obtained using the continuum plate theory. The continuum mechanics results in general tend to overestimate the natural frequencies. The two approach agree more for the cantilevered boundary condition compared to the bridged boundary condition. The results obtained in the paper may be useful for the design and analysis of vibrating SLGS based NEMS and sensor devices. Dependence of the molecular structure on the fundamental frequency of SLGS. It is found that, natural frequencies of zigzag SLGS (18,0) are higher compared to armchair SLGS [START_REF] Rangel | Graphene Terahertz Generators for Molecular Circuits and Sensors[END_REF][START_REF] Rangel | Graphene Terahertz Generators for Molecular Circuits and Sensors[END_REF] for bridged case, whereas for cantilevered case, it is almost comparable. Trust for the award of the Philip Leverhulme Prize.

  Continuum plate model for graphene.

Fig. 1 .

 1 Fig. 1.A rectangular single layer graphene sheet (SLGS) and its mathematical idealization using a thin continuum plate. The SLGS is assumed to be of dimension a × b.

Fig. 2 .

 2 Fig. 2. Bridged boundary condition -First four vibrational frequencies of zigzag SLGS as a function of the length of SLGS. The widths are -(10,0): 9.317 Å; (14,0): 13.554 Å; (18,0): 17.803 Å.

Fig. 3 .

 3 Fig. 3. Cantilevered boundary condition -First four vibrational frequencies of zigzag SLGS as a function of the length of SLGS. The widths are -(10,0): 9.317 Å; (14,0): 13.554 Å; (18,0): 17.803 Å.

Fig. 4 .

 4 Fig. 4. Bridged boundary condition -First four vibrational frequencies of armchair SLGS as a function of the length of SLGS. The widths are -(11,11): 12.31 Å; (15,15): 17.23 Å; (19,19): 122.15 Å.

Fig. 5 .

 5 Fig. 5.Cantilevered boundary condition -First four vibrational frequencies of armchair SLGS as a function of the length of SLGS. The widths are -[START_REF] Scarpa | The bending of single layer graphene sheets: Lattice versus continuum approach[END_REF][START_REF] Scarpa | The bending of single layer graphene sheets: Lattice versus continuum approach[END_REF]: 12.31 Å; (15,15): 17.23 Å; (19,19): 122.15 Å.

  Fig. 6.Dependence of the molecular structure on the fundamental frequency of SLGS. It is found that, natural frequencies of zigzag SLGS (18,0) are higher compared to armchair SLGS[START_REF] Rangel | Graphene Terahertz Generators for Molecular Circuits and Sensors[END_REF][START_REF] Rangel | Graphene Terahertz Generators for Molecular Circuits and Sensors[END_REF] for bridged case, whereas for cantilevered case, it is almost comparable.

  Frequency of the second bending mode.

  Frequency of the third bending mode.

Fig. 7 .

 7 Fig. 7. Bridged boundary condition -vibrational frequencies of the zigzag SLGS (10,0) for the second and third bending mode. Molecular mechanics results and continuum mechanics results are compared for different values of the length.

  Frequency of the second bending mode.

  Frequency of the third bending mode.

Fig. 8 .

 8 Fig. 8.Cantilevered boundary condition -vibrational frequencies of the zigzag SLGS (10,0) for the second and third bending mode. Molecular mechanics results and continuum mechanics results are compared for different values of the length.

Table 1 .

 1 Coefficients for the free-free boundary conditions of the edges of length (side

	a)
		Mode index (j) G y H y	J y
	.	1 2	0 0	0 0 1.216 0
		3	1.506 1.248 5.017
	conditions on the edges of width (side b) are given in Table

Table 2 .

 2 Coefficients for the two boundary conditions on the edges of width (side b)

			Clamped-clamped	Clamped-free
		Mode index (i) G x H x	J x	G x	H x	J x
	.	1	1.506 1.248 1.248 0.597 -0.0870 0.471
		2	2.5 4.658 4.658 1.494 1.347 3.284
		3	3.5 10.02 10.02 2.5 4.658 7.842

Table 3

 3

) and

Table 3 )

 3 with width 9.317 Å the frequency decreases form 784.10 to 101.88 GHz when the length increases from 23.391 Å to 121.822 Å, while for certain length of 23.391 Å it only varies from 784.10 GHz to 752.30 GHz, when width increases from 9.317 Å to 17.803 Å. Similar behavior is also observed for armchair SLGS in Table

Table 3 .

 3 Vibrational frequencies of zigzag SLGS in GHz -Bridged boundary condition.

	Index and	Length ( Å) ω 1	ω 2	ω 3	ω 4	ω 5
	Width ( Å)
		23.391 784.10 1178.43 1842.16 2511.38 2960.40
		47.995 293.30 504.48 623.61 967.17 1006.04
	(10,0)	72.603 180.15 321.84 370.86 474.31 580.16
	9.317	97.213 130.12 236.47 264.56 285.31 406.86
		121.822 101.88 186.91 194.26 205.92 313.98
		23.391 763.00 1014.94 1799.25 1850.99 2220.43
		47.995 281.72 420.08 600.04 861.22 983.66
	(14,0)	72.603 172.47 265.86 355.04 536.70 558.29
	13.554	97.213 124.42 194.62 252.85 361.47 387.81
		121.822 97.34 153.53 196.67 244.08 299.77
		23.391 752.30 929.62 1430.21 1777.36 2059.21
		47.995 275.09 374.56 587.23 772.02 897.41
	(18,0)	72.603 168.02 235.48 345.94 477.20 544.22
	17.803	97.213 121.12 171.99 246.01 346.27 379.20
		121.822 94.72 135.45 191.26 272.03 291.82

Table 4 .

 4 Vibrational frequencies of zigzag SLGS in GHz -Cantilevered boundary condition.

	Index and	Length ( Å) ω 1	ω 2	ω 3	ω 4	ω 5
	Width ( Å)
		23.391 76.71 437.63 483.34 659.82 1323.45
		47.995 19.18 115.03 166.98 206.05 321.91
	(10,0)	72.603	8.58 50.46 73.69 134.18 140.80
	9.317	97.213	5.20 28.57 41.25 78.61 99.77
		121.822 3.75 18.70 26.34 50.53 79.32
		23.391 75.35 339.71 471.07 854.05 1035.23
		47.995 18.41 109.56 158.49 230.04 311.10
	(14,0)	72.603	8.34 47.75 102.90 103.10 133.99
	13.554	97.213	5.84 27.91 58.08 74.40 77.05
		121.822 3.74 17.52 37.11 47.57 60.62
		23.391 75.00 279.65 467.10 857.12 883.47
		47.995 17.80 104.88 130.68 287.79 305.18
	(18,0)	72.603	8.15 45.04 84.69 128.84 131.38
	17.803	97.213	4.68 24.89 62.71 70.16 74.43
		121.822 3.67 16.34 44.52 47.76 49.86

Table 5 .

 5 Vibrational frequencies of armchair SLGS in GHz -Bridged boundary condition.

	Index and	Length ( Å) ω 1	ω 2	ω 3	ω 4	ω 5
	Width ( Å)
		41.21	282.84 411.22 628.37 858.54 1072.30
		83.83	117.86 185.65 243.91 375.70 385.14
	(11,11)	126.46	74.47 120.09 151.31 197.16 232.77
	12.31	169.08	63.05 96.48 135.33 170.34 176.12
		211.71	42.94 70.46 80.55 86.31 130.72
		41.21	287.54 359.84 636.08 765.41 953.01
		83.83	120.58 158.87 248.98 323.45 392.11
	(15,15)	126.46	76.32 102.08 154.90 205.84 237.93
	17.23	169.08	55.85 75.30 112.58 151.24 151.67
		211.71	43.43 59.43 87.37 101.74 119.02
		41.21	290.01 335.32 640.24 681.45 720.46
		83.83	122.14 145.19 251.80 297.89 395.86
	(19,19)	126.46	77.40 93.09 156.92 188.19 240.78
	22.15	169.08	55.91 68.18 112.83 136.68 171.43
		211.71	44.70 54.26 89.80 108.84 123.21

corresponding to i = 2, 3 for the clamped-clamped case in Table

2

and j = 1 in Table

Table 6 .

 6 Vibrational frequencies of armchair SLGS in GHz -Cantilevered boundary condition.

	Index and	Length ( Å) ω 1	ω 2	ω 3	ω 4	ω 5
	Width ( Å)
		41.21	24.18 149.23 156.26 267.23 416.07
		83.83	6.72 37.14 67.35 74.23 102.89
	(11,11)	126.46	3.63 17.02 29.88 46.14 48.69
	12.31	169.08	2.50 10.08 16.82 26.42 36.25
		211.71	2.04 7.06 10.81 17.54 28.85
		41.21	24.06 118.14 149.10 353.65 387.94
		83.83	6.82 37.24 55.23 92.43 102.87
	(15,15)	126.46	3.62 16.98 36.00 41.29 46.08
	17.23	169.08	2.50 10.02 23.26 26.35 26.76
		211.71	1.29 6.03 14.86 16.47 21.08
		41.21	23.96 96.16 148.95 330.70 413.55
		83.83	6.73 37.20 43.58 102.76 116.43
	(19,19)	126.46	3.59 16.90 28.47 45.93 52.51
	22.15	169.08	1.91 9.04 21.27 25.46 29.60
		211.71	1.69 6.82 16.53 17.27 19.03

Acknowledgements

RC acknowledges the support of Royal Society through the award of Newton International Fellowship. SA gratefully acknowledges the support of The Leverhulme