Higher-order contributions and non-perturbative effects in the non-degenerate nonlinear optical absorption of direct-gap semiconductors

TL;DR

This paper uses the semiconductor Bloch equations to analyze nonlinear optical absorption in direct-gap semiconductors, revealing higher-order effects and non-perturbative phenomena that influence multi-photon absorption and pulse interactions.

Contribution

It provides a comprehensive analysis of higher-order and non-perturbative effects in nonlinear optical absorption, including analytical and numerical results for multi-photon processes.

Findings

Strong enhancement of two-photon absorption under non-degenerate conditions

Rich frequency dependence of higher-order multi-photon absorption

Non-monotonous intensity-dependent absorption behavior

Abstract

The semiconductor Bloch equations for a two-band model including inter- and intraband excitation are used to study the nonlinear absorption of single and multiple light pulses by direct-gap semiconductors. For a consistent analysis the contributions to the absorption originating from both the interband polarization and the intraband current need to be included. In the Bloch equation approach theses contributions as well as different excitation pathways in terms of sequences of inter- and intraband excitations can be evaluated separately which allows for a transparent analysis, the identification of the dominant terms, and analyzing their dependence on the excitation conditions. In the perturbative regime, we obtain analytical expressions for the multi-photon absorption coefficients for continuous-wave excitation. These results are shown to agree well with numerical results for short pulses and/or finite dephasing and relaxation times and we confirm the previously predicted strong enhancement of two-photon absorption for non-degenerate conditions for pulsed excitation. We discuss the dependencies on the light frequencies, initial band populations, and the time delay between the pulses. The frequency dependence of the two-photon absorption coefficient for non-degenerate excitation is evaluated perturbatively in third-order. The higher-order contributions to the optical absorption include three- and four-photon absorption and show a rich frequency dependence including negative regions and dispersive lineshapes. Non-perturbative solutions of the Bloch equations demonstrate a strongly non-monotonous behavior of the intensity-dependent optical absorption for a single incident pulse and in a pump-probe set-up.

Figures

Equations

$${\hat{H}}_e$$

$${\hat{H}}_i$$

$$\frac{\partial }{\partial t}{p}_{\mathbf{k}}$$

$$\frac{\partial }{\partial t}{n}_{q\mathbf{k}}$$

$$\frac{\partial }{\partial t}{p}_{\mathbf{k}}^{(1)}$$

$$\frac{\partial }{\partial t}{n}_{q\mathbf{k}}^{(1)}$$

$$\frac{\partial }{\partial t}{p}_{\mathbf{k}}^{(m)}$$

$$-i({\omega }_{\mathbf{k}}-i{\gamma }_p)p_{\mathbf{k}}^{(m)},$$

$$\frac{\partial }{\partial t}{n}_{q\mathbf{k}}^{(m)}$$

$$-{\gamma }_n{n}_{q\mathbf{k}}^{(m)}.$$

$${\mathbf{j}}_p$$

$${\mathbf{j}}_n$$

$${\omega }_k={\omega }_{\mathrm{g}}+\frac{{\omega }_{\mathrm{b}}}{2}[1-\cos (ka)],$$

$${\omega }_{\mathbf{k}}={\omega }_{\mathrm{g}}+\hslash {\mathbf{k}}^2/2m^{\ast },$$

$${\alpha }^{(l)}({\omega }_1,{\omega }_2)$$

$${\hat{\mathbf{e}}}_1\cdot {\mathbf{P}}^{(3)(1|0)}$$

$${\hat{\mathbf{e}}}_1\cdot {\mathbf{j}}_n^{(3)(1|0)}$$

$${\alpha }^{(l)}({\omega }_1,{\omega }_2)=-\frac{1}{2}{\left(\frac{2}{{ϵ}_{0}c}\right)}^l$$

$$\frac{\partial }{\partial t}{p}_{\mathbf{k}}$$

$$\frac{\partial }{\partial t}{n}_{q\mathbf{k}}$$

$$\frac{\partial }{\partial t}{p}_{\mathbf{k}}^{(1)}$$

$$\frac{\partial }{\partial t}{\bar{p}}_{\mathbf{k}}$$

$$+\frac{e}{\hslash }{\mathbf{E}}_1\cdot {\nabla }_{\mathbf{k}}{p}_{\mathbf{k}}+\frac{e}{\hslash }{\mathbf{E}}_2\cdot {\nabla }_{\mathbf{k}}(p_{\mathbf{k}}^{(1)}+{\bar{p}}_{\mathbf{k}}),$$

$$\frac{\partial }{\partial t}{\bar{n}}_{q\mathbf{k}}$$

$$+\frac{e}{\hslash }{\mathbf{E}}_2\cdot {\nabla }_{\mathbf{k}}{\bar{n}}_{q\mathbf{k}}+\frac{i}{\hslash }\mathbf{d}\cdot {\mathbf{E}}_2({\bar{p}}_{\mathbf{k}}^{\ast }-{\bar{p}}_{\mathbf{k}}+p_{\mathbf{k}}^{\ast (1)}-p_{\mathbf{k}}^{(1)}),$$

$${\chi }^{(3)}({\omega }_1,{\omega }_2,-{\omega }_2)$$

$${\alpha }^{(2)}({\omega }_1,{\omega }_2)$$

$$|v_{k_{\mathrm{r}}}|=\sqrt{2\hslash ({\omega }_{k_{\mathrm{r}}}-{\omega }_{\mathrm{g}})/m^{\ast }}.$$

$$(2/3)\sqrt{2\hslash /m^{\ast }}({\omega }_1+{\omega }_2-{\omega }_{\mathrm{g}}{)}^{3/2},$$

$${\alpha }^{(3)}({\omega }_1,{\omega }_2)$$

$$\displaystyle\quad\cdot\big{[}\omega_{1}\omega_{2}\left(\omega_{1}^{2}+\omega_{1}\omega_{2}+4\omega_{2}^{2}\right)\omega_{k_{\mathrm{r}}}^{\prime\prime}$$

$$\displaystyle\quad\ \;-2\left(\omega_{1}+\omega_{2}\right)\left(\omega_{1}^{2}+2\omega_{2}^{2}\right)\omega_{k_{\mathrm{r}}}^{\prime}{}^{2}\big{]}{}^{2}.$$

$${\alpha }^{(l+1)}({\omega }_1,{\omega }_2)=G({\omega }_1,{\omega }_2)\frac{\partial }{\partial {\omega }_1}{\alpha }^{(l)}({\omega }_1,{\omega }_2),$$

$${\alpha }^{(l)}({\omega }_1,{\omega }_2)=-[2^{2l-1}({\omega }_1+{\omega }_2+i{\gamma }_n{)}^{l-1}({\omega }_1-{\omega }_0+i{\gamma }_p{)}^l{]}^{-1}$$

$$\exp \left(-\frac{2(l-1)(t_1-t_2{)}^2}{(l-1){\tau }_1^2+{\tau }_2^2}\right).$$

$$\frac{\partial }{\partial t}{p}_{\mathbf{k}}^{(1)[e]}$$

Full text

I Introduction

