1. Introduction

Strong terahertz (THz) few-cycle field may change drastically properties of a bulk medium [1,2]. Such phenomena as ionization [3], bleaching [4], the Kerr effect [5], etc. were observed at the THz field strength below 1 MV/cm. A relatively strong THz pulse can induce symmetry breaking in a centrosymmetric medium, thus allowing for the THz pulse characterization using second harmonic (SH) generation of an optical probe [6,7,8]. Here the SH field E_2ω of the probing optical field E_ω is generated due to the third order non-linearity: ${E_{2\omega }}{ \propto {\chi^{(3 )}}}{E_{THz}}E_\omega ^2$. The same process can be observed in any medium that is transparent for the THz field, optical field, and its SH [9]. In the non-perturbative mode this provides for the generation of even-order high harmonics for imaging of operating electronic circuits [10].

Recently a few cycle THz pulses (central frequency ω_THz ∼1.5-1.8 THz) with the pulse energy of 0.9 mJ and field strength as high as E_THz>40 MV/cm were reported [11]. This radiation was used for studying a thin metallic layer breakdown [12], n-Si self-induced transparency [13], and other phenomena. Preliminary data on the non-linear optical probing of p-Si under the action of such an extreme THz pulse were presented recently in [14]. The optical SH yield exhibits complex dependence on the time delay between the pumping THz and probing optical pulses, but this issue was not touched in [14].

In this paper starting from the elaborated analysis of the THz pump-optical SH probe temporal data we suggest physical mechanisms for the extreme THz field induced non-linearity in the optical domain. We proposed that the delayed SH yield is due to the deformation-dependent susceptibility, i.e., lattice strain caused by the ponderomotive force of the THz pulse. This breaks the central symmetry of a medium, thus inducing field dependent χ⁽²⁾. Such a mechanism comes into play if the field strength E_THz is above 1 MV/cm (hereafter we reference the E_THz quantity as the amplitude of the strongest half period of the THz field, second one in our case, see Fig. 1(a)). Lattice deformation increases up to the Lindemann limit [15], long range order vanishes, the medium becomes centrosymmetric again, and the associated SH signal drops down. We also introduced the ionization induced photocurrent mechanism of SH production that comes into play at field strengths above ∼6–8 MV/cm. It determines the instantaneous SH yield from the first strong half-period of the THz field, while the delayed deformation-dependent susceptibility establishes later, during the second half period. This ionization nonlinearity arises from a weak modulation by the optical probing field of the impact ionization rate by free carriers moving in an extreme THz field.

 

Fig. 1. Output THz field waveforms (red curves), the input THz field waveform (E_THz∼22.8 MV/cm, x0.25, black dashed curve in Fig. 1(a)) and SH yield temporal shapes (blue curves) at the THz field strength E_out ∼5.7 (a), 12.3 (b) and 14.5 (c) MV/cm.

Download Full Size | PDF

2. Time-dependent experimental data analysis

The detailed description of the experimental setup was presented in [13,14]. Briefly, a few cycle THz pulses were generated through the optical rectification of 120 fs, 1.24 µm, 6 mJ pulses in the DSTMS crystal [16]. The maximum energy of the THz pulse reaches 100 µJ corresponding to the E_THz ∼22 MV/cm in the focal spot. The field waveform of the THz pulse was characterized by the electro-optical sampling (EOS) using 100 µm GaP crystal after an appropriate attenuation of the signal before the crystal (see black dashed curve in Fig. 1(a)). The energy of the THz pulse measured with the Golay cell was varied by changing the optical pulse energy. It was checked that the THz pulse waveform does not depend on the optical pulse energy. The small part of the optical pulse (1 µJ) was used as a probe. The delay between the THz and optical probe pulses was changed within 3 ps window covering THz pulse field envelope. Polarization of both fields were linear and parallel to each other. The SH of the optical probe at 612 nm measured after the sample was used as a nonlinear probing signal. The sample under study was a p-Si (100) wafer 250 µm thick (carrier density 1.6 × 10¹⁵ cm⁻³). Since the p-Si sample is transparent at 1240 nm and opaque at 612 nm the SH signal comes from the thin, 2 µm, layer at the exit surface of the sample.

