Abstract
Implicit summation is a technique for the conversion of sums over intermediate states in multiphoton absorption and the highorder susceptibility in hydrogen into simple integrals. Here, we derive the equivalent technique for hydrogenic impurities in multivalley semiconductors. While the absorption has useful applications, it is primarily a loss process; conversely, the nonlinear susceptibility is a crucial parameter for active photonic devices. For Si:P, we predict the hyperpolarizability ranges from χ^{(3)}/n_{3D} = 2.9 to 580 × 10^{−38} m^{5}/V^{2} depending on the frequency, even while avoiding resonance. Using samples of a reasonable density, n_{3D}, and thickness, L, to produce thirdharmonic generation at 9 THz, a frequency that is difficult to produce with existing solidstate sources, we predict that χ^{(3)} should exceed that of bulk InSb and χ^{(3)}L should exceed that of graphene and resonantly enhanced quantum wells.
Introduction
Multiphoton absorption requires a high intensity, and was first observed shortly after the invention of the laser using impurities in solids^{1} and alkali vapor^{2}. Although multiphoton absorption is useful for metrology and modulators, and can be enhanced where there is nearresonance of an intermediate state as in the case of Rb^{3}, it is essentially a loss process contributing an imaginary part to the nonlinear susceptibility. The corresponding real part is responsible for a great variety of wavelength conversion processes such as harmonic generation, first observed in quartz^{4} and later in atomic vapors^{5} including alkalies^{6}. THz multiphoton absorption has been shown to be very large in hydrogenic shallow impurities in semiconductors, even without intermediate state resonances^{7}, due to the large dielectric screening and low effective mass. Here, we predict giant values for the real part of the THz nonlinear susceptibility for doped silicon and germanium. This finding opens access to novel applications for these materials in THz photonics. For example, tripling the output of a 2–4 THz quantum cascade laser through thirdharmonic generation would fill the frequency gap currently only filled by larger, more expensive systems. We show that a good efficiency can be obtained for thirdharmonic generation with doped silicon and germanium. Our theory can be readily applied to any donor in any semiconductor host where the effective mass approximation is valid, and our discussion makes it clear that a giant value of χ^{(3)} is expected for donors with a small binding energy in a host with a large dielectric constant and small effective mass.
The theory developed in this paper is appropriate for frequencies both near to and far from lossinducing resonances, including the effects of effective mass anisotropy, multivalley interactions and the central cell correction. The method could easily be applied to other systems with complicated potentials, such as multiquantum wells. Although this work focuses on perturbative harmonic generation, we anticipate that shallow impurities may also be useful for nonperturbative highharmonic generation (HHG)^{8,9} taking advantage of the excellent control over the carrierenvelope phase of fewcycle pulses in this THz regime, which can be used to enhance HHG^{10}.
Results
The implicit summation technique
From Nthorder perturbation theory^{7,11} the Nphoton absorption (NPA) transition rate may be written as
where \(I_a = E_H^2/\hbar a_B^2\), a_{B} is the Bohr radius, E_{H} the Hartree energy, and α_{fs} the fine structure constant. M^{(N)} is a dimensionless transition matrix element, and I_{m} is the intensity of the light in the medium with relative dielectric permittivity ε_{r}. The lineshape function Γ^{(N)}(ω) has unit area. For silicon and germanium donors, the factors inside the bracket are renormalized, and of particular importance here I_{a} is ten orders of magnitude smaller for silicon than it is for hydrogen. This is apparent from the formulae of the Hartree energy and Bohr radius for donors in these materials: \(E_H = m_t(e^2/4\pi \epsilon _0\varepsilon_r \hbar )^2\), and \(a_B = 4\pi \epsilon _0\varepsilon_r \hbar ^2/m_te^2\), where m_{t} is the transverse effective mass and \(\varepsilon_r\) the dielectric constant^{12}. Both germanium and silicon have a small m_{t} and large \(\varepsilon_r\), raising the Bohr radius and lowering the binding energy. The wavefunction is therefore significantly larger than that of alkali atoms, leading to an enhanced dipole matrix element and hence a substantially stronger interaction with light.
The details of the spectrum given by Eq. (1) are controlled by M^{(N)}, which is influenced in silicon by the indirect valley structure, the anisotropic effective mass, and the donor central cell correction potential. Our main aim here is to calculate these effects. For singlephoton absorption (N = 1) between states ψ_{g}〉 (the ground state) and ψ_{e}〉 (the excited state), \(M^{(1)} = \left\langle {\psi _e\left {{\boldsymbol{\epsilon }}.{\boldsymbol{r}}} \right\psi _g} \right\rangle /a_B\), where \({\bf{\epsilon }}\) is a unit vector in the polarization direction, and Eq. (1) reduces to Fermi’s golden rule. For twophoton absorption,
in the E.r gauge, which may be written as M^{(2)} = 〈ψ_{e}ζG_{1}ζψ_{g}〉 where \(\zeta = {\boldsymbol{\epsilon }}.{\boldsymbol{r}}/a_B\),
and ω = ω_{eg}/N. The states j〉 are intermediate states, and along with ψ_{e}〉 & ψ_{g}〉 they are eigenstates of \(H\left j \right\rangle = \hbar \omega _j\left j \right\rangle\), where H is the Hamiltonian in the dark. For general multiphoton absorption,
The summation in Eq. (2) can be avoided^{11} by noticing that (H − W_{n})G_{n} = E_{H}, where W_{n} = ℏω_{g} + nℏω, and ω = ω_{eg}/N as already mentioned, and by using the completeness relation \(\mathop {\sum}\nolimits_j {\left j \right\rangle \left\langle j \right} = 1\). In other words,
Rewriting Eq. (3), M^{(N)} = 〈ψ_{e}ζψ_{N−1}〉 where ψ_{0}〉 = ψ_{g}〉 and ψ_{n}〉 is the solution of the partial differential equation (PDE) \(G_n^{  1}\left {\psi _n} \right\rangle = \zeta \left {\psi _{n  1}} \right\rangle\). Instead of finding M^{(N)} by repeated application of Eq. (2), which requires infinite sums (that might be reduced down to a few terms if there are obvious resonances), we may now use Eq. (4) and the PDE at each stage, which can be simpler.
The Nthorder susceptibility far from any multiphoton resonances may also be calculated using the Nthorder perturbation theory^{13}. For example, the “resonant” term in the thirdorder susceptibility, χ^{(3)}(3ω), is
where e is the electron charge, and n_{3D} is the concentration. χ^{(3)} may be written in a similar form to Eqs (1) and (3), and for N^{th} order,
where C^{(N)} = 〈ψ_{g}ζG_{N}…G_{2}ζG_{1}ζψ_{g}〉 is a dimensionless matrix element that may be found in a similar way to M^{(N)}, either by repeated application of Eq. (2)—as has been done previously for alkali metal vapors^{6}—or by using the implicit summation method of Eq. (4) with the only difference being ω ≠ ω_{eg}/N. The antiresonant terms^{13} and other nonlinear processes, such as sumfrequency generation, can be calculated with simple modifications to W_{n} at each step.
Multivalley theory for donors in silicon and germanium
In this section, we develop the multivalley theory for the nonlinear optical processes of donors based on the effective mass approximation (EMA). For simplicity of presentation, we describe the derivation for silicon; the case of germanium is discussed in the Supplementary Materials. It will become apparent that our theory is readily applicable to any donor in any host as long as the EMA is reliable.
To apply the method to donors, we require ψ_{g}〉, ω_{g}, ψ_{e}〉, ω_{e} and Hψ_{n}〉. Silicon and germanium are indirect with equivalent conduction band minima (valleys) near the Brillouin zone edge; each minimum is characterized by a Fermi surface that is a prolate ellipsoid with transverse & longitudinal effective masses, m_{t,l}. According to the KohnLuttinger effective mass approximation^{14}, the state ψ_{j}〉 of a shallow donor can be decomposed into slowly varying hydrogenic envelope functions, one for each valley, modulated by planewave functions corresponding to the crystal momenta at the minima, k_{μ} (and a lattice periodic function that is unimportant here). We write \(\psi _j({\boldsymbol{r}}) = \mathop {\sum}\nolimits_\mu {e^{i{\boldsymbol{k}}_\mu .{\boldsymbol{r}}}F_{j,\mu }({\boldsymbol{r}})}\) where F_{j,μ}(r) is the slowly varying envelope function. We have neglected the lattice periodic part, u_{μ}(r), of the Bloch functions for the simplicity of presentation. A rigorous derivation with u_{μ}(r) included is provided in the Supplementary Materials, but it does not lead to any change in the final equations for the envelope functions (Eqs (7) and (8) below).
We separate the potential into the slowly varying Coulomb term of the donor V(r), and a rapidly varying term due to the quantum defect that is short range, U(r), referred to as the central cell correction (CCC). Within the EMA, the kinetic energy term in the Hamiltonian operates only on the envelope function, and the EMA Schrodinger equation may be written as
where H_{0} includes the Coulomb potential V(r): \(E_H^{  1}H_0 =  \frac{1}{2}a_B^2\left[ {\partial _x^2 + \partial _y^2 + \gamma \partial _z^2} \right]  a_Br^{  1}\) using a valleyspecific coordinate system (x, y, z where z is the valley axis, i.e., the valleyframe is rotated relative to the labframe of x_{1}, x_{2}, x_{3}). The kinetic energy has cylindrical symmetry because γ = m_{t}/m_{l} ≠ 1, and V(r) and U(r) are spherical and tetrahedral respectively. H_{0} produces wave functions that are approximately hydrogenlike, and U(r) mixes them to produce states that transform as the A_{1}, E and T_{2} components of the T_{d} point group.
We take U(r) to be very short range, and we neglect the small change in the envelope functions over the short length scale 2π/k_{μ}. Premultiplying Eq. (6) by \(e^{  i{\boldsymbol{k}}_{\mu \prime }.{\boldsymbol{r}}}\) and averaging over a volume (2π/k_{μ})^{3} around r, the Schrodinger eqn now reads \(\left[ {H_0  \hbar \omega _j} \right]F_{j,\mu }({\boldsymbol{r}}) + \mathop {\sum}\nolimits_{\mu \prime } {U_{\mu \mu \prime }\delta ({\boldsymbol{r}})F_{j,\mu \prime }({\boldsymbol{r}})} = 0\), where δ(r) is the Dirac delta function, and \(U_{\mu \mu \prime } = {\int} {d{\boldsymbol{r}}{\kern 1pt} e^{i({\boldsymbol{k}}_{\mu \prime }  {\boldsymbol{k}}_\mu ).{\boldsymbol{r}}}U({\boldsymbol{r}})}\). For an A_{1} state, all the envelope functions have the same amplitude at r = 0, hence, \(\mathop {\sum}\nolimits_{\mu \prime } {U_{\mu \mu \prime }\delta ({\boldsymbol{r}})F_{j,\mu \prime }({\boldsymbol{r}})} =  U_{cc}\delta ({\boldsymbol{r}})F_{j,\mu }({\boldsymbol{r}})\), where \(U_{cc} =  \mathop {\sum}\nolimits_{\mu \prime } {U_{\mu \mu \prime }}\). It is found experimentally that for E and T_{2} states, the CCC has a rather small effect, and so we neglect it. Since H_{0} has cylindrical symmetry, the component of angular momentum about the valley axis is a conserved quantity, i.e., F_{j,μ}(r) = e^{imϕ}f_{j,m,μ}(r, θ), where m is a good quantum number, and now f_{j,m,μ} is a 2D function only. Substituting into the Schrodinger eqn, premultiplying by e^{−im′ϕ} and finally integrating over ϕ, the eigenproblems are
where \(H_0^{(m)} = H_0 + E_Ha_B^2m^2\)/2(r sin θ)^{2}. We solve Eq. (7) using a 2D finite element method (FEM) (see Supplementary Materials).
We focus on silicon, in which case the valley index, μ, runs over (±1, ±2, ±3), where 1, 2, 3 are the three crystal axes, and we let the light be polarized along a crystal axis, x_{1}, by way of illustration; the calculation for germanium and other polarization directions is described in the Supplementary Materials. For the μ = ±1, ±2, ±3 valleys, a_{B}ζ_{μ} = z, x, y = r cos θ, r sin θ cos ϕ, r sin θ sin ϕ, respectively, because each has its coordinate rotated so that z is the valley axis. Following the expansion of ψ_{j} in terms of the f_{j,m,μ}, we write the intermediate state functions as \(\psi _n({\boldsymbol{r}}) = \mathop {\sum}\nolimits_{m,\mu } {e^{im\phi }e^{i{\boldsymbol{k}}_\mu .{\boldsymbol{r}}}f_{n,m,\mu }(r,\theta )}\), substitute them into \(G_n^{  1}\psi _n = \zeta \psi _{n  1}\), premultiply by \(e^{  i{\boldsymbol{k}}_{\mu \prime }.{\boldsymbol{r}}}\), average over a volume of (2π/k_{μ})^{3}, premultiply by e^{−im′ϕ}, and finally, integrate over ϕ. Since f_{0,0,μ} = f_{g,0,μ} for all μ, we find that f_{n,m,3} = i^{−m}f_{n,m,2} and f_{n,m,−μ} = f_{n,m,μ}, and
where \({\cal{D}} = U_{cc}\delta ({\boldsymbol{r}})\delta _{m,0}/3\) and δ_{m,0} is the Kronecker delta. In the above equations we drop the valleyspecific coordinates in f_{n,m,μ} for notational simplicity, and the coordinates in \(H_0^{(m)}\) and the right hand side are understood to belong to the valley of the envelope function that they act on.
It is evident that Eq. (8) are not coupled by U_{cc} when the envelope function is zero at the origin. The ground state ψ_{0}〉 = ψ_{g}〉 has only m = 0 components, and it has even parity. Therefore, ψ_{1}〉 has odd parity according to Eq. (8), so the U_{cc} coupling term is suppressed. By the same logic, the U_{cc} coupling is only nonzero for even n and m = 0. In the case of \(\left {f_{n,m,1}} \right\rangle\), there is only dipole coupling to the functions with the same m, while for \(\left {f_{n,m,2}} \right\rangle\) the dipole coupling is to states with Δm = ±1. The latter couplings are identical, so f_{n,−m,μ} = f_{n,m,μ}. Figure 1 shows how the intermediate states are coupled by dipole excitation and the CCC.
Equation (8) can be solved by sequential application of the 2D FEM^{15}. To test our numerical calculation we first compute C^{(3)} for hydrogen, and each of the resonant and antiresonant terms is shown in Fig. 2. Their sum is shown in Fig. 3, and we find excellent agreement within 0.2% of the previous result obtained from a Sturmian coulomb Green function in ref. ^{16}.
Discussion
Giant thirdorder nonlinear susceptibility
Since silicon and germanium donors have an isotropic potential in an isotropic dielectric, the lowestorder nonlinear response is determined by χ^{(3)}. The χ^{(3)} spectrum for each (including the antiresonant terms) is shown in Fig. 3. We took the parameters for silicon obtained from spectroscopic^{17} and magnetooptical measurements^{12,18}, which are γ ≈ 0.208, a_{B} ≈ 3.17 nm and E_{H} ≈ 39.9 meV. The parameters for germanium are γ ≈ 0.0513, a_{B} ≈ 9.97 nm and E_{H} ≈ 9.40 meV^{19}. Resonances occur when 3ω = ω_{eg}, labeled according to ψ_{e}〉, and there are also signchanges at which χ^{(3)} goes to zero. In the range of frequency shown, we also observe a twophoton resonance for 1sA_{1} → 1sE, which is an obvious illustration of the need for a multivalley theory. There is no 3ω resonance with 1sT_{2} within the approximations made above in which there is no intervalley dipole coupling. The effect of U_{cc} on χ^{(3)} and the NPA matrix element is shown in Fig. 4. The lowfrequency response of C^{(3)} is illustrated at 100 GHz. Two higherfrequency curves are included, with both far from 3ω resonances, half way between the 2p_{0} and 2p_{±} resonances, and between the 3p_{0} and 3p_{±}. We choose these average frequencies since χ^{(3)} for Si:P varies slowly around them (see Fig. 3) and hence would not be sensitive to small experimental variations in the light frequency. For the 2paverage frequency, the 2ω resonance with the 1sE produces a coincidental zerocrossing for Si:Bi. Example results for the intermediate state wave functions produced in the calculation are shown in Fig. 5. The state ψ_{2}〉 is much larger in extent (and in magnitude) than ψ_{0}〉, and the extra node in the radial dependence due to the contribution of 2s is visible at about 5 nm. Similarly, the state ψ_{3}〉 is much larger in extent (and in magnitude) than ψ_{1}〉.
The square bracket in Eq. (5) gives the scaling of χ^{(N)} from hydrogenic atoms in vacuum to hydrogenic impurities in semiconductors, just as that in Eq. (1) does for w^{(N)}, and as before, the much smaller I_{a} greatly increases the strength of the nonlinearity. For example, the lowfrequency limit of the hyperpolarizability χ^{(3)}/n_{3D} for Si:P is much larger than that for hydrogen or alkali metal vapors such as Rb^{6}, as shown in Fig. 3.
Some of the highest values of χ^{(3)} have been reported for solids, e.g., 2.8 × 10^{−15} m^{2}/V^{2} for InSb^{20} and 2 × 10^{−16} m^{2}/V^{2} for GaTe^{21}. To convert the hyperpolarizability to a bulk χ^{(3)} value requires the concentration. To match InSb with Si:P at low frequency where C^{(3)} ≈ 1 (Fig. 4) (and χ^{(3)}/n_{3D} = 2.9 × 10^{−38} m^{5}/V^{2}) requires a donor density of n_{3D} = 10^{17} cm^{−3} (where the donor–donor distance is 10a_{B}). At high frequency, the hyperpolarizability is much higher, but the density should be lower to avoid inhomogeneous concentration broadening of the nearby excited levels. For example, C^{(3)} ≈ 20 between the 2p_{0} and 2p_{±} resonances at \(\omega = \bar \omega _{2p}/3 = 2\pi \times 3.2\,{\text{THz}}\) (Fig. 4), and we match InSb at a density of n_{3D} = 5 × 10^{15} cm^{−3} at which concentration the 2p lines are well resolved^{22}. If 3ω is moved even closer to the 2p_{±} resonance (or if the resonance is tuned with a magnetic field^{18}), then χ^{(3)} could easily exceed InSb. Losses due to dephasing by phonon scattering may become important if the time spent in the intermediate states exceeds the phonon lifetime. Since the inverse of the former is given approximately by the detuning (ΔfΔt ≥ 1/2π) and the inverse phononlimited width (1/πT_{2} = 1 GHz^{23,24}), then this loss is negligible for much of the spectrum. At 50 GHz below the 2p_{±} line so that such losses may be ignored, C^{(3)} ≈ 200, and χ^{(3)} is an order of magnitude above InSb.
We are not aware of any larger values for bulk media, but higher “bulk” values have been reported for 2D systems such as graphene and MoS_{2} for which χ^{(3)}L data are divided by an interaction thickness L to obtain χ^{(3)}; in particular, reports for graphene range from 10^{−19} ^{25,26} to 10^{−15} m^{2}/V^{2} ^{27} for nearIR excitation and up to 10^{−10} m^{2}/V^{2} in the THz region under resonant enhancement by landau levels in a magnetic field^{28}. A recent experiment with singlelayer graphene at room temperature reports a remarkably high value of 1.7 × 10^{−9} m^{2}/V^{2} for the THz thirdorder nonlinear susceptibility^{29}. In the case of coupled quantum wells (QW), large values of χ^{(3)} may be engineered through resonances, as demonstrated up to 10^{−14} m^{2}/V^{2} ^{30}. However, since the nonlinear effect is limited by the interaction length, the 2D χ^{(3)}L is probably a better figure of merit in these cases. For THz fieldenhanced graphene with 50 layers, χ^{(3)}L = 9 × 10^{−20} m^{3}/V^{2 }^{28}, and for singlelayer graphene χ^{(3)}L = 5.1 × 10^{−19} m^{3}/V^{2} ^{29}, or χ^{(3)}L = 1.4 × 10^{−18} m^{3}/V^{2} for resonant coupled QWs^{30}. Even higher values are predicted for doped QWs up to χ^{(3)}L = 5 × 10^{−17} m^{3}/V^{2 }^{31}. To match this value with Si:P at \(\omega = \bar \omega _{2p}/3 = 2\pi \times 3.2\,{\text{THz}}\) and n_{3D} = 5 × 10^{15} cm^{−3} (see above) would require a sample thickness of L = 2 cm. Obviously, the required thickness can be significantly reduced when close to resonance, or for germanium.
Efficient thirdharmonic generation
The nonlinear susceptibility is important for predicting the strength of frequency conversion processes such as thirdharmonic generation (3HG), and we use this as an example application to investigate the utility of the medium. A solution for the amplitude of the generated wave produced by 3HG, neglecting absorption, is given by^{32}. Converting to irradiance in MKS units,
where I_{in} is the irradiance of the input pump wave at frequency f_{in}, and n is the geometric mean of the refractive indexes for the input and output waves, and n_{2D} = n_{3D}L. Note that the isotropy mentioned earlier means that the polarization of the input and output waves must be parallel. We ignored a factor for the phase matching, which is unity if the length of the sample L ≪ L_{c}, where the coherence length L_{c} = πc/(3ω_{in}[n_{out} − n_{in}]). Si:P at room temperature has a nearly constant n = 3.4153 in the range from 1 THz to 12 THz^{33}, leading to typical values of L_{c} ≈ 10 cm. The factor x = 6.9 × 10^{23} W/cm^{2} × THz × cm^{−2} for silicon. For comparison, germanium has x = 9.2 × 10^{19} W/cm^{2} × THz × cm^{−2}.
To illustrate the possible applications of this high χ^{(N)}, we note that two types of THz diode lasers are available, the quantum cascade laser (QCL) from 0.74 THz^{34} to 5.4 THz^{35} with output powers of up to a few W^{36,37}, and the hot hole (pGe) laser^{38,39} with a similar range and power. However, there is a large gap in the availability of solidstate sources from about 5 THz to about 12 THz^{40}, where the GaAs Reststrahlen band renders laser operation impossible. This is an important region for quantum qubit applications^{41,42,43,44}. Currently, the gap is only filled by larger, more expensive systems (difference frequency generators and free electron lasers). Tripling the output of 2–4 THz QCLs would fill the gap, but their output powers are far smaller than those typical for a pump laser in standard tripling applications, so a giant nonlinearity is critical. At \(\omega = \bar \omega _{2p}/3 = 2\pi \times 3.2\,{\text{THz}}\), C^{(3)} ≈ 20, so for n_{2D} = 10^{16} cm^{−2} (see above), a 1% predicted conversion may be obtained with 100 kW/cm^{2}, and by moving to 50 GHz below the 2p_{±} resonance this value could be brought down to 10 kW/cm^{2}, which is just about achievable with a well focussed QCL, and would thus provide enough output for spectroscopy applications. A nonlinear process that may possibly reduce the 3HG efficiency is multiphoton ionization^{45} since it reduces the population of the donors in the ground state. When \(\omega = \bar \omega _{2p}/3\), for example, a fourphoton absorption takes the electron to the continuum. We estimate this ionization in Si:P using the implicit summation method and find that the rate is w = 3.17 s^{−1} for I_{in} = 10 kW/cm^{2}. This result simply means that the pulses must be kept significantly shorter than a second to avoid significant ionization.
In summary, we calculated the absolute values of the THz nonlinear coefficients for the most common semiconductor materials, lightly doped silicon and germanium, which are available in the largest, purest and most regular single crystals known. The values we obtain for offresonance rival the highest values obtained in any other material even when resonantly enhanced, and the material could gain new applications in THz photonics. We also predict the highly efficient thirdharmonic generation of THz light in doped silicon and germanium. Our multivalley theory for nonlinear optical processes of donors in silicon and germanium can be readily applied to any donor in any semiconductor host in which the effective mass approximation is reliable.
Materials and methods
Details of the finite element computation used for solving the coupled partial differential equations (Eq. (8)) are provided in the Supplementary Material.
Data availability
Data for Nguyen Le et al. Giant nonlinear susceptibility of hydrogenic donors in silicon and germanium, https://doi.org/10.5281/zenodo.3269481. The data underlying this work is available without restriction.
References
 1.