Two- and multi-photon absorption (2PA/MPA) are among the most fundamental effects Göppert‐-Mayer (1931); Braunstein and Ockman (1964); Loudon (1962); Hopfield and Worlock (1965); Mahan (1968) in nonlinear optics that has been proposed and studied for several decades considering various physical systems. Applications include infrared detection,Schneider et al. (2005, 2009); Fishman et al. (2011); Pattanaik et al. (2016); Piccardo et al. (2018) optical switching,Venkataraman et al. (2011); Yavuz (2006) and optical control of electrical currents.Atanasov et al. (1996) The time-reversed process of stimulated two-photon emissionHayat et al. (2008); Reichert et al. (2016) is expected to eventually lead to the realization of two-photon lasing. For solid-state systems this topic is currently studied intensively, also with regards to applications in the ultrafast and strongly nonlinear regimes. However, presently still several fundamental aspect for MPA in solids are not well understood, in particular, with regards to higher-order processes and the dynamical behavior.

Traditionally 2PA has been studied mainly by transition rate approaches.Braunstein and Ockman (1964); Hopfield and Worlock (1965); Loudon (1962); Basov et al. (1966); Mahan (1968) Even simplistic band structure models lead to reasonable predictions of 2PA spectra for various semiconductors.Bolger et al. (1993); Said et al. (1992); Sheik-Bahae et al. (1991) However, transition rate approaches are not easy to evaluate for higher numbers of simultaneously absorbed photons and also do not account for dephasing or relaxation effects. A more versatile approach describing the nonlinear-optical response of semiconductors is based on susceptibilities. Typically this approach is applied in a semiclassical and/or an adiabatic approximation. In Refs. Sipe and Ghahramani, 1993; Aversa et al., 1994; Aversa and Sipe, 1995 it has been evaluated in the framework of a semiclassical density matrix formalism for the steady-state nonlinear-optical response.

In this study we use a more general approach based on the semiconductor Bloch equations (SBE) extended by the intraband terms,Haug and Koch (2009); Meier et al. (1994) thereby neglecting electron-electron interactions. Our approach follows Refs. Sipe and Ghahramani, 1993; Aversa et al., 1994; Aversa and Sipe, 1995 in using the length gauge and separating the intraband and interband motion from the start. This avoids unphysical divergencies that plagued some earlier calculations.Sheik-Bahae et al. (1990, 1991, 1994); Moss et al. (1990) The perturbative expansion, where intraband and interband motion are treated equivalently, yields macroscopic susceptibilities and ac conductivities, from which in turn one obtains nonlinear absorption coefficients. However, a main aspect of this work is the non-perturbative evaluation of the total absorption of (one or several) strong pulses or the absorption of a weak probe pulse in the presence of a strong pump pulse. Our method can be applied to arbitrary band structures. To demonstrate the basic principles we focus on a two-band model with $k$-independent dipole matrix element. Due to their self-consistent calculation of the interband dipole matrix element, Aversa et al.Aversa et al. (1994) obtained a different frequency dependence of the 2PA coefficient in an analogous model, but the general tendency of a steep increase towards strongly non-degenerate 2PACirloganu et al. (2011) is fully reflected by our work. In other aspects we find good agreement with the literature on 2PA. More unexpected results are found for higher order MPA, which we evaluate analogously. In particular, the 3PA coefficient is found to vanish in the one-dimensional model if the frequency sum assumes a certain value not far above the band gap. The transients can be generally explained by pulse area overlaps as long as no saturation effects occur. Dynamic effects are discussed in the regime of high pump intensities.

Our paper is organized as follows. The general formalism for obtaining nonlinear absorption rates from the SBE is outlined in Sec. II. We are interested not only in absorption coefficients describing a perturbative steady state response, but also in the non-perturbative absorption of short and possibly strong pulses. While the absorption of a single pulse is anticipated in Sec. II.3, the results for the pump-probe scheme are discussed in Sec. III. Analytical 2PA coefficients are presented in Sec. III.1 explicitly for the two-band model and compared with previous works. In Sec. III.2 the SBE are evaluated to higher orders in the intensity of the pump pulse. We discuss both the change of third order effects (2PA) as well as the appearance of higher order effects (MPA) for more intense pump pulses. In principle this expansion can be taken to arbitrary order. In Sec. III.3 we evaluate the approach including the pump pulse to infinite order. In Sec. IV we summarize our main results and discuss further directions.

II Theoretical formalism

In the following, we present the approach which we use to analyze the nonlinear optical response of bulk semiconductor and extended semiconductor nanostructures (of dimension $d$) within a two-band model. Following previous works we use the length gauge to describe the light-matter interaction including inter- and intraband excitations. The resulting SBE have been employed previously to study, e.g., coherent effects induced by static and time-dependent electric fields in semiconductor superlattices Meier et al. (1994, 1995a), the generation of photocurrents by two-color fields Duc et al. (2005, 2006), high-harmonic generation resulting from the excitation with intense optical Golde et al. (2008) and THz fields Golde et al. (2011), as well as transient Wannier-Stark localization Schmidt et al. (2018).

Here, we consider the case of a homogeneous excitation with several radiation fields of arbitrary frequencies. Propagation effects as well as the many-body Coulomb interaction are neglected for simplicity since we focus on the interplay between inter- and intraband excitations in particular for strongly non-resonant excitation conditions. Dephasing and relaxation processes are treated phenomenologically and sometimes disregarded in order to obtain analytical expressions.

We start by presenting the SBE and then describe the general perturbative expansion in powers of the optical field. Next we consider a pump-probe set-up and include the propagation directions of the optical pulses into the analysis. Furthermore, we present equations which allow to include the pump field non-perturbatively.

II.1 Non-perturbative extended Bloch equations

The system is described microscopically by a two-band model Hamiltonian $\hat{H}=\hat{H}_{0}+\hat{H}_{\mathrm{LM}}$, where $\hat{H}_{0}$ is the single-particle Hamiltonian of the valence and conduction band electrons, $\hat{H}_{0}=\hbar\sum_{\mathbf{k}}(\omega_{v\mathbf{k}}\hat{a}_{v\mathbf{k}}^{\dagger}\hat{a}_{v\mathbf{k}}+\omega_{c\mathbf{k}}\hat{a}_{c\mathbf{k}}^{\dagger}\hat{a}_{c\mathbf{k}})$. The semiclassical light-matter interaction is given by $\hat{H}_{\mathrm{LM}}=\hat{H}_{e}+\hat{H}_{i}$, where in dipole approximation the interband part $\hat{H}_{e}$ and the intraband part $\hat{H}_{i}$ are given by

[TABLE]

where $\mathbf{d}$ is the interband transition dipole matrix element which is taken as $k$-independent.

The Bloch equations are obtained from the Heisenberg equations of motion for the hole occupation $n_{h\mathbf{k}}=1-\langle\hat{a}_{v\mathbf{k}}^{\dagger}\hat{a}_{v\mathbf{k}}\rangle$, the electron occupation $n_{e\mathbf{k}}=\langle\hat{a}_{c\mathbf{k}}^{\dagger}\hat{a}_{c\mathbf{k}}\rangle$, and the interband polarization $p_{\mathbf{k}}=\langle\hat{a}_{v\mathbf{k}}^{\dagger}\hat{a}_{c\mathbf{k}}\rangle$ and read Meier et al. (1994, 1995b); Duc et al. (2005)

