December 1, 2011 / Vol. 36, No. 23 / OPTICS LETTERS

4569

Full-band quantum-dynamical theory of saturation and four-wave mixing in graphene Zheshen Zhang1,2 and Paul L. Voss1,2,* 1

2

Georgia Tech Lorraine, Georgia Tech-C.N.R.S. UMI 2958, 2 rue Marconi, Metz, France School of Electrical and Computer Engineering, Georgia Institute of Technology, 777 Atlantic Drive NW Atlanta, Georgia 30318, USA *Corresponding author: [email protected] Received August 5, 2011; revised October 17, 2011; accepted October 18, 2011; posted October 18, 2011 (Doc. ID 152431); published November 25, 2011 The linear and nonlinear optical response of graphene are studied within a quantum-mechanical, full-band, steady-state density-matrix model. This nonpurtabative method predicts the saturatable absorption and saturable four-wave mixing of graphene. The model includes τ1 and τ2 time constants that denote carrier relaxation and quantum decoherence, respectively. Fits to existing experimental data yield τ2 < 1 fs due to carrier–carrier scattering. τ1 is found to be on the timescale from 250 fs to 550 fs, showing agreement with experimental data obtained by differential transmission measurements. © 2011 Optical Society of America OCIS codes: 190.4380, 190.4720.

Graphene possesses unusual optical properties. Over a broad range of wavelengths, measurement and theory result in 2.3% absorption per atomic layer [1]. This uniform absorption is due to its zero bandgap and linear energy dispersion near the Dirac points. Nonlinear optical properties have also been studied. Graphene was shown to be an excellent saturable absorber with very low saturation threshold [2,3]. Studies of higher-harmonic generation [4,5] and four-wave mixing (FWM) [6,7], show theoretically and experimentally, that graphene possesses extraordinarily high nonlinearity, eight orders of magnitude greater than silica for a comparable interaction length. This extremely large nonlinearity can potentially be used in conjuction with micro- and nanolithography techniques to achieve novel nonlinear-optical devices for applications such as frequency conversion and χ ð3Þ -based quantum entanglement generation [8]. From a theoretical point of view, however, the physical model of the nonlinearity has not yet been completed. The saturation effect of graphene is investigated in [2] by a semiclassical method, resulting in a decay constant much lower than that measured in differential transmission (DT) experiments [9–15]. The nonlinear response of graphene up to the THz frequency range was calculated in [16]. [6] theoretically studied FWM in graphene based on a dipole-interaction Hamiltonian used for calculating the optical response of typical semiconductors. We have recently proposed an analytical quantumdynamical FWM perturbation theory for graphene based on the minimal substitution [7]. This theory uses phenomenological time constants τ1 and τ2 to approximate carrier relaxation (CR) and quantum decoherence (QD) in graphene. This Letter, on the other hand, studies saturation effects by numerically solving the density-matrix equations for strong continuous wave excitation where perturbation theory fails. This is important due to graphene’s low saturation power. This Letter differs also in including full-band calculations without resorting to the ideal Dirac fermion assumption. We evaluate the linear optical response, saturable absorption, and saturable FWM as a function of the time constants. The results agree broadly with the τ1 , τ2 , and saturation intensity from femtosecond experiments. 0146-9592/11/234569-03$15.00/0

By use of the tight-binding model with nearestneighbor interaction approximation, the Hamilton for an electron with momentum p in graphene is written in the first quantized language as [4]
H0 ¼

0 h
p

 hp ; 0

ð1Þ