Kaiser, W. & Garrett, C. G. B. Twophoton excitation in CaF2: Eu2+. Phys. Rev. Lett. 7, 229–231 (1961).
 2.
Abella, I. D. Optical doublephoton absorption in cesium vapor. Phys. Rev. Lett. 9, 453–455 (1962).
 3.
Saha, K. et al. Enhanced twophoton absorption in a hollowcore photonicbandgap fiber. Phys. Rev. A 83, 033833 (2011).
 4.
Franken, P. A. et al. Generation of optical harmonics. Phys. Rev. Lett. 7, 118–119 (1961).
 5.
Ward, J. F. & New, G. H. C. Optical third harmonic generation in gases by a focused laser beam. Phys. Rev. 185, 57–72 (1969).
 6.
Miles, R. & Harris, S. Optical thirdharmonic generation in alkali metal vapors. IEEE J. Quantum Electron. 9, 470–484 (1973).
 7.
Van Loon, M. A. W. et al. Giant multiphoton absorption for THz resonances in silicon hydrogenic donors. Nat. Photonics 12, 179–184 (2018).
 8.
Vampa, G. et al. Theoretical analysis of highharmonic generation in solids. Phys. Rev. Lett. 113, 073901 (2014).
 9.
Beaulieu, S. et al. Role of excited states in highorder harmonic generation. Phys. Rev. Lett. 117, 203001 (2016).
 10.
