Unconventional double-bended saturation of carrier occupation in optically excited graphene due to many-particle interactions

Saturation of carrier occupation in optically excited materials is a well-established phenomenon. However, so far, the observed saturation effects have always occurred in the strong-excitation regime and have been explained by Pauli blocking of the optically filled quantum states. On the basis of microscopic theory combined with ultrafast pump-probe experiments, we reveal a new low-intensity saturation regime in graphene that is purely based on many-particle scattering and not Pauli blocking. This results in an unconventional double-bended saturation behaviour: both bendings separately follow the standard saturation model exhibiting two saturation fluences; however, the corresponding fluences differ by three orders of magnitude and have different physical origin. Our results demonstrate that this new and unexpected behaviour can be ascribed to an interplay between time-dependent many-particle scattering and phase-space filling effects.

S aturation of light absorption is a long-known central phenomenon in nonlinear optics. It results from the fermionic character of optically driven electrons exhibiting Pauli blocking in the excited states. As a result, with increasing light fluence, the probability to increase the carrier occupation in the excited state due to light absorption is reduced. The saturation behaviour can be exploited in optical devices for the modulation of ultrashort laser pulses. In particular, mono-or multilayer graphene can be used as passive mode-locking medium [1][2][3][4][5][6] . Here the graphene sample is placed into the laser cavity 1,3 or into the output coupler 5,6 , which affects the relative amplification of the frequency with highest intensity. As a result, the modulated pulse becomes sharpened in bandwidth. The advantage of graphene as a saturable absorber is a relatively low saturation intensity and an extremely short recovery time (reflecting the ultrafast carrier dynamics [7][8][9] ) enabling high repetition rates 1 . Furthermore, the broadband absorption of graphene allows for the modulation at different wavelengths, where the intensity can be controlled through the number of graphene layers 1,2 .
The widely used approach for modelling the saturation behaviour is the two-level saturable absorber model [10][11][12] . It describes the saturation of the carrier occupation after continuous wave excitation accounting for the basic interplay of excitation strength, recombination, dephasing and the Pauli principle. The latter governs the high-intensity regime and explains the observed saturation behaviour. In the low-intensity regime, the Pauli blocking is negligible and the carrier occupation is proportional to the light intensity. In a strict sense, the standard saturation model can only be applied for two-level systems and continuous wave excitations. Therefore, for solid-state absorbers, such as graphene, exhibiting many electronic degrees of freedom and for ultrashort pulses produced by mode-locked solid-state lasers, the described scenario can only be of limited value.
In this work, we present a joint theory-experiment study on the saturation behaviour of the maximal optically induced carrier occupation in graphene. We apply a sophisticated theoretical approach based on graphene Bloch equations 8,13,14 , which significantly go beyond the standard saturation model. The approach includes optical pumping and Coulomb-and phononinduced many-particle carrier scattering processes on the same microscopic footing. We provide a microscopic view on the carrier dynamics in optically excited graphene and reveal a new and surprising saturation regime in the weak excitation regime that cannot be explained by the standard Pauli-blocking-based saturable absorber model. This results in a double-bended saturation of the maximal optically induced carrier occupation. The theoretical prediction is experimentally confirmed through differential transmission measurements over a large range of excitation strengths.

