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.


Introduction
The lowest-order non-linearity in centrosymmetric materials is χ (3) , which describes that part of the response that is third order in the amplitude of the driving beams 1 . It is responsible for effects like degenerate four-wave mixing (DFWM), in which all four photons have the same energy and two are excited and two are emitted, Fig. 1. A substantial degenerate (or near-degenerate) FWM response is a prerequisite for applications of media in active optical systems ranging from modulators 2,3 to quantum repeaters [4][5][6] .
Although many non-linear effects have been demonstrated in the THz domain 7 , there are no quantitative measurements of susceptibilities for transparent bulk materials-indeed there are very few values of χ (3) available for any material in this part of the spectrum [8][9][10][11] . Doped silicon at low temperature has already been shown to produce giant values for the imaginary part of the non-linear refractive index (via multi-photon absorption) 12 , and there have been theoretical predictions that the real part of the non-linear refractive index is also very large 13 and of large experimental nonlinearities 14 but there have been no experimental reports of χ (3) till now, largely because of the challenge of quantitative non-linear THz metrology.

Experiment
We performed non-collinear DFWM experiments as illustrated in Fig. 1, using THz pulses from the free electron laser FELIX, both on and off resonance with the 1s → 2p transitions in Si:P and Si:Bi at 10 K. We chose this geometry because it enables the measurement of dynamical relaxation and dephasing times needed to make detailed theoretical comparisons, under identical experimental conditions. It is very difficult to obtain clean beam profiles with low diffraction in the THz regime, and great care was taken in avoiding apertures and optical imperfections in order to obtain them, as shown in Fig. 1. Care was also taken to accurately calibrate absolute pulse energies. It may be seen immediately from the relative strength of the output beam (k 3 ) in Fig. 1 that the DFWM process is very efficient.
In the plane-wave limit (i.e., for infinitely long pulses and infinitely broad beams), the complex polarisation amplitude of the generated beam (P 3 ) is related to the complex field amplitudes of the input beams (F 1;2 ) inside the material by i.e. the intensity of the output is determined by χ (3) . The definition of χ (3) in Eq. (1) suggests that, for a pulsed experiment, the internal pulse energies E i of the three beams k i (Fig. 1) are related by where E c is a constant that is inversely proportional to χ (3) L and L is the sample thickness. E c defined by Eq. (2) is a critical pulse energy at which the output would become equal to the inputs, and we generally stay well below this limit so as to avoid the need to consider higher-order non-linear effects.
We varied E 1 keeping the ratio E 2 /E 1 fixed, as shown in  Table 1. At high intensity, a saturation occurs for resonant cases, due to an intensity-dependent reduction in dephasing time 15 , which reduces χ (3) .
Conversion of E c to χ (3) Away from resonance and in the limit of long pulses, the relationship between E c (given on Fig. 2) and χ (3) has a straightforward dependence on the geometry and pulse duration. For short pulses, the dynamics are important, and on resonance there is loss that attenuates the pumps and the output, which must also be taken into account. We integrated the equations describing propagation of light through a lossy non-linear medium for the case of inhomogeneous broadening and finite pulse durations to find χ (3) from E c (see Supplementary Materials Section IV). In this case, the conversion from the experimental E c of Fig. 2 to the value of χ (3) is, for a beam with a Gaussian spatial profile, where n is the refractive index (which we took to be n = 3.4), λ 0 is the free-space wavelength, Z 0 is the  Fig. 1 The degenerate four-wave mixing geometry. a The energylevel scheme for DFWM with two excitation photons from beam 2, a stimulated emission from beam 1 and an output photon in beam 3. The left hand process involves two virtual excited states, and the right hand permutation is the strongest near to a resonance with the ground state g j i and excited state e j i. b A camera image of the beam at the sample position is shown superimposed on the sample. Farfield images were taken by scanning an iris after collimating, which requires careful conversion from space to angle. Each image has been normalised to the peak power density, and the scale factors for the far-field images are indicated relative to beam 1. Note that the far-field image of the beam 3 has only been scaled by a very small factor in this example (×3), i.e. the DFWM efficiency is very high. The phase matching condition is also shown , and the solid lines are fits to the low intensity portion. The fitted values of E c are also indicated. For the high density Si:P sample, only one intensity was measured at each laser frequency and a cubic dependence (dashed lines) is shown for comparison with the other measurements characteristic impedance of free space, and r 0 and t 0 are the root mean square (r.m.s.) beam radius and pulse duration, respectively (at which the intensity has fallen by 1= ffiffi e p ). The dimensionless factor f appearing in Eq. (3) depends on the loss and also the pulse shape and duration relative to the dynamical timescales of the system. The equation defines f in such a way that f = 1 when the loss is negligible (which is our case when far from resonance) and in the monochromatic limit of very long pulses with Gaussian temporal profile (t 0 ≫ T 1,2 , i.e. pulse duration much larger than the population decay, T 1 , and dephasing time, T 2 , of the system). For negligible loss but with pulses that are very short compared with the inverse line-width, then f becomes of order T 1 /t 0 , which can evidently be larger than unity (effectively replacing t 0 in Eq. (3) with T 1 because now the atomic polarisation P 3 lasts much longer than the drive pulses). For significant loss, f becomes very large and sample thickness dependent.
Using perturbation theory for temporally overlapping, weak beams within the two-level approximation 1 , and averaging over the distribution for a Gaussian (fully inhomogeneously broadened) line, we calculated values of f for our experimental circumstances. See Supplementary Materials for more details. The results are shown in Table 1. As expected, the off-resonant values of f in Table 1 are of order unity and are not significantly affected by the details of the model chosen. They are slightly greater than unity primarily because of the short pulses. The on-resonance values of f in Table 1 are large primarily because of the loss. The two-level model is expected to give a reasonably good estimate of f in resonant cases because there is one dominant transition: the one shown in Fig. 1a 1 .
The experimental values of E c from Fig. 2 along with the calculated f have been converted to values of χ ð3Þ expt in Table 1.

