Terahertz strong-field physics in light-emitting diodes for terahertz detection and imaging

Intense terahertz (THz) electromagnetic fields have been utilized to reveal a variety of extremely nonlinear optical effects in many materials through nonperturbative driving of elementary and collective excitations. However, such nonlinear photoresponses have not yet been obeserved in light-emitting diodes (LEDs), let alone employing them as fast, cost-effective, compact, and room-temperature-operating THz detectors and cameras. Here, we report ubiquitously available LEDs exhibiting photovoltaic signals of ~0.8 V and ~2 ns response time with signal-to-noise ratios of ~1300 when being illuminated by THz field strengths ~240 kV/cm. We also demonstrated THz-LED detectors and camera prototypes. These unorthodox THz detectors exhibited high responsivities (>1 kV/W) with response time four orders of magnitude shorter than those of pyroelectric detectors. The mechanism was attributed to THz-field-induced impact ionization and Schottky contact. These findings not only help deepen our understanding of strong THz field-matter interactions but also contribute to the applications of strong-field THz diagnosis. Interest in the exploration of non-perturbative nonlinear optical phenomena driven by intense terahertz fields has seen a leap forwards with the recent development in femtosecond laser-based table-top sources for strong THz radiation. The authors present intense THz-field-induced effects in ubiquitously available LEDs illuminated by strong THz pulses, paving the way to their use in detecting and imaging intense THz radiation.

C urrently, there is much interest in exploring nonperturbative nonlinear optical phenomena in various types of matter driven by intense terahertz (THz) electric and magnetic fields [1][2][3][4][5][6][7][8][9][10] . This is fueled by the impressive advancements made during the last two decades in the development of femtosecond laser-based table-top sources for strong THz radiation [11][12][13] . These intense THz electromagnetic fields have been utilized [14][15][16][17][18] to efficiently accelerate and manipulate electrons to extreme degrees, advancing the frontiers of both the technology and science of nonlinear THz-matter interactions. For example, acceleration of relativistic electron bunches in free space is promising for the realization of next-generation table-top particle accelerators and compact X-ray sources. In solids, electrons have been excited and accelerated by strong THz fields, leading to an ultrafast phase transition 19 , interband tunneling and impact ionization 20 , exciton generation, and fluorescence 21 .
These groundbreaking studies have stimulated much excitement in diverse fields of photonics, materials science, and condensed matter physics, requiring further sophistication and advancement of intense-field THz technologies. However, what is strikingly absent in this context is a convenient detector for intense THz radiation. Currently available detectors are bulky and/or expensive, often requiring liquid nitrogen or liquid helium, and are usually very slow in response. Fast, sensitive, small, inexpensive, and room-temperature-operating detectors and cameras for intense THz radiation are being sought.
Here, we report that one of the most commonly available semiconductor devices around us, light-emitting diodes (LEDs), serve this purpose. That is, these inexpensive, ordinary devices work very well for detecting and imaging THz radiation. We observed a gigantic photovoltaic signal when such LED devices were illuminated by strong-field THz pulses. This observation presents us a rare opportunity to explore the fundamental science of the ultrafast response of electron-hole pairs in a semiconductor device induced by strong-field THz pulses. More importantly, this work demonstrates that these roomtemperature operating, inexpensive LEDs are promising for developing next-generation THz detectors, cameras, and functional devices.

