Predicting Graphene's Nonlinear‐Optical Refractive Response for Propagating Pulses

Nonlinear‐optical refraction is typically described by means of perturbation theory near the material's equilibrium state. Graphene, however, can easily move far away from its equilibrium state upon optical pumping, yielding strong nonlinear responses that cannot be modeled as mere perturbations. So far, one is still lacking the required theoretical expressions to make predictions for these complex nonlinear effects and to account for their evolution in time and space. Here, this long‐standing issue is solved by the derivation of population‐recipe‐based expressions for graphene's nonperturbative nonlinearities. The presented framework successfully predicts and explains the various nonlinearity magnitudes and signs observed for graphene over the past decade, while also being compatible with the nonlinear pulse propagation formalism commonly used for waveguides.


Introduction
Over the past decade there has been a rapidly growing interest in the theoretical [1,2] and experimental [3][4][5][6][7][8][9][10][11][12][13] investigation of nonlinear-optical refraction in graphene using free-space and waveguided excitation configurations. This new research field was pioneered by the experiments of Hendry et al. showing an exceptionally large effective graphene nonlinearity | (3) gr | ≈ 10 −7 esu [3] or, equivalently, an effective nonlinear index |n 2,gr | ≈ 10 −13 m 2 W −1 . The research results reported since then seem to point in different directions, and have made it a major challenge to fully understand graphene's nonlinear-optical behavior. The experimental data include both positive- [4,9] and negativevalued [6][7][8]12] effective nonlinearities with a magnitude compatible with that of Hendry's experiments, [3] as well as much smaller nonlinearities. [10] From the theory point of view, calculations for the perturbative nonlinearity (3) gr [1,2] give rise to nonlinearities DOI: 10.1002/lpor.201900402 that are two orders of magnitude smaller than | (3) gr | ≈ 10 −7 esu measured in the aforementioned experiments. Strictly speaking, using perturbation theory implies that the graphene remains close to its initial equilibrium state, but this is not necessarily the case. [14] Most of the experiments carried out so far were conducted using an exfoliated or chemical-vapor-deposition (CVD)-grown graphene sample without intentional doping. As such, when using optical excitation wavelengths, onephoton absorption (1PA) in the graphene layer is unavoidable, and the resulting free-carrier generation can give rise to nonlinearities outside the perturbative regime with strong changes not only in temperature [11,15] but also in chemical potential. [13] These changes will be time-and spacedependent when dealing with propagating excitation pulses. What is more, our recent investigations on the evolution of the spectral bandwidth of pulses propagating in graphene-covered waveguides [12] have shown the need for a nonperturbative treatment of graphene's nonlinear-optical refraction, with saturation playing an important role in the carrier dynamics. For telecom excitation wavelengths, we observed that this "saturable photoexcited carrier refraction" yields a strong, negative nonlinearity from a phenomenological point of view, [12] but we did not yet present a mathematical description for the refraction efficiency as a function of more fundamental properties of graphene.
Microscopic models used for calculating the ultrafast carrier dynamics and the resulting nonlinear-optical response of graphene are generally built upon the equation of motion for the density matrix. [14] Over the years, these models have been applied for describing the nonlinear-optical physics in both large-area graphene and graphene nanoribbons. [1,2,[16][17][18][19][20][21][22][23] Although manyparticle interactions such as carrier-carrier and carrier-phonon scattering can be rigorously included in this framework [16][17][18][19] and provide new insights in, for example, time-resolved optical pump-probe experiments, [17,18] the complexity and computational intensiveness of such modeling motivate a simpler treatment of this microscopic physics based on the relaxationtime approximation. Indeed, such simplification was employed in several works: on one hand, in, for example, refs. [2,20], where graphene's perturbative optical nonlinearity was calculated by solving the semiconductor Bloch equations [2] or by solving the equation of motion for the density matrix with the electric field described by means of a scalar potential [20] ; on the other hand, in, for example, refs. [21][22][23], where graphene's nonperturbative www.advancedsciencenews.com www.lpr-journal.org response was calculated within a steady-state density-matrix formalism [21] or by means of the generalized Bloch equations for massless Dirac fermions, [22,23] originally introduced in ref. [24]. In microscopic theories like these, one calculates the surface current in the graphene induced by a given time-dependent electric field going through the infinitesimally thin graphene sheet. However, when studying graphene's nonperturbative response for pulses propagating along the graphene sheet, yielding electric fields evolving in time and space as in, for example, graphene-covered waveguides, one needs to simultaneously solve the time-and space-dependent Maxwell equations. Hence, when adhering to microscopic calculations as outlined above, it would become necessary then to adopt a Maxwell-Bloch framework which could only be solved by means of complex numerical calculations. Therefore, another, more accessible, approach is needed to link graphene's nonlinear response for propagating pulses to its fundamental properties and to calculate the nonlinear response by means of a closedform expression of the electric field as is commonly done for waveguides. [25] In this letter, we show that graphene's nonlinear refraction efficiency for time-and space-dependent excitations outside the perturbative regime can be predicted using an easily accessible formalism based on the so-called population recipe. [14] More specifically, the refraction efficiency can be obtained from the excitationinduced instantaneous change in both temperature and chemical potential, and from the resulting instantaneous change in the linear conductivity of graphene, in line with the population recipe. We derive an expression that links the conductivity change with the commonly used effective nonlinear index n 2,gr and with the free-carrier refraction (FCR) coefficient FCR introduced in ref. [12]. Finally, we showcase the validity of our formalism for a wide variety of experiments in free space and in waveguides.

Rate Equation and Population Recipe
Semiconductor nonlinearities resulting from the optical excitation of free carriers can often be described by means of the so-called population recipe, where the equilibrium carrier concentrations that enter into the linear conductivity are replaced by the photoexcited carrier concentrations. [14] In this section, we investigate whether this general ansatz for semiconductors could also hold for graphene. Let us assume a singlemode electric field propagating in the z-direction, , where 0 is the mean spectral frequency of the light beam, N is a normalization constant such that |A(z, t)| 2 equals the power, e denotes the transverse distribution of the electric field, and corresponds to the propagation constant. If a current density J is produced due to the interaction between the electric field and a material such as graphene, exhibiting 1PA characterized by its conductivity ij ( 0 ), then the electromagnetic energy converted into free carriers per unit volume and unit time is given by , where the last equation assumes a cycle averaging of the field. [26] In view of this equation, the dependence on space and time coordinates in w can be separated as w(x, y, z, t) = u(z, t)v(x, y). Since u can show complex dependences on A as several carrier-related processes may exist, effective models are often employed to deal with these phenomena. [12,25] Along the same lines, we will study the temporal dynamics of the photoexcited carrier concentration based on the following rate equation for u: where u sat and c are phenomenological parameters that account for saturation and decay mechanisms, respectively. The 1PAinduced carrier concentrations evolving in z and t can then be computed by means of a spatial average (over (x, y)) of w∕(ℏ 0 ), where ℏ is the reduced Planck constant. In the Supporting Information, we show that Equation (1) can be derived from the generalized Bloch equations for massless Dirac fermions [24] on condition that the optical pulse width T 0 is sufficiently long (T 0 ≳ 100 fs) and that the light intensity I 0 satisfies I 0 × 4 0 ≪ 1.8 × 10 −9 W × m 2 [22,23] with 0 = 2 c∕ 0 and c the speed of light in vacuum (e.g., I 0 ≪ 3 × 10 14 Wm −2 for 0 = 1550 nm). The implications of these conditions are further discussed in Section 4. Now we explore how u behaves far from t u = 0 solving Equation (1) by the integrating factor method.
where |A(z, t)| 2 = P 0 U(z, = t∕T 0 ), with P 0 being the input peak power. To facilitate interpretation, a saturation time sat = u sat ∕P 0 and a saturation power P sat = u sat ∕ c have been defined. The last step in Equation (2) assumes T 0 ≫ sat and T 0 ≈ c . [12] Important here is that Equation (2) can be recovered from the steady-state solution of Equation (1) ( t u = 0) if the timeindependent steady-state power |A(z)| 2 is replaced by the instantaneous power |A(z, t)| 2 . Hence, the photoexcited carrier densities in graphene could behave as carrier densities in "instantaneous" steady states. As such, graphene's non-perturbative response could still be well described by means of the steady-state linear conductivity but with the carrier density at equilibrium replaced by the photoexcited densities calculated along Equation (1). This carrier-induced instantaneous change in the linear conductivity then translates into the nonlinear response of the material.

Nonlinear Index and FCR Coefficient
To implement the approach outlined above for modeling graphene's nonlinear response for propagating pulses in the non-perturbative regime, we need to find explicit functional dependences of the chemical potential and temperature T on the carrier concentration. In principle, one should consider Fermi-Dirac distributions with different temperatures [18] and chemical potentials [27] for the electrons and the holes. However, supported by the conclusions of refs. [13,18], we here use a single temperature and chemical potential for both electrons and Laser Photonics Rev. 2020, 14,1900402 www.advancedsciencenews.com www.lpr-journal.org holes. To determine and T, we consider that the probability of occupation of a state with energy E within the time scale of the excitation pulse is given by instantaneous Fermi-Dirac distributions. The density of states in graphene near the Dirac point is approximately gr = 2|E|∕( (ℏv F ) 2 ), where v F ≈ 10 6 m s −1 is the Fermi velocity. The 2D densities of electrons n and holes p are calculated as n( , T) = 2∕( (ℏv F ) 2 ) 2 J 1 ( )∕ 2 and p( , T) = 2∕( (ℏv F ) 2 ) 2 J 1 (− )∕ 2 , where = ∕(k B T), J 1 ( ) is a Fermi-Dirac integral and k B is the Boltzmann constant. [28] We aim at deriving closed expressions for = (n) and T = T(n) taking into account the electroneutrality condition, which relates n and p. To illustrate our approach, we assume an initial hole density p 0 , [7,8,12] so p = p 0 + n.
To obtain explicit solutions = (n) and T = T(n), we use two fitting functions as specified in the Supporting Information and after some algebra derive with = n∕(n + p 0 ), = 2.682 and T = 0.957. Note that for n = 0, T = 0 and = −ℏv F ( p 0 ) 1∕2 are retrieved. Now we turn to the imaginary part of graphene's conductivity at finite temperature, which is given by ref. [2] Im Here, 0 = e 2 ∕4ℏ is the universal conductivity, with −e being the electron charge, and Γ intra and Γ inter are scattering parameters for intraband and interband transitions, respectively. This result accounts for the interband and intraband motion of electrons around the Dirac cone at the single-particle level. [1,2] Equation (5) with and T as defined in Equations (3)  It should be noted that the quantity often measured in nonlinear experiments is the temporal average of the refractive in-dex change ∫ +∞ −∞ Δn(t) I(t) dt∕ ∫ +∞ −∞ I(t) dt where I(t) represents the instantaneous light intensity, and subsequently, the nonlinear index defined via n 2 = ∫ +∞ −∞ Δn(t) I(t)dt∕ ∫ +∞ −∞ I 2 (t)dt. Modeling graphene as a thin layer of 3D material as done in many experiments, [3][4][5][7][8][9][10] we derive the following expression to determine the magnitude and sign of graphene's n 2 (see Supporting Information): where = e 2 ∕(4 0 ℏc) ≈ 1∕137 denotes the fine-structure constant, d gr = 0.33 nm is the effective graphene thickness, n sat is the saturation value of n, and T FWHM indicates the full-width-athalf-maximum of the pulse duration.
In contrast to the works reporting a nonlinear index for graphene, we already analyzed our experiments with graphenecovered waveguides in ref. [12] in terms of a nonlinear carrierrefraction response of the form − FCR N c embedded in the commonly used nonlinear pulse propagation formalism. Here, N c is a 1D carrier density defined at each z-distance along the waveguide. Using ref. [25] as our starting point, we derive an expression for our FCR coefficient: where e y represents the y−component, parallel to the graphene sheet of width w graph , of the waveguide's quasi-transverse electricfield mode,Ñ sat = (v F ∕ 0 )n 1∕2 sat , and n eff denotes the waveguide effective index. We point out that FCR does not depend on d gr nor on the pulse duration, making it a robust measure for graphene's nonlinear response.

Comparison with Experiments in Literature
Most of the experiments to measure graphene's nonlinearity at optical wavelengths have been performed using exfoliated or CVD graphene without intentional doping. As such, for most of those graphene samples the exact p 0 value has not been measured. However, typical ranges of unintentional doping values for exfoliated and CVD graphene can be found in refs. [29,30], respectively. Based on these works we will consider an average doping level of p 0 = 9 × 10 12 cm −2 for exfoliated graphene and an average doping level of p 0 = 4 × 10 12 cm −2 for CVD graphene. To illustrate the behavior of n 2,gr , Figure 1a-d shows the calculated Im(Δ (1) yy ∕ 0 ) versusÑ c = (v F ∕ 0 )n 1∕2 at two commonly used wavelengths 0 = 0.8 µm and 0 = 1.55 µm and at the doping values for exfoliated and CVD graphene. Note that these graphs are robust against changes in the scattering rates Γ intra,inter . According to these results and Equation (7), n 2,gr is expected to switch sign between 0.8 µm and 1.55 µm with positive (negative) values at the shorter (longer) wavelengths (see Supporting Information for a discussion in terms of interband and intraband contributions). This prediction is in full agreement with the www.advancedsciencenews.com www.lpr-journal.org  [2]) as a function of the normalized square-root of the carrier density ranging up to n sat = 10 17 m −2 . [12] experiments reported so far (see Table 1) and shows the tunability of graphene's nonlinearity sign in the near-infrared. Our theory also puts forward the so far undetermined (negative) sign of n 2,gr in Hendry's seminal experiments. [3] Besides predicting correct signs for n 2,gr , Equation (7) also provides orders of magnitude compatible with most n exp 2,gr values in Table 1 (including the highest values), within its one-order-ofmagnitude precision. This correspondence is remarkable as our calculations were carried out with estimated parameter values n sat = 10 17 m −2 and c = 1 ps [12,31] (see Table 1). Equation (7) also recovers the scaling of n 2,gr with pulse duration in the sub-ps excitation regime as experimentally observed in refs. [8,9]. Furthermore, as discussed in the Supporting Information, according to our theory, the nonlinear response is mainly introduced by the dynamical behavior of the temperature rather than chemical po-tential and its switching in sign is mainly produced by changes in the interband contribution.
When verifying again the conditions defining the validity area of our theory (namely T 0 ≳ 100 fs and I 0 × 4 0 ≪ 1.8 × 10 −9 W × m 2 ) we find that these conditions are essentially fulfilled for all works enlisted in Table 1 except for that in ref. [5], where 15fs-long pulses at high optical intensities were used. The latter explains why our theory underestimates n exp 2,gr reported in ref. [5] by several orders of magnitude (see result marked with an asterisk in Table 1). In experiments with very short pulses such as in ref. [5], the intraband carrier thermalization that follows right after the interband excitation might still be ongoing at the end of the pulse duration. Under such circumstances, the intraband contribution to the conductivity change might not be important, in which case a much higher positive nonlinearity value could be produced by theory (see Supporting Information). Further theoretical developments will be required to check this hypothesis. This being said, the accessible theory presented here does yield a good match with all other data in Table 1 observed for ps and sub-ps pulse lengths. This underlines the general applicability of our model for a wide range of excitation conditions. Finally, we want to emphasize that FCR as specified in Equation (8) is a more robust parameter than n 2,gr to quantify graphene's nonlinear response outside the perturbative regime. For the FCR value that we measured in graphene-covered waveguides, [12] we find an excellent agreement with our theoretical prediction (see bottom of Table 1).
In conclusion, with our population-recipe-based theory we can adequately predict graphene's nonlinear refractive response for propagating pulses outside the perturbative regime. With our expressions for n 2,gr and FCR , one can readily calculate how the nonlinear effects evolve over time and distance. Our work explains the differences in sign and magnitude of graphene's nonlinear coefficients measured over the past decade, and also provides the long-needed framework for conceptualizing, designing, and optimizing graphene-enhanced nonlinear-optical (on-chip) devices. Table 1. Comparative table of the magnitudes and signs of graphene's experimental n 2,gr and FCR values and the theoretical values calculated by means of Equation (7) with estimated parameter values n sat = 10 17 m −2 and c = 1 ps [12,31] and Equation (8). Works in refs. [3,5,10] were carried out with exfoliated graphene, and all others with CVD graphene. The symbol ± indicates that the sign of n exp 2,gr could not be determined in those experiments. * See text for a discussion of this case.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.