Giant non-linear susceptibility of hydrogenic donors in silicon and germanium

Implicit summation is a technique for the conversion of sums over intermediate states in multiphoton absorption and the high-order susceptibility in hydrogen into simple integrals. Here, we derive the equivalent technique for hydrogenic impurities in multi-valley semiconductors. While the absorption has useful applications, it is primarily a loss process; conversely, the non-linear susceptibility is a crucial parameter for active photonic devices. For Si:P, we predict the hyperpolarizability ranges from χ(3)/n3D = 2.9 to 580 × 10−38 m5/V2 depending on the frequency, even while avoiding resonance. Using samples of a reasonable density, n3D, and thickness, L, to produce third-harmonic generation at 9 THz, a frequency that is difficult to produce with existing solid-state 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 non-linear 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 non-linear 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 third-harmonic generation would fill the frequency gap currently only filled by larger, more expensive systems. We show that a good efficiency can be obtained for third-harmonic 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 loss-inducing resonances, including the effects of effective mass anisotropy, multi-valley 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 non-perturbative highharmonic generation (HHG) 8,9 taking advantage of the excellent control over the carrier-envelope phase of few-cycle pulses in this THz regime, which can be used to enhance HHG 10 .

The implicit summation technique
From Nth-order perturbation theory 7,11 the N-photon absorption (NPA) transition rate may be written as where , 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: where m t is the transverse effective mass and ε r the dielectric constant 12 . Both germanium and silicon have a small m t and large ε 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 single-photon absorption (N = 1) between states |ψ g 〉 (the ground state) and |ψ e 〉 (the excited state), M ð1Þ ¼ ψ e ϵ:r j jψ g D E =a B , where is a unit vector in the polarization direction, and Eq. (1) reduces to Fermi's golden rule. For two-photon absorption, in the E.r gauge, which may be written as M (2) = 〈ψ e | ζG 1 ζ|ψ g 〉 where ζ ¼ ϵ:r=a B , and ω = ω eg /N. The states |j〉 are intermediate states, and along with |ψ e 〉 & |ψ g 〉 they are eigenstates of H j j i ¼ hω j j j i, 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 P j j j i j h j ¼ 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 . 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 Nth-order susceptibility far from any multiphoton resonances may also be calculated using the Nth-order perturbation theory 13 . For example, the "resonant" term in the third-order susceptibility, 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 non-linear processes, such as sum-frequency generation, can be calculated with simple modifications to W n at each step.

Multi-valley theory for donors in silicon and germanium
In this section, we develop the multi-valley 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 Kohn-Luttinger 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 plane-wave functions corresponding to the crystal momenta at the minima, k μ (and a lattice periodic function that is unimportant here). We write ψ j ðrÞ ¼ P μ e ik μ :r F j;μ ð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 X μ e ik μ : where H 0 includes the Coulomb potential V(r): E À1 using a valley-specific coordinate system (x, y, z where z is the valley axis, i.e., the valley-frame is rotated relative to the lab-frame 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 hydrogen-like, 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 Àik μ′ :r and averaging over a volume (2π/|k μ |) 3 around r, the Schrodinger eqn now reads is the Dirac delta function, and U μμ′ ¼ R dr e iðk μ′ Àk μ Þ:r UðrÞ. For an A 1 state, all the envelope functions have the same amplitude at r = 0, hence, P μ′ U μμ′ δðrÞF j;μ′ ðrÞ ¼ ÀU cc δðrÞF j;μ ðrÞ, where U cc ¼ À P μ′ U μμ′ . 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 ðmÞ 0 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 ψ n ðrÞ ¼ P m;μ e imϕ e ik μ :r f n;m;μ ðr; θÞ, substitute them into G À1 n ψ n ¼ ζψ nÀ1 , premultiply by e Àik μ′ :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 D ¼ U cc δðrÞδ m;0 =3 and δ m,0 is the Kronecker delta.
In the above equations we drop the valley-specific coordinates in f n,m,μ for notational simplicity, and the coordinates in H ðmÞ 0 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 non-zero for even n and m = 0. In the case of f n;m; 1 , there is only dipole coupling to the functions with the same m, while for f n;m;2 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 .

Giant third-order nonlinear susceptibility
Since silicon and germanium donors have an isotropic potential in an isotropic dielectric, the lowest-order 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 magneto-optical 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 sign-changes at which |χ (3) | goes to zero. In the range of frequency shown, we also observe a two-photon 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 low-frequency response of C (3) is illustrated at 100 GHz. Two higher-frequency 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 2p-average frequency, the 2ω resonance with the 1sE produces a coincidental zero-crossing 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 non-linearity. 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 ω ¼ ω 2p =3 ¼ 2π 3:2 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 phonon-limited 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 near-IR 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 single-layer graphene at room temperature reports a remarkably high value of 1.7 × 10 −9 m 2 /V 2 for the THz third-order nonlinear susceptibility 29  To match this value with Si:P at ω ¼ ω 2p =3 ¼ 2π 3:2 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 third-harmonic generation
The non-linear susceptibility is important for predicting the strength of frequency conversion processes such as third-harmonic 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 (p-Ge) laser 38,39 with a similar range and power. However, there is a large gap in the availability of solid-state 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 non-linearity is critical. At ω ¼ ω 2p =3 ¼ 2π 3:2 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 ω ¼ ω 2p =3, for example, a four-photon 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 non-linear 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 off-resonance 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 third-harmonic generation of THz light in doped silicon and germanium. Our multi-valley 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.