Highly scalable multicycle THz production with a homemade periodically poled macrocrystal

The THz regime is widely appealing across many disciplines including solid-state physics, life sciences, and increasingly in particle acceleration. Multicycle THz pulses are typically formed via optical rectification in periodically poled crystals. However the manufacturing procedures of these crystals limit their apertures to below ~1 cm, which from damage limitations of the crystal, limits the total pump power which can be employed, and ultimately, the total THz power which can be produced. Here we report on the simple in-house fabrication of a periodically poled crystal using ~300 μm thick wafers. Each wafer is consecutively rotated by 180∘ to support quasi-phase matching. We validate the concept with a Joule-class laser system operating at 10 Hz and measure up to 1.3 mJ of energy at 160 GHz, corresponding to an average peak power of approximately 35 MW and a conversion efficiency of 0.14%. In addition, a redshifting of the pump spectrum of ~50 nm is measured. Our results indicate that high-power THz radiation can be produced with existing and future high-power lasers in a scalable way, setting a course toward multi-gigawatt multicycle THz pulses. Multicycle pulsed Terahertz radiation has appealing applications in medicine, applications in medicine, biology, material science, and could power emerging laser-based table-top accelerators and diagnostics. However a simple way to generate such pulses with high energy is not straightforward. The authors report on generation of high-peak-intensity THz pulses using a Lithium-Niobate wafer-chain and achieving a 35 megawatt in peak intensity.

H igh-power multicycle Terahertz (THz) radiation is highly sought after, with applications in medicine 1 , imaging 2 , spectroscopy 3 , characterization, and manipulation of condensed matter 4,5 , and could support the development of next-generation compact laser-based accelerators with applications in electron microscopy, ultrafast X-ray sources 6,7 , and subfemtosecond longitudinal diagnostics 8 .
Two prevailing paths toward powering THz-based accelerators have emerged. Cavity-based approaches seek to fill accelerating cavities with narrowband pulses [9][10][11] ; over the fill time, the accelerating field strength increases until a properly timed electron bunch is accelerated. This approach requires appropriate bandwidths (<1%) to match the quality factor of the cavity.
Alternatively, traveling wave approaches in, e.g., dielectric-lined waveguides, do not rely on filling times, rather they use the fields directly [6][7][8]12,13 . Here the duration of the THz pulse and group velocity in the waveguide limit the interaction length and energy gain; consequently, high THz peak-powers are advantageous to produce large accelerating gradients.
High-power single-cycle THz generation has flourished with the introduction of the tilted pulse front technique in lithium niobate (LN) [14][15][16] . Multicycle THz radiation has also been extensively studied, and can be produced with quasi-phase matching approaches in a periodically poled lithium niobate (PPLN) [17][18][19][20][21][22][23][24][25][26][27] . Alternative techniques have also been proposed and demonstrated with, e.g., semiconductor-based materials [28][29][30] . Multicycle pulses have also been generated in parametric oscillators 31 . Unfortunately, the manufacturing processes of PPLNs requires strong electric fields Oð10 kV m À1 Þ across the crystal width to locally reverse the polarization domains to overcome the coercive field strength of the material; this limits the crystal apertures to below~1 cm. The damage threshold of lithium niobate thereby limits the laser power which can be handled by the crystal, which inherently limits the production of high-power THz pulses. For high-power applications, the limited aperture sizes have motivated the development of wafer-bonded quasiphase matching crystals 32,33 . However, these fabrication techniques are challenging for large diameter wafers and require delicate in-vacuum procedures, limiting their widespread use.
Here we show that in the THz regime, a PPLN crystal can be effectively constructed in air by stacking lithium niobate wafers together; each wafer is rotated by 180 ∘ to the previous wafer to support quasi-phase matching. The relatively long (mm) wavelengths of the generated THz radiation compared to the small gaps (~10 μm) between wafers support a near-ideal THz transmission between wafers. To overcome the pump laser reflection at each boundary, each wafer surface is applied an anti-reflective coating. We demonstrate the concept using a Joule-class laser system with~50-mm diameter wafers and measure up to 1.3 mJ of energy at a central frequency of~160 GHz, corresponding to an average peak power of~35 MW, a 50 times increase in THz power compared to previous demonstrations 25 . Our results indicate that high-power THz radiation can be produced with existing and future high-power lasers in a scalable way, setting a course toward multi-gigawatt THz pulses. Moreover, the simplicity of the scheme provides a simple way to synthesize waveforms 21,34 for a variety of applications. We expect our results to have a broad range of potential use, including nonlinear optics, high-power pump-probe spectroscopy 35 , and to illuminate a path toward laser-based table-top high-brightness particle accelerators and diagnostics.