Haworth, C. A. et al. Halfcycle cutoffs in harmonic spectra and robust carrierenvelope phase retrieval. Nat. Phys. 3, 52–57 (2007).
 11.
Gontier, Y. & Trahin, M. On the multiphoton absorption in atomic hydrogen. Phys. Lett. A 36, 463–464 (1971).
 12.
Li, J. et al. Radii of Rydberg states of isolated silicon donors. Phys. Rev. B 98, 085423 (2018).
 13.
Boyd, R. W. Nonlinear Optics. 3rd edn, 640. Academic Press, New York (2008).
 14.
Kohn, W. & Luttinger, J. M. Theory of donor states in silicon. Phys. Rev. 98, 915–922 (1955).
 15.
The Mathematica FEM code used in this paper is available at, https://github.com/lehnqt/chi3.git.
 16.
Mizuno, J. Use of the sturmian function for the calculation of the third harmonic generation coefficient of the hydrogen atom. J. Phys. B: At. Mol. Phys. 5, 1149–1154 (1972).
 17.
Ramdas, A. K. & Rodriguez, S. Review article: spectroscopy of the solidstate analogues of the hydrogen atom: donors and acceptors in semiconductors. Rep. Prog. Phys. 44, 1297–1387 (1981).
 18.
Murdin, B. N. et al. Si:P as a laboratory analogue for hydrogen on high magnetic field white dwarf stars. Nat. Commun. 4, 1469 (2013).
 19.