Theory
We now obtain theoretical estimates for χ (3) to compare with the experimental results. Silicon donors at low temperature are hydrogen-like, with a series of levels and orbitals closely resembling the Rydberg series 1s, 2p 0 , etc. 16 . The energies are scaled down and the orbital sizes scaled up, thanks to the small effective mass and large dielectric constant. The large orbitals give a commensurately large dipole moment, and this has a very large influence on non-linear optical coefficients.
Using the same two-level model mentioned above, the following limits may be found (see Supplementary Materials) for the contribution per bound electron in the vicinity of its resonance: where n D is the donor concentration, hΔ is the detuning from resonance in energy and μ eg ¼ ejhψ e jrjψ g i:ϵj is the component of the dipole moment transition matrix element between ground and excited states along the polarisation direction, ϵ. The total dephasing time T Ã 2 is defined by the total absorption line half-width in energy, h=T Ã 2 , which was obtained from the small-signal absorption spectrum. The population relaxation time, T 1 , was obtained by performing a pump-probe experiment 17 , and the homogeneous dephasing time, T 2 , was obtained using a photon echo experiment 1,15 . The results are shown in Table 1. These time-resolved experiments were performed with the same set-up that was used for the main theory =n D 0.0024 100 0.015 3100 23 18 hω is the photon energy, and labels R and T refer to resonant and transparent excitations. Values of μ eg are all taken from ref. 29 . All values for T 1,2 were found from photon echo and pump-probe performed under the DFWM conditions, except: a taken from ref. 21 . All values of the half-width, h=T Ã 2 , were found from the small-signal absorption spectrum, except: b assumed equal to the 2p 0 halfwidth. x is the ratio of the intensities of the pump pulses from Fig. 2. L is the sample thickness and r 0 is the spot size. The dimensionless factor f, which is unity for zero loss and infinitely long pulses, appearing in Eq.

