Highly efficient THz four-wave mixing in doped silicon

Third-order non-linearities are important because they allow control over light pulses in ubiquitous high-quality centro-symmetric materials like silicon and silica. Degenerate four-wave mixing provides a direct measure of the third-order non-linear sheet susceptibility χ(3)L (where L represents the material thickness) as well as technological possibilities such as optically gated detection and emission of photons. Using picosecond pulses from a free electron laser, we show that silicon doped with P or Bi has a value of χ(3)L in the THz domain that is higher than that reported for any other material in any wavelength band. The immediate implication of our results is the efficient generation of intense coherent THz light via upconversion (also a χ(3) process), and they open the door to exploitation of non-degenerate mixing and optical nonlinearities beyond the perturbative regime.

In these Supplementary Materials we provide some additional experimental details about the spectral dependence of the signal obtained in the experiment in section I. We provide details of the calculation of χ (3) for a two level model in section II. In section III this two-level result is extended to include Gaussian inhomogeneous broadening and produce the approximate results for the non-linear susceptibility in Eq.  Table I  The FTIR absorption spectra are shown with black lines and the grey shaded curves indicate the laser spectrum used. Absorption is defined as 1 − Iout/I bckgnd where Iout is the transmitted intensity normalised by the system response, and I bckgnd is the value of Iout between resonances. In (a) there are four different laser frequencies indicated ( ω)and four corresponding beam profiles. In (b) a single laser frequency is used an the beam profiles are shown for four different laser intensities indicated by the resulting output pulse energies (E3). Noise caused by scatter from the pumps is evident when the efficiency is low, at low intensity and also at high intensity when the DFWM saturates (while scatter does not). The color scale for each beam image has been normalized. The FTIR absorption spectrum of the nD = 10 16 cm −3 Si:P sample is shown in red (right axis). The thick black curve is the DFWM output for this sample (left axis), which follows notionally the product of the FTIR and air transmission from (b). For clarity two portions of the DFWM signal have been scaled by factors of 0.5 and 10 as indicated, and in one section of the scan the laser input beams were attenuated by a factor of 0.1 as indicated.

I. SPECTRA OF THE SMALL SIGNAL ABSORPTION, THE DFWM SIGNAL AND THE LASER
We monitored the DFWM output beam throughout the experiment in order to ensure that it was clearly identifiable and separated from scatter from the pump beams. Example beam images are shown in Fig. S1, which also shows the laser spectrum and sample absorption spectrum for comparison.
An example calibration of the ratio between E 2 and E 1 and the absolute pulse energy E 1 as a function of laser frequency is given in Fig. S2. The figure also shows the transmission of air, and the resulting laser frequency dependence of the output pulse energy, which clearly scales with the energy of the pump pulse and the absorption cross-section, as expected.

II. TWO-LEVEL MODEL FOR THIRD ORDER SUSCEPTIBILITY
The Hamiltonian for the interaction of light with the donor is where H 0 is the Hamiltonian of the donor's electron with eigenstates H 0 |j = ω j |j and V (t) = −µF (t) is the dipole potential with electric field F (t). The density matrix of the donor's electron evolves according to the quantum Liouville equation where W is responsible for relaxation and dephasing processes. In the basis of unperturbed states |j we take where Γ jk is the damping rate of ρ jk ≡ j| ρ |k back to its value at thermal equilibrium in the dark, ρ jk . The rate Γ jj is the population relaxation rate of the jth level and the Γ j,k , with j = k, are the dephasing rates of the off-diagonal elements. At low temperature we can assume that ρ (eq) gg = 1 while the equilibrium values of the other elements are zero.
From Eq. (2) we have where the last step follows from inserting |l l| = 1 between the operators µ and ρ, we have One can express the solution in terms of the Green function which is given by where δ(t) is the Dirac delta function, which yields where and Θ(t) is the step function. The function G jk (t) can be thought of as the envelope of the Green function. One can verify that G jk (t) is the solution of Eq. (6) by using dΘ(t)/dt = δ(t). For jk = gg, the equilibrium value is zero at low temperature, hence where we have assumed the initial condition ρ jk (−∞) = 0 for jk = gg. Although this integral equation does not provide a direct evaluation of ρ jk (t) since ρ lk (t ) on the RHS is unknown, it is useful for obtaining the perturbative solution for weak field. We expand ρ jk (t) = n ρ (n) jk (t), then Eq. (9) gives ρ (n) which can be applied successively to find ρ jk up to the nth order [1]. When the carrier frequency of the pulse is tuned close to the transition frequency between the ground state |g and excited state |e , one can use the two level approximation. Now there are only two damping rates: The population relaxation rate, Γ ee = 1/T 1 and the dephasing rate, Γ eg = Γ ge = 1/T 2 . The time-dependent polarization of a sample with donor density n D is so all we need to find is ρ eg (t). The zero order elements of ρ(t) are ρ jk =gg = 0. Substituting this into Eq. (10) and using µ gg = µ ee = 0, one obtain the first order off-diagonal elements and then the second order population and finally the third order off diagonal element where the last step follows from ρ (2) gg (t) + ρ (2) ee (t) = 0 as the total population is conserved, hence (15) The polarization up to the third order is then P (t) = P (1) Note that the second order term P (2) (t) = 0 since ρ (2) eg (t) = 0. We write the electric fields of the pump and rephasing pulses as where ω 0 is the carrier frequency, F 1,2 (z, t) the pulse envelopes, and the z axis is parallel to the light beams and hence perpendicular to the sample surface. The total electric field is To find the 3rd-order polarization we substitute the electric field of Eq. (18) into Eq. (15). Each electric field is a sum of 4 terms, and an expansion yields 4 3 = 64 terms in ρ (3) eg (t). However, we need to keep only the terms that propagate in the 2k 2 − k 1 direction of the ouput field, and we can also neglect the counter rotating terms. Upon inspection there are only two terms left, and the third order polarization is where the envelope is where ∆ = ω 0 − ω eg is the detuning andF 1,2 (z, t) = F 1,2 (z, t)e −i∆t .
In the monochromatic limit of infinitely long pulses, and when there is no loss due to absorption, F 1,2 (z, t) and hence P (3) (z, t) are constants. It's now straightforward to carry out the integrals in Eq. (20), which yields where the third order susceptibility is On resonance ∆ 1/T 2 and while sufficiently far from resonance ∆ 1/T 2 and in this case