Results
Theory. The relatively long (mm) wavelengths of THz radiation supports a simple fabrication of a wafer stack in ordinary laboratory conditions. Just as a large wheel hardly feels a small bump on the road, here the generated THz wavelengths are much longer than the separation between wafers, supporting very small reflection coefficients at the wafer interfaces. For normal incidence, the THz transmission coefficient between wafers can be calculated 36 , where the THz wavelength, λ THz ¼ c f THz , f THz is the corresponding THz frequency and c is the speed of light in vacuum. The transmission as a function of the ratio between the wafer spacing and generated wavelength, d/λ, is illustrated in Fig. 1a. As an example, for d/λ = 0.001, T~0.999. Figure 1b depicts the assembled wafer stack; see "Methods" for a detailed mathematical description of the THz generation process. The domain period Λ (which is two times a wafer thickness), along with the refractive indices of lithium niobate for THz (n THz~5 .05) and IR (n IR2 .2) determine the central frequency of the THz radiation via, The wafer thickness and walk off (n THz − n IR ) limit the range of the accessible THz frequencies. However, from a practical point of view, the ratio between the wafer separation distance and wavelength (d/λ) strongly limits good transmission through the stack, see Eq. (1). For example, considering 300-μm-thick ZnTe wafers instead of LN, where n IR = 3.13 and n THz = 3.17 37 , phase matching to~16 THz would be possible. The corresponding wavelength of~19 μm would experience a small transmission coefficient, falling randomly in the oscillatory behavior of Eq. (1), which would generate a fast decaying signal. Alternatively, for example, LN, requiring a minimum transmission of T = 0.9 per wafer interface, for wafer separations of~10 μm would limit the generated frequencies to~1 THz. Pressing the wafers together could also help reduce the wafer separation to achieve better transmission, especially at higher frequencies.
Experiment. The experiment was carried out at the LASERIX facility of Paris-Saclay University at Orsay 38 , the laser parameters of the facility are presented in Table 2 in "Methods". The LASERIX laser driver is a chirped pulse amplification 39 based Ti: Saphire system which provides different femtosecond beamlines Fig. 1 Overview of the THz generation scheme. The THz transmission between wafers as a function of the ratio between the wafer separation and THz wavelength (d/λ) is illustrated in (a) corresponding to Eq. (1). A schematic of the homemade periodically poled lithium-niobate (PPLN) macrocrystal is shown in (b); the wafers are successively rotated by π for quasi-phase matching. Here the wafer separation, d, is approximately given by the total thickness variation (TTV).
from the mJ-level to~1.2 J after compression. All beams operate at 10 Hz and the main beam can be compressed down to 50 fs with the in-vacuum movable grating compressor. While conventional lasers generally produce transverse Gaussian mode profiles, high-energy multi-amplification systems tend to super-Gaussian (SG) profiles due to saturation in high-energy amplifiers. SG beams are more challenging to transport, their increased transverse uniformity permits a larger portion of the beam to operate near the damage threshold limitation of LN, thereby improving the conversion efficiency which scales linearly with the pump intensity. Moreover, the generated THz shares a similar mode profile to the pump laser and the propagation dynamics of the generated THz beam depend strongly on the beam size of the pump-laser beam. The half-angle of the THz divergence goes as θ ≈ λ THz /πn THz w 0 , where w 0 is the beam waist. In this respect, larger laser beam size in the crystal are advantageous for efficient waveform generation along the wafer stack. Figure 2 displays a diagram of the experimental setup. The main beam was rerouted into air by placing a converging lens with focal length f = 3 m in the vacuum chamber to reduce the beam size to pass through the 2 inch MgF 2 window. MgF 2 is chosen for its low dispersion and low nonlinear index n 2 which is suitable for femtosecond applications. A diverging lens was appropriately placed in air to provide a magnification of M = 0.5. A zero order half-wave plate was used to adjust the polarization of the IR beam with respect to the extraordinary axis of the wafers. The wafer stack was assembled by hand in air (see "Methods"). After propagation through the wafer stack, the IR beam was dumped onto a white diffuser supported on a rectangular Teflon (polytetrafluoroethylene (PTFE)) plate. The plate was tilted to spread the IR pump across a larger surface area. Teflon has excellent optical properties in the sub-THz range with a small absorption coefficient (see "Methods"). A fiber spectrometer was pointed at the diffuser and captured the depleted IR spectrum. The THz radiation propagated through the Teflon and went to our THz-diagnostics setup.
Two THz diagnostics were used (see Fig. 2). The first was dedicated to measuring the THz energy and consisted of two offaxis parabolic mirrors (OAP) placed to provide 4F imaging onto a THz pyroelectric detector. A black 2-mm-thick polyethylene (PE) filter was present on each THz detector to ensure that any remaining scattered infrared light did not contribute to an artificial THz signal; no THz signal was present without the wafers in place. A second diagnostic was constructed to measure an autocorrelation of the THz signal. We were unable to use the electro-optic sampling (EOS) technique to measure the THz waveforms directly in the allotted beamtime. Alternatively, an autocorrelation splits an original signal into two and recombines the two pulses onto a detector. By delaying one path with respect to the other, an interference signal between the two pulses is produced at the detector location, thereby providing some information on the temporal and spectral content of the coherent part of the signal. We note that the frequency content of the autocorrelation is the same as the original signal, i.e., they are Fourier pairs, see "Methods" for more details.
THz production with varying wafer counts. Initial wafer tests indicated differences between individual wafers. The wafers were accordingly sorted (see "Methods"). Subsequently, we investigated the output energy and evolution of the THz waveforms for stacks with different numbers of wafers with a pump-laser energy of E = 263 mJ, rms pulse duration of τ~700 fs and a FWHM beam diameter of 18 mm. The resulting energy measurements are illustrated in Fig. 3a with black markers and illustrate the quasiphase matching process. Here each additional wafer is rotated by 180 ∘ , and adds an approximately half-cycle to the waveform. Without absorption and other unwanted effects, each additional wafer would produce a linear increase in THz energy production. However, the generated THz is absorbed while propagating through the crystal. The absorption coefficient of the THz in the LN stack g can be calculated from the fit of the data (magenta trace), and Eq. (3) (see "Methods"), to g = 1.86 × 10 −3 THz −1 . A drop in THz energy production compared to the anticipated production from theory is noticeable beyond the 10th wafer. The reduced THz output likely arises from redshifting effects which leads to larger reflection coefficients at the wafer surfaces. In addition, dispersive effects can lead to lower intensities and phase mismatching. The fit was performed with the least-squares regression technique which minimizes the square of the sums diffuser Fig. 2 Diagram of the experimental setup. Here, HWP = half-wave plate, GM = gold mirror, OAP = off-axis parabolic mirror, PEF = polyethylene filter, and pyro = pyroelectric detector. The setup could be easily modified for energy measurements (using dashed-GM, no Si wafer) and for an autocorrelation measurement (no dashed-GM, with Si wafer beam splitter). Here, the red beam represents the Ti:Saphire pump laser and the blue represents the generated THz beam. between the data values and fitted values over the entire data set (see "Methods"). Generally, lower absorption coefficients are favorable as less radiation will be absorbed while propagating through a material; in this case, the relatively small absorption leads to the generation of semi-flattop temporal waveforms. Large absorption coefficients would result in the formation of asymmetric, rapidly decaying waveforms since the THz produced at the entrance face of the crystal would be absorbed while propagating through the crystal compared to THz generated close to the output face. The absorption factor is frequency and temperature dependent, generally lower frequencies have smaller absorption factors, and lower temperatures reduce the absorption which primarily arises from coupling to phonon modes in the LN crystal 19 . We have assumed g to be constant, but we note that the employed highpower laser could contribute to local temperature increases in the crystal via IR absorption, thereby increasing the THz absorption coefficient. Finally, autocorrelation signals for 1, 2, 6, and 12 wafers are shown in Fig. 3b.
THz production with 12 wafers. We further investigated THz production with 12 wafers. Figure 4a illustrates the comparison between measurement (black) and simulation (magenta) using the calculated absorption factor g and Eq. (3) (see "Methods"). The power spectrum is shown in Fig. 4b for both the measured and simulated signals from (a). The central frequency occurs at f THz~1 60 GHz, in good agreement with Eq. (2). Figure 4c shows the simulated and normalized electric field used in the autocorrelation comparison (b). Finally, Fig. 4d shows the THz energy scaling with pump intensity using the FWHM pulse duration as for comparison with Bach et al. 40 . Here two data sets are used, one with constant pulse duration and one with a constant pump energy which is shown in Fig. 5 (263 mJ). THz energy production decreases at large pump intensities and, as discussed above, is mainly attributed to redshifting effects which reduce pump beam transmission through the AR-coated surfaces and can also produce phase mismatching from dispersive effects. The drop off behavior for the constant pump energy data set could reasonably be attributed to other nonlinear processes like self-phase modulation in the stack. THz conversion efficiencies generally scale linearly with the pump intensity up to the saturation level. However, as is well discussed and derived in ref. 21 , the optimal FWHM pulse duration of the pump laser is given by the ratio between the domain length d to the walkoff length d w ¼ cτ FWHM n THz À n IR = 2.2. This criterion minimizes coupling to higher order modes of the PPLN and maximizes the coupling to the fundamental mode. The ideal pulse duration would be 1.3 ps for these conditions.
We explored the IR spectrum after 12 wafers for various pulse durations tuned by changing the gratings separation in the pulse compressor (which changes mainly the 2nd order spectral phasei.e., chirp-of the pulse). We did this at low and high pump energies (16 and 263 mJ) with a FWHM beam diameter of 18 mm (see Fig. 5). We were limited by SPM in the negative lens for pulse durations below τ ≈ 700 fs. For the case of 16 mJ, we were able to explore down to the best pulse compression of τ ≈ 50 fs. At large pump intensities, a significant redshift was observed along with some smaller amounts of blueshift for negative chirps. The simple intensity-based model from Eq. (3) (see "Methods") does not include cascading effects and is not suited to model this behavior. More elaborate theoretical models from Ravi et al. 26,27 or Vodopyanov 22 study these intricate behaviors in detail. We note that such behaviors can also come from other nonlinear effects in the wafers such as self-phase modulation (SPM); LN has an n 2~4 0 times larger than fused silica, which is favorable for SPM, see ref. 41 .
With no sign of damage from the LN wafers, we developed a final experimental configuration to explore the maximum producible THz signal. With the pump energy limited to E = 910 mJ, we removed the negative lens from our experimental setup to produce a slowly converging beam downstream of the 3 m positive lens located in vacuum. With a single wafer, we determined a safe operating location where the wafer took no damage; we then replaced the single wafer with the 12 wafer stack used in previous experimental setups. In addition, we also removed the 4F diagnostic setup and manually increased the collection efficiency by varying the position of the OAP and detector. We performed two measurements with different pulse durations, τ = 1.65 ps and τ = 1.2 ps; the corresponding measured THz energies were, respectively, E THz = 1.07 and 1.3 mJ, the latter corresponding to an average peak power of~35 MW using the calculated FWHM pulse duration of 37.5 ps, see Fig. 4c.
For the latter case of 1.3 mJ, the maximum field strength in Fig. 4c can be calculated for different spot sizes (see "Methods"). Focusing a Gaussian beam down to the diffraction limited beam waist of w 0 = 2λ/π~1.194 mm (with a numerical aperture of 0.5) with this energy would produce a maximum electric field strength of~1.6 MV cm −1 . However, this would require a focal length of 18 mm from our experimental conditions. By focusing the beam with a more realistic focal length of 100 mm, the beam waist would increase by a factor of 100/18~5.6, reducing the electric field tõ 285 kV cm −1 (because E ∝ I 1/2 ). Moreover, the anticipated quasiflattop or high-order super-Gaussian would induce an energy spread at focus into Airy rings which would further reduce the intensity, where 60% of the encircled energy could reasonably be considered in the central spot (assuming the same spot radius w 0 ); this would effectively reduce the maximum field to~170 kV cm −1 . This estimation is based on the originally employed beam size of 18 mm in the near-field. Such large field strengths at low frequencies are especially attractive considering the simplicity of the overall experimental setup. Moreover, by employing thinner wafers, and scaling to higher frequencies, significantly larger electric fields could be achieved.
The associated redshift of the pump beam associated with the production of 1.3 mJ, which was the maximum obtained at the end of the beamtime on LASERIX, is shown in Fig. 6. Here a center-of-mass shift of the pump spectrum by~50 nm is exhibited. From energy conservation, an upper limit on the IR-to-THz conversion efficiency can be calculated from the spectrum shift via η IR ¼ where ν i and ν o are the frequency barycenter of the IR spectrum before (λ i = 808 nm, ν i = 3.713 × 1015 Hz) and after (λ o = 858 nm, ν o = 3.496 × 10 15 Hz) the stack, respectively. For Fig. 6, the redshift would correspond to a maximum efficiency of η = 5.8%.