[TABLE]

where the second equation holds for $q=e,h$ and we defined the transition frequency ${\omega}_{\mathbf{k}}\equiv\omega_{c\mathbf{k}}-\omega_{v\mathbf{k}}$ and introduced $n_{\mathbf{k}}\equiv n_{e\mathbf{k}}+n_{h\mathbf{k}}$, i.e., the sum of the electron and hole carrier occupations. The relaxation and dephasing constants are denoted by $\gamma_{n}=1/T_{1}$ and $\gamma_{p}=1/T_{2}$, respectively. Both evolution equations contain two types of source terms, an interband term (labeled $e$) and an intraband term (labeled $i$) which are treated equivalently in the following analysis.

II.2 Perturbative expansion

Macroscopic response functions (susceptibilities and conductivities) are obtained by expanding the quantities $X_{\mathbf{k}}=\{p_{\mathbf{k}},n_{e\mathbf{k}},n_{h\mathbf{k}}\}$ in powers of the exciting field $\mathbf{E}(t)$ as $X_{\mathbf{k}}(t)=\sum_{m=0}^{\infty}X_{\mathbf{k}}^{(m)}(t)$ with $X_{\mathbf{k}}^{(m)}\propto\mathbf{E}^{m}$. In our approach $n_{q\mathbf{k}}^{(0)}=n_{q\mathbf{k}}(t\to\pm\infty)$ acts at the same time as the initial density and as the equilibrium distribution, cf. Eq. (4). Due to the intraband excitation, a finite initial density $n_{q\mathbf{k}}^{(0)}$ serves as a source for a linear electronic response which contributes to the optical response in third order.

In first-order in the optical field ($m=1$) the evolution equations are given by

[TABLE]

The $e$ and $i$ source terms lead to ${p}_{\mathbf{k}}^{(1)}$ and $n_{q\mathbf{k}}^{(1)}$, respectively. Also for higher orders ($m>1$) both types of source terms appear in each equation

[TABLE]

Due to the two source terms in each equation, the total result in $m$th order is given by a sum of 2m terms involving all combinations of the different source termsAversa and Sipe (1995), cf. Figure 1. The evolution equations for the individual pathways are given exemplary for up to $m=2$ in the Appendix A. Since only some of these pathways are relevant for a particular nonlinear optical process, it is instructive to study each contribution separately by zeroing all others. Below we label the pathways by a series of $e$’s and $i$’s, where the lowest order is given on the right. For example, in second order there are two paths $ie$ and $ei$ leading to $p_{\mathbf{k}}^{(2)}$, as well as two paths $ee$ and $ii$ leading to $n_{\mathbf{k}}^{(2)}$. A path with an odd (even) number of $e$ terms leads to a polarization (density). All paths starting from $i$ in first order vanish for a clean, cold semiconductor, i.e., without any electron and hole occupations present prior to the optical excitation.Aversa and Sipe (1995)

A truncation of the perturbative expansion of $X_{\mathbf{k}}(t)$, i.e., Eqs. (5)-(8), requires that the incident electric fields are sufficiently weak such that $|X_{\mathbf{k}}(t)^{(m)}|\ll 1$ for all $m>0$.

II.3 Absorption rate - general formalism

The rate at which the energy of a light beam is absorbed per unit volume is generally given by a macroscopic spatial average of the product of the current density $\mathbf{j}(\mathbf{r},t)$ and the field amplitude $\mathbf{E}(\mathbf{r},t)$, i.e., $\langle\mathbf{j}\cdot\mathbf{E}\rangle$.Mahr (1975) This definition still holds in a non-perturbative regime, where the incident and transmitted pulse shapes generally differ. Since propagation effects are disregarded here, i.e., we assume that variations of the light intensity within the sample are small, this averaging does not need to be carried out explicitly and we can set $\mathbf{r}=0$. For the case of two or more non-copropagating incident light beams, the spatial averaging is carried out as a directional expansion as described in Sec. II.4.

The electric current density $\mathbf{j}$ is induced in the medium by $\mathbf{E}$ either via a change of the interband polarization or via carrier occupations, $\mathbf{j}=\mathbf{j}_{p}+\mathbf{j}_{n}$.Sipe and Ghahramani (1993) For an isotropic band structure the part $\mathbf{j}_{n}$ vanishes if the densities are symmetric distributions in $k$-space. A contribution $\langle\mathbf{j}_{n}\cdot\mathbf{E}\rangle$ arises from density distributions which are asymmetric along the $k$-direction parallel to $\mathbf{E}$. For symmetric $n_{q\mathbf{k}}^{(0)}$ these contributions correspond to conductivity paths with an odd number of intraband terms $i$. Transforming the $k$-sums into integrals, $\sum_{\mathbf{k}}\to(L/2\pi)^{d}\int_{\,\mathrm{BZ}}\!\!\mathrm{d}\mathbf{k}$, and omitting the prefactor, the two contributions to the current density are given by

[TABLE]

Substituting perturbative quantities $X_{\mathbf{k}}^{(m)}$ into these equations one analogously obtains perturbative results for the current density $\mathbf{j}^{(m)}$.

Inspecting the perturbative expansion of the Bloch equations we see that the excited densities are equal in both bands, $n_{e\mathbf{k}}-n_{e\mathbf{k}}^{(0)}=n_{h\mathbf{k}}-n_{h\mathbf{k}}^{(0)}$, if the optically irrelevant contributions from paths without interband excitation are neglected. Thus, for our optical studies instead of Eq. (10) we use

[TABLE]

and have to solve the $n$-equations not individually for each band but only as a sum.

Here we describe the two-band models used in this work. As long as no initial densities $n_{q\mathbf{k}}^{(0)}$ are present, only the valence to conduction interband transition frequency $\omega_{\mathbf{k}}$ enters the Bloch equations along with Eqs. (9),(11) and the individual band dispersions are not relevant. We either use a one-dimensional tight-binding band structure with lattice constant $a$,

[TABLE]

where $\hbar\omega_{\mathrm{g}}$ is the band gap and $\hbar\omega_{\mathrm{b}}$ the joint band width. This model is particularly suitable for very strong fields where the electronic excitations are driven through the entire Brillouin zone. The other model used for analyzing the optical response near the band gap is a parabolic dispersion

[TABLE]

where $1/m^{*}=1/m_{e}^{*}+1/m_{h}^{*}$ is the inverse reduced effective mass. When considering initial quasiequilibrium carrier distributions $n_{q\mathbf{k}}^{(0)}$, the individual effective masses ($m_{e}^{*},m_{h}^{*}$) need to be specified. Here we choose the values for the lowest conduction and the heavy hole bands of GaAs.