The EOS traces of the THz field after the sample provide us with both clear waveform and strength of the acting THz field E_out (see the red solid curves in Fig. 1(a)-(c)). The latter value refers to the maximal THz field inside the sample near its exit surface, calculated following [17] and considering the exit surface THz transmittivity. It is even more essential since the sample transmission by energy for the THz pulse under study changes with its input amplitude (see section 3.2 and Fig. 3(a)). Propagation effects and self-phase modulation of the THz pulse (observed for the thin but opaque Si sample in [18]) as well as the phase mismatch did not play a role in our study. Traces of the THz field before the sample (at the highest field strength of 22 MV/cm, dashed black line in Fig. 1(a)) and after it (at the “low” output field strength E_out ∼5.7 MV/cm, solid red line in Fig. 1(a)) are nearly identical except for the small ∼10% difference in the amplitudes of the second negative field oscillation and small changes at t >1.7 ps. There are certain changes in the relative amplitudes of the THz field oscillations after the sample with the THz field strength (compare red lines in Figs. 1(a)-(c)), and we will address this issue below (Fig. 3(b)). Note also that no surface or bulk damages were observed for the sample even under the highest field strength used.

Let us now present first the detailed temporal dynamics of the experimentally measured SH yield. Figure 1 shows how the SH yield depends on time at different THz field strengths along with the output THz field waveform. The SH yield is strongly non-monotonic in time. It has three main maxima at the “low” field strength E_out ∼5.7 MV/cm (Fig. 1(a)). Their temporal positions coincide approximately with positions of THz field’s maxima, but the first and second maxima are slightly delayed. The maximal SH yield occurs from the second positive half-period of the THz field ($T_2^ + $ in Fig. 1(a)), that is almost twice less than the first positive half-period ($T_1^ + $). Note, that the first strong (negative) half-period of the THz field ($T_1^ - $) induces a very tiny SH signal. Besides, the SH yield for the second negative half-period ($T_2^ - $) is dumped with respect to the SH yields from the $T_1^ + $ and $T_2^ + $ half-periods of the field.

The SH signal changes at the higher field of E_out ∼12.3 MV/cm (Fig. 1(b)). Now one can see only one strong maximum corresponding to the $T_1^ + $ half-period, while the subsequent half-periods induce much weaker SH yield with almost equal amplitudes. There is also a peak at the left shoulder of the main maximum of the SH yield. And it is this peak that dominates at the highest field strength of E_out ∼14.5 MV/cm (Fig. 1(c)). The maximum corresponding to the $T_1^ + $ half-period becomes weaker, while other maxima have exceedingly small amplitudes.

Figure 2 shows how the amplitudes of those main maxima (labelled A_i where i is the maximum number, see Fig. 1(b)) in the SH yield temporal shape depend on the THz field strength E_out. There are a few remarkable features in this plot. Firstly, the A₁ amplitude is almost zero at “low” field strengths E_out < 6 MV/cm. Second, the A₂ value experiences fast growth with power dependence ${A_2} \propto E_{out}^{2.7}$ if E_out is less than 6–7 MV/cm, and A₃ and A₄ have the similar scalings, ${A_{3(4 )}} \propto E_{out}^{2.3({2.4} )}$, below E_out ∼6 MV/cm. This power is slightly higher than that characteristic of the usual χ⁽³⁾ mechanism, where the SH signal should be proportional to $E_{out}^2$ (or should be linear with a THz field intensity). Increase in the A₃ and A₄ amplitudes becomes weaker within 6–8 MV/cm range with average scaling of ${A_{3,4}} \propto E_{out}^{1.7}$ within the 2–8 MV/cm range. The A₁ value starts growing if E_out>6 MV/cm. Next, increase in the A₂ and A₃ values saturates at E_out∼8 MV/cm. The A₄ drops down here. The A₁ growth is exceptionally fast, ${A_1} \propto E_{out}^{4.2}$, i.e. it is proportional to the THz intensity in power ∼2.1. Finally, the A₂ starts decreasing at E_out>12.5 MV/cm and at E_out >13 MV/cm the A₁ value becomes higher than the A₂ one.

 

Fig. 2. The SH yield in the four main maxima A₁-A₄ labelled in Fig. 1(b). Solid lines show A_i∝ $E_{\textrm{out}}^{{b_i}}$ fits with b₁= 4.2 ± 0.4, b₂=2.7 ± 0.1, b₃=2.3 ± 0.3, b₄=2.4 ± 0.2. The black dashed line represents the fit obtained from the numerical simulations (see section 3.4). The vertical dashed lines depict ranges of the THz field strength with remarkably different behaviour of the A₁-A₄ amplitudes.

Download Full Size | PDF