Discussion
In the transparent regions, away from resonance (labelled T in Table 1), we obtain very good agreement between experiment and theory; there is also suitable but not perfect scaling with impurity density, n D . The approximate theory in Eq. (4) consistently underestimates the experiment by about an order of magnitude, implying that the neglected terms due to higher intermediate states and other permutations are additive. On resonance (R), the agreement is almost as good, with a similar magnitude of discrepancy but this time in either direction (particularly notable when we compare 2p ± transitions for P and Bi), presumably because of the strong sensitivity to the effect of the loss. It is obvious that resonance significantly reduces E c in Fig. 2 and enhanced χ (3) relative to the nonresonant cases. Figure 3 shows a survey of coherent χ (3) measurements in other materials, systems and frequency bands. The figure shows χ (3) L since this is the quantity that has actually been measured in each case, and it is the quantity that is relevant for frequency mixing applications. In Fig. 3, experiments in which the pump transition is virtual have been labelled as transparent, and those where there is a real absorption process at the pump frequency have been labelled as resonantly enhanced. For example, free carrier processes can produce not only a very large χ (3) L 8,9 but also very significant absorption; Dirac materials like graphene produce large χ (3) L 11,18 but have resonant interband or free-carrier processes depending on the chemical potential; and resonant enhancements by quantum well design 19 or Landau levels 10,20 also naturally induce absorption pathways. In such cases (where absorption loss is present), the volume susceptibility χ (3) is not an especially useful figure of merit for the material, because the output varies in a non-trivial way with sample thickness thanks to the loss. We note that very large apparent values of χ (3) have been reported in twodimensional and quantum well systems 9,11,19 . In all these cases, the measured output is normalised by the (very small) thickness, and the sheet susceptibility, χ (3) L, is very small by comparison with the values reported here and would remain so even for stacks of very many layers. It is interesting to note the general trend in Fig. 3 to increased susceptibility as the frequency is reduced. It can be seen from Eq. (4) that there is no intrinsic frequency dependence, so this increase is likely to be due to the difficulty of observing all but the strongest effects at THz frequency. It happens that the material used here has particularly large dipole moments 12 , which enter Eq. (4) with the fourth power, and an advantageous combination of long dephasing and decay times 15,21,17 . Our measured χ (3) L far from resonance is a record for any transparent material, and the only such measurement for THz pumping. The χ (3) values are all larger even than low temperature bulk InSb close to its band edge frequency 22 , meaning that here the contribution per electron (the hyperpolarisability) is far, far larger. This material can easily be produced in macroscopic thicknesses relevant for devices. An obvious immediate application is for metrology of a weak (k 1 ) THz beam with an arrival time clocked by a strong coupling pulse (k 2 ). Further prospects arise because compact and efficient semiconductor sources 23 now cover the entire electromagnetic spectrum from radio waves to the ultraviolet with just one gap between about 5-15 THz  (Table  4.6.1) and references therein (thanks to phonons in common polar solids): doped silicon could fill the gap by tripling the emission wavenumber for existing semiconductor THz lasers. For perspective, we point out that generation of THz light by downconversion from near-infrared (near-IR) [24][25][26] and mid-IR 27 sources is well established, but the efficiency is rather low, typically parts per thousand 24 . Our experiments were performed at cryogenic temperature, but many THz sources and detectors already require cryogenic environments 23 , and the operating temperature might be raised in future by exploiting deeper donors.

Samples
Samples used were single-crystal silicon doped with either bismuth or phosphorous and kept at a temperature of 5-10 K during the experiment.
The donor densities were determined by four-point resistivity measurements. In all cases, the surfaces were chemically and mechanically polished with a wedge of about 1°. The small-signal absorption was measured with Fourier transform spectroscopy with the samples mounted in liquid helium at 2.2 K (see Supplementary Materials for absorption spectra), and the half-width of each inhomogeneously broadened transition, h=T Ã 2 , was obtained from Gaussian fits. One transition, the 1s → 2p ± line in Pdoped samples, was saturated, and we took the line-width to be the same as for the 1s → 2p 0 line in this case. The concentrations, sample thicknesses and line-width values are given in Table 1.