Further simplifications arise since we consider linearly and parallel polarized incident pulses here. In this case the transition dipole $\mathbf{d}$ enters only through its projection $\mu$ onto the field direction. Here we use a constant, i.e., $k$-independent, interband dipole matrix element suitable for bulk GaAs of $\mu=3ea/5.65$, that is $\mu=3e$Å  after substituting the lattice constant of GaAs, $a\approx$ 5.65Å. The lattice constant $a$ does not enter our calculations explicitly, but is transferred to other quantities such as the field strength in Figure 2.

Before we describe the generalization of the light absorption formalism to multiple propagation directions in the next Section, we evaluate here the approach for a single pulse. Besides serving as a simple example, this is experimentally relevant since degenerate MPA can be measured this way. We assume an initially unexcited, one-dimensional sample which is described by a two-band tight-binding band structure (12). Figure 2(a) and (b) show the non-perturbative absorption, calculated as the time integral of $jE$, of a single pulse with a Gaussian envelope $e^{-t^{2}/\tau_{1}^{2}}$ and different central frequencies.

The one-photon absorption (1PA) resonant region, $\omega_{1}>\omega_{\mathrm{g}}$, shows regular Rabi oscillations of period $\Delta\hat{E}_{1}=\hbar/(\sqrt{\pi}\mu\tau_{1})\approx 0.014$V$\cdot 2\pi/a$. For weak fields, $\hat{E}_{1}\ll\Delta\hat{E}_{1}$, the 1PA increases quadratically with the field amplitude $\hat{E}_{1}$ as expected from perturbation theory. Figure 2(b) shows explicitly for $\hbar\omega_{1}\approx 1.7$eV that the absorption increases roughly linearly with each Rabi cycle in the case $T_{2}\to\infty$. For finite $T_{2}$ the increase is stronger and follows a $\hat{E}_{1}^{1/2}$-dependence when the excitation frequency is within the parabolic region of the band structure. As expected, the effect of the intraband motion is negligible for 1PA case in the considered field amplitude range.

In the frequency region $\omega_{\mathrm{g}}/2<\omega_{1}<(\omega_{\mathrm{g}}+\omega_{\mathrm{b}})/2$, which partly overlaps with other MPA regions, 2PA can be clearly recognized. Since 2PA scales as $\propto\hat{E}_{1}^{4}$ at field amplitudes well below saturation effects, it is negligible compared to 1PA for weak fields. Saturation sets in at $\hat{E}_{1}\approx 0.05$V$\cdot 2\pi/a$, corresponding to a peak intensity of 41 GW/cm2 (using the lattice constant of GaAs). This agrees in order of magnitude with the 2PA saturation observed in bulkLami et al. (1996) and monolayerDong et al. (2018) semiconductors. We note that Rabi-like oscillations caused by 2PA appear in two different ways. For higher frequencies, e.g. $\hbar\omega_{1}\approx 1.2$eV, weak modulations of large period are visible. On the other hand, more regular Rabi-modulations with shorter period appear beyond a certain resonance signal which starts slightly above $\omega_{\mathrm{g}}/2$ (the 2PA resonance with the lower band edge) and is Stark shifted to higher frequencies. The shift is first quadratic (at about $\hat{E}_{1}\approx 0.06$V$\cdot 2\pi/a$) and later linear in $\hat{E}_{1}$. Below of this resonance line further lines appear regularly and the entire absorption edge remains roughly at $\omega_{\mathrm{g}}/2$ (hard to distinguish due to the overlap of 2PA and 3PA). Generally, Rabi-like modulations for 2PA appear to be weaker than for 1PA (due to the effect of intraband motion). This may be the reason why two-photon Rabi-cycles have been observed so far only in quasi zero-dimensional systems, such as atoms,Gentile et al. (1989) molecules,Linskens et al. (1996) and semiconductor quantum dots.Stufler et al. (2006)

For 3PA the behavior is somewhat analogous and even more easily distinguishable. The most dominant signal is the band gap resonance $\omega_{\mathrm{g}}/3$, which is again Stark shifted upwards. In the perturbative analysis in Sec. III.2 we will discuss in more detail that this dominant signal is separated from a broader 3PA range further above the band gap, which is also visible in Figure 2(a). Relatively pronounced Rabi oscillations appear beyond this upwards shifting band gap resonance. The field amplitude range where perturbative analysis is valid for 3PA is very small since the weak field regime for 3PA with a scaling $\propto\hat{E}_{1}^{6}$ has only a narrow validity range before also higher orders contribute. Note also that all higher MPA sets in at roughly comparable field amplitudes. Another strong 3PA signal is the resonance with the upper band edge $\omega_{\mathrm{g}}+\omega_{\mathrm{b}}$.

For higher MPA basically only narrow resonance lines appear, which are resonant with frequencies slightly above the band gap. Here we mainly see the linear frequency shift of these resonances. We confirm in the perturbative analysis in Sec. III.2 that 4PA is indeed restricted to a rather narrow range above the band gap. In our model, for certain MPA, e.g. 6PA, again another signal appears corresponding to a resonance with the upper band edge and Stark-shifting to lower frequencies for strong fields.

II.4 Perturbative treatment of pump-probe experiments

We consider MPA as measured via the differential absorption of a probe pulse (’1’) in the presence of a non-copropagating pump pulse (’2’). The probe pulse is assumed to be sufficiently weak so that material excitation of higher than linear order in its field amplitude can be neglected. Depending on the pump intensity $I_{2}$, nonlinear-optical processes consisting of the simultaneous absorption of one probe photon and $(l-1)$ pump photons are relevant, which are described by a $l$PA-coefficient $\alpha^{(l)}$. The total absorption coefficient is thus expanded in the pump intensity $I_{2}$ as $\alpha=I_{1}\sum_{l=1}^{\infty}\alpha^{(l)}I_{2}^{l-1}$. The formalism presented here is also applicable without the restriction to weak probe pulses or to the number of pulses, but for such cases the number of required equations in the directional expansion, see below, grows strongly.

The different propagation directions of pump and probe are taken into account by including the dependence of the field on the macroscopic position vector $\mathbf{r}$.Binder et al. (1991) The individual fields are plane waves with a temporal envelope and a central frequency $\mathbf{E}_{\nu}=\hat{\mathbf{e}}_{\nu}\hat{E}_{\nu}(t)\cos(\mathbf{k}_{\nu}\cdot\mathbf{r}-\omega_{\nu}t)$, for $\nu=1,2$. We assume the two fields to have the same linear polarization, i.e., $\hat{\mathbf{e}}_{1}=\hat{\mathbf{e}}_{2}$. As mentioned above, the transition dipole matrix element $\mathbf{d}$ enters only through its projection onto the field direction, $\mu=\mathbf{d}\cdot\hat{\mathbf{e}}_{1}$. The intensities are obtained from the field amplitudes as $I_{\nu}(t)=(\epsilon_{0}c/2)|\hat{E}_{\nu}(t)|^{2}$. In this work we consider either \hat{E}_{\nu}(t)=\hat{E}_{\nu}=\const. (continuous wave limit) or finite pulses with a Gaussian envelope function $\hat{E}_{\nu}(t)$. Using a standard perturbative expansion in the two field amplitudes $\mathbf{E}_{\nu}$ (and abbreviating $\mathbf{j}^{(m)}\equiv\mathbf{j}^{(m)(1|m-1)}$ where in the latter expression the superscripts $1$ and $m-1$ denote the orders in the probe and pump field amplitude, respectively), we find the MPA coefficientMahr (1975) in the steady-state response as