Hence, there are four ranges of E_out in the time domain with different non-linear response (see vertical dashed lines in Fig. 2): (i) the 2–6 MV/cm range where the response looks like the usual ${E_{2\omega }}{ \propto {\chi^{(3 )}}}{E_{THz}}E_\omega ^2$ one, but delayed (the $T_1^ - $ half-period of the THz field produces a very weak SH yield); (ii) the transient narrow 6–8 MV/cm range where the A₂-A₄ yields grow slowly, while the A₁ signal is still weak but grows fast; (iii) the 8–12 MV/cm range, where the A₂-A₄ SH yields saturate or even decrease, while the A₁ experiences even faster growth; (iv) above 12–13 MV/cm range where even the A₂ yield drops down and the main impact to the SH yield comes from the $T_1^ - $ half-period of the field. The non-linearity is high and has distinctly different physical origin in the latter range.

3. Modeling and discussion

There are two key mechanisms for the SH generation from a centrosymmetric medium in optical and IR domains under action of an external DC or slowly varying electromagnetic field: non-linearity from the symmetry breaking by this field due to valence electrons’ polarization [19] and non-linearity from the symmetry breaking caused by the current of free carriers arose due to the applied DC field [20,21]. Our data shows that the $T_1^ - $ half period produces very weak SH signal at “low” driving field strengths. This signal comes from the ${E_{2\omega }}{ \propto {\chi^{(3 )}}}{E_{THz}}E_\omega ^2$ process in the transparent semiconductors as in [22,23]. Note that the large Kerr non-linearity observed in the thin Si wafer opaque to the THz pulse [18] should also lead to the efficient SH generation, but we consider the transparent Si here. Hence the induced non-linearity is rather slow in our case, while the above mentioned non-linearities are almost instantaneous (see also recent data at the much lesser field strength in [22,23]).

Below we are considering two additional mechanisms under the action of an extreme THz pulse with the maximal field strength well above the breakdown of Si by a DC field [24]. The first one is the local lattice strain by the ponderomotive force of a THz pulse resulted in deformation-dependent susceptibility $\chi^{(2 )}({E_{THz}^2} )$ and generation of the SH signal ${E_{2\omega }}{ \propto \chi^{(2 )}}({E_{THz}^2} )E_\omega ^2$ [25]. Secondly, a strong THz field causes ionization in semiconductors [26,27], the induced ionization current is strongly non-linear and could provide a source of the SH of the probing optical field. Below we will consider these two mechanisms and show how impacts of different non-linearities might give an explanation to the data observed.

3.1. Local deformation by the ponderomotive force (E_THz>2 MV/cm)

The SH generation observed in our experiment can be partially explained using the notion of the dielectric susceptibility ɛ dependent on the local deformation tensor $\widehat {\sigma}$. Indeed, the electric field E_THz contributes to the energy density of a crystal as ${\raise0.7ex\hbox{${{\varepsilon}E_{THz}^2}$} \!\mathord{\left/ {\vphantom {{E_{THz}^2} {8\pi }}} \right.}\!\lower0.7ex\hbox{${8\pi }$}}$. In the linear approximation, a part of the energy associated with a strain reads ${\raise0.7ex\hbox{${E_{THz}^2}$} \!\mathord{\left/ {\vphantom {{E_{THz}^2} {8\pi }}} \right.}\!\lower0.7ex\hbox{${8\pi }$}}{\frac{\partial{\varepsilon}}{{{\partial {\sigma}_{ij}}}}{\sigma}_{ij}}$. Thus, the electric field makes lattice strain energetically favorable. Since the scalar part of the deformation tensor is a relative variation of the elementary lattice volume, the THz electric field induces a local pressure gradient proportional to the ${\raise0.7ex\hbox{${E_{THz}^2}$} \!\mathord{\left/ {\vphantom {{E_{THz}^2} {8\pi }}} \right.}\!\lower0.7ex\hbox{${8\pi }$}}\textrm{Tr}\left( {\frac{\partial{\varepsilon}}{{{\partial{\sigma}_{ij}}}}} \right)$ and directed along the wavevector of the THz field. This pressure gradient can be treated as the ponderomotive pressure of a short THz pulse [28].

The SH generation from a strained surface layers of Si was first discussed for implanted layers [25], and then for the Si-SiO₂ interfaces [29,30]. The huge few orders of magnitude increase in the SH yield was observed. The THz field experiences prominent growth near the exit surface of the sample in our case [27] thus enhancing its pressure exactly inside the SH emitting layer. Besides, the ponderomotive force also contains the term proportional to the $E_{THz}^2\nabla{\varepsilon}$, that might also strongly impact near the exit surface.