Discussion
There is a drastic difference between the measured conversion efficiency, η = 0.14% and the spectrum-calculated maximum efficiency η IR = 5.8%. Several factors can account for this discrepancy. The limited aperture of the THz detector and imaging optics of the THz limited the collection efficiency. Second, the large refractive index for the THz in LN (n THz = 5.05) generates a large reflection at the output surface into air, R = 0.45, see "Methods". THz absorption in the wafers (as visible in Fig. 3a) also reduces the measured energy by~31%. Finally, as discussed in ref. 26 , the generated THz is also inadvertently phase matched with the IR, leading to some THz energy being absorbed and blueshifting the IR spectrum.
Another experiment is being planned to incorporate significant improvements to the experimental setup, including a better THz transport, cryogenic cooling of the wafers to minimize absorption losses 19 , output coupling using for example, a fused silica wafers to decrease the reflection loss at the end of the wafer stack, and an investigation of longer wafer-stacks. The future experiment will also be used in preparation for an experiment aiming at injecting the THz produced with LASERIX into a circular dielectric-lined waveguide and postaccelerate an electron bunch produced at the laboratory, Irène Joliot-Curie (IJC-Lab) with the photoinjector, Photo-Injector at LAL (PHIL) 42 . Finally, we note that the super-Gaussian profiles are suspected to be beneficial to improve mode coupling into dielectric-lined waveguides, where the transverse field profiles are uniform for phase velocities v p = c (see refs. 6-8 ).  We have proposed and demonstrated a simple scheme to produce THz radiation with high-peak powers and energies in a highly scalable way. In this scheme, AR-coated lithium-niobate wafers were stacked in air to build a periodically poled macrocrystal. We investigated the scaling with the number of wafers in the stack through energy and spectroscopy measurements. Furthermore, we explored the redshift of the IR pump laser with varying pulse durations. We described a final measurement aiming to maximize the conversion efficiency and measured an energy of 1.3 mJ corresponding to an average peak power of~35 MW. To our knowledge this is the largest multicycle THz energy ever generated and the corresponding peak power is around 50 times larger than previous demonstrations 25 . The simplicity and scalability of the scheme could be attractive for THz waveform synthesis 21 and to power next-generation THz-based particle accelerators 6,7 and diagnostics 8 . The large peak powers may also be attractive to scientific endeavors for pump-probe spectroscopy and other applications for user platforms 35 , where field strengths could exceed 1 MV cm −1 with adequate magnification optics.