[TABLE]

The factor $l$ in the denominator ensures that for $\omega_{1}=\omega_{2}$ the obtained absorption coefficient agrees with the single beam setup, i.e., $\alpha^{(l)}(\omega_{1},\omega_{1})=\alpha^{(l)}(\omega_{1})$. For pulses of finite duration we have accordingly

[TABLE]

Due to the spectral width of the pulses, this expression yields approximately the frequency dependence of the absorption coefficient.

Here we replace the spatial averaging appearing in Eqs. (14) and (15) by a directional expansion. The directional expansion differs from the power expansion only in the splitting of each field into its two Fourier components, $\cos(\mathbf{k}_{\nu}\cdot\mathbf{r}-\omega_{\nu}t)\to\frac{1}{2}(\exp[i(\mathbf{k}_{\nu}\cdot\mathbf{r}-\omega_{\nu}t)]+\exp[-i(\mathbf{k}_{\nu}\cdot\mathbf{r}-\omega_{\nu}t)])$, corresponding to forward/backward propagating complex fields. The current density $\mathbf{j}^{(m)(a|b)}$ is obtained from the complex quantities $p_{\mathbf{k}}^{(m)(a|b)}$ and $n_{q\mathbf{k}}^{(m)(a|b)}$, where as before the superscript $(m)$ refers to the total order and the additional superscript $(a|b)$ denotes that the term is proportional to $\exp[i(a\mathbf{k}_{1}+b\mathbf{k}_{2})\cdot\mathbf{r}]$, i.e., $(1|0)$ corresponds to the probe direction.Lindberg et al. (1992); Sieh et al. (1999); Weiser et al. (2000); Binder et al. (1991) In the steady-state response the probe-directed signals are related to the field amplitudes, e.g., in third order, by

[TABLE]

where $\chi$ is the susceptibility and $\sigma$ the ac-conductivity. Here, we do not study the tensor character of the response functions since we consider linearly and parallel polarized fields and an isotropic medium, thus only the diagonal elements are relevant. Using relations of the type (16)-(17) in Eq. (14) we obtain

[TABLE]

where single (double) prime stands for the real (imaginary) part, and the arguments of $\chi$ and $\sigma$ are $(\omega_{1};\omega_{1},\omega_{2},-\omega_{2},\ldots)$ according to the order $l$.

Note that a directional expansion may not be necessary for certain choices of parameters. For instance, in the case of non-degenerate 2PA the time integration in Eq. (15) automatically singles out the probe-directed contributions for sufficiently long pulses. The absorption coefficient for 2PA (18) was given by Bredikhin and Genkin,Bredikhin and Genkin (1970) but without the conductivity term, which plays an important role for $l>1$ even for the case of filled bands. Third-order expressions for the 2PA coefficient for varying excitation frequencies were obtained using a transition rate approachSheik-Bahae et al. (1994) and by a susceptibility approach Aversa et al. (1994) and we compare to these results in Sec. III.1. Furthermore, the optical injection of photocurrents by two light beams having frequencies $\omega$ and $2\omega$ with $2\hbar\omega$ larger than the gap energy was studied by Atanasov et al.Atanasov et al. (1996).

II.5 Bloch equations non-perturbative in the pump field

So far we have treated both probe and pump field perturbatively, yielding absorption coefficients for each order of the pump intensity. When a distinction between these orders is not required or when the pump field is so intense that an order expansion cannot be applied, it is possible to include the pump field non-perturbatively. We still consider a weak probe pulse so that its absorption remains linear in the probe intensity. If both pulses are treated non-perturbatively, the directional expansion is inevitable and immensely involved. As discussed below and in Sec. III.3, a directional expansion may be omitted for the weak probe scheme as long as the pump pulse intensity does not reach very high values.

The differential probe absorption in the presence of an arbitrarily strong pump field is extracted from the mixed quantities $\bar{p}_{\mathbf{k}},\bar{n}_{q\mathbf{k}}$ in the following set of SBE:Schmidt et al. (2018)

[TABLE]

with $\bar{\omega}_{\mathbf{k}}\equiv\omega_{\mathbf{k}}-i\gamma_{p}$. The quantities $p_{\mathbf{k}},n_{q\mathbf{k}}$ are purely due to the pump, while $p_{\mathbf{k}}^{(1)}$ is the linear probe polarization. A possibly present initial density $n_{q\mathbf{k}}^{\mathrm{(0)}}$ needs also to be included by the initial condition $n_{q\mathbf{k}}(t\rightarrow-\infty)=n_{q\mathbf{k}}^{\mathrm{(0)}}$ whereas the mixed densities that are induced by both pulses start with $\bar{n}_{q\mathbf{k}}(t\rightarrow-\infty)=0$. For sufficiently strong pump pulses, frequency components will appear near the probe frequency and other directions than the probe may contribute to the absorption obtained analogously to Eq. (15). For weaker or longer pulses, such that only a maximum of $m$ orders is relevant, the response described by Eqs. (19)-(23) agrees with the perturbative expansion.

III Results for pump-probe scheme

In the following we present and discuss our analytical and numerical results. Various regimes of pump intensities are discussed successively, starting in Sec. III.1 from sufficiently weak pump intensities and an off-resonant probe such that a third-order analysis is valid. Higher-order effects are evaluated for $m=5$ and $m=7$ in Sec. III.2. We analyze both the corresponding MPA coefficients ($\alpha^{(3)}$ and $\alpha^{(4)}$) and the corrections to the lower orders, e.g. $\chi^{(5)}$-corrections for a resonant sum frequency, i.e., $\omega_{1}+\omega_{2}$ larger than the gap frequency. Finally, the regime of very high pump intensities leading to dynamic saturation effects is investigated in Sec. III.3.

III.1 Two-photon absorption in third order

Here we start by deriving analytical expressions for the 2PA coefficients in third order in the continuous wave limit and compare them with numerical results obtained for Gaussian pulses. The susceptibility $\chi^{(3)}$ and the conductivity $\sigma^{(3)}$ both have the following form

[TABLE]

where the paths leading to $\chi^{(3)}$ are $\{iie$, $eee$, $iei$, $eii\}$, while those for $\sigma^{(3)}$ are $\{eie$, $iee$, $eei$, $iii\}$. The two sums over subpaths ($\ell$ and $\eta$) are due to the order in which the fields $E_{1},\ E_{2},\ E_{2}^{*}$ appear in the perturbation expansion. $\ell$ labels the order in which field $E_{1}$ appears, while $\eta$ labels whether $E_{2}$ or $E_{2}^{*}$ appears first in the series.