ap ap where hp ¼ −η½expði pffiffi3yℏÞ þ 2 cosðak2x Þ expð−i 2 pffiffi3yℏÞ, η ≈ 2:7 eV is the hopping energy for electrons in graphene, and a ¼ 3:3Å is the lattice constant of graphene. Given the momentum p of an electron, the Hamiltonian in Eq. (1) gives two energy eigenstates jC p i and jV p i with eigenenergies jhp j, corresponding to the conductionband and valence-band states. When an electron is exposed to an external x-polarized electromagnetic field described by its vector potential ~ AðtÞ, the Hamiltonian of the electron can be obtained ~ by the minimal substitution p → p þ eAðtÞ on the freeelectron Hamiltonian in Eq. (1). We then expand the obtained Hamiltonian around p and separate it into the sum of the free-electron Hamiltonian H 0 and a timevarying part V ðtÞ as HðtÞ ¼ H 0 þ V ðtÞ, where V ðtÞ ¼ ~ To study both the linear and nonlineare∇p H 0 · AðtÞ. optical response, in particular pump-degenerate FWM, ~ ¼ Ap e−iωp t þ As e−iωs t þ c:c: ¼ of graphene, we let AðtÞ P~ j Aj ðtÞ, where Ap is the amplitude of the pump mode and As is the amplitude of the signal mode. Once electrons are excited into the conductionband from the valence-band by an external field, they undergo ultrafast relaxation processes, caused by a fast electron–electron scattering process followed by a slower electron–phonon relaxation process [9–15]. These processes play an important role for the quantumdynamics of electrons in graphene since they not only change the state-occupation probability, but also destroy quantum coherence. A complete treatment of the electron-relaxation processes requires complicated fully quantum-mechanical calculations [17]. However, since we restrict ourselves to steady-state response where the state-occupation probability does not change in time, it is

© 2011 Optical Society of America

4570

OPTICS LETTERS / Vol. 36, No. 23 / December 1, 2011

possible to simply introduce two phenomenological decay constants Γ1 ¼ 1=τ1 and Γ2 ¼ 1=τ2 , denoting the CR rate and QD rate. The quantum-dynamics of electrons in graphene are derived within the density-matrix formalism: ϱ_ ¼ −Γ1 ðϱ − ϱeq Þ þ

2ie X ~ A ðtÞV VC ðρVC þ ρCV Þ; ℏ j j

ρ_ VC ¼ ðiωCV − Γ2 ÞρVC −

ie X ~ A ðtÞð2V VV ρVC þ V VC ϱÞ; ℏ j j

ρ_ CV ¼ ð−iωCV − Γ2 ÞρCV ie X ~ A ðtÞð−2V VV ρCV þ V VC ϱÞ; − ℏ j j

ð2Þ

P~ where hmp jV ðtÞjnp i ¼ e j A j ðtÞV mn with m; n ∈ fV ; Cg. This approach is valid for pulses longer than the relaxation timescale, chemical potential close to zero, and ℏωpðsÞ ≫ kT. In Eq. (2), ϱ ¼ ρCC − ρVV is the population inversion, ρVC is the quantum coherence, and ρeq is the population inversion in thermal equilibrium. ωCV is the bandgap frequency. ϱ, ρVC , and ρCV are composed of several frequency components with frequencies mωp þ nωs . P We let ρðtÞ ¼ m;n ρðm;nÞ ðtÞe−iðmωp þnωs Þt þ c:c:, where ρ ∈ fϱ; ρVC ; ρCV g. Together with Eq. (2) and by use of the steaðm;nÞ ðm;nÞ dy-state condition ϱ_ ðm;nÞ ðtÞ ¼ ρ_ CV ðtÞ ¼ ρ_ VC ðtÞ ¼ 0, we obtain the solution for the density-matrix of a single electron with momentum p. To evaluate the optical conductivity, we first derive the expectation velocity of an electron on polarization x, which is related to the velocity operator by hvx i ¼ Trðvx ρÞ, where the velocity operator vx ¼ ∂H=∂px . The surface-current density of graphene is then obtained by integrating hvx i over the whole Brillouin zone: Z 1 ~ JðtÞ ¼ 2 −ehvx idkx dky : 4π BZ X ¼ J ðm;nÞ e−iðmωp þnωs Þt þ c:c: ð3Þ m;n

We define σ ¼ 2 × J ð1;0Þ =E p as the linear conductivity composed of both a real part and an imaginary part, where the factor two accounts for the spin degeneracy.