Faulkner, R. A. Higher donor excited states for prolatespheroid conduction bands: a reevaluation of silicon and germanium. Phys. Rev. 184, 713–721 (1969).
 20.
Yuen, S. Y. & Wolff, P. A. Differencefrequency variation of the freecarrierinduced, thirdorder nonlinear susceptibility in nInSb. Appl. Phys. Lett. 40, 457–459 (1982).
 21.
Susoma, J. et al. Second and third harmonic generation in fewlayer gallium telluride characterized by multiphoton microscopy. Appl. Phys. Lett. 108, 073103 (2016).
 22.
Thomas, G. A. et al. Optical study of interacting donors in semiconductors. Phys. Rev. B 23, 5472–5494 (1981).
 23.
Steger, M. et al. Shallow impurity absorption spectroscopy in isotopically enriched silicon. Phys. Rev. B 79, 205210 (2009).
 24.
Greenland, P. T. et al. Coherent control of Rydberg states in silicon. Nature 465, 1057–1061 (2010).
 25.
Woodward, R. I. et al. Characterization of the second and thirdorder nonlinear optical susceptibilities of monolayer MoS2 using multiphoton microscopy. 2D Mater. 4, 11006 (2017).
 26.
Karvonen, L. et al. Rapid visualization of grain boundaries in monolayer MoS2 by multiphoton microscopy. Nat. Commun. 8, 15714 (2017).
 27.