Dynamical measurements, beam imaging and overlap
The optical set-up was a standard, time-resolved, forwards DFWM arrangement 10 , Fig. 1. All beams were focussed into a cryostat, recollimated and refocused using off-axis parabolic mirrors onto a high sensitivity liquid He-cooled Ge:Ga photo-conductive detector with a response time of about 100 ns.
A mechanical moving stage controlled the delay, and for the photon echo experiment (used to measure T 2 ), we measured the k 3 beam pulse energy as a function of the delay between k 1 and k 2 beams, while for the pump-probe experiment we simply moved the iris to detect the transmitted k 1 beam, which then functions as a weak probe. To ensure optimal overlap, we imaged the beams at the sample position with a pyroelectric camera with an effective pixel pitch of 80 μm (Spiricon Pyrocam IV). To obtain beam selection and optimal discrimination of the far-field beams after the sample, a motorised iris with a controllable aperture was mounted on a x-y scanning stage between the collimating mirror and the detector. The dependence of the output DFWM pulse energy shown on Fig. 2, E 3 , was measured with the iris open wide enough to capture the whole beam (while still excluding the pumps). The resulting T 1 and T 2 data are given in Table 1.

Pulse energy calibration
For metrology of the pulses in each beam for the data of Fig. 2, we calibrated the photon energy-dependent responsivity of the detector before each measurement. As a reference standard, we used a calibrated pyro sensor (SLT PEM 34 IR) with an accuracy of 2%.
For each set of measurements, we determined the ratio x = E 2 /E 1 by simultaneously recording both beams with the pyroelectric camera just before the sample, while scanning the undulator of the FEL.
The cryostat window transmission was calibrated by measuring the laser transmission through the empty cryostat (i.e. without the sample), referenced to the signal with the cryostat removed. We assumed both windows had the same transmission. The reflection loss at the sample surface was estimated using the Fresnel transmission coefficient T % 0:7 for the interface of the sample, which approximately agrees with the laser transmission signal when very far from resonance.
A polariser pair (Infraspecs P03) before the beam splitter was used to adjust the total laser pulse energy in fine steps for the intensity dependence of Fig. 2.
The values of r 0 used are given in Table 1 and came from the measurement with the Pyrocam mentioned above.

Pulse shape
In order to make the conversion from critical pulse energy to χ (3) (see below), the pulse duration is required. This autocorrelation signal appears Gaussian to a good approximation (red line) with r.m.s. width σ t = 7.5 ps, so the inferred r.m.s. duration of the pump beams is t 0 ¼ ffiffiffiffiffiffiffi 2=3 p σ t ¼ 6.1 ps The FELIX laser produces trains of intense, tuneable, short pulses. The train is emitted at 10 Hz in so called macropulse bunches, which contain approximately 200 micropulses each, at a repetition rate of 25 MHz. The pulse duration can be estimated from the spectrum since the pulses are approximately bandwidth limited. For a Gaussian pulse, its r.m.s. intensity duration t 0 = 1/4πσ ν where the r.m.s. intensity bandwidth in frequency ν averaged over the macropulse is typically about σ ν /ν ≈ 0.3% corresponding to a pulse duration of a few picoseconds, and there was little variation throughout the experiments. For this work, we made use of the fact that, when off-resonance, the measured E 3 as a function of the delay between the inputs k 1,2 gives a third-order autocorrelation (3AC), shown in (Fig. 4). For a Gaussian pulse, the r.m.s. duration is t 0 ¼ σ t ffiffiffiffiffiffiffi ffi 2=3 p ¼ 6.1 ps, where σ t is the r.m.s. duration of the 3AC given on the figure. For the data in Table 1, we took the value of t 0 from the 3AC measurement of Fig. 4, and we assumed that it was constant for all experiments.