The real part contributes to absorption and the imaginary part results in a phase shift. In Fig. 1, we fit our theory to the experimental transmittance data obtained in [1], giving τ2 ≈ 0:64 fs. In Fig. 1 the wavelength dependence of the real part is plotted in the left insert and the imaginary part is plotted in the right insert. They are compared to the analytical solution in [7]. When the wavelength decreases to <450 nm, the optical conductivity based on the full-band calculation starts to diverge from the analytical solution, resulting in part from the fact that at lower wavelengths, the linear energy dispersion relation becomes invalid. The imaginary part of σ mostly depends on the QD rate τ2 [7]. Given the complex refractive index data in [18–20] at 550 nm, extracted values of τ2 range from 0:25 fs to 0:72 fs, showing agreement with the previous absorption fitting. The τ2 value of ∼1 fs for graphene is shorter than that of similar materials such as carbon nanotubes (τ2 ∼ 10 fs to 4 ps) [21], perhaps due to the rippled surface of graphene and interactions with substrate or other graphene layers. τ1 for both graphene and other materials can be extracted from DT experiments. The τ1 value of ∼100 fs for carbon nanotubes remains comparable to that of graphene. With QD time less than 0:5 fs, τ1 begins to decrease because, with short QD time, the optical linear conductivity becomes significantly different from the universal value σ 0 ¼ e2 =4ℏ. The saturation threshold at 800 nm is reported to be 4  1 GW=cm2 [2]. In the insert of Fig. 2, the real part of σ is plotted with different decay time constants versus intensity. We also compare the saturation curve obtained by quantum calculations adopted in this Letter with classical saturation curve given by the nonsaturating optical linear conductivity multiplied by intensity-dependent coefficient 1=ð1 þ I=I th Þ. At saturation threshold intensity I th , σ is decreased by half. Knowledge of I th and τ2 allow determination of the CR rate τ1 . In Fig. 3(a) we plot the relation between τ1 and τ2 given saturation powers of 3, 4, and 5 GW=cm2 . With τ2 experimentally obtained from [18–20] and I th obtained by [2], the resulting value of τ1 ranges from 250 fs to 550 fs, within the range of previously reported DT experimental data [9–15]. The reported experimental parameters in [6] give a peak power density of about 2 GW=cm2 , suitable for use of [7] or this Letter. The QD and CR rate seem to differ from sample to sample. Factors that have impact on the two

2

x 10

5

1 1

0.5

0

1

1

Re( )/

(3) (3) / eff silica

1.5

= 250 fs, = 125 fs,

2

= 1 fs, Ith = 1.53 GW/cm 2

2

= 0.5 fs, Ith = 5.27 GW/cm 2

2

= 0.25 fs, Ith = 11.8 GW/cm 2

1 0.8 0.6 0.4 0.2 10

10

= 500 fs,

6

6

7

8

9

10

10 10 10 10 10 Light intensity (W/cm2)

10

7

10

8

11

10

9

10

10

10

11

2

Light intensity (W/cm )

Fig. 1. (Color online) (Outer) Fit to the experimental data in [1], yielding τ2 ¼ 0:64 fs. (Inner) Optical linear conductivity for different τ1 and τ2 . Solid lines: full-band calculation. Hollow squares: the analytical solution in [7].

ð3Þ

Fig. 2. (Color online) Saturation of χ eff . Pump at 775 nm and signal at 1000 nm. (Inset) Saturable absorption at 800 nm. Solid lines: from quantum calculations in this Letter. ×-marks: fit to classical saturation curve.

December 1, 2011 / Vol. 36, No. 23 / OPTICS LETTERS 600

7

χ(3) /χ(3) eff silica

1

400 300

4571

detuning and at large detunings. We have also proposed the use of FWM to probe CR and QD rates.

6

500

τ (fs)

x 10

7

5

We gratefully acknowledge support from National Science Foundation (NSF) grant DMR-0820382.

4 3 2

200

1 0.5

1

1.5 τ2 (fs) 2

Ith = 3 GW/cm

2

Ith = 4 GW/cm

2

Ith = 5 GW/cm

2.5 2

−4

−2 0 2 Detuning (THz)

4

τ1 = 500 fs, τ2 = 1 fs

τ1 = 125 fs, τ2 = 0.25 fs

τ1 = 250 fs, τ2 = 0.5 fs

Standard Lorentzian

Fig. 3. (Color online) (a) Dependence of τ1 on τ2 (at 800 nm). ð3Þ (b) χ eff compared to silica. Solid lines: full-band calculation. Hollow squares: the analytical solution in [7]. ×-marks: a Lorentzian fit to the top blue curve.