Säynätjoki, A. et al. Rapid largearea multiphoton microscopy for characterization of graphene. ACS Nano 7, 8441–8446 (2013).
 28.
KönigOtto, J. C. et al. Fourwave mixing in landauquantized graphene. Nano Lett. 17, 2184–2188 (2017).
 29.
Hafez, H. A. et al. Extremely efficient terahertz highharmonic generation in graphene by hot Dirac fermions. Nature 561, 507–511 (2018).
 30.
Sirtori, C. et al. Giant, triply resonant, thirdorder nonlinear susceptibility \(\chi _{3\omega }^{(3)}\) in coupled quantum wells. Phys. Rev. Lett. 68, 1010–1013 (1992).
 31.
Yildirim, H. & Aslan, B. Donorrelated thirdorder optical nonlinearites in GaAs/AlGaAs quantum wells at the THz region. Semicond. Sci. Technol. 26, 085017 (2011).
 32.
Shen, Y. R. The Principles of Nonlinear Optics. 3rd edn, 576. WileyInterscience, New York, 2002.
 33.
Chick, S. et al. Metrology of complex refractive index for solids in the terahertz regime using frequency domain spectroscopy. Metrologia 55, 771 (2018).
 34.
Scalari, G. et al. Magnetically assisted quantum cascade laser emitting from 740GHz to 1.4THz. Appl. Phys. Lett. 97, 081110 (2010).
 35.