Results
Saturation behaviour. The saturation behaviour has been mostly studied by exploiting the two-level saturable absorber model [10][11][12]   (a) Saturation of the carrier occupation r c in a two-level system as a function of the continuous wave intensity I (in units of the saturation fluence I s ). The maximal occupation is limited (dashed line) by Pauli blocking according to the standard saturation model in equation (1). (b) In the low-intensity regime of the standard saturation model the occupation is proportional to the intensity (with slope 1/gG). Theoretical prediction of the saturation of the maximal optically excited carrier occupation r c k0;max in graphene in the (c) high-and (d) low-fluence regime, respectively. The green and blue line represent the extrapolation of the high-and the low-fluence data with equation (1) yielding saturation fluences of F ðhÞ s ¼2 mJ cm À 2 and F ðlÞ s ¼7:8 mJ cm À 2 , respectively (upper abscissa). Experimental data demonstrating the saturation behaviour of the maximal differential transmission DT max /T 0 measured in the (e) high-and (f) low-excitation regime. The saturation fluences (denoted by the vertical lines) are 10 mJ cm À 2 and 7 mJ cm À 2 , respectively. (g-j) Schematic illustration of different saturation regimes in solid-state absorbers: (g) linear regime with low-scattering efficiency. (h) The first surprising saturation already at low intensities can be ascribed to Coulomb-induced out-scattering processes (blue arrows). describing the saturation of the carrier occupation r c after continuous wave excitation with the intensity I. The saturation intensity I s is determined by I s /I ¼ gG/4O 2 , with the Rabi frequency O (p ffiffi I p ), a constant recombination G and dephasing rate g. Obviously, equation (1) accounts for optical excitation (O), carrier recombination (G) and dephasing (g). Here Pauli blocking governs the high-intensity regime I=I s ) 1, where an increase of carrier occupation to 0.5 results in optically induced transparency, see Fig. 1a. In the low-intensity regime, the carrier occupation r c is proportional to the intensity I/gG, where the slope is given by the intrinsic time scales g À 1 and G À 1 of the two-level system (Fig. 1b).
The standard saturation model from equation (1) cannot be applied for solid-state absorbers with many electronic degrees of freedom nor for optical excitations with short pulses. There are fundamental differences in the electron dynamics in solids compared to the simplest saturation model in equation (1): first, for pulsed excitations, the transmission depends on time and, therefore, the maximal transmission at the energy of the optical excitation e 0 and at a specific time t max during an applied pulse is used to characterize the saturation behaviour. This has been applied in several recent studies investigating the transmission saturation behaviour in graphene 2,[5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] . Second, the more nontrivial electronic band structure e l k of solids (here we discuss graphene) makes the decay channels more complex, where carrier-carrier and carrier-phonon-induced in-and outscattering for electronic states k (band l) plays the crucial role rather than radiative recombination. In particular, many-particleinduced scattering channels sensitively depend on the filling of the surrounding phase space that is driven by the intensity of the excitation pulse 8,20 . This results in different regimes characterized by an interplay of time-and fluence-dependent scattering processes and phase-space filling effects.
Here we investigate the saturation behaviour of the maximal optically induced carrier occupation r c k0;max in graphene, which approximately determines the maximal differential transmission DT max /T 0 measured in pump-probe experiments 8 . Here, k 0 describes the optically excited momentum state in the conduction band. We demonstrate a new and surprising saturation regime of the carrier occupation in the excited state. The observed scenario goes beyond the standard Pauli-blockingbased saturable absorber model. We theoretically predict and experimentally confirm double-bended saturation of the maximal optically induced carrier occupation by measuring the differential transmission of a probe pulse. A first illustration of our results is shown in Fig. 1c-f. We find that the Pauli-blocking saturation is valid in the strong-excitation regime and can be fitted by the standard model from equation (1), see Fig. 1c,e for theory and experiment, respectively. For low fluences, however, where the standard model would predict a linear behaviour as shown in Fig. 1b, we find a clearly sub-linear relation between pump fluence and optically induced carrier occupation, see Fig. 1d,f for theory and experiment, respectively. Again, the curvature at low fluences follows saturation-like bending, equation (1), but the corresponding saturation fluence lies three orders of magnitude lower than in the high-fluence regime of Pauli blocking. Note that this behaviour cannot be explained by the standard Pauliblocking saturation, since at the considered low fluences the carrier occupations r c k are below 10 À 2 , suggesting a linear relation between optically induced carrier occupation and pump fluence. In contrast, our calculations reveal that the unconventional saturation-like behaviour can be ascribed to Coulomb-induced out-scattering processes that efficiently balance the pump-pulse-driven occupation in the state k 0 (Fig. 1h). The depopulation due to Coulomb scattering at low intensity occurs at a comparable rate as the optical pumping increases the occupation resulting in a deviation from the initial linear increase of r c;max k 0 with the excitation intensity. This saturationlike behaviour ceases with the increasing importance of Coulomb-induced in-scattering processes (into state k 0 ) that balance out the impact of out-scattering (Fig. 1i) resulting in the regular Pauli-blocking-based saturation behaviour at high pump fluences (Fig. 1j).
Theoretical and experimental approach. Before we analyse the double-bended saturation and its microscopic origin in more detail, we briefly introduce our theoretical and experimental methods. The theoretical approach is based on Heisenberg equations of motion formalism combined with tight-binding wave functions 8,[20][21][22][23] . To accurately model the dynamics of optically excited carriers, we account for the light-carrier interaction as well as carrier-carrier and carrier-phonon scattering on a consistent microscopic footing by solving the many-electron graphene Bloch equations. They constitute a coupled set of differential equations for the occupation probability r l k in the state k in the conduction and the valence band (l ¼ c, v), the microscopic two-band polarization p k and the phonon occupation n j q (not explicitly shown) with the momentum q for different optical and acoustic phonon modes j (refs 7,8,13,21): where e 0 (m 0 ) is the free electron charge (mass), M k the carrier-light coupling strength and A(t) is the vector potential of the applied laser pulse. The in-and outscattering rates G in k;l ðtÞ and G out k;l ðtÞ contain contributions from the carrier-carrier as well as carrier-phonon interaction and depend explicitly on time and momentum, that is, implicitly on the applied pump fluence. They have been calculated microscopically within the second-order Born-Markov approximation. The total dephasing contains non-diagonal U k and diagonal dephasing g k , which for symmetric conduction and valence bands reads g k ðtÞ¼G in k ðtÞ þ G out k ðtÞ. More details on these equations can be found in refs 7,8,13,20,21,24. The modelling is performed for an intrinsic graphene monolayer; nevertheless, the results are still valid for multilayer graphene samples, since the pump fluence variation for different layers is negligibly small compared to the studied pump fluence ranges 25 .
In the following, we consider pulsed excitation and probe at 1.5 eV where the initial thermal carrier occupation can be neglected, so that the saturation of differential transmission DT max /T 0 of the probe pulse is proportional to the maximal occupation r c k;max in the excited state. We briefly note that compared to a previous study, where we have shown that for ultrashort excitation pulses with a duration of just 10 fs, the saturation behaviour of graphene follows the standard model (equation (1)) over eight orders of magnitude in excitation strength 9 , we increase here the excitation pulse duration to B30 fs Using these parameters, we focus in this work on the most interesting regime, where the internal time scales 1=G in k and 1=G out k are on the same order of magnitude as the pulse duration s. In this regime, the competition between in-scattering, outscattering and the optical pump is most important.
Degenerate pump-probe experiments were performed at room temperature on multilayer epitaxial graphene (50 layers) grown by thermal decomposition on the C-terminated face of SiC 26 . A pulsed near-infrared beam of photon energy 1.55 eV was split into a strong pump and a weaker probe beam. We measure the pump-induced change in transmission for varied time delays between pump and probe beam by monitoring the intensity of the probe beam with a silicon photodiode. To this end, lock-in amplification with a mechanical chopper modulating the pump beam was employed 27 . Two laser systems were used including an oscillator and an amplifier for the low and high-fluence range, respectively. The oscillator (amplifier) delivered nJ (mJ) pulses of 30 fs (40 fs) duration. In all experiments, the polarization of pump and probe beam was parallel. The maximal differential transmission as a function of delay between pump and probe beam was extracted for each fluence (Fig. 1e,f). We have verified that the substrate does not contribute to the transmission via nonlinear effects in the investigated fluence range by performing experiments on the bare substrate.
Double-bended saturation in graphene. Now we turn to the microscopic explanation of the double-bended saturation observed in Fig. 1c-f and discuss its many-particle-induced origin as well as the resulting implications for saturation experiments. To get the basic physical picture of the elementary processes we restrict our analytical analysis in good approximation to the diagonal dephasing and derive the counterpart of equation (1) for the solid-state two-band model based on equations (2 and 3). Similar to the derivation of equation (1) for a two-level system, we apply a parametric time dependence of r k (t) and p k (t) by solving the graphene Bloch equations quasi-stationary ( _ r k % 0, _ p k % 0) in rotating wave approximation, yielding with the time-dependent pulse intensity I(t) and also the formal counterpart of the saturation intensity, denoted with I k (t). Their ratio reads I k ðtÞ=IðtÞ¼g 2 k ðtÞ=4OðtÞ 2 . Provided that the full-time dependence of the decay channels G in k (t), G out k (t) and g k (t) ¼ G in k (t) þ G out k (t) is considered, the quasi-stationary solution from equation (4) exhibits excellent agreement with the numerical solution of the graphene Bloch equations (2 and 3), including the saturation behaviour, see Fig. 2 illustrating the saturation according to equation (4) and to the full numerical solution (red and blue triangles).
The saturation dynamics for the solid-state continuum, determined by equation (4), obviously differs from the two-level case in equation (1): in a solid state, in-scattering À G in k 0 ðtÞ Á and out-scattering processes À G out k0 ðtÞ Á populating or depopulating the optically excited electronic state k 0 are possible. For graphene, there is a linear relation between energy and momentum with e 0 ¼ :n F |k 0 |. The first term in equation (4) closely resembles the standard carrier occupation saturation from equation (1), however a total rate g k (t) appears instead of a pure out-scattering contribution Gg in the two-level case. Most importantly, an additional purely many-particle term occurs. This new contribution scaling with the ratio of the in-scattering rate G in k (t) and the total scattering rate g k (t) turns out to be of crucial importance to reproduce the experimentally observed differential transmission saturation behaviour (proportional to the pump-induced carrier occupation) including the qualitative double bending as well as the quantitative saturation intensities. The neglect of this manyparticle term leads to only a single saturation in the range of 3 mJ cm À 2 at very low occupations (green line in Fig. 2). This means that a standard saturation model according to the first term in equation (4), that is without considering carrier inscattering, would fail rigorously. It cannot even reproduce the standard Pauli-blocking-based saturation behaviour in the strong-excitation regime with pump fluences above 1 mJ cm À 2 . Therefore, the next step is to analyse the interplay of the in-and out-scattering rates occurring in equation (4) to get a microscopic understanding of the saturation behaviour in a solid-state absorber.
The microscopically determined scattering rates G in and G out are plotted in Fig. 3 as a function of fluence. Focusing first on the low-excitation limit, where G out ) G in , we find for the outscattering rate a constant offset G 0 (owing to electron-phonon scattering) and a dominating contribution proportional to accumulated carriers with F ðtÞ¼ R t À 1 IðtÞdt (owing to Coulomb scattering), where F ðtÞ is the fraction of the pump fluence until time t (red lines in Fig. 3a). For very low pump fluence, we find G out EG 0 , that is independent of the fluence and therefore r c k0;max ðIÞ $ I, as expected from the standard model,  equation (1). However, with increasing pump fluence, the maximal carrier occupation r c k0;max saturates. It turns out that due to the strong Coulomb interaction, the out-scattering rate in the Bloch equations for r c k0 rises for increasing intensity faster than the contribution of the optical drive (O 2 (t)pI(t)) itself, finally balancing the optical excitation. To understand this in more detail, we consider equation (4): for the first bending, where all results can be explained without Pauli blocking, we find r c k0;max ðIÞ % IðtÞ=I k ðtÞ % 4OðtÞ 2 G out ðtÞ 2 ¼4OðtÞ 2 = G 0 þ cF ðtÞ ð Þ 2 , where c is a numerical factor depending on the material properties. That is, r c k 0 ;max (I) saturates due to the denominator. The interplay between the optically driven increase of the carrier occupation r c k0 and an increasing intensity-sensitive depopulation of r c k 0 due to the efficient Coulomb-induced out-scattering (Fig. 1h) results in the first bending of r c k 0 ;max in the low-fluence regime, see Fig. 1d. This means that already at low intensities, where Pauli blocking, that is the denominator in equation (4), is neglected, the fluence dependence of G out k (I) starts to become important and causes the deviation from the linear dependence and thus leads to a bending. Note that this bending has nothing to do with the regular Pauliinduced saturation, since at the considered very low fluences, the carrier occupation of the excited state is below 0.01 (Fig. 1d). In contrast, our calculations reveal that the observed saturation-like behaviour is a signature of Coulomb-induced redistribution (outscattering) at the expense of the simultaneous optical pumping in a solid-state absorber.
As shown in Fig. 2, the additional many-particle in-scattering dominated contribution (second term in equation (4)) is essential to understand the high-intensity saturation, where Pauli blocking is dominant. The dynamics of the Coulomb-induced in-scattering rate G in ðccÞ is shown in Fig. 3b. We find a linear and quadratic dependence on the fluence F ðtÞ with G in ðccÞ ðtÞ¼c 1 F ðtÞ þ c 2 F 2 ðtÞ, where c 1 and c 2 are fitting parameters, see red lines in Fig. 3b. As a result, at sufficiently high fluences, the second, in-scattering dominated term in equation (4) becomes as important as the out-scattering from r c k 0 for the saturation behaviour. At elevated intensities, the in-scattering finally balances out the loss of carrier occupation in the optically excited state r c k0 through out-scattering processes, see Fig. 1i (detailed balance limit). In this regime, the carrier occupation increases again due to the optical pump, that is linear with the pump fluence, until a strong-excitation regime (mJ cm À 2 ), the Pauliblocking regime, is reached. Here in-and out-scattering are in balance and the conventional Pauli blocking of the excited state results in a standard saturation behaviour.