time constants may also be due to differences in substrate interactions, temperatures, impurities, and/or excitation wavelengths. However, we do not expect significant variations in τ2 in different samples given the similarity in complex refractive indices of two different samples measured by the same experimental setup [19]. On the other hand, τ1 ranges from tens of fs to ps timescale in the DT experimental measurements [9–15]. The FWM conductivity at ωi ¼ 2ωp − ωs is given by σ FWM ¼ 2 × 83 × J ð2;−1Þ =E 2p E
s . We define the effective thirdð3Þ order susceptibility of graphene as χ eff ¼ jσ FWM j=ϵ0 ωi d ð3Þ [6]. In Fig. 3(b), we plot the ratio between χ eff with ð3Þ the third-order susceptibility of silica χ silica ¼1:84 × 10−22 m2 =V2 . The full-band calculation is compared with the analytical solution in [7], confirming its validity under low-excitation power, and with a Lorentzian fit that illustrates the qualitative difference between graphene and a two-level atom, notably graphene’s stronger nonlinearity at large detunings. Since different time constants lead to significantly different nonlinearity, FWM is a good probe of the time constants. Longer time constants imply less impurity scattering, longer ballistic times, and higher nonlinearity. With increasing pump power, the FWM surfacecurrent density also saturates near 4 GW=cm2 . In Fig. 2, saturation of FWM conductivity is plotted. Such knowledge will be useful for future microfabricated FWM devices using graphene for its large nonlinearity. In conclusion, we have examined FWM and saturable absorption in graphene and described their dependence on CR and QD rates and the strength of FWM at zero

References 1. R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, Science 320, 1308 (2008). 2. G. C. Xing, H. C. Guo, X. H. Zhang, T. C. Sum, and C. H. A. Huan, Opt. Express 18, 4564 (2010). 3. D. Popa, Z. Sun, F. Torrisi, T. Hasan F. Wang, and A. C. Ferrari, Appl. Phys. Lett. 97, 203106 (2010). 4. B. Rosenstein, M. Lewkowicz, H. C. Kao, and Y. Korniyenko, Phys. Rev. B 81, 041416(R) (2010). 5. K. L. Ishikawa, Phys. Rev. B 82, 201402(R) (2010). 6. E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, and S. A. Mikhailov, Phys. Rev. Lett. 105, 097401 (2010). 7. Z. S. Zhang and P. L. Voss, arXiv:1106.4838v1 [cond-mat. mes-hall]. 8. X. Y. Li, P. L. Voss, J. E. Sharping, and P. Kumar, Phys. Rev. Lett. 94, 053601 (2005). 9. D. Sun, Z-K. Wu, C. Divin, X. B. Li, C. Berger, W. A. de Heer, P. N. First, and T. B. Norris, Phys. Rev. Lett. 101, 157402 (2008). 10. J. M. Dawlaty, S. Shivaraman, M. Chandrashekhar, F. Rana, and M. G. Spencer, Appl. Phys. Lett. 92, 042116 (2008). 11. X. Zhao, Z.-B. Liu, W.-B. Yan, Y. P. Wu, X.-L. Zhang, Y. S. Chen, and J.-G. Tian, Appl. Phys. Lett. 98, 121905 (2011). 12. M. Breusing, C. Ropers, and T. Elsaesser, Phys. Rev. Lett. 102, 086809 (2009). 13. M. Breusing, S. Kuehn, T. Winzer, E. Malić, F. Milde, N. Severin, J. P. Rabe, C. Ropers, A. Knorr, and T. Elsaesser, Phys. Rev. B 83, 153410 (2011). 14. P. J. Hale, S. M. Hornett, J. Moger, D. W. Horsell, and E. Hendry, Phys. Rev. B 83, 121404(R) (2011). 15. F. Carbone, G. Aubock, A. Cannizzo, F. van Mourik, R. R. Nair, A. K. Geim, K. S. Novoselov, M. Chergui, Chem. Phys. Lett. 504, 37 (2011). 16. S. A. Mikhailov and K. Ziegler, J. Phys. Condens. Matter 20, 384204 (2008). 17. M. Kira and S. W. Koch, Prog. Quantum Electron. 30, 155 (2006). 18. M. Bruna and S. Borini, Appl. Phys. Lett. 94, 031901 (2009). 19. X. F. Wang, Y. P. Chen, and D. D. Nolte, Opt. Express 16, 22105 (2008). 20. Z. H. Ni, H. M. Wang, J. Kasim, H. M. Fan, T. Yu, Y. H. Wu, Y. P. Feng, and Z. X. Shen, Nano Lett. 7, 2758 (2007). 21. S. Roche, J. Jiang, F. Triozon, and R. Saito, Phys. Rev. Lett. 95, 076803 (2005).