Wienold, M. et al. Frequency dependence of the maximum operating temperature for quantumcascade lasers up to 5.4 THz. Appl. Phys. Lett. 107, 202101 (2015).
 36.
Li, L. et al. Terahertz quantum cascade lasers with >1 W output powers. Electron. Lett. 50, 4309 (2014).
 37.
Li, L. H. et al. MultiWatt highpower THz frequency quantum cascade lasers. Electron. Lett. 53, 799 (2017).
 38.
Pfeffer, P. et al. ptype Ge cyclotronresonance laser: Theory and experiment. Phys. Rev. B 47, 4522–4531 (1993).
 39.
Hübers, H. W., Pavlov, S. G. & Shastin, V. N. Terahertz lasers based on germanium and silicon. Semicond. Sci. Technol. 20, S211–S221 (2005).
 40.
Ohtani, K., Beck, M. & Faist, J. Double metal waveguide InGaAs/AlInAs quantum cascade lasers emitting at 24 μm. Appl. Phys. Lett. 105, 121115 (2014).
 41.
Saeedi, K. et al. Optical pumping and readout of bismuth hyperfine states in silicon for atomic clock applications. Sci. Rep. 5, 10493 (2015).
 42.
Litvinenko, K. L. et al. Coherent creation and destruction of orbital wavepackets in Si:P with electrical and optical readout. Nat. Commun. 6, 6549 (2015).
 43.