Discussion
Now, we can summarize our results from a microscopic point of view by distinguishing four distinct intensity regimes (Fig. 4a): (A) linear regime at very low pump fluences (few mJ cm À 2 ), where the maximal optically induced carrier occupation scales linearly with the intensity r c k 0 ;max (I)pI/G 0 with nearly intensityindependent scattering rate G 0 (owing to electron-phonon scattering; Fig. 1g). (B) First saturation-like behaviour driven by Coulomb-induced carrier out-scattering (B10 mJ cm À 2 ) that efficiently depopulates the optically excited state causing a deviation from the linear scaling (Fig. 1h). Here equation (4) simplifies to r c k 0 ;max (I)E4(O(I)/G out k (I)) 2 and the bending results from the strong increase of the intensity-dependent out-scattering rate. (C) Intermediate regime (few 100 mJ cm À 2 ), where the Coulomb-induced in-scattering balances out the out-scattering resulting again in an approximately linear increase of the carrier occupation due to the optical pump (Fig. 1i). (D) Second, Pauli blocking caused saturation in the high-excitation regime (few 10 mJ cm À 2 ; Fig. 1j). The different regimes are also recognizable regarding the time t max maximizing r c k 0 ;max (t) (with respect to the centre of the excitation pulse). The discrete steps of t max arise from our calculation beyond the rotating wave approximation and corresponds to a single-cycle duration of the exiting intensity, Fig. 4b. In the out-scattering dominated regime, t max shifts to earlier times due to the faster increasing out-scattering rate compared to the optical drive. Consequently, the time where outscattering overbalances the excitation shifts towards t max E0 fs. As soon as in-scattering becomes significant the delay between pulse centre and t max rises again. In the Pauli-blocking regime, t max is constant.
In summary, based on a joint theory-experiment study, we have found a yet undiscovered double-bended saturation of the optically excited carrier occupation in graphene for optical pulse excitations having a duration in the range of the predominant scattering channels. This unconventional behaviour beyond the standard saturable absorber model is explained by the universal structure of Coulomb-induced many-particle scattering processes, clarified via in-and out-scattering and their fundamental scaling with the pump fluence. In particular, we have identified both in theory and experiment a surprising low-intensity saturation-like behaviour of the optically excited carrier occupation that is not due to Pauli blocking, but can be ascribed to the out-scatteringinduced depopulation of electronic levels. Even if our accompanying analytical approach is an approximation, it covers important aspects of the dynamics. Most importantly, we find that in solids, a sub-linear relation between the optically excited carrier occupation and the pump fluence is not a sufficient criterion for an absorption saturation by means of fermionic Pauli blocking. We note that due to the universal structure of Bloch equations for solids, the observed effects are not restricted to graphene. The gained insights will be also valuable for the technological application of graphene as saturable absorber in lasers generating ultrashort pulses 28 . Here the optically induced transmission of graphene at comparably low fluences is a specific advantage over conventional materials. Specifically, the fluences are typically in the 10 mJ cm À 2 range, which corresponds to the regime of the predicted unconventional saturation behaviour.
Data availability. The data that support the findings of this study are available from the corresponding author upon request.