1
Strain dependent elastic modulus of graphene
Supplementary Information
Authors: Guillermo López-Polín†, Miriam Jaafar†, Francisco Guinea, Rafael Roldán, Cristina
Gómez-Navarro* and Julio Gómez-Herrero
*Correspondence to: [email protected]
Index.
SI1. Substrate fabrication.
SI2. Raman spectroscopy.
SI3. AFM tip characteristics.
SI4. Leak rate of pressurized graphene blisters.
SI5. Finite element numerical simulations.
SI6. Additional data: statistics on membranes.
SI7. Suppression of anharmonicities in strained graphene and its effect on the elastic modulus: A
continuum elastic theory analysis.
2
SI1. Substrate fabrication.
The substrates were fabricated from 6 inch silicon wafers (As-doped, resistivity 1 – 3.5 mΩ·cm). A 300
nm thick silicon oxide was grown by dry oxidation. The back of the wafer was coated with a 0.8 μm
thick layer of AlSi(1%)Cu(0,5%).
The silicon oxide was patterned using projection optical lithography and reactive ion etching (RIE).
Figure S 1 shows a SEM image of one part of the chip after all the processes and a zoom in one of the
holes to appreciate the verticality of the walls.
SI2. Raman spectroscopy.
Raman spectra were performed using a WITEC/ALPHA 300AR Raman confocal microscope in
ambient conditions. The laser wavelength and power were 532nm and 0.7mW respectively.
Graphene monolayers were first identified by optical microscopy and then corroborated by Raman
spectroscopy as described in ref [1] (fig. S3)
Figure S 2: Example of a graphene structure displaying mono- bi- and tri- layer regions.
3
SI3. AFM tip characteristics.
In order to have a constant and well defined contact geometry we use commercial tips from
NanoScience Instruments with hemispherical geometry and low wear coating of Tungsten carbide with
nominal final tip radius of 60 nm (fig. S2). Cantilever spring constants range between 25 and 35 N/m.
Each cantilever was individually calibrated using Sader’s method [2].
Figure S 3: (Scale bar 200 nm). SEM image of the AFM tip used for the measurements.
SI4. Leak rate of pressurized graphene blisters.
As described in the main text strain is induced in graphene membranes by applying a pressure difference
(ΔP) across them. When subjected to these conditions graphene membranes deflect in a spherical cup
geometry as dictated by ΔP. The pressure inside and outside the graphene blisters tends to equal with
time, relaxing the strain of the membrane. This is mainly due to small leakage between the graphene
flake and the substrate [3]. Here we present representative graphs of the dependence of the strain
induced by pressure with time. In our samples, the characteristic relaxation time of the membrane is in
the order of few hundreds of minutes. This is long enough to induce a steady strain during a set of
measurements (~2 minutes) but it is also fast enough to change the strain by waiting for relaxation
instead of changing the chamber pressure.
4
0,6
0
0 (nm)
· Drum 1
· Drum 2
-20
-30
-40
a)
(a)
-50
Strain (%)
0,5
-10
0,4
· Drum 1
· Drum 2
0,3
0,2
b)
(b)
0,1
0,0
-60
0
200
400
600
800
Time (min)
0
200
400
600
800
Time(min)
Figure S 4: Leak rate. a) Evolution of the maximum deflection (δ0) with time for two different
representative drumheads with initial ΔP=4 atm. b) Corresponding geometrical strain vs. time for the
same two drumheads.
SI5. Finite element numerical simulations.
As described in the main manuscript there is not any analytical equation to obtain the elastic modulus of
a pressurized membrane by nano-indentation experiments. Curves acquired at non zero ΔP do not
follow the Schwerin expression (eq.1 in main text) and required a deeper analysis for a correct
interpretation.
A plausible possibility is to consider the transversal deformation caused by ΔP as an initial indentation.
Then, we can reformulate the Schwerin equation as:
F ( )   0    0  
E2 D
   0 3