Chick, S. et al. Coherent superpositions of three states for phosphorous donors in silicon prepared using THz radiation. Nat. Commun. 8, 16038 (2017).
 44.
Stoneham, A. M. et al. Letter to the editor: optically driven siliconbased quantum gates with potential for hightemperature operation. J. Phys.: Condens. Matter 15, L447–L451 (2003).
 45.
Bebb, H. B. & Gold, A. Multiphoton ionization of hydrogen and raregas atoms. Phys. Rev. 143, 1–24 (1966).
Acknowledgements
We acknowledge financial support from the UK Engineering and Physical Sciences Research Council [ADDRFSS, Grant No. EP/M009564/1] and EPSRC strategic equipment grant no. EP/L02263X/1.
Author information
Affiliations
Contributions
N.H. Le and B.N. Murdin worked on the multivalley theory and the finite element calculation of the thirdorder susceptibility. B.N. Murdin and G.V. Lanskii calculated the thirdharmonic generation efficiency. B.N. Murdin, N.H. Le and G. Aeppli wrote the manuscript. All authors contributed to the discussion of the results.
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Le, N.H., Lanskii, G.V., Aeppli, G. et al. Giant nonlinear susceptibility of hydrogenic donors in silicon and germanium. Light Sci Appl 8, 64 (2019). https://doi.org/10.1038/s4137701901746
Received:
Revised:
Accepted:
Published:
Further reading

Highly efficient THz fourwave mixing in doped silicon
Light: Science & Applications (2021)

The multiphoton induced Fano effect
Nature Communications (2021)