Results and discussion
THz-field-induced photovoltaic signal in LEDs. Intense THz pulses were generated from lithium niobate (LN) crystals via optical rectification 13,22 , and a schematic diagram of the experimental setup is shown in Fig. 1a. The generated THz radiation was vertically polarized, and its strength was adjusted by a pair of THz polarizers. We focused the THz beam to achieve a fluence of, at best,~4.1 μJ/mm 2 ; the LED detection size was only~200 × 200 μm 2 , and the focused THz beam diameter was~3.3 mm in full width half maximum (FWHM). See more details in Supplementary Fig. 1f. Under these excitation conditions, we were readily able to obtain the polarization-dependence responses via rotating the azimuthal angle of the printed circuit board (PCB) (see Fig. 1b) and monitor photosignals in an oscilloscope (Fig. 1c) in all LEDs at the same time. Typical THz temporal waveforms (see, e.g., Fig. 1d) extracted by a single-shot spectrum-encoding method 23 indicated that the generated THz pulses were near single-cycle with a frequency bandwidth of~0.8 THz (see Supplementary Fig. 1c). The measured LEDs in our experiment with different colors (red, yellow, green, blue, and white in Fig. 1e) which correspond to the bandgaps (central emission wavelength) of 1.9, 2.1, 2.3, 2.8 eV (except for the white LEDs which includes three types of LEDs of red, yellow and blue). More detailed information of the LEDs can be found in "Methods".
THz field strength and polarization dependence. To understand the origin of the photovoltaic effect in the LEDs, we examined the THz pump fluence and polarization dependence of the response of the blue LED. In the pump-fluence-dependent experiments, we used six different peak fields: 81. 4,119,155,190,223, and 241 kV/ cm (see Fig. 2a, the mentioned THz field strength value was the maximum peak electric field over the whole temporal duration. The details of the field strength calculation can be found in Supplementary Note 4). Observed time-dependent photoresponse curves are summarized in Fig. 2b. From this figure, we can see that all these signals are negative, i.e., the THz-pulse-induced current flew in the reverse direction, from the n-side to the p-side of the device. When the THz electric field was 81.4 kV/cm, the peak-to-peak photoresponse signal was~41 mV, while it increased to~600 mV when the THz field strength was increased to 241 kV/cm. The minimum detected signal was under the incident field strength of~50 kV/cm. Figure 2c plots the photovoltaic signals of the LED when it was illuminated by THz radiation with three different pulse energies: 80.0, 52.8, and 24.7 μJ. The response curves exhibit a two-fold symmetric period of 180°, which is similar to the previous work on the nonlinear photoresponse of type-II Weyl semimetals 24 , manifesting an anisotropic photoresponse related to crystallographic symmetry. We interpret this polarization-dependent photoresponse as a result of different threshold energy due to the orientation-dependent electron effective mass of GaN (see "Methods"). Figure 2d shows two photoresponse curves as a function of THz pump fluence for specific THz polarizations (marked (1) and (2) in Fig. 2c). For pump fluences lower than 1.2 μJ/mm 2 , the response tendency exhibited a quasi-quadratic behavior, while it increased linearly when the fluence was larger than 1.2 μJ/mm 2 . As an example, a response signal of a blue LED is illustrated in Fig. 2e. The response time of the detected photovoltaic signal was four orders of magnitude shorter than that of a commercial pyroelectric detector (SDX-1152, Gentec, sensor diameter = 8.0 mm, see Fig. 2e). Furthermore, the responsivity of these LEDs was as high as~10 times that of the pyroelectric detector without an amplifier, reaching 750 V/W.
THz field induced impact ionization. When a strong THz field interacts with a semiconductor, various nonperturbative nonlinear optical phenomena can occur, including high-harmonic and sideband generation 25 , the dynamic Franz-Keldysh effect 26,27 , Zener tunneling 28 , metallization 29 , and impact ionization 6,30 . Carrier generation can result in unconventional ways through some of these processes, even though the THz photon energies used are much smaller than the bandgap. For example, Zener tunneling-induced photocarrier generation has been demonstrated, although the density achieved remained relatively low 28 ; metallization can also instantaneously produce carriers but requires much higher (>100 MV/cm) field strengths. Hence, in order to explain our observed huge photovoltaic signals in the LEDs, we propose impact ionization to be the dominant mechanism. A simplified Monte Carlo simulation based on this mechanism was performed, which reproduced all salient features of our experimental results, as detailed below.
Impact ionization depicts a three-particle down-conversion process in which a high-energy (hot) conduction band electron (or valance band hole) interacts with a valence band electron, and excites the electron across the bandgap, leaving behind a hole 31 . Over the entire interaction process, the hot carrier is enforced by the THz electric field according to hdkðtÞ=dt ¼ qEðt; NÞ þ γ tot (q = e for holes, and q = −e for electrons). Here k(t) denotes the wavevector of carrier, h is the reduced Planck constant, e is elementary charge, and E(t, N) is the screened THz electric field which concerns the carriers density N. γ tot presents the total scattering rate which involves impact ionization (see "Methods"). Getting use of the impact ionization transition rate of GaN from the first principle calculation 32 , we found that the hole-initiated impact ionization transition is much more efficient than that of electrons (see "Methods"), when the incident THz pulse approached the p-doped region prior to n-doped region in the blue LED (see Fig. 3a). Accordingly, we only consider the holeinitiated events in our modeling.
From the elaboration above, a carrier multiplication process occurred in this system can be described by Q = Q 0 × MF, where Q 0 is the initial quantity of charge, and MF denotes the multiplication factor. Q 0 is determined by the initial carrier density which multiplied by elementary charge e and the volume around the positive electrode (see "Methods"). Q is found through integrating the photocurrent through the whole response temporal profile obtained by the oscilloscope. Then, MF under different field strength can be deduced. (Hollow red circles in Fig. 3b). On the other hand, MF can be evaluated by Monte Carlo simulation for different carrier relaxation time τ. Screening due to the multiplied carriers must be considered, for it changed the dielectric property of material, and affected the transmitted THz field through p-doped region (see "Methods"). Figure 3c exhibits the THz pulses with an initial field strength 200 kV/cm inside the material with (red curve) and without (gray curve) screening effect respectively. Apparently, the generated carriers greatly wakened the THz pulse, particularly at the tail.
We found that MF obtained from our experiment was sandwiched between two MF curves obtained by Monte Carlo simulation with τ = 75 fs and τ = 100 fs, respectively (green dotdot line for τ = 75 fs and blue dot line for τ = 100 fs in Fig. 3b), and believe this was a reasonable result which is consistent with the experimental result of~100 fs relaxation time for holes 33 .
So far, it can be concluded that the simplified Monte Carlo model provides a consistent physical interpretation and is in good semiquantitative agreement with the experimental results. In order to present the simulation more vividly, we depicted the carrier density multiplication dynamics during the THz-induced acceleration under peak field of 100 kV/cm (blue curve in Fig. 3d) and 200 kV/cm (red curve in Fig. 3d) respectively, so are the corresponding energy evolutions of carriers (Fig. 3e). By comparison, it is distinct that there happened one impact ionization event under 100 kV/cm incident THz field, while four times under 200 kV/cm. This model further tells us that this linear trend will persist till the built-in field is offset by the potential difference induced by THz excited photocarriers, such that the saturation voltage will reach more than 3 V based on the built-in potential inside the blue LED. In short, we fortunately observed such gigantic photovoltaic signals in blue LEDs that stimulated much interest in investigating the impact ionization process. Although this process is common in some typical devices like photomultiplier, the large and effective photocarrier multiplication is often obscured by other effects, such as phonon absorption 34 , valley scattering 35 , and exciton dissociations 36 . Huge photovoltaic signals observed here show that the carrier multiplication process induced by intense THz pulses is efficient and not obscured by other effects. We believe studying this field-induced high efficiency of photocarrier multiplication is a very important driving force for future material physics and device optimization in the THz frequency range.
THz photovoltaic response in different color LEDs. We next studied four other LEDs with different colors (green, white, yellow, and red, see Fig. 4a-d) and observed negative signals in all LEDs. Data were taken for θ = 0°(see Supplementary Fig. 4) with THz field strengths of 160, 126, and 54 kV/cm. The optimal responsivity of the green LED reached 1250 V/W. Since all the materials for these LEDs are in the III-V group (GaN-based for blue and green LEDs, GaP-based for yellow and red LEDs), the carrier lifetimes were roughly identical. Therefore, the decay time of all examined LEDs were in a similar time scale (~1 ns), as extracted in the experiments. In terms of the comparison of amplitude for distinct LEDs, early works focused on the nanosecond far infrared laser pumped photodiodes showed lager photoresponse when the bandgap of material is smaller 37 . Nonetheless, we didn't observe such simple law in the current system that many ingredients needed to be taken into account which may impact on the performance of LEDs profoundly, such as the junction depth and the doping level 38 . Thus, the complicities of LED detection system truly enriched the physical significance in the THz strong-field region.
The most striking aspect is that the sign of the THzphotoresponse of all these LEDs is opposite to (see Fig. 2b) that of conventional photovoltaic signals. To explain this, we propose a Schottky contact based mechanism, as shown in Fig. 4e. Under a simple approximation, there will be a single abrupt junction between n + like metal and p-type semiconductor (Schottky contact), which leads to the generation of an electric field E 1 with direction opposite to the electric field E 2 in the p-n junction. THz induced photocarriers move towards the n region, and positive signals are also detectable; see Fig. 4c, d. Since the THz pulses first interact with the surface region, the induced photocarriers have a much higher density in the p region than in the n region; thus, the Schottky-based negative photovoltaic signals must be stronger than those positive signals flowing towards the n region. Moreover, because of the parasitic capacitance, the negative signals are obviously faster Corresponding ultrafast photovoltaic signals monitored by the oscilloscope (load = 50 Ω). The maximum photovoltaic signal is~600 mV. c Anisotropic photovoltaic signals obtained under different THz energies of 80.0, 52.8, 24.7 μJ, respectively. The photovoltaic response has strong correlation to the crystal symmetry. However, it has no strong correlation to the THz energy. Two specific THz polarization dependent responses labeled as Perpendicular (1) and Parallel (2) in c are systematically measured by illuminating various THz pump fluences, as exhibited in (d). When the THz pump fluence is <1.2 μJ/mm 2 , both polarization dependent photovoltaic signals obey quasi-quadratic relationship with respect to the THz pump fluence, while a linear behavior is predominant for the higher fluence region. e Photoresponse comparison between the blue LED and the commercial pyroelectric detector (SDX-1152, Gentec). The photoresponse signal and response time from the LED is 11 times larger and 4 orders of magnitudes faster than those obtained in the pyroelectric detector, respectively. The error bars are defined as standard deviations which are calculated automatically via an oscilloscope after 32 samplings. than positive signals, further implying the existence of a Schottky contact, which fortunately improves the rising time up to a few hundred picoseconds. Finally, we even observed a small positive signal in red LEDs earlier than abovementioned two signals. Among these LEDs, because the barrier height of the Schottky junction is relatively low in the red LED, hot holes within the valence band are easier to emit from the boundary of p side, resulting in a small positive signal. Such small positive signals among other LEDs were not observed because the Schottky barriers are higher and the hot holes cannot generate. Since this positive signal is not used in the intense THz detection, it does not affect the performance of the LED detectors. More experiments are needed to fully understand the intensive THz-matter interactions in such devices.
According to the proposed model in Fig. 4e, it is possible to observe a positive photovoltaic signal from LEDs if their structures did not contain Schottky junctions. We observed such a phenomenon from LEDs produced by the same company with 850 nm and 940 nm emission wavelengths. As shown in Supplementary Note 8, appreciable positive photovoltaic signals were observed in the 850 nm LED, which can also be well understood by the proposed model. In this case, there is no Schottky junction replacing with the ohmic contact. We also observed a saturation behavior in the 940 nm LED, which has a lower saturation voltage than the 850 nm LED as expected because the saturation value is determined by the built-in potential which is proportional to the bandgap.
As a detector working in the strong-field region, LEDs have some advantages over other commercially available THz detectors. Both Golay cells and pyroelectric detectors suffer from a slow response time with approximately several hundred microseconds, while bolometers usually work under low temperature conditions. Schottky barrier diodes (SBDs), which are widely used in the radio-and microwave frequency ranges, are high-speed devices (<100 ps) but require advanced material growth and device fabrication techniques. Field-effect transistors usually need a direct-current bias between the gate and source. Thus, LEDs have the advantages of high speed, broadband response, small size, low cost, no bias requirement, easy fabrication, and room-temperature operation. With regard to the Schottky contact discussed in the LEDs, it mainly affects the polarity of the photovoltaic signal and response time. Therefore, we cannot name it as a Schottky barrier diode because SBDs usually respond to a sub-cycle electric field component with a different detection mechanism. Moreover, the surface contact does not affect the THz-LED detection mechanism because we observe a positive photovoltaic signal in GaAs LEDs (see Supplementary Fig. 7). THz-LED camera prototypes. Since LEDs can be used for THz detection, one can straightforwardly think about fabricating a THz-LED camera. We demonstrated two prototypes. One of them is a scanning THz-LED camera. The other is a onedimensional sing-shot THz-LED array. In the context of scanning prototype, a blue LED was mounted onto an automatically controlled three-dimensional translation stage with 25 μm spatial resolution. We exploited this prototype to scan the LED at the focal plane, obtaining the profile of a THz beam, as shown in Fig. 5a. The obtained circular beam profile has a diameter of~2 mm (FWHM), which agrees well with that imaged by the commercial camera (Spiricon-IV, Ophire) (see Fig. 5b). LED displays and various large screen devices are widely available at low costs, indicating the feasibility of a large-area, high-resolution THz-LED real-time camera. Furthermore, we extended the THz-LED camera to a one-dimensional 1 × 6 array. With this LED prototype, we achieved a real-time THz camera, and the recorded THz beam profile is depicted in Fig. 5c. In this first-generation device, due to the large size of the market available LEDs, we could not directly image a focused THz beam. However, we could measure the THz beam profile in a non-focal plane, which in turn proved the high responsivity of the THz-LED detectors. The photographs of these prototypes are exhibited in Fig. 5d, e, respectively. With these results, we believe that LED-based THz detectors and cameras may have some valuable applications in strong-field THz science and technology.
Conclusions. In summary, we showed that gigantic photovoltaic signals were produced in market available LEDs through ultrafast detection of picosecond THz transients with moderate electric field strengths. Compared with conventional photovoltaic signals, an unexpected abnormal negative signal was observed for LEDs with different colors examined in our experiments. We systematically tested several kinds of LEDs in THz polarization and pump fluence dependent studies. Experimental results can be well understood by a simplified Monte Carlo model. Based on the high responsivity, ultrafast response time, low cost, and easy integration, we further demonstrated scanning and integrated THz-LED camera prototypes. The LED-based THz detection not only help us further deeply understand intense THz matter interaction physics, but also enables other applications such as THz real-time video imaging.