From the analytical expressions for $\phi_{\mathbf{k}}^{(\mathrm{path};\ell)}(\omega_{1},\omega_{2})$, see Appendix B, all contributing third-order effects, i.e., 2PA, electronic Raman, and optical Stark effect, can be deduced. Compact results are obtained in the limit of vanishing relaxation and dephasing constants, i.e., $\gamma_{n}\to 0$ and $\gamma_{p}\to 0$. Strictly speaking, the results depend on the order of the limits and we take $\gamma_{n}\to 0$ first since the relaxation time is typically longer than the dephasing time.

We restrict the analytical evaluation to the continuous wave limit to 2PA, hence we consider $0<\omega_{\nu}<\omega_{\mathbf{k}}$ for $\nu=1,2$. In third-order four out of eight paths (indicated by solid circles in Figure 1) show a resonance factor $(\omega_{1}+\omega_{2}-\bar{\omega}_{\mathbf{k}})^{-1}$ and are thus relevant for 2PA. The other four paths, for which the second-order quantity is a density $n_{\mathbf{k}}^{(2)}$, actually lead to divergent absorption contributions, arising from a forbidden truncation of the perturbation series for continuous-wave fields and $\gamma_{n}=0$ and $\gamma_{p}\to 0$. Such contributions are thus omitted in our evaluations. For sufficiently short pulses, such these $n_{\mathbf{k}}^{(2)}$ remain small, the truncation of the perturbation series is allowed and these paths do not significantly contribute as confirmed by our numerical solutions. The 2PA coefficient scales as $\mu^{2}$ if only ,,contributing” paths are considered - this is in contrast to the 3PA coefficient, where different powers of $\mu$ appear as indicated in Figure 3. It is noteworthy that the path $iei$ ($eei$), which vanishes for initially filled bands, is completely canceled by other terms of the path $iie$ ($eie$). As a consequence the derivatives of the initial density, e.g., $\nabla_{\mathbf{k}}n_{k}^{\mathrm{(0)}}$, do not enter the 2PA coefficient.

Using the response functions and the recipe outlined in Appendix B, the 2PA coefficient for a one-dimensional system is found from Eqs. (18) and (24) in the limit $\gamma_{n}\to 0$ and $\gamma_{p}\to 0$ as

[TABLE]

where $k_{\mathrm{r}}$ is the positive resonant point in the first Brillouin zone, i.e., $\omega_{k_{\mathrm{r}}}=\omega_{1}+\omega_{2}$ and $v_{k_{\mathrm{r}}}=\mathrm{d}\omega_{k}/\mathrm{d}k|_{k=k_{\mathrm{r}}}$ is the group velocity at this point, where for simplicity we assume that the band structure is isotropic. The dependencies of $\alpha^{(2)}$ on the frequency ratio $\omega_{1}/\omega_{2}$, the band structure, and initial band fillings are independent of each other and will be discussed separately in the following. Numerical results using Gaussian pulses of duration $\tau_{\nu}$ (solid lines) are compared with analytical continuous-wave expressions (dashed lines) in Figures 3-5.

Figure 3 shows the probe frequency dependence of the absorption for a fixed sum frequency $\omega_{1}+\omega_{2}$. Since the sample is initially unexcited, only one susceptibility path ($iie$) and one conductivity path ($eie$) contribute (with 2 of 6 subpaths for each path). All other paths are negligible, i.e., they vanish for the investigated analytical limit and are very small in the numerical evaluations with Gaussian pulses. Nearly perfect agreement between numerical and analytical results is obtained as long as both dephasing and relaxation times are longer than the chosen pulse durations $\tau_{\nu}$.

We note that since each beam delivers one photon to the 2PA process, $\alpha^{(2)}(\omega_{1},\omega_{2})/\omega_{1}$ should be symmetric with respect to the two frequencies. For weak relaxation this symmetry is fulfilled automatically. When we consider very short dephasing times in our numerical evaluations the absorption is reduced but still remains symmetric. On the other hand, since the relaxation time enters only in the conductivity path $eie$ the symmetry is lost for $T_{1}<\tau_{\nu}$. The absorption coefficient diverges if one of the light frequencies vanishes and the other is resonant. This drastic enhancement for non-degenerate excitation was also obtained previously.Sheik-Bahae et al. (1994); Aversa et al. (1994)

The general behavior of the 2PA coefficient for varying excitation frequencies is in good qualitative agreement with previous studies. However, regarding the exact dependence differing results exist in the literature. The frequency dependence obtained by the transition rate approach in Ref. Sheik-Bahae et al., 1994 is $\alpha^{(2)}(\omega_{1},\omega_{2})\propto(\omega_{1}+\omega_{2})^{2}/\omega_{1}^{3}\omega_{2}^{4}$. Aversa et al.Aversa et al. (1994) applied a susceptibility approach to a two-band model and obtained that $\alpha^{(2)}(\omega_{1},\omega_{2})\propto(\omega_{1}+\omega_{2})^{3}/\omega_{1}^{3}\omega_{2}^{4}$. For the case of a $k$-independent dipole matrix element, i.e., assuming that the interband matrix element of the position operator does not depend on $k$, the dependence obtained here is, however, in agreement with the result of Ref. Aversa et al., 1994, which is based on limiting the general $\chi^{(3)}$ analysis of Ref. Aversa and Sipe, 1995 to a two-band model.

The band structure enters $\alpha^{(2)}$ explicitly only via the group velocity in field direction. In the parabolic region near the band gap, Eq. (13), the group velocity is given by

[TABLE]

The results of Eq. (25) can be easily generalized to three dimensions in the parabolic region of the band structure, since the two directions perpendicular to the field can be included by a continuum of energetically shifted one-dimensional parabolic bands. These correspond to a constant two-dimensional density of states in the directions perpendicular to the linear polarization of the incident fields. In this case the factor $|v_{k_{\mathrm{r}}}|$ in Eq. (25) is replaced by

[TABLE]

while all other factors remain unchanged. The proportionality of 2PA coefficients to $(\omega_{1}+\omega_{2}-\omega_{\mathrm{g}})^{3/2}$ is well-established from the scaling rules by Wherrett.Wherrett (1984) The dependence of $\alpha^{(2)}$ on the band gap of the three-dimensional model is shown in Figure 4 and again compared with our numerical results. We find a very good agreement apart from a very small tail below the band gap which simply arises from the finite pulse durations.

The band-filling factor $(1-n_{\mathbf{k}}^{\mathrm{(0)}})$ appearing in $\alpha^{(2)}$ is the same as for the 1PA coefficient and implies a corresponding change from 2PA to stimulated emission when population inversion is present. Also for this case we find good agreement between our analytical and numerical results as is shown in Figure 5 for a three-dimensional two-band GaAs model with initial quasi-equilibrium Fermi-Dirac carrier distributions corresponding to a carrier density of $N_{eh}$, which is equal for electrons and holes, with a temperature of $T=300$ K. Due to the strongly different effective masses of the conduction and the heavy-hole valence bands, $\alpha^{(2)}$ drops steeply to 0 due to the filling of the conduction band near the $\Gamma$ point until $N_{eh}\lesssim 5\cdot 10^{18}$/cm3. Then $\alpha^{(2)}$ slowly becomes negative and the 2PA coefficient approaches the negative of the value without initial carriers present. The plot also highlights the cancellation of the terms depending on the $k$-gradient of the carrier density in the two susceptibility paths $iie$ and $iei$ as well as in the two conductivity paths $eei$ and $eie$.