2
a
If we expand this equation we obtain a full polynomial cubic function as described in eq. 3 in the main
manuscript. Interestingly, the coefficient in 3 still remains as E2D/a2 as described by Lin et al.[4]
According to this expression we could tentatively fit our experimental curves obtained on pressurized
membranes to a full polynomial cubic function where the elastic modulus of the membrane would be
given by the third order coefficient.
In order to check this assumption we performed finite element numerical simulations to model
indentations on membranes under pressure. To this end we used Comsol Multiphysics 4.3, modeling a
conventional circular clamped membrane (in this context conventional means without entropic
anharmonic effects) under hydrostatic pressure, with radius a. The elastic modulus of the membrane was
pre-fixed to 340 N/m as an input of the simulations. The topographic profiles obtained from the
numerical simulations agreed with the profiles obtained by imaging drumheads under pressure previous
to indentations. As illustrated in Figure S5, the in plane strains yielded by numerical simulation also
fitted those estimated according to equation 2 in the main text.
5
a)
0
20
40
0.27
60
80
0.26
100
0.25
120
140
0.24
160
R 
COMSOL
1st principal strain (%)
a)
R
R-0
a
0
S
b)
b)
(nm)
0.28
-5000
-800 -600 -400 -200
-1000
200 400 600 800
0
Radius (nm)
X(nm)
500
1000
1400
400
E2D
(N/m)
E2d
(N/m)
E2D E2d
(N/m)
(N/m)
Figure S 5: Strain in pressurized membranes. Panel a) is a sketch of the membrane geometry under
d) numerical simulations across the center of
ΔP. Panel bc)
is the in plane strain calculated by finite1200
element
380 membrane with 1.5 μm diameter and a ΔP=3 atm. Following the scheme shown in figure a)
a circular
1000
δ= 0 in good agreement
COMSOL yields a strain of 0.25%
with figure b).
δ= δ 0
360
δ=0
800
The indentation force exerted by the tip is modeled by inducing a high pressure Pic within an inner circle
with radius ain (a/ain=20) at the center of the drumhead,
600 as sketched in figure S6. ain was calculated
340 [5]. The force exerted by the indenter is calculated as the pressure multiplied by the surface
following
area of the inner circle. This model is more realistic400
than a point load and it is easier to converge for
finite 320
element calculations. The topographic profiles200
of a representative pressurized membrane under
4
3
2
1
0
4 S6.
3 in figure
1 forces 2are shown
0
different indentation
Pressure [atm]
Pressure [atm]
a)
b)
ain
calculated E2D
1200
a
c)
F=c1d+c3d3
1000
800
 0