Methods
THz generation in PPLN. The generated field in a PPLN can be described well 18,19 via the partial differential equation, ∂ 2 ∂z 2 þ n 2 ðωÞ ω 2 where z is the longitudinal propagation coordinate, ω is the angular frequency, c the speed of light, μ 0 the magnetic vacuum permeability, E(z, ω) the electric field, and P(z, ω) is the nonlinear polarization density of the medium. The electric field can be calculated as þ g c ðL À z 0 Þ, τ is the IR rms pulse duration, L is the stack length, N is the number of wafers, A 0 is the field amplitude, and g is the decay time associated to losses.
Electric field strengths of the generated THz. The electric field strengths can be inferred from the measured THz pulse energies which can be useful to provide an estimate, for example, THz-driven condensed matter physics. The standard definition of intensity is given by the beam power over the area of the beam, where P is the peak power, A is the area of the beam, c is the speed of light in vacuum, ϵ 0 , is the vacuum permittivity in free space, and E is the instantaneous electric field. From the agreement between simulation and experiment, e.g., Fig. 4a, c, the instantaneous power in the pulse can be calculated from the total measured pulse energy, i.e., with the universal relation P ¼ E=τ. Subsequently by attributing a beam size, the electric field E can be calculated from Eq. (4) numerically.
Wafer specifications. The wafers were provided by Precision Micro Optics, see Table 1 of "Methods" for specifications. When stacking wafers together, good surface contacts were evident from adhesion forces between the wafers, i.e., the wafers could not be separated by pulling, but rather only by sliding the wafers along the surface faces. The wafers had alignment markers to distinguish the polarization axes and were accordingly stacked and placed into an optical mount. The total thickness variation (TTV) < 10 μm of the wafers indicated an expected transmission coefficient from wafer to wafer of~99.4% (Eq. (1)). To overcome the reflection of the pump laser operating at significantly smaller wavelengths (~800 nm), each wafer face had an antireflective coating which provided a transmission above 99.75% per face, see Fig. 1b for a schematic of the wafer stack Table 2.
Wafer sorting. In early tests, we observed different responses from individual wafers and sought to identify the reason for this behavior. We sorted the wafers by measuring the polarization rotation of the pump laser by using an energy detector placed after a polarizing cube. X-cut wafers should not have any birefringence for polarizations aligned in the z-plane, according to our purchase order; however, this was not the case. We distinguished approximately three distinct families of wafers from our order; there were a small amount of outliers which we discarded. Each wafer was tagged with a non-permanent marker. The remainder of all measurements were conducted with the 'purple' family, which exhibited the expected response. We also sorted the wafers in terms of thickness to determine if there was any contribution to our observations using a micrometer-caliper. We confirmed that the thickness variation was not responsible for the observed differences by mixing two families with the same thickness. Moreover, the thickness variation was within specifications. We suspect the wafers we received from the manufacturer were not properly cut.
Wafer stacking. To minimize THz reflections at wafer-wafer boundaries, the wafer separation distance should also be minimized according to Eq. (1). Each wafer was thoroughly cleaned with acetone and optical-wipes. We explored various ways to minimize the contact between wafers. Our most successful approach was to slide the wafers onto each other in parallel. This avoided the formation of airpockets by contacting two wafers directly. Good wafer contacts were evident from the adhesion forces between the wafers. We tried using various optical mounts to support the wafers. We used a typical rotation mount (FMP2) from Thorlabs, but the rotation clamp led to the rotation of some wafers. This was solved by applying some tape to the sides of the wafer stack. We also used a simple c-clamp to hold the wafers together and did not notice any change in the signal strength.
Pump and THz separation. The high-power laser pump and the THz need to be separated for THz diagnostics and future applications. We tested various approaches but found that a white diffusing thin sheet, usually used for laser dumping mounted onto a 6-mm-thick block of Teflon (PTFE) block worked well. The diffuser scattered the infrared pump beam with a good transmission of the THz; this was tested by measuring the THz signal with and without the diffuser present at low infrared pump intensities. Without the diffuser, the Teflon quickly damaged at high pump energy. The Teflon does, however, also reduce the measured THz energy, as it absorbs some of the generated THz and produces reflections at both interfaces. At f ≈ 160 GHz, the THz absorption in Teflon is very low, α ≈ 0.01 cm −1 43 , according to the Beer-Lambert law where the initial intensity I 0 decays with distance through the material as IðzÞ ¼ I Àαz 0 . The reflection coefficient of a wavelength at an interface between two boundaries (n 1 , n 2 ) with different refractive indices can be calculated from the Fresnel equations, The refractive index of Teflon at our frequency is n THz = 1.43 43 which corresponds to a reflection coefficient at each surface of~3%. Altogether, the total associated losses due to the Teflon is~6%.
Autocorrelation. The autocorrelator was assembled with two gold off-axis parabolic mirrors (OAPs), two planar gold mirrors, a silicon wafer and a THz detector, see Fig. 2. The input THz beam is directed onto a silicon wafer, which serves as a beam splitter. The reflected portion is subsequently reflected by a gold mirror mounted on a linear translation stage, the transmitted portion of the THz beam is also reflected on another gold mirror. Both the transmitted and reflected parts are recombined on the beam splitter and transported to the pyroelectric detector. The interference between both paths leads to the autocorrelation signal via a scan of the delay stage. The autocorrelation signal strength is limited by the transmission and reflection coefficients on the silicon wafer. By exchanging the silicon beam splitter with a gold mirror, we were able to measure the THz energy directly, see Fig. 2.  Mathematically, the autocorrelation is given by the convolution of the electric field with its complex conjugate, From the Wiener-Khinchin theorem, the Fourier transform of A(t) is equivalent to the Fourier transform of E(t), i.e., they are Fourier pairs. See below for a brief discussion the Fourier transformation.
THz detectors. All three THz detectors used were from GentecEO. A QE9SP-B-MT-L-BNC-D0 was used for low-energy measurements and is capable of measuring as low as~100 nJ of THz radiation. A second detector, QE8SP-B-MT-D0 is much less sensitive but capable of measuring up to the mJ range. This was used for the majority of the measurements. Finally, a third calibrated THz detector was used, THZ5I-MT-BNC to cross-calibrate the other two detectors for our energy measurements. A black 2-mm-thick polyethylene (PE) filter was present on each THz detector to ensure that the remaining diffused infrared radiation after the diffuser and Teflon (PTFE) did not contribute to an artificial THz signal.
IR pulse measurements and calibration. Laser pulse durations τ < 00 < τ < 1000 fs, a single-shot autocorrelator (TiPA, Light Conversion) was used. Pulse durations longer than 1 ps were estimated using the forward projection of the measured τ as a function of the compressor grating separation. Changing the distance D between the two gratings of the laser compressor introduces a spectral dephasing which increases the pulse duration from the optimal Fourier limited duration. The main contribution is called chirp, group velocity dispersion (ψ″) which enlarges the duration from T 0 (Fourier limitation) to T with the following approximation (considering a spectro-temporal Gaussian profile): T % 4ln ð2Þψ 00 =T 0 , and ψ″ varies linearly with D.
Data acquisition and error analysis. The IR laser pulses were measured with an Ophir power meter with 2% rms uncertainty. The employed THz detectors had a measurement uncertainty of 3% rms. The shot-to-shot energy fluctuations of the IR laser was 0.5% rms corresponding to 1% rms shot-to-shot energy fluctuations of the measured THz energies. The measurements were averaged over 100 shots with an oscilloscope. The error bars on our figures reflect the square root of the sum of the squares for these uncorrelated errors.
Processing techniques for acquired data. The Fourier transform is essential to understand the bandwidth or frequency content of a signal. In essence, a signal in time, s(t) can be transformed into the frequency domain (F ðωÞ) via, note here, the angular frequency used here is defined as ω = 2πf. Generally, however, detectors cannot measure signals continuously, i.e., they have a limited sampling rate, this detail leads to the requirement of a discrete Fourier transform (DFT) or the fast Fourier transform (FFT) to analyze the spectral content of the signal. Figure 4b shows the resulting power spectrum of the FFT of Fig. 4a, and is equivalent to the FFT of Fig. 4c. For details behind the FFT algorithm, see ref. 44 . The THz absorption factor g is determined by fitting the measured pulse energy vs. the number of wafers against the theoretical model as described above. The leastsquares method is employed for this purpose, where the sum of squared residuals with r i = y sim,i − y mes,i is minimized. Here N refers to the number of measured data points and sim and mes correspond to the simulated and measured values, respectively.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.