In general, we see that the dynamic response for quite short pulses is in the considered frequency range and resolution nearly identical to the continuous wave response, apart from minor deviations arising from the finite spectral width of the pulses.

III.2 Higher-order response and multi-photon absorption

We now discuss the probe absorption in the presence of more intense pump pulses where higher-order processes are relevant. The absorption is again calculated using the perturbative equations (5)-(8) including also an expansion in the relevant directions. This procedure is carried out both numerically for short pulses and analytically for continuous-wave excitation. Figure 6 shows spectra separately for the three orders $m=3$, $5$, and $7$, as function of the pump frequency $\omega_{2}$ for a fixed probe frequency of $\hbar\omega_{1}=1.2$ eV for the one-dimensional tight-binding model given in Eq. (12). We use the sum frequency $\omega_{1}+\omega_{2}$ as a common axis and indicate MPA resonances with the band gap by vertical thin dashed lines.

In third order, see Figure 6(a), the TPA vanishes below the band gap except the tail arising from the finite spectral width of the pulses. Above the gap the numerically calculated spectrum agrees well with the analytical result of Eq. (25) shown by the dashed line (scaled). After an initial increase which is determined by the group velocity $|v_{k_{\mathrm{r}}}|$ in the parabolic regime, Eq. (26), the absorption decreases with further increasing sum frequency $\omega_{1}+\omega_{2}$. This decrease arises from the decrease of the ratio $\omega_{1}/\omega_{2}$, i.e., we move towards a more degenerate situation for which according to Fig. 3 the TPA is minimal.

In fifth order, see Figure 6(b), 3PA starts at $\omega_{1}+\omega_{2}=1.35$ eV$/\hbar$, i.e., when $\hbar\omega_{1}+2\hbar\omega_{2}$ is equal to the band gap. Significant 3PA is present mainly in two separate frequency ranges and nearly vanishes in between (here for $\omega_{1}+2\omega_{2}\approx 1.54$ eV$/\hbar$). For the used value of the dipole matrix element $\mu=3e\,$Å the paths with the lowest power of $\mu$ are dominant. According to Figure 1 three paths are quadratic in $\mu$ when the sample is initially unexcited, while all other paths are either of higher order in $\mu$ or are irrelevant. Solving the perturbative expansion for the three dominant paths analytically, we find the following contribution to the 3PA coefficient in the continuous-wave limit for an initially unexcited one-dimensional system

[TABLE]

In the derivation of these expressions we assumed that the sum frequency $\omega_{1}+\omega_{2}$ is off-resonant and therefore the fifth-order corrections to 2PA are not included in this result. The frequency tail of the solid line at $\omega_{1}+2\omega_{2}\lesssim\omega_{\mathrm{g}}$ which represent the numerical results for pulsed excitation is not included either, since it is again due to the finite spectral width of the pulses. Otherwise there is again very good agreement with the numerical result. Indeed, the analytical solution (III.2) has a zero of order 6 (for a parabolic band) in the relevant frequency range. In the degenerate case $\omega_{1}=\omega_{2}$, the zero is found at $\omega_{1}=2\omega_{k}^{\prime}{}^{2}/\omega_{k}^{\prime\prime}$.

The difference between the numerical and the analytical results at $\omega_{1}+\omega_{2}\gtrsim\omega_{\mathrm{g}}$ is due to the abovementioned fifth-order corrections to 2PA. Even though its shape is not dispersive it results, when added to the TPA, in a reduction and an optical Stark shift of the 2PA near the band gap to higher frequencies with increasing pump intensity. In the next order $m=7$, see Figure 6(c), further corrections occur to both 2PA and 3PA. For the chosen frequency values, the Stark-shift correction of the 3PA peak at $\omega_{1}+2\omega_{2}\gtrsim\omega_{\mathrm{g}}$ to higher frequencies is dominant, while there is no significant correction to the broad 3PA range above. We would like to point out that the shape of the higher-order corrections to 2PA and 3PA roughly correspond to a derivative of the lower-order absorption coefficient with respect to the probe frequency. In other words, the correction can be written as

[TABLE]

where $G(\omega_{1},\omega_{2})$ varies slowly over the frequency range of the feature, presumably since it contains nonresonant factors only. Here we do not discuss analytical Stark shift expressions to confirm this relation for a crystal. However, consider a two-level system with transition energy $\omega_{0}$ excited by two beams where $\omega_{1}$ is near-resonant and $\omega_{2}$ is detuned. Then

[TABLE]

and (29) holds with $G(\omega_{1},\omega_{2})=-[4l(\omega_{1}+\omega_{2}+i\gamma_{n})]^{-1}$ (the term $\propto\partial G/\partial\omega_{1}$ is negligible). As long as the features are confined to a narrow frequency range the derivative alone fits our numerical results surprisingly well. The corresponding curves, scaled separately for each feature, are shown by the dotted lines in Figure 6(b)-(c).

For resonant excitation the absorption changes as measured in pump-probe experiments typically alternate with the order, e.g., as demonstrated in Ref. Meier et al., 2000 the fifth-order contribution is proportional to the negative of the third-order result. Interestingly, the relation given in (29), i.e., that the next higher order is proportional to the derivative of the lower order qualitatively also holds for the results shown in Ref. Meier et al., 2000 when exciting with co-circularly polarized pulses below the exciton resonance since the third-order Stark shift becomes a broadening in fifth order.

In Figure 6(c) we also recognize the genuine 4PA with a maximum slightly above $\omega_{1}+3\omega_{2}=\omega_{\mathrm{g}}$, i.e., $\omega_{1}+\omega_{2}=$1.3 eV, which is restricted to a quite narrow frequency range. Finally, we would like to note that further calculations have shown that the MPA coefficients have a similar strong enhancement for strongly non-degenerate excitation as the 2PA coefficient shown in Figure 3. However, the exact dependence on the frequency ratio obviously differs, in particular, the minima occur for $\omega_{1}<\omega_{2}$.

Let us now discuss the influence of a time delay between the pump and the probe pulse on the MPA. In bulk GaAs 2PA transients have been measured in Ref. Penzkofer and Bugayev, 1989. Evaluating the denominator of Eq. (15) for two Gaussian pulses, one finds that the $l$PA depends on the time delay $t_{1}-t_{2}$ as

[TABLE]

Note that for $l\gg 1$ the temporal width is mainly determined by the duration of the probe pulse. Figure 7(a) shows that this scaling is well fulfilled for $\tau_{1}=\tau_{2}=50$ fs for the three cases $l=2,3,4$.

Another influence to be discussed for MPA is that of an initial carrier density $n_{q\mathbf{k}}^{\mathrm{(0)}}$. Using quasi-equilibrium distributions we find that the entire MPA spectra of Figure 6 scale very well with the band filling factor $(1-n_{\mathbf{k}}^{\mathrm{(0)}})$. It is thus excepted that in the continuous wave limit of MPA (using $\gamma_{p}\to 0$ and $\gamma_{n}\to 0$) this factor describes the entire density dependence as for 1PA and 2PA, at least as long as a perturbative treatment is valid.