III. EFFECT OF INHOMOGENEITY IN THE TWO LEVEL MODEL FOR THIRD ORDER SUSCEPTIBILITY
For an inhomogeneously broadened sample we model the distribution of the transition frequency, ω eg , by a Gaussian where 2 ln(2)Γ = 1/T inh is the half width at half maximum of the distribution in angular frequency, which defines T inh , the inhomogeneous contribution to the dephasing time, and ω eg the peak frequency of the broadened transition.
For such a sample the third order susceptibility given in Eq. (22) has to be averaged over the distribution of ω eg : where and δ = ω eg − ω eg . ∆ = ω 0 − ω eg as before, which now means the detuning from the centre of mass. If ΓT 2 1 then we have a homogeneous line, and g(δ) is only large near δ = 0, and since g has unit area η ≈ 1/(1 + ∆ 2 T 2 2 )(i + ∆T 2 ) and we recover Eq. (22).
Far from resonance, ∆ 1/T 2 and g(δ) is only large when δ is small compared with ∆ so and we recover Eq (24), which evidently holds for both homogeneously and inhomogeneously broadened transitions. When the laser is on resonance ∆ Γ, where Erfc is the complementary error function and on the first line we used the fact that the real part of the integrand is odd. For a homogeneously broadened line T inh T 2 so ΓT 2 1 and (as it must be comparing Eqs (23) and (26)). For an inhomogenously broadened line T inh T 2 , so ΓT 2 1 and The linear absorption line shape in the presence of both homogeneous and inhomogeneous broadening is given by a Voigt profile (Main Text Ref [17]) which is the convolution of the Lorentzian of width 1/T 2 in angular frequency for the homogeneous contribution, and the Gaussian g(δ) inhomogeneous contribution: where x = T 2 ∆. The half width of the Voigt lineshape in angular frequency, 1/T * 2 , is given by T * 2 must be found numerically in general, though the asymptotic limits are T * 2 ≈ T inh for T inh T 2 and T * 2 ≈ T 2 for T inh T 2 , and a common analytical approximation based on these limits is (Main Text Ref Similarly, a useful analytical approximation for η based on its asymptotic behaviour Eq.s (30) & (31) is Using 2 √ π ln 2 ∼ 1 in Eq. (35) and substituting (34), or by inspection of Fig S3, a further simplifying approximation for η is These approximations are shown on Fig S3, and it may be seen that they are satisfactory over a wide range of T 2 /T inh . The maximum discrepancy in Eq.
This approximation now covers the resonant situation for both homogeneous and inhomogenous broadening, as did Eq (24) for the off-resonant situation, given in the Main Text as Eq. (4). For self-consistency within Main Text we use the approximate values resulting from Eq. (37) in Table I.

A. Wave propagation
To get the output field from the nonlinear polarisation we need to solve the wave propagation equation. The electric field satisfies the Maxwell equation [2] The total polarization in a doped semiconductor is the sum of the host's polarization (silicon in our case), P host = 0 ( r − 1)E ≈ 0 (n 2 − 1)F where n is the refractive index of the host, and the polarization component P of the donors that propagates in the k = 2k 2 − k 1 direction. Substituting P tot = 0 (n 2 − 1)F + P yields For general non-monochromatic pulses with a finite frequency broadening this equation is best solved with Fourier transform. Starting with where the envelope functions have the Fourier transforms where ω is the deviation from ω 0 . Substituting Eqs. (41) and (40) into Eq. (39) we obtain the propagation equation for each Fourier component When the first spatial derivative of the envelope functions changes very little over a wavelength, which is true for our experiment because the wavelength is of the order of 0.01 mm while the length of the Gaussian pulse is of the order of 2 mm, we have Neglecting the 2nd order spatial derivative in Eq. (42) and using k = ωn/c and µ 0 = 1/( 0 c 2 ), we obtain where Z = Z 0 /n and Z 0 = 1/c 0 is the impedance of the vacuum. In our experiment the frequency bandwidth is much smaller than the carrier frequency, ω ω 0 , hence ω + ω 0 ≈ ω 0 and In the monochromatic limit of infinitely long pulses there is only a single mode, ω = 0, in the spectrum. Assuming that there is no loss, substituting the third order polarisation from Eq. (21) and integrating, we obtain the four wave mixing output field where L is the thickness of the sample.

B. Loss of the input due to absorption
In our experiment the pump and the rephasing pulses are approximately Gaussian before entering the sample where the time envelope at the beginning of the sample, z = 0, is where τ is the r.m.s. duration of the electric field, and is related to the r.m.s duration of the intensity in the Main Text, t 0 , by τ = √ 2t 0 . As the input pulses travel into the sample they are attenuated by loss due to absorption. This is predominantly a linear optical process due to the 1st order polarisation, which is where k in is the wave vector of the input field, either the pump or rephasing. Substituting the electric field from Eq. (17) into Eq. (12) and then the resulting into Eq. (16), we obtain The input field inside the sample is given by the wave propagation equation, Eq. (45), which can be solved exactly by first taking the Fourier transform then Substituting Eq.(53) into Eq. (45) we obtain which is the usual equation for absorption in the frequency domain where the frequency dependent absorption coefficient is The Fourier component of the input pulses inside the sample is thus given by a simple exponential Therefore, the time envelope of the field inside the sample has the form where the Fourier component isŜ

C. Output field
Substituting Eq. (47) into Eq. (20) the third order polarisation is now where whereS(∆, z, t) = S(z, t)e −i∆t . We now use the Fourier transform of S(z, t) in the integral of , Similarly, .
The Fourier transform of Q 1 (∆, z, t) iŝ Using the Dirac delta function's identity we obtainQ Similarly,Q For the output field the wave propagation equation Eq. (45) becomes where P (1) is the linear polarisation propagating in the same direction as the output field, 2k 2 − k 1 , and is responsible for output loss due to absorption. Therefore, whereQ = (1/2)(Q 1 +Q 2 ) and The solution for the output field at the end of the sample, z = L, iŝ Usinĝ and hence we haveF whereĴ = (1/2)(Ĵ 1 +Ĵ 2 ) and

D. Effect of inhomogeneous broadening
To take into account the effect of inhomogeneous broadening the output field in Eq. (75) and thusĴ of Eqs. (76) and (77) has to be averaged over the Gaussian distribution of Eq. (25) and the output field of Eq. (75) now becomeŝ Before calculating the input and output energy, we first note that the pump and the rephasing beam in our experiment has a Gaussian radial profile, i.e., the field F 1,2 decreases with the distance from the center of the beam, ρ, through the factor where ρ 0 is the r.m.s. radius of the electric field profile of the beam, and is related to the r.m.s radius of the intensity profile in the Main Text, r 0 , by ρ 0 = √ 2r 0 . Since the four wave mixing output field is proportional to F * 1 (0, 0)F 2 2 (0, 0), its radial profile is W 3 (ρ).
To calculate the energy one has to integrate the intensity over time and then calculate the spatial integration over the radial profile. The input energy at z = 0 is and the output energy at z = L is Using Eq. (79) we obtain The critical energy discussed in the Main Text is From Eq. (26) we have Substituting this into the formula for E c , we obtain where the critical intensity is defined as with λ 0 = 2πc/ω 0 , the carrier wavelength of the laser, and E. Analytical approximations for short pulses

Far off resonance
Far away from resonance the loss is negligible because α(∆, ω)L 1. In addition, the detuning is much larger than the width of the inhomogeneous broadening, ∆ Γ, as well as the spectral width of the laser, ∆ 1/τ , so the effect of the inhomogeneous broadening can also be ignored. Moreover, the detuning of the central laser frequency is much larger than the width of the line, ∆ 1/T 2 , 1/T 1 , and for short pulses the duration satisfies τ T 1 , T 2 . This leads to simplification in the expression of the output field.
We use the following approximation: ∞ 0 e −i∆t−t/T2 S(t) To prove this consider the integration from t = 0 to t = 2π/∆ which is the very short period of the oscillating factor.
has a simple analytical result. Using ∆T 2 1 and the fact that the duration of S is much smaller than T 2 one can show that this analytical result can be approximated by and Eq. (98) follows from applying this repeatedly until infinity. Now we evaluate Q 1,2 (∆, z, t) in Eq. (61) (we drop the z variable because when the loss is negligible S and hence Q 1,2 are independent of z) Note that S(t 3 − t 2 ) is non-negligible only when t 2 is within 3τ around t 3 and since τ T 1 the exponential factor e −t 2 /T1 changes very little over this interval, hence And finally,