Methods
THz generation and characterization. In our experiment, a commercial Ti:sapphire laser amplifier system (Pulsar 20, Amplitude Technology, Institute of Physics, Chinese Academy of Sciences) was employed to produce strong-field THz pulses in a LN crystal based tilted pulse front technique 39 . This laser can deliver maximum single pulse energy of~500 mJ with central wavelength of 800 nm, pulse duration of 30 fs at 10 Hz repetition rate. In order to obtain high optical-to-THz efficiency, we optimized the pulse duration by altering the group velocity dispersion of an acousto-optic programmable dispersive filter 13 (AOPDF, Dazzler, Fastlite). In the tilted pulse front setup, a grating with density of 1480 lines/mm and of 140 × 140 × 20 mm 3 (length×width×height) was used for tilting the intensity of the pump pulses. A half-wave plate and a plano-convex lens with f = 85 mm focal length were inserted between the grating and the LN crystal. The LN was a z-cut congruent crystal with 6.0 mol% MgO doped to improve its damage threshold. The crystal was a triangle prism with dimensions of 68.1 × 68.1 × 64 mm 3 in x-y plane and the height in z-axis was 40 mm. For the LED detection experiment, the maximum pump energy used was~60 mJ and the LN crystal was not cryogenically cooled, resulting in a maximum THz electric field of~240 kV/cm with a spectrum from 0.1 to 0.8 THz. The radiated THz signal was first collimated by a 90°off-axis parabolic mirror (OAP1) and then went through two THz polarizers (Tydex) which were used for tuning the THz pump energy. After the second parabolic mirror (OAP2), the THz pulses were focused onto the LED which was connected to an oscilloscope (DPO 4104, Tektronix), in which a photovoltaic signal was recorded.
For THz temporal waveform and spectrum characterization, we used a singleshot spectrum coding method. As illustrated in Fig. 1a, when there was no LED in the setup, the THz waves were collimated again by OAP3 and focused together with the probing beam which came from the zero-order reflection of the grating into a 1 mm thick ZnTe detector. The probing beam was first stretched by a grating pair (Grating reticle density of 1200 lines/mm, 50 × 50 × 10 mm 3 ) and then propagated through a delay line. With this method, the THz temporal electric field vector was encoded onto the chirped spectrum of the probing beam in the ZnTe crystal through electro-optic effect. The modulated probing beam spectrum was recorded by an analysis system including a lens, a BBO, a Glen-prism (GP), a spectrometer and a CCD camera, and the THz temporal waveform can be decoded and obtained. More details can be found in Supplementary Notes 1-4. The THz single pulse energy was detected by a pyroelectric detector (SDX-1152, Gentec), and its beam profile was extracted by a commercial THz camera (Pyrocam IV, Spiricon). Thus, we can estimate the error of the THz peak field strength with the information from the measured THz energy, beam size, and temporal duration (see Supplementary Note 4). All these measurements were conducted at room temperature and the THz system was not purged.
Screening effect. For a given carrier density, the dielectric permittivity of a conductive medium in the Drude model is given by, where ω p ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Ne 2 =ε 0 ε 1 m* p is the plasma frequency, τ is the carrier momentum relaxation time, and N is the density of free carriers, ε 0 is the vacuum dielectric constant, ε ∞ is the high frequency limit of the dielectric permittivity, and m* presents the carrier effective mass. The square root of the dielectric permittivity is the complex refractive index. Thus, the carrier density in the medium determines its refractive index. From Eq. (1), it is easy to see that the complex dielectric permittivity is frequency dependent, and the lower frequency components affect the dielectric response greater than that of higher frequency components. Thus, the variation of frequency-dependent refractive index due to impact ionization is crucial to be involved in the modeling of the ultrafast carrier generation.
By introducing the screening effect into the GaN-based system, we set the parameters m* = 0.2m (considering the light hole of GaN in reference 40 ), where m is electron rest mass, ε ∞ = 5.35 (in reference 41 ), N = 5 × 10 16 cm −3 , which is a typical hole density in p-doped GaN-based blue LED 42 . The initial refractive index n 0 ∼3 was referred to the experimental results in reference 43 such that the initial transmitted THz waveform is E(t) = 2E inc (t)/(1 + n 0 ), where E inc (t) is the incident THz pulse, and E(t) is the transmitted THz pulse at the interface between GaN and air. Furthermore, as to obtain the full frequency-dependent refractive index, we implemented the Fast Fourier Transform (FFT) against the THz temporal profile and calculated frequency-dependent refractive index every time when the carrier density changed. After calculating the frequency-dependent Fresnel transmission e EðωÞ ¼ 2 e E inc ðωÞ=ð1 þ e nðωÞÞ, inverse FFT was exploited to, in turn, obtain the subsequent THz temporal profile, as exemplified in Fig. 3c. Eventually, only one free parameter i.e., carrier momentum relaxation time τ was left to reproduce the experimental results.
Monte Carlo simulation. To simulate the physical process of the carrier multiplication, a Monte Carlo simulation is established. In the simulation, a carrier is propagated in a linearly polarized THz electric field with screening (see the "Screening effect" in "Methods"). The electric field strength is adjusted such that the generated carrier density and the impact times under different incident field are recorded to obtain the multiplication factor. A carrier is accelerated in the THz field, and the evolution of its energyε ca ðtÞ is calculated under the parabolic approximation ε ca ðtÞ ¼ h 2 k 2 ðtÞ=ð2m * Þ in steps of 20 fs. At every time step, the scattering rates of different scattering mechanisms are calculated considering the energy of free carriers and are referred from literature. The total scattering rate can be written as γ tot ðε ca Þ ¼ γ ph ðε ca Þ þ γ imi ðε ca Þ, where γ ph ðε ca Þ is the energy dependent phonon scattering rate 33 , and γ imi ðε ca Þ is the energy dependent impact ionization rate 32 . From the first principle calculation in reference 32 , the impact ionization rate of holes in GaN-based system is about two orders of magnitude lager than that of electrons, such that we only consider the hole-initiated impact ionization rate.
The scattering rates of above two channels are utilized to determine the scattering probabilities at every step, and in turn determine the scattering event which the hole proceeds. If the hole scatters with the optical phonon, it will lose 100 meV energy 33 , and the wavevector will change accordingly to the value of parabolic approximation. But as seen from the energy evolution in Fig. 3e, phonon does not reduce the carrier energy significantly. If impact ionization occurs, the number of holes doubles as in the process new electron and hole generate. The carrier will lose its energy. The carrier density is dynamically upgraded through the simulation by taking the Auger recombination and impact ionization according to the following equation, where ξ is a stochastic number determined by whether impact ionization happens or not. If impact ionization happens, ξ will receive a value of 2 or 1 otherwise. The Auger recombination rate is given by γ Auger ðNÞ ¼ ÀCN 3 , where C is Auger coefficient 44 . In the Monte Carlo simulation, 1000 times are simulated at every given THz field strength. The average MF is taken as the final output. A carrier momentum relaxation time τ is tuned as a parameter to reproduce the experimental results.
Initial quantity of charge estimation. What we actually measured was the photovoltaic signal with an external load of R = 50 Ω. In order to correlate the simulation and experimental results, the quantity of charge was calculated such that the photovoltaic signal was divided by the load value, following the integration through the whole photoresponse temporal profile to get Q = ∫S(t)/R · dt, where S(t) is the photovoltaic signal. For the geometry of device, because of the negative signal feature, the detected carriers should be close to the surface at the positive electrode. We believe the carriers near the negative electrode can be ignored due to the recombination of carriers and longer diffusion distance. It can also be verified experimentally that no positive signals were observed in the blue LED. Thus, we can roughly estimate the initially generated quantity of charge near the positive electrode, choosing the electrode size to be~50 μm, doping concentration 5 × 10 16 cm −3 , diffusion length~1 μm in reference 45 , and the thickness of p-doped region 0.2 μm (Fig. 3a). Combining the ingredients above, we obtained 8 × 10 −14~1 0 −13 C quantity of charge. According to the experimental integration, the detected quantity of charge under 95 kV/cm field strength was 4.6 × 10 −13 C which approaches the estimated value the most. Therefore, we choose integrated value under 95 kV/cm field as initial quantity of charge which corresponding to the MF = 1 in Fig. 3b.
Applicability of the simplified monte carlo model. Inspection of the anisotropic curves in Fig. 2c, which is consistent with Baraff's prediction that the impact ionization process is mainly due to lucky carriers at low applied electric fields 46 . In other words, these carriers will reach threshold in some certain directions more easily than others. Thus, this phenomenon reflects that the LED certainly worked at the low-field region. Only Shockley's 'lucky carriers' 47 can be accelerated to high energy level over threshold and initiated impact ionization. Others will remain below threshold. Nonetheless, how many 'lucky carriers' are lucky enough to reach this goal is unknown. It is difficult to be solved and beyond the model in this work.
In terms of other factors that may affect the applicability of our model, we need to conduct more discussions. The assumptions of this model base on a uniform system which does not include information about defects and other channels, such as an impact ionization event occurring between a free carrier and a confined carrier within an impurity state, an impurity band or confined quantum-well state. For defects, we actually measured three regular and symmetric anisotropic two-fold photoresponse curves and plotted them in Fig. 2c. However, the defects can prevent carriers from attaining high multiplication 48 and render the photoresponse irregular while the anisotropic photoresponse of blue LEDs is relatively regular. Therefore, we can conclude that the defects affect little and attest to the applicability of the model in this sense.
With regard to the second factor, we pick up the confined quantum-well state as an example, the impact ionization between the free carriers and confined carriers will result in the space-charge neutrality no longer being maintained 49 . This will alter the electric field profile within the device and subsequently the device performance. However, the measured Volt-ampere characteristics before and after the illumination of THz pulses exhibit no hysteresis phenomenon (see Supplementary Note 6), implying the confined quantum-well state affect little on this model as well. This is a reasonable result because the activation layer (MQWs part) is much thinner than other layers. The confined states are very few such that the second effect plays a small role in the total carrier multiplication process. In addition to the influence from the confined states, we eliminated the possibility arise from resonance excitations as well, which is detailed in Supplementary Note 7.
Recovery process in blue and green LEDs. After the carrier multiplication process, combining the previous theory with the large negative photovoltaic signal we have observed, we can infer the carrier transport process in the blue LED and green LEDs. Because the strong-field THz pulse excited carriers at 15 ps, it was quite short compared with the carrier lifetime of a few nanoseconds. The carrier multiplication can be regarded as an instantaneous process. In the blue and green LEDs, since there is no positive signal at the tail of the photovoltaic signal, it can be deduced that most of the carriers inside the device were mostly generated in the p-type region. Due to the presence of the surface electric field (Fig. 4e), the generated electrons drifted outwards. The holes diffused into the bulk. Electrons accumulated near the surface of positive electrode side, while holes accumulated in the middle of the p-type region. Subsequently, the accumulated electrons and holes formed an electric field opposite to the surface electric field direction.
These fields partially canceled out the original surface electric field. In the external circuit, the electrons on the surface reduced the potential at the positive electrode, forming a large negative signal. The internal carrier recombination process gradually reduced the surface potential, and finally decayed to zero while all excess carriers recombined.
LED material and structure. In our experiment, the LEDs we examined were the most common LED lamps in the market for children's toys, TV sets, monitors, telephones, computers, and circuit warnings. We mainly selected several kinds LEDs with different colors of blue, green, white, yellow and red for systematical investigations. The materials of blue and green LED lamps are InGaN with 3 mm round type (4204-10SUBC/C470/S400-X9-L, 04-10SUGC/S400-A5, EVERLIGHT), while the white, orange, and red ones are AlGaInP (204-15/FNC2-2TVA, 204-10UYOC/S530-A3, 4204-10SURC/S530-A3, EVERLIGHT). The blue LED was exemplified, and its parameters were referred (https://www.powerwaywafer.com/ epitaxial-wafer.html?gclid=EAIaIQobChMIr8qD2ZGd4wIVGIfVCh0wYwA0EAA YASAAEgJi-vD_BwE) in Fig. 3a. From top to bottom, the structure is composed by different layers with various thicknesses of 0.2 μm thick p-GaN, 0.03 μm thick p-AlGaN, the activation layer of 0.2 μm thick InGaN/GaN, 2.5 μm thick n-GaN, 3.5 μm thick u-GaN, and a 430 μm thick Al 2 O 3 substrate. For this LED, it emitted blue light when a voltage of 3 V was applied. When no voltage was applied and the THz pulse illuminates onto the LED, we can probe a photovoltaic signal. We also saw similar phenomena in LEDs with other colors, but the responses were different.
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.

Code availability
The code that supports the findings of this study are available from the corresponding author upon reasonable request.