III.3 Non-perturbative regime with high pump intensities

In this section we present results for the probe absorption obtained from the set of SBE (19)-(23) which is applicable for arbitrarily strong pump fields. In this regime, the linear probe absorption shows a complex dependence on both the time delay $t_{1}-t_{2}$, see Fig. 7(b), and the pump field amplitude $\hat{E}_{2}$, see Fig. 8). As we include the probe field to first order, the absorption of probe photons by the simultaneous absorption of one (or more) pump photons does not cause saturation effects itself. The irregular characteristics reflects the population dynamics induced by the strong pump pulse, here mainly in the region of the Brillouin zone that is resonant with $\omega_{1}+\omega_{2}$. Therefore, the pump absorption is shown on a different scale in Figure 8 (analogous to Figure 2). Roughly-speaking, the increase of the pump absorption starting at $\hat{E}_{2}\approx 0.15$V$\cdot 2\pi/a$ results in a decrease of the probe absorption. More precisely, the dips in the probe absorption are caused by high pump-induced carrier densities $n_{\mathbf{k}}(t)$, where $\mathbf{k}$ is in the $\omega_{1}+(l-1)\omega_{2}$-resonant region of the Brillouin zone (where $l\geq 2$) and $t$ is near the peak time $t_{1}$. For instance, the curve for $\hat{E}_{2}=0.06$V$\cdot 2\pi/a$ in Figure 7 is flattened compared to a Gaussian, and correspondingly a dip is visible near this $\hat{E}_{2}$-value in Figure 8.

The data in Fig. 8 is again obtained for the 1D tight-binding model, (12), since the entire Brillouin zone becomes important for strong fields. For weak fields, the absorption completely agrees with the perturbative regime (the seventh order expansion is accurate up to roughly $\hat{E}_{2}\approx 0.04$V$\cdot 2\pi/a$). For very strong pump intensities (near the right edge of Figure 8), Eqs. (19)-(23) are not accurate anymore, since signals in directions different from the probe direction are mixed into the calculated absorption. The proper directional expansion is shown in a small region of the plot in Figure 8 only since it is computationally expensive. Note that the probe absorption may take positive or negative values in the regime of very intense pump fields. Overall, the non-monotonous intensity dependence is related to Rabi oscillations which occur in regions of the Brillouin zone which are resonantly excited by MPA and the negative regions can be understood to arise from (transient) population inversion. In fact, the origin and the non-monotonous behavior and possible sign changes of the absorption are similar to effects that have been predicted for two-color photocurrents.Pasenow et al. (2008)

IV Summary

Using a semiclassical approach based on the semiconductor Bloch equations, we have studied nonlinear optical absorption from one or two light beams within a two-band model. In general, the presence of far off-resonant light waves makes it necessary to consider both inter- and intraband excitation effects. In the length gauge of the electromagnetic field both types take a simple form which facilitates perturbative analytical solutions of the Bloch equations. The imaginary (real) part of the obtained susceptibility (conductivity) together determine the MPA coefficient for the steady state response. As long as the order of the perturbation expansion is sufficiently high, the MPA coefficients are in good agreement with absorption spectra obtained numerically for pulsed excitation.

Our approach, however, is not limited to the perturbative regime. We also demonstrate how to obtain the absorption of a single strong pulse or a weak probe pulse in the presence of an intense pump pulse in a non-perturbative fashion and demonstrate, e.g., a strongly non-monotonous behavior of the intensity-dependent optical absorption which also appears in the transient response.

The simplicity of the employed model enables us to analyze and point out several novel fundamental aspects. Future extensions of this model might lead to a more realistic description. The evaluation of nonlinear optical susceptibilities in more sophisticated band structure models including more than two bands and the $k$-dependence as well as the vector character of the matrix elements would lead to more quantitative predictions and would also allow to study the complete tensor character of the nonlinear susceptibilities. Furthermore, the inclusion of electron-electron-interactions is expected to lead to an excitonic enhancement of MPA appearing spectrally below resonances originating from the fundamental band gap.

Acknowledgements.

We thank John E. Sipe for several useful discussions. Support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 231447078 – TRR 142 (project A07) is gratefully acknowledged. We also thank the PC2 (Paderborn Center for Parallel Computing) for providing computing time.

Appendix A Path evolution equations

Splitting the evolution equations (5)-(8) into individual paths of interband and intraband terms, the series takes the following form in the first two orders:

[TABLE]

where the path labels are given in square brackets. The higher orders can be obtained analogously.

Appendix B Third-order response functions

Here we give explicit expressions for the functions $\phi_{\mathbf{k}}^{(\mathrm{path};\ell)}(\omega_{1},\omega_{2})$ defining the third-order susceptibility and conductivity in the continuous-wave limit according to Eq. (24). We define primed quantities as derivatives with respect to $k$ along the field direction, e.g., $\omega^{\prime}_{\mathbf{k}}=\hat{\mathbf{e}}_{1}\cdot\nabla_{\mathbf{k}}\,\omega_{\mathbf{k}}$, and similarly for $n_{\mathbf{k}}$. The subpath index $\ell$ runs from 1 to 3 and we define $\omega_{3}\equiv\omega_{2}$, $\Omega_{1}=\Omega_{2}\equiv\omega_{1}+\omega_{2}$, and $\Omega_{3}\equiv 0$. Note that the occurrences of these frequency expressions are not invariable, but change according to the arguments of $\phi_{\mathbf{k}}^{(\mathrm{path};\ell)}(\omega_{1},\omega_{2})$. For simplicity we neglect all terms which do not feature a resonance denominator $\omega_{1}-\bar{\omega}_{\mathbf{k}}$ or $\omega_{1}+\omega_{2}-\bar{\omega}_{\mathbf{k}}$. This approximation is sufficient for considering 2PA, electronic Raman effect (where $\omega_{2}\to-\omega_{2}$ according to $\eta=-$ in Eq. (24)), and ac Stark shift.

For the four susceptibility paths we find:

[TABLE]

The solutions for the four conductivity paths read:

[TABLE]

The paths contributing to 2PA (or Raman for $\eta=-$) are easily identified by the occurrences of the resonance denominator $\omega_{1}+\omega_{2}-\bar{\omega}_{\mathbf{k}}$. All paths lacking this term do not contribute to $\alpha^{(2)}$ (despite unphysical divergences for $\gamma_{n}\to 0$). The contributing paths are further evaluated by carrying out the $k$-integral from Eq. (24) and taking $\gamma_{n}\to 0$ and $\gamma_{p}\to 0$. This limit may be immediately carried out in all cases except at the resonance poles. According to Eq. (18), $\alpha^{(2)}$ is determined by the imaginary part of the poles,

[TABLE]

As discussed in the main text, we restrict the further evaluation to the one-dimensional case. We transform the $k$-integral in Eq. (24) into frequency space:

[TABLE]

where the sum is over the resonant parts of the spectrum. Then we use the asymptotic expressions

[TABLE]

which leads to the result given in Eq. (25).