A crude estimation might be done by comparing the binding energy W≈1.8 eV in the c-Si lattice and the maximal energy difference Δ of two electrons separated by d∼0.13 nm (interatomic distance in Si) and subjected to the 10 MV/cm THz field:

Lattice strain becomes larger if the THz field strength increases and may approach the Lindemann limit (∼10% of interatomic distance), and this should lead to a subsequent loss of a long-range order [15]. Hence at a high enough field strength the medium might become centrosymmetric again after a certain time and the SH signal vanishes. This correlates well with the observed behaviour of the A₃ and A₄ amplitudes (Fig. 1(b) and (c), Fig. 2). The drop in A₂ at E_out>12 MV/cm might be linked with the fact that the higher the THz field strength the earlier in time lattice deformation approaches the Lindemann limit.

3.2. Impact ionization by the extreme THz field

A high amplitude THz field might give rise to ionization in semiconductors. There are two possible mechanisms of ionization – Zener effect (or tunnel ionization) [31] and impact ionization by free carriers gaining energy from the THz field [27]. The Zener effect is negligible for the Si being an indirect bandgap semiconductor [32]. In the following we are using impact ionization model described in detail in [27]:

Here ${n_\textrm{e}},v,q,m$ are carrier density, velocity, charge and effective mass respectively, $E_{\textrm{THz}}^{(0 )}(t )$ is the THz pulse envelope. This model treats a carrier movement in the THz field considering dependencies of both electron-phonon $\eta_{\textrm{ph}}(v )\; $ and ionizing $\eta_{\textrm{ion}}(v )\; $ collisions rates on the instant carriers’ velocity $v\; $, as well as carriers’ density ${n_\textrm{e}}\; $ dynamics. Initial free carriers are holes for a p-Si. The impact ionization rate for holes in p-Si is order of magnitude less than for electrons [33,34], but any ionization event creates a new electron-hole pair, thus free electrons come into play very fast and govern the ionization process later on [35]. Note that the inelastic scattering rate is quite high, ∼10¹³s⁻¹, in our study.

Carrier density can be estimated from the experimentally measured and calculated drops in the Si sample’s integrated energy transmission measured with the Golay cell (Fig. 3(a)). In this figure we also plotted calculated carrier density achieved right after the THz pulse termination. One can see good coincidence of the experimental and calculated data. Maximum carrier density amounts to ∼2 × 10¹⁶ cm⁻³.

 

Fig. 3. Experimental (open black squares) and calculated (full red squares) data on the Si sample’s normalized integrated energy transmission by the THz field along with calculated carrier density inside the sample (full blue circles) (a) and experimental ratios R of amplitudes of the $T_1^ - $, $T_2^ + $ and $T_2^ - $ half-periods to the amplitude of the $T_1^ + $ half-period (b) in dependence on the E_out.

Download Full Size | PDF

Note, that the background doping plays significant role setting up the initial conditions for the ionization process. The electron density increases 14 times at the highest field strength considered in our study. Hence, the ionization is negligible if the initial density is 10 times smaller than that in our sample, and vice versa it could be catastrophic (i.e., creates opaque plasma for the THz pulse) at higher doping levels.

Calculations also show that ionization is transient, and the main impact comes from the $T_1^ + $ half-period of the field. The THz field experiences damping within $T_2^ - $ and $T_2^ + $ half-periods since ionization changes both real and imaginary parts of the refractive index of the sample. This is shown in Fig. 3(b), where we plotted dependencies on the E_out of amplitudes of $T_1^ - $, $T_2^ - $ and $T_2^ + $ half-periods normalized to the $T_1^ + $ amplitude. One can see that the $T_2^ \pm $ half periods experience damping while the $T_1^ - $ - even increases. This is most likely due to the damping of the $T_1^ + $ half period since the $T_1^ - $ half period is almost unaffected by ionization. In total these amplitudes’ behaviour describes well the experimentally observed drop in the Si sample transmission.

Comparing data in Figs. 2 and 3 one can conclude that there exists a strong correlation between changes in the carrier density and SH yields A₂-A₄. Ionization changes transmission if the field strength is above 5–6 MV/cm, and this correlates well with changes in the slopes of SH yields A₂-A₄ (see Fig. 2). Further, strong saturation of the yield A₂ and drop in the A₃ and A₄ yields is linked to increase in the carrier density by a few times from the initial one and decrease of amplitudes of $T_2^ + $ and $T_2^ - $ - half periods for the field strength of 8–12 MV/cm. Finally, even the A₂ yield declines and this obviously correlates with the prominent drop in the sample transmission and fast rise in the carrier density at the field strength above 12 MV/cm.