600
F=c0+c1d+c2d2+c3d3
400
200
0
1
2
P(atm)
3
4
Figure S 6: Finite element modeling of a membrane under pressure. a) Cartoon of a circular
membrane. The red inner area illustrates the region where the force is exerted by the AFM tip. b)
Topographic profiles obtained in Comsol Multiphysics across the center of a pressurized membranes
under different indentation forces. The upper profile corresponds to the deformation of the membrane
just by hydrostatic pressure with a pressure difference of 4 atm. c) E2D obtained from the coefficient C3
in F()= c0+c1+c22+c33(green) and F()= c1+c33 (red) for ΔP from 0 to 4 atm. The predefined E2D
for the finite elements calculation is 340 N/m.
6
In order to test our suggested analysis, the F () curves generated by the finite element analysis were
fitted both to equation 1 and to a full 3rd order polynomial (eq. 3 in main text). For both cases the elastic
modulus is calculated from the third order coefficient as c3= E2D/a2. The values obtained from such
fittings for membranes under different pressures are displayed in figure S6c. Here we can appreciate that
while the coefficient of the complete third order polynomial yields the correct (pre-fixed in simulations)
value of E2D even in the case of high pressure difference, fitting to a Schwering type equation is only
applicable in the case of zero pressure difference. Based in these simulations, in what follows we will
obtain E2D from c3= E2D/a2 using a complete third order polynomial fitting, in agreement with ref. [4].
In order to check the robustness of the results obtained by finite element simulations and with the aim of
comparing with our experimental data, the simulated indentation curves where also fitted to a complete
third order polynomial function up to different loading forces. Results are displayed in Figure S7 panel
a. In this plot we observe that the numerically simulated indentation curves yield robust values of elastic
modulus above ~500 nN, where the curves converge to the pre-fixed EM of 340 N/m. In figure S7b we
represent three experimental curves under the same analysis. Here we appreciate that the curves
acquired at lower ΔP agree with those obtained from numerical simulation, further supporting the
applicability of eq.3 in the main text to analyze our experimental results. In contrast, the experimental
curve acquired at ΔP=3.5 atm shows a marked different behavior; instead of converging to the ~300
N/m value, it presents a growing and eventual saturation tendency at a value of ~700 N/m. This result
indicates a change in the elastic constant of the membrane for the situation where ΔP=3.5 atm.
1200
a)
P= 0 atm
P= 1 atm
P= 2 atm
P= 3 atm
P= 4 atm
E2D (N/m)
1000
800
600
400
200
1400
1200
b)
P= 0 atm
P~ 1.5 atm
P~ 3.5 atm
1000
E2D (N/m)
1400
800
600
400
200
0
0
0
300
600
900
1200
Indentation force (nN)
1500
0
300
600
900
1200
1500
Indentation force (nN)
Figure S 7: Elastic modulus as a function of pressure difference across the membrane. a) Elastic
modulus estimated from c3 when simulated curves are fitted up to different indentation forces. This
fitting provides robust values for forces higher than 500 nN. The over estimation observed at low forces
is reflect that considering the transversal deformation caused by ΔP as an initial indentation fails for low
indentations. b) Same analysis than in panel a), now performed on experiential curves.
7
SI6. Additional data: statistics on membranes.
Similar results than those displayed in figure 3 in the main text were obtained in more than 10 graphene
drumheads. Some of them are shown in figure S8.
a)
b)
Figure S 8: Elastic modulus as a function of total prestrain for several graphene drumheads.
8
SI7. Suppression of anharmonicities in strained graphene and its effect on the
elastic modulus: a continuum elastic theory analysis
In this Section we use the statistical mechanics of flexible membranes, beyond the harmonic
regime, to study the effect of an external strain in the sample. Our main focus here is to model
the effects of tension on the elastic (Young) modulus of 2D membranes, giving a theoretical
analysis that can explain qualitatively the experimental observations shown in the main text.
The flat phase of a 2D membrane as graphene at sufficiently long scales, much larger than the
interatomic separation, is well described by a free energy that is a sum of a bending and a
stretching part [6]
Z
h
i
2
1
F[u, h] =
d2 r κ ∇2 h + 2µε2αβ + λε2αα
(1)
2
where κ is the bending rigidity, λ and µ are the first Lam´e constant and the shear modulus,
respectively [7], and εαβ is the internal strain tensor expanded to the lowest order in the in-plane
displacements u(r) and the out-of-plane displacements h(r),
1
εαβ ≈ (∂α uβ + ∂β uα + ∂α h∂β h).
2
(2)
The most commonly used harmonic approximation neglects the interaction between the bending
and stretching modes, and therefore the u(r) and h(r) fields are decoupled. If an external strain
is applied to the membrane, its effect can be modeled by the inclusion of a new term in the free
ext
energy (1), that couple to the internal strain tensor, ταβ εαβ , where ταβ = λδαβ uext
αβ + 2µuαβ is
expressed in terms of the external strain tensor uext
αβ .
We start by considering the harmonic correlation function for the out-of-plane (flexural)
modes, G0 (q) = h|h(q)|2 i0 , where the subscript 0 indicates that we neglect anharmonic couplings between in-plane and out-of-plane modes. In the presence of an isotropic expansion, the
external strain tensor takes the form uext
¯δαβ , where u¯ = δS/2S accounts for the uniform
αβ = u
dilation of the membrane, S is its surface, and δS is the change in area due to the application
of strain. It is straightforward to show that G0 (q) can be expressed as [8]:
G0 (q) =
q 2 [κq 2
kB T
+ 2(λ + µ)¯
u]
(3)
where kB is the Boltzman constant and T is the temperature. It is illustrative to notice that
a purely harmonic theory predicts, in the absence of strain, a divergence in the mean-square
amplitude
displacements hh2 i ∝ L2 , where L is the size of the sample and
P of the out-of-plane
2
2
hh i = q h|h(q)| i. This is the well known result that harmonic theory of membranes predicts not a flat but a crumpled membrane. In particular, the harmonic theory predicts that
graphene should not be stable. This contradiction is solved by including the anharmonic coupling between bending and stretching modes in the calculation. Anharmonic effects have been
considered from different approaches. Among them, the self-consistent screening approxima-
9
tion (SCSA) is especially useful because it is relatively simple, and their results agree well with
more sophisticated methods, like atomistic Monte Carlo simulations [9, 10]. The exact correlation function obtained from the SCSA is well approximated as an interpolation between the
high-q and the low-q limits, by considering the effective self-energy Σa (q) that accounts for the
anharmonic interaction between bending and stretching modes:
Geff (q) =
κq 4
kB T
+ 2(λ + µ)¯
uq 2 + kB T Σa (q)
where [10]
Σa (q) = A(T )q
4
q0
q
(4)
η
(5)
where ηpis a characteristic exponent, which within the SCSA takes the value η ≈ 0.82 [9],
0
/κ, in terms of the bending rigidity κ and the non-interacting Young modulus
q0 = 2π E2D
0
=
E2D
4µ(µ + λ)
2µ + λ
(6)
which we assume to take the value measured by the present experiments when the graphene
0
∼ 700N/m ∼ 44 eV˚
A−2 , as observed by the plateau in Fig.
membrane is highly stretched, E2D
3a) of the main text. On the other hand, A(T ) is a numerical factor whose value is estimated to
be of the order of A(T ) ≈ 5.9 T [K]η/2−1 , where T [K] is the temperature expressed in Kelvins.
As it is well understood [6], anharmonicities modify the elastic constants of the theory and
make them scale dependent, for example κR (q) ∼ q −η , and λR (q), µR (q) ∼ q 2+ηu , where the
anomalous exponents η and ηu are related by ηu = 2 − 2η [9]. As a consequence, the Young
modulus is scale dependent as well, and one can define an effective Young modulus (considering
anharmonicities) as [11]
0
E2D
E2D (k) =
(7)
0
1 + E2D
b(k)
where the interaction between the h(r) and u(r) fields is considered using a Feynman-diagram
technique, leading to
Z
1
d2 q ˆ
b(k) =
|k × q|4 Geff (q)Geff (k − q).
(8)
2
2
(2π)
It is useful to take into account that
Z
ˆ×q
Γ 1+
ˆ |4
d2 k m |k
3 Γ 1 − n+m
2 I(m, n) =
k
=
n+m
2
4−n
ˆ|
(2π)
|k − q
16π Γ 2 + 2 Γ 2 −
m
Γ
2
m
Γ
2
1 + n2
2 − n2
(9)
Based on the above description, we have calculated the Young modulus at room temperature
10
as a function of wave-vector for different values of strain and the Young modulus as a function
of strain for different wave-vectors, and the results are shown in Fig. 9. The figures have been
obtained numerically from Eq. (7) and (8). Depending on the characteristic wavelength of the
considered mode, we can distinguish three different regimes, which are indicated in Fig. 9 by
Roman characters I-III:
• I. Long wavelengths: Strain dominated regime.
We can define a long wavelength region for very small wave-vectors for which k qs ,
where qs is defined from the solution of [8]
η q0
κ + kB T A(T )
qs2 = 2u(λ + µ).
(10)
qs
In this regime, the Young modulus is dominated by the strain term and can be approximated by
0
E2D
u
E2D
(k) ≈
(11)
2 .
kB T
3
0
2
1 + E2D 64π 2(λ+µ)¯u qs
• II. Intermediate wavelengths: Anharmonic regime.
This regime can be defined for qs k qc , where qc is a characteristic wave-vector that
determines the onset of anharmonic effects, which can be estimated perturbatively from
the Ginzburg criterium [8]
r
0
3kB T E2D
.
(12)
qc =
16πκ2
In this regime the Young modulus can be approximated by
anh
E2D
(k) ≈
1
k 2−2η
= A2 (T )q02η 1 ,
b(k)
I
2 η
(13)
where Iη = I(η, η) can be determined from Eq. (9).
• III. Short wavelengths: Harmonic regime.
Finally, for qc k we are in the harmonic dominated region, for which the Young
modulus is not affected by interactions between bending and stretching modes, and it is
well approximated by
har
0
(k) ≈ E2D
.
(14)
E2D
The different regimes described in the previous section can be easily identified in Fig. 9(a),
in which we observe that for a given value of applied strain u, the Young modulus presents
three different behaviors. First we have one region, for small wave-vectors, in which E2D is
almost independent on the wave-vector, taking an approximate value given by Eq. (11). This
long wavelength region, as expected [8], is dominated by the effect of strain. Then, there is a
11
��
��
III
�∼�� %
�%
��
�� -� %
II
�� -� %
��
(�)
�� -� %
I
��
�� -� %
��
-�
�[Å ]
��
��
��
� � � [�/�]
� � � [�/�]
� = � × �� -� Å-�
��
��
� × �� -� Å-�
��
� × �� -� Å-�
��
� × �� -� Å-�
��
� × �� -� Å-�
� × �� -� Å-�
� × �� -� Å-�
(�)
��
��
��
� × �� -� Å-�
��
��
�[%]
Figure S 9: (a) Young modulus as a function of wave-vector E2D (q) for different values of applied
strain. (b) Young modulus as a function of applied strain, for different wave-vectors. The dashed
0 . The dashed vertical line indicates the crossover wave-vector q
black horizontal line is the bare E2D
c
below which anharmonic effects are relevant, as given by (12).
transition to the anharmonic region of the spectrum, in which the Young modulus is highly
dependent on the wave-vector k, following a power law E2D ∼ k 2−2η as given by Eq. (13).
Finally, for even larger wave-vectors, we enter the region where the anharmonic effects are not
important (regardless the value of applied strain), and the Young modulus is well reproduced
0
by its harmonic value E2D
(shown by the dashed black line in the figure).
The effect of strain for a given wave-vector is more clearly seen in Fig. 9(b) in which we
represent the strain dependence of E2D for different wave-vectors. First we observe that, for
a given wave-vector, the Young modulus has a plateau at low values of strain. This is due to
the fact that the strain is still too weak that cannot suppress the anharmonic effects at that
wavelength. If we keep increasing the strain, we reach a threshold for which anharmonic effects
are being suppressed, and the Young modulus is affected by the external strain. This growing
0
of E2D with u happens until we reach the non-interacting value of the Young modulus E2D
,
which is fully 2D. Notice that this behavior compare qualitatively well with the experimental
results shown in Fig. 2 of the main text.
In order to compare more quantitatively to experiments, in Fig. 10 we use three wavevectors associated to the expected characteristic wavelengths which are relevant in graphene
(between l ∼ 10nm − 1µm). Again, we observe for those wave-lengths the same qualitative
behavior as compared to Fig. 2 of the experiments. Notice that, in agreement with previous
experiments [12], our theoretical results quantitatively agree better for a wavelength of the
order of l ∼ 100 nm. Indeed, the Young modulus for shorter wavelengths like l ∼ 10 nm
is almost in the harmonic regime, so that only a small difference is observed between the
anharmonic and harmonic regimes. On the other hand, for long wavelengths like l ∼ 1µm, the
transition between anharmonic and harmonic regimes, driven by the strain, is rather abrupt
with a factor ×3 difference between the anharmonic and harmonic regimes. However, according
to the present theory the strain at which the anharmonicities are suppressed for l ∼ 1µm is
12
��
��
� � � [�/�]
� � � [�/�]
��
��
��
��
��
��
� ~ ��
(�)
��
��
��
��
��
��
��
� ~ ��
(�)
��
�[%]
��
��
��
��
��
��
�[%]
� � � [�/�]
��
��
��
��
� ~ � μ�
(�)
��
��
��
��
��
�[%]
Figure S 10: Young modulus as a function of applied strain, for three wave-lengths l. We show the
0 and the dotted black line
results for l ∼ 10nm, l ∼ 100nm and l ∼ 1µm. The dashed black line is E2D
anh
is E2D (k) from Eq. (13).
too small, compared to the experimental results. Therefore, we the present results based
on a continuum theory of elastic membranes qualitatively capture the effects observed in the
experiments. A more quantitative theoretical description will require of more sophisticated
numerical methods, as atomistic Monte Carlo simulations.
13
REFERENCES
[1] Sader, J. E., Chon, J. W. M. & Mulvaney, P. Calibration of rectangular atomic force microscope
cantilevers. Review of Scientific Instruments 70, 3967-3969, (1999).
[2] Ferrari, A. C. et al. Raman spectrum of graphene and graphene layers. Physical Review Letters 97,
187401, (2006).
[3] Lin, Q.-Y.; Jing, G.; Zhou, Y.-B.; Wang, Y.-F.; Meng, J.; Bie, Y.-Q.; Yu, D.-P.; Liao, Z.-M. Acs
Nano, 7 (2), 1171-1177, (2013).
[4] Koenig, S. P., Wang, L. D., Pellegrino, J. & Bunch, J. S. Selective molecular sieving through porous
graphene. Nature Nanotechnology 7, 728-732, (2012).
[5] Begley, M.R. & Mackin, T.J. Spherical indentation of freestanding circular thin films in the
membrane regime. Journal of the Mechanics and Physics of Solids 52, 2005-2023 (2004).
[6] Nelson, D. R., Piran, T. and Weinberg, S., Statistical Mechanics of Membranes and Surfaces. World
Scientific (Singapore, 2004)
[7] Typical parameters for graphene at room temperature are κ ≈ 1.1eV, λ ≈ 2.4eVÅ-2 and μ ≈9.95eVÅ-2
[8] Roldán, R, Fasolino, A., Zakharchenko, K. V. and Katsnelson, M. I., Suppression of anharmonicities
in crystalline membranes by external strain. Phys. Rev. B 83, 174104(2011)
[9] Le Doussal, P. and Radzihovsky, L. Self-consistent theory of polymerized membranes. Physical
Review Letters. 69, 1209 (1992).
[10] Zakharchenko, K.V., Roldán, R., Fasolino, A., Katsnelson, M.I., Phys. Rev. B 82, 125435 (2010)
[11] Katsnelson, M. I. Graphene: carbon in two dimensions. Cambridge University Press (2012)
[12] López-Polín, G., Gómez-Navarro, C., Parente, V., Guinea, F., Katsnelson, M. I., Pérez-Murano, F.
and Gómez-Herrero, J. Increasing the elastic modulus of graphene by controlled defect creation. Nature
Physics 11, 26 (2015)