Nevertheless, changes in the transmission with the THz field strength are arbitrary small and cannot explain strong changes in the A₂-A₄ amplitudes, especially for the intermediate THz field strengths. In [14] the non-linearity from symmetry breaking by the current of free carriers under action of an extreme THz pulse was considered as one of the causes of optical SH generation in the initially centrosymmetric medium. Authors [14] concluded that the ionization current mechanism is feasible even at the highest THz field used but the induced effective χ⁽²⁾ susceptibility is few times weaker than the ordinary χ⁽³⁾E_THz one. Moreover, more accurate estimations give even much lower effective χ⁽²⁾: the ionization current is determined by the carriers’ mobility that decreases at least by one order under action of the extreme THz field due to the huge increase in the velocity dependent elastic collision rate [27].

3.3. Ionization induced SH generation

The most intriguing feature of the data obtained is emerging of the SH signal from the very first half period of the THz field followed by its extremely fast increase at the E_out∼8–10 MV/cm. Certainly, this is manifestation of an alternative mechanism of the Si non-linear response under action of an extreme THz field, and we attributed it to the ionization induced photocurrent mechanism. Note again, that the field strength inside the sample is well above the breakdown limit of the silicon by a DC electric field (< 1 MV/cm).

This mechanism’s essence can be understood from the following. A photocurrent J inside the Si is proportional to ${\raise0.7ex\hbox{${{q^2}}$} \!\mathord{\left/ {\vphantom {{{q^2}} m}} \right.}\!\lower0.7ex\hbox{$m$}}{n_\textrm{e}}(t){\varepsilon}(t )$ with the optical pulse ${\varepsilon}(t)= E_{\textrm{opt}}^{}(t )\textrm{exp}({ - i{\omega_{\textrm{opt}}}t} )$. The impact ionization rate and therefore carrier density ${n_e}$ in Eq. (1b) depend on the full kinetic energy ${K_f} = \frac{{m{v^2}}}{2}$ of an electron. This energy should oscillate at the optical frequency ${\omega _{\textrm{opt}}}$ (and its second harmonic) hence providing for the $2{\omega _{\textrm{opt}}}$ photocurrent source. To calculate the K_f we used Eq. (1a) with an additional weak optical driving force in its righthand side (here we assumed linear and parallel polarizations of both fields):

Figure 4(a) shows the normalized spectral power density of ${n_\textrm{e}}$ (t) obtained from the numerical solution of Eqs. (1b) and (2) with experimental envelopes $E_{\textrm{THz}}^{(0 )}(t )$ and $E_{\textrm{opt}}^{}(t )$ and maximal field strengths of 14.5 and 20 MV/cm respectively. The latter value corresponds to the field strength of the probe optical pulse in our experiments, but it is “weak” if the kinetic energy of an oscillating charged particle is considered. One can see the first and second harmonics of the optical field here. The spectral amplitude A_opt at the fundamental frequency ${\omega _{opt}}$ scales as $E_{\textrm{opt}}^{1.95 \pm 0.05}$ with the optical field strength (Fig. 4(b), black squares and line) and strongly depends on the field strength of the THz field (red squares and lines). The fit gives a power of 2.8 ± 0.4 for the 7<E_out<12 MV/cm range and 4 ± 0.2 for the 12<E_out<15 MV/cm range for the latter dependence.

 

Fig. 4. Spectral power density S of the ${n_e}(t )$ from numerical solution of Eqs. (1b) and (2) with E_out=14.5 MV/cm and E_opt=20 MV/cm (a) and dependence of the A_opt maximum on the E_out (red squares at E_opt=20 MV/cm) and on the E_opt (black squares at E_out=14.5 MV/cm) (b). Lines in Fig. 4(b) show A_opt=a(E_out(opt))^b fits.

Download Full Size | PDF

The SHG signal scaling with the optical field strength (${\sim} E_{\omega}^{1.9}$) looks like the normal intensity squared dependence, characteristic for the standard χ⁽²⁾ process. This is not surprising since we assume that the optical field is weak (in terms of kinetic energy gained from the field). Hence the SH source looks like $\chi^{(2)}({E_{THz}^2} )E_\omega ^2$. Much faster growth of the SH signal with the THz field strength comes from the complex dependence of the ionization rate on the kinetic energy (and hence the field strength). This is also confirmed by the fact that the deduced electron density scalings from Fig. 3(a) (not shown) are remarkably close to the scalings from Fig. 4(b).

Hence the non-linear ionization induced current at 2 ${\omega _{opt}}$ is proportional to $E_{\textrm{out}}^{2.8 \ldots 4}E_{\textrm{opt}}^{1.9}$. Substituting this current source into the Maxwell equations one should obtain the second harmonic yield proportional to $E_{\textrm{out}}^{5.6 \ldots 8}E_{\textrm{opt}}^{3.8}$. The scaling obtained corresponds well to the observed scaling for the A₁ yield in Fig. 2. Here dashed black line was calculated using fits from Fig. 4(b) within 7–12 and 12–15 MV/cm ranges. This line was linked to the vertical scale in Fig. 2 at E_out=9.5 MV/cm only. Note that the experimental curve A₁(E_out) is obviously more complicated, but our simplified model describes the main features well (the “threshold” for the SH generation, overall dependence on the E_out).

4. Conclusions

Hence our study opens up a new intriguing area of research, where an extreme ultrashort THz pulse changes medium’s properties drastically. Non-linear optical probing of interaction of an extreme few cycle THz pulses with the p-doped Si reveals an essential role of (i) lattice deformation by the ponderomotive force and (ii) impact ionization. The central symmetry of the medium (p-Si) vanishes, and SH is generated due to these two processes.

The extremely strong few cycle THz pulse causes lattice deformation after a single half-oscillation of the THz field (∼300–500 fs). The strain induced second order optical non-linearity plays major role in the SH generation by a weak ultrashort optical probe even at the lowest THz fields ∼2–3 MV/cm used in our study. This optical SH signal appears at the second strongest half-period of the THz field and following half-periods but is negligible for the first half-period that is only twice less than the second one. The SH yield grows fast, ∼$E_{\textrm{out}}^{2.3 - 2.7}$ for all the 2^nd, 3^rd and 4^th half-periods of the THz field and then saturates at ∼7–8 MV/cm. SH yields corresponding to the 3^rd and 4^th half periods of the THz field decrease with further rise in the E_out. This is due to (i) increase in the lattice deformation up to the Lindemann limit (∼10% of interatomic distance) then medium becomes centrosymmetric and (ii) onset of impact ionization that increases absorption mainly for the 3^rd and 4^th half periods of the THz field. Estimations showed that the carrier density is doubled at E_out ∼8 MV/cm. Both effects strengthen with the E_out increase to 10 MV/cm and higher, and this is the first cause of the SH yield vanishing from the 3^rd and 4^th half periods of the THz field and the SH yield decreasing from the 2^nd (strongest) half-period.

Impact ionization is responsible for the new SH generation mechanism that manifests itself at E_out ∼6–8 MV/cm as the extremely fast-growing SH yield from the 1^st half period of the THz field. The SH yield here rises up as $E_{\textrm{out}}^5$ and faster. We attributed this effect to the ionization induced photocurrent mechanism that can be reproduced only with an extreme few cycle THz pulses. Here the kinetic energy of a free electron moving in the THz field is modulated by the weak optical field. This leads to modulation in the impact ionization rate at the fundamental frequency ${\omega _{\textrm{opt}}}$ and in its turn to the modulation in the free carrier density. Hence non-linear ionization current gained modulation at 2 ${\omega _{\textrm{opt}}}$ and this is a new source for the SH generation. Our simplified model shows good qualitative coincidence both in the threshold and growth rate for this effect. It should play an essential role also for the 2^nd half period of the field.

We considered the both mechanisms of the THz-induced non-linearity within frameworks of phenomenological semi-quantitative models. Much more insight can be obtained from the Monte-Carlo and/or Maxwell-Bloch approaches accounting for the real p-Si properties with its complex energy band structure, different channels of inelastic scattering, phono-electron coupling, etc. It is of great interest to study the same effect in direct bandgap semiconductors, where the Zener effect will play a role instead of impact ionization [22,36]. The interplay between deformation and ionization related mechanisms depends badly on the mutual polarization of the THz and optical fields, crystal orientation, relative amplitudes of the THz field half-periods and its number. It is also interesting to measure the SH yield from the samples transparent both at the fundamental optical frequency and its SH, as well as for the polar semiconductors. Interestingly the THz field itself should also be enriched by the harmonics of fundamental THz frequency (without an optical field) under ionization.