---
author:
- 'I. Babushkin'
- 'A. J. Galan'
- 'V. Vaičaitis'
- 'A. Husakou'
- 'F. Morales'
- 'A. Demircan'
- 'J. R. C. Andrade'
- 'U. Morgner'
- 'M. Ivanov'
title: 'All-optical attoclock: accessing exahertz dynamics of optical tunnelling through terahertz emission '
---

**The debate regarding attosecond dynamics of optical tunneling has so far been focused on time delays associated with electron motion through the potential barrier created by intense ionizing laser fields and the atomic core. Compelling theoretical and experimental arguments have been put forward to advocate the polar opposite views, confirming or refuting the presence of tunnelling time delays. Yet, such delay, whether present or not, is but a single quantity characterizing the tunnelling wavepacket; the underlying dynamics are richer. Here we propose to complement photo-electron detection with detecting light, focusing on the so-called Brunel radiation – the near-instantaneous nonlinear optical response triggered by the tunnelling event. Using the combination of single-color and two-color driving fields, we determine not only the ionization delays, but also the re-shaping of the tunnelling wavepacket as it emerges from the classically forbidden region. Our work introduces a new type of attoclock for optical tunnelling, one that is based on measuring light rather than photo-electrons. All-optical detection paves the way to time-resolving multiphoton transitions across bandgaps in solids, on the attosecond time-scale.**

Main
====

Ionization induced by intense low-frequency laser fields is the key resource in attosecond science [@corkum07; @kim14]. The so-called optical tunnelling of the electron through the potential barrier created by the oscillating electric field and the binding potential of the core lies at the heart of high harmonic generation and high harmonic spectroscopy [@baker06; @smirnova09; @bruner15; @torres10; @pedatzur15] generation and characterization of attosecond pulses [@dudovich06], sub-femtosecond laser-driven electron diffraction and holography [@huismans11; @spanner04; @meckel08; @blaga12; @schell18; @calegari14], ionizaton-triggered charge migration in molecules, etc (see e.g. [@blaga09; @smirnova09; @yan10; @haessler10; @huismans11; @lepine14]). Most of these phenomena, such as high harmonic generation, responsible for converting intense infrared driving field into extreme ultraviolet and soft X-ray radiation, are based on laser-driven electron return to the parent ion.

Interestingly, even when the electron does not return back to the core, where it can recombine with the hole left by ionization and release the energy gained from the laser field, radiation can still be emitted. It is often referred to as the “Brunel radiation” or the ”Brunel harmonics” [@brunel90; @babushkin17; @balciunas13]. The Brunel radiation is associated with the bursts of electron current triggered by ionization. These ionization bursts occur every time the modulus of the oscillating electric field goes through a maximum, i.e. twice per laser cycle. Overall, the bursts of ionization current are balanced: the electrons librerated every half-cycle are driven by the field in the opposite directions. The emitted radiation is associated with the time-derivative of this current, which is sharply modulated by the ionization bursts and contains odd harmonics of the fundamental field [@brunel90].

An important new frequency component of the Brunel emission emerges when the incident laser field combines fundamental radiation (frequency $\omega$) with its second harmonic (frequency $2\omega$). Even a small admixture of the $2\omega$ field breaks the symmetry between the two subsequent half-cycles of the fundamental field, leading to the emergence of even harmonics. These even harmonics also include the ’0-th’ order Brunel harmonic. This harmonic corresponds to the emission frequencies $\sim 1/T$, $T\sim 10^2-10^4$ fs, which fall within the terahertz (THz) frequency range for (e.g [@babushkin11; @kim08b]).

Since the Brunel radiation directly encodes ionization dynamics, this radiation can be used to resolve and reconstruct it. Here we demonstrate this possibility and show that different frequency components of the Brunel radiation generated in single-color and two-color laser fields provide detailed dynamical picture of optical tunnelling. The ionization dynamics is encoded in the polarization ellipse of the Brunel radiation. Measurement of the polarization ellipse allows one to get access to several parameters characterizing the electron wavepacket as it emerges into continuum. In analogy with characterizing optical pulses, we gain access not only to the “ionization time-delay”, the analogue of the group delay for optical pulses, but also to the distortions of the wavepacket – the analogue of the higher order dispersion for optical pulses.

Our approach to imaging the ionization dynamics is distinctly different from those based on detecting photo-electrons, be it in the multi-photon (strong field) or one-photon ionization regime [@eckle08a; @eckle08; @huismans11; @shafir12; @kienberger04; @pfeiffer11; @landsman14; @torlina15; @schultze10; @kluender11; @gruson16; @heuser16]. These photo-electron based approaches, which currently dominate, face an important challenge: the electrons are measured well after the ionization step. Thus, the ionization dynamics has to be deconvolved from the additional effects imparted on the outgoing electron as it flies to the detector.

In the strong field regime, the attoclock approach to measuring ionization dynamics [@eckle08a; @eckle08] uses the mapping between the moment at which the electron “escapes the tunnelling barrier” and its final velocity at the detector. The rotating electric field of the nearly circularly polarized pulse “sprays” the electrons, which are escaping the atom at different times, into different detection angles (see [Fig. \[fig:base\]]{}a). Interpretation of the photoelectron spectra in terms of time-delays associated with the ionization step is a subject of active discussions [@yakaboylu14; @landsman14; @landsman15; @ni16; @teeny16; @camus17; @kaushal15; @torlina15]. One of the key difficulties in the interpretation stems from the need to deconvolve ionization and post-ionization dynamics.

All optical techniques based on measuring high harmonic emission [@blaga12; @zhao13; @landsman14; @larsen16; @pedatzur15; @bruner15] face similar challenge: the need to decouple the ionization and the photo-recombination steps. Thus, the Brunel radiation has an important potential advantage, since it is produced at the ionization step. Moreover, the available precision in measuring angle-resolved photoelectron spectra $\sim$ 1$^{\circ}$ is orders of magnitude lower than that achievable when measuring the light polarization ellipse, not only for optical but also for the THz frequencies [@ikebe08; @nemoto14]. As a result, the lowest (THz) frequencies in the Brunel radiation can provide immediate, [*in-situ*]{} [@kim14] attosecond-resolved information, even though the period of the associated light oscillations is orders of magnitude longer. In particular, the ionization delay is imprinted in the “0-th” order Brunel harmonic, whereas the reshaping of the ionization distribution is encoded in the 3rd harmonic. Detecting Brunel radiation instead of detecting electrons also paves the way to extending attosecond measurements of multiphoton transitions across bandgaps in solids, while boosting the measurement accuracy to the sub-attosecond range by taking advantage of the extraordinary accuracy available in polarization measurements in THz, IR, and optical frequency range.

We consider atomic hydrogen, since it captures all relevant physics while being simple and amenable to numerically exact calculations. An H-atom interacts with the laser field made of an intense circularly polarized infrared (IR) pump with frequency $\omega_0$ and its co-rotating second harmonic $2\omega_0$ of equal amplitude, see [Fig. \[fig:base\]]{}b. The nonlinear response contains harmonics $n\omega_0$, with $n=0,1,2\ldots$. All harmonics with $n\geq 1$ rotate in the same direction, while the signal near $n=0$th harmonic is almost linearly polarized, with ellipticity increasing as one moves to higher frequencies within the broad “0-th harmonic” line (see Methods, including the details of the analytical model). For short pulses, the “$n=0$ harmonic” lies in the terahertz (THz) range, reaching into the mid-infrared for few-cycle pulses. Below we show that the orientation of the corresponding polarization ellipse identifies the ionization delay $\tau$: the rotation angle $\phi$ with respect to the $y$-axis of [Fig. \[fig:base\]]{}b is $\phi=-\omega_0\tau$. This property is preserved throughout the 0th harmonic signal, extending up to the optical frequencies.

Further details about the ionization dynamics can be obtained if higher-order Brunel harmonics are considered. We show that, for a single-color elliptically polarized field (see [Fig. \[fig:base\]]{}c), the polarization state (ellipticity and polarization direction) of the 3rd Brunel harmonic encodes the attosecond reshaping of the ionization distribution with respect to the field maxima.

In general, harmonics of intense infrared driving fields can be generated via three different mechanisms: the bound-free transitions associated with sub-cycle ionization bursts and yielding the Brunel radiation, the bound-bound transitions (non-resonant frequency mixing) which manifest in the perturbative regime, and recollision-based harmonics [@corkum07] responsible for high orders. Their relative contributions vary for different harmonics orders. For our approach, we seek harmonics dominated by the Brunel radiation. [Fig. \[fig:sy\]]{}a,c shows the contribution of each of these mechanisms as a function of the harmonic order for a co-rotating 800-400 nm bi-circular driver and an 800 nm single-color elliptical driver (ellipticity $\epsilon$=0.6), respectively.

The black circles show ab-initio simulations of the time-dependent Schrödinger equation (TDSE) for the hydrogen atom, using the method of [@patchkovskii16]. Here all three mechanisms are naturally included. The red squares show the results of ab-initio simulations for a short-range (Yukawa) potential, with a single bound state having the same ionization potential as the hydrogen atom (see Methods). In particular, this calculation allows us to evaluate the contribution of bound-bound transitions, which are absent in the Yukawa potential, and calibrate out the effects of the core potential on the optical attoclock in a way similar to that of Ref.  [@torlina15]. The blue triangles are analytical results for the Brunel emission computed using the Drude model. There, the response ${\mathrm{\mathbf{E}}}_r$ is calculated as the product of the electric field ${\mathrm{\mathbf{E}}}(t)$ and the free electron density $\rho(t)$: ${\mathrm{\mathbf{E}}}_r(t) = {\mathrm{\mathbf{E}}}(t) \rho(t)$ (atomic units are used; see Methods for details of calculations). Finally, the grey crosses in [Fig. \[fig:sy\]]{}a,c show the strength of the conventional bound-bound contribution: the response is calculated as ${\mathrm{\mathbf{E}}}_r = \chi^{(3)} |{\mathrm{\mathbf{E}}}|^2{\mathrm{\mathbf{E}}}$, with the coresponding susceptibility $\chi^{(3)}=0.8\times10^{-20}$ W/cm$^2$ obtained from Ref. [@koehler13].

These results allow us to separate the ab-initio harmonic spectra into three regions, depending on the dominant mechanism. Region I ([Fig. \[fig:sy\]]{}a,c blue shading) corresponds to linear response. In region II ([Fig. \[fig:sy\]]{}a,c green shading), the Brunel mechanism, evaluated using the Drude model, dominates. This region contains the 0th and 3rd harmonics in both the two-color and the single-color case. The influence of bound-bound transitions in region II, as well as in the other two regions, evaluated directly from the $\chi^{(3)}$ model, is several orders of magnitude smaller than the Brunel mechanism, see [Fig. \[fig:sy\]]{}a,c. Similarly, the influence of the recollision-based contribution for 0th harmonic is negligible, and is few orders of magnitude below the Brunel mechanism for the third harmonic. Conservatively, we have placed the beginning of the region III ([Fig. \[fig:sy\]]{}a,c orange shading) beyond $n=3$. For $n\geq 5$ the Drude model yields rapidly decreasing response [@brunel90; @babushkin11], while the ab-initio calculations for both the Coulomb and the Yukawa potential are developing plateau-like structure, suggesting that recollision-based contribution is important for $n\geq 5$. We note that the difference in the harmonic intensity between the Coulomb and the Yukawa potentials comes from the fact that the tunnelling barrier is higher in the latter case, reducing the ionization probability by orders of magnitude compared to the Coulomb potential.

![ \[fig:base\] Ionization dynamics and delay observed by detecting electrons (a) and the Brunel harmonics (b,c) using two- (b) and single-color (c) fields. In (a), angle ($\theta$)-resolved photo-electron spectra generated by elliptically polarized laser field (blue line) are detected, yielding the electron density $\rho(\theta)$ (blue gradient) at a remote detector (black line). It has the maxima shifted from the extrema of the electric field due to a nonzero ionization delay time $\tau$, providing the “attoclock” measurement of the ionization delay. In (b), this delay is probed, instead of electrons, by the 0th-order Brunel radiation ${\mathrm{\mathbf{E}}}_r$ in a two-color field ${\mathrm{\mathbf{E}}}$. The Brunel radiation is linearly polarized (red line) and, in the presence of a nonzero ionization delay $\tau$, rotated by the angle $-\tau\omega_0$. In (c), a single-color elliptically polarized pump ${\mathrm{\mathbf{E}}}$ produces the 3rd Brunel harmonic ${\mathrm{\mathbf{E}}}_r$, with the polarization state (ellipticity and polarization direction) encoding the sub-cycle dynamics of the ionization process. The spectra below show schematically the harmonic structure of the atomic response. The blue lines show the pump frequencies, the red lines indicate the harmonics of interest. ](base){width="\columnwidth"}

Two-color field
---------------

![image](fig2_latest.pdf){width="17cm"}

We now focus on 0th and 3rd order harmonic (region II), where the Brunel mechanism dominates. For the 2-color case, the 3rd harmonic is nearly circularly polarized in all calculations, owing to the angular momentum selection rules. In contrast, the polarization of the broad emission line corresponding to the “0-th” harmonic is much closer to linear. For the Yukawa potential, the ab-initio results agree very well with the analytical predictions based on the Drude model. The polarization ellipse of the emission is strictly parallel to the polarization ellipse of the driving laser field, in agreement with the analytical model based on zero ionization delay.

For the Coulomb potential, the polarization ellipse of the “0-th harmonic” is rotated by $-\phi\approx 10.0$ degree, which corresponds perfectly to the so-called ionization delay $\tau=-\phi/\omega_0 \approx 74$ as calculated for the hydrogen atom for sufficiently strong field using the standard attoclock technique [@torlina15]. Such ionization delay is caused by electron deflection in the Coulomb field, disappearing in the short-range Yukawa potential. For sufficiently strong fields, the rotation of the polarization ellipse becomes intensity- independent, just like in the conventional attoclock measurements.

To confirm that the rotation of the polarization ellipse is indeed related to the ionization delay, we introduce ionization delay $\tau$ into the Drude model of the Brunel radiation: $$\label{eq:delayed_drude}
 {\mathrm{\mathbf{E}}}_r(t) = {\mathrm{\mathbf{E}}}(t)\rho(t-\tau)$$ where $\tau$ is the ionization delay (see Methods for the derivation).

Equation \[eq:delayed\_drude\] can be rewritten in the frequency domain as: $$\label{eq:14}
   {\mathrm{\mathbf{E}}}_{r,n} = \sum_j
 {\mathrm{\mathbf{E}}}_j\rho_{n-j}e^{i(n-j)\omega_0 \tau},$$ where $\rho_j$ and ${\mathrm{\mathbf{E}}}_j$ are the $j$-th Fourier components of the the free electron density and of the pump field, respectively. In Jones vector notation and for the two-color field, the latter can be expressed as, up to a constant factor (see Methods): ${\mathrm{\mathbf{E}}}_1 \sim (1,i)$, ${\mathrm{\mathbf{E}}}_{-1} \sim (1,-i)$, ${\mathrm{\mathbf{E}}}_2 \sim (1,i)$, ${\mathrm{\mathbf{E}}}_{-2} \sim (1,-i)$. For the 0th-order we have thus: ${\mathrm{\mathbf{E}}}_{r,0} = {\mathrm{\mathbf{E}}}_1 \rho_{-1} + {\mathrm{\mathbf{E}}}_{-1} \rho_{1} + {\mathrm{\mathbf{E}}}_2 \rho_{-2} + {\mathrm{\mathbf{E}}}_{-2} \rho_{2}$, which leads to (see Methods): $$\label{eq:11}
  {\mathrm{\mathbf{E}}}_{r,0} \sim (-\sin(\omega_0 \tau),\cos(\omega_0\tau)).$$ Thus, the ionization delay $\tau$ does indeed translate into the rotation of the polarization. For $\tau=0$ the 0th harmonic is linearly polarized in the $y$-direction, while for nonzero $\tau$ the polarization vector is rotated by $\omega_0\tau$ along the field rotation direction (see [Fig. \[fig:base\]]{}b), just as observed in our TDSE simulations, [Fig. \[fig:sy\]]{}b.

The deviation from zero ellipticity seen in the TDSE and the Drude model calculations is due to the non-zero frequency at which we monitor the response (around 100 THz), with excellent agreement between the Drude model and the TDSE simulations for the Yukawa potential; in fact, the ellipticity monotonously increases with frequency.

Another important observation is the much more narrow polarization ellipse of the Brunel radiation from the Coulomb potential, compared to the Yukawa potential. We see that the optical response encodes the Coulomb focusing of the departing electron wavepacket [@brabec96; @yudin01b]: the core potential clamps the transverse impulse and thus the transverse component of the current driven by the rotating laser field.

As one can see from [Eq. (\[eq:11\])]{}, apart from the ionization delay, the polarization state of the 0th harmonic [Eq. (\[eq:11\])]{} does not depend on the particular temporal shape of $\rho(t)$. As we shall see, for higher Brunel harmonics this is no longer the case.

Single-color field
------------------

![image](fig3_new.png){width="\textwidth"}

We now show that the 3rd harmonic provides further information about the ionization dynamics. For the 2-color driver, the 3rd harmonic is circularly polarized, and thus the polarization direction is undefined. However, the situation is different for the single-color elliptical pump ${\mathrm{\mathbf{E}}}_1\exp(-i\omega t)+{\mathrm{\mathbf{E}}}_{-1}\exp(-i\omega t)$, [Fig. \[fig:base\]]{}c. Here ${\mathrm{\mathbf{E}}}_1= (1,i\epsilon)$, ${\mathrm{\mathbf{E}}}_{-1}= (1,-i\epsilon)$ and $\epsilon$ is the ellipticity. In this case, the 0th harmonic is absent, but the 3rd harmonic lies in region II (see [Fig. \[fig:sy\]]{}b), that is, the Brunel radiation dominates the atomic response. According to Eq. [Eq. (\[eq:14\])]{}, the 3rd Brunel harmonic is ${\mathrm{\mathbf{E}}}_{r,3} = {\mathrm{\mathbf{E}}}_1 \rho_{2} + {\mathrm{\mathbf{E}}}_{-1} \rho_{4}$, which maps on an electric field given by (see Methods for details): $$\label{eq:1}
    {\mathrm{\mathbf{E}}}_{r,3} \sim \left(1+r/2,
      i\epsilon(1-r/2)\right), \, r=2\rho_4e^{2i\tau\omega_0}/\rho_2.$$

[Eq. (\[eq:1\])]{} shows that the polarization state of the 3rd harmonic depends critically on the 2nd and 4th harmonics of the ionization burst $\rho(t)$. This explains different results for the ellipticity and the polarization direction in the Yukawa and the Coulomb potentials (both analytical and the TDSE), since the subcycle ionization dynamics in these cases is different.

The map shown in [Fig. \[fig:del-simp\]]{} allows one to reconstruct the complex-valued parameter $r$ defined by [Eq. (\[eq:1\])]{} via the polarization direction and the ellipticity of the 3rd harmonic. We are interested in the dynamics of the ionization rate $W(t)$. In the frequency domain, and for small $\rho$, its harmonics $W_n$, are given by $W_n=in\omega_0 \rho_n$. Thus, [Eq. (\[eq:1\])]{} implies that $W_4/W_2=re^{2i\tau\omega_0}$. Comparing to the results of the TDSE simulations, we see that not only $r$ but also $W_4/W_2$ is complex. Its absolute value $|r|$ determines how sharp the ionization spike is. For instance, for $\delta$-function-like $W(t)$, $|r|=1$. Smoother ionization process corresponds to smaller $|r|$. As for the phase of $r$, we note that delaying $W(t)\to W(t-\tau)$ induces delay-dependent phase shift $W_n \to W_ne^{in\omega_0 \tau}$ in the frequency domain.

Our result shows that in addition to this spectral-phase shift, there must be an additional phase modification which is only present in higher harmonics. Both effects modify the ionization rate (represented by $\{W_n\}$) as follows: $$\label{eq:2}
    W_0\to W_0,\, W_2 \to W_2e^{2i\omega_0 \tau},\, W_4 \to W_4e^{4i
      \omega_0 \tau - i\delta}, \, \ldots$$ where $\delta$, according to the definitions above, is related to $r$ as: $\delta=-\arg{(r)} + 4\omega_0\tau$. It is easy to see that $\delta$ leads to an additional reshaping of $W(t)$ making it asymmetric in time with respect to the oscillations of the electric field, leading to non-instantaneous harmonic response.

Although the observation of 3rd harmonic provides only limited information about the ionization dynamics, it is already sufficient to reconstruct the dynamics with a fairly good level of precision. Indeed, we observe from our reconstruction that $|r|$ is quite small, namely $|r|=0.14$ for the Yukawa and $|r|=0.27$ for the Coulomb potential. It is reasonable to expect, that the next term in ionization is of the order of $|r|^2$ which gives us 7% uncertainty from unknown $W_6$ for the Coulomb potential and 2% for the Yukawa one. Thus, we can neglect $W_6$ and the higher terms in the present case. Finally, the value of $W_0$ can be estimated from the ionization rate at the minimum of our elliptical field. Having this in mind and also using the known ionization delay time $\tau$, we are able to construct the ionization rate as a sum of harmonics $W_0$, $W_2$, $W_4$ with the corresponding amplitudes and phases. The result of the reconstruction is shown in [Fig. \[fig:del-simp\]]{}c. One can see that the ionization for the Coulomb potential is not only time-shifted but also shows a certain asymmetry in respect to the ionization maximum. Reconstruction from the Yukawa potential still displays a reshaping of the ionization rate, albeit to an almost negligible extent. Thus, we were able to demonstrate that the capabilities of the proposed technique go beyond the mere delay measurement, and include the reconstruction of the deep-subcycle features of the ionization rate.

![ \[fig:exp\] The power of 3rd harmonic (solid lines) as well as of the fundamental (dashed line) as a function of the polarizer angle for the ellipticity of the pump pulse $\epsilon=0.6$ and various pump powers (encoded in the line color). The inset shows the polarization angle of the 3rd harmonic relative to the pump in dependence on the pump power. In the yellow-shaded region the Brunel radiation dominates the $\chi^{(3)}$-based one yet the ionization degree is still small. ](exp-fig.pdf){width="50.00000%"}

Our theoretical analysis predicts an interesting new phenomenon – rotation of the polarization ellipse of the nonlinear optical response in a gas of atoms or randomly oriented, achiral molecules. This rotation is confirmed by our experimental measurements shown in [Fig. \[fig:exp\]]{}, where we observed the third harmonic generated in a plasma spot created by a strong elliptically polarized pulse focused in air. The experimental method is described in detail in the Methods section. After the plasma onset (around 0.1 mJ pulse energy) the polarization angle of 3rd harmonics varies around 10-13 degrees (see yellow box in the inset to [Fig. \[fig:exp\]]{}). Deviation for still higher energies above $\sim$0.28 mJ is attributed to the saturation of ionization. Even though the theory analyzed the response of a Hydrogen atom, which is amenable to exact numerical simulations, and the experiment used ambient air, the predicted effect is thus confirmed not only qualitatively but also quantitatively. We note that the rotation of (higher order) harmonics were discussed before [@weihe96; @strelkov12; @becker97]. Here we not only explain the physics of this effect, but also show how to use it to recover the underlying electron dynamics.

Conclusion and Discussions
==========================

In conclusion, we have shown that the Brunel radiation, including the radiation in the THz range, is able to provide crucial information about attosecond-scale details of strong-field ionization dynamics. In particular, this means the ability to measure attosecond dynamics using a much slower, picosecond, clock. For the 0th-order Brunel harmonic in the case of a 2-color field, one attosecond of ionization time delay translates into about 0.1 degree of polarization emission ellipse rotation. Taking into account that the polarization can be determined with better than milliradian precision, this is a promising mechanism for sub-attosecond time resolutions. The polarization state of the 0th harmonic demonstrates some indications of the Coulomb focusing in its polarization ellipsis. Furthermore, we have shown that in the case of single-color driver the polarization and the ellipticity of the 3rd Brunel harmonic is very sensitive to the temporal shape of the ionization step and allows for its reconstruction. Reconstruction from the TDSE results yields both the noticeable temporal asymmetry and the delay of the ionization burst.

These features demonstrate the promising capabilities of the proposed new tool, which allows one to explore attosecond ionization dynamics — and even beyond it — in a simpler and more general setup than the electron-detection based one. The key potential of this method is its applicability beyond the gas-phase ionization dynamics, e.g. within the solid state context, where detection of light is much simpler than the detection of electrons and where attosecond-scale delays were recently discussed [@hofmann15; @sommer16; @hassan16].

Methods
=======

Experimental setup
------------------

For the experiments we have used a 1 kHz repetition rate femtosecond Ti:sapphire chirped pulse amplification laser system (Legend elite duo HE+, Coherent Inc.), delivering 35-40 fs (FWHM) light pulses centered at 790 nm with maximal pulse energy of 7 mJ, which was used to pump an optical parametric amplifier (OPA, Light Conversion, Inc.). Though the OPA output was tunable within the broad wavelength range, for the experiments we have used the infrared (IR) radiation with the central wavelength of about 1.1 $\mu$m. The duration and single-pulse energy of the OPA pulses were about 40 fs and 0.25-0.3 mJ, respectively. The ellipticity of the IR wave was controlled by rotating the quarter-wave plate around the optical beam axis. Then the IR pulses were focused in air at the atmospheric pressure by a lens with a focal length of 2 cm and a plasma filament of about 0.1-1 mm length was produced, where the third harmonic generation (THG) took place. The THG signal was filtered from the pump radiation, collimated and through the computer-controlled broadband rotary Glan polarizer was sent to the photodiode. By rotating the Glan polarizer we were able to record the pump or TH power as a function of their polarization angle. The angle between the incident and generated TH polarization ellipses can be obtained from the phase shift of these nearly sinusoidal functions. Besides, we also compared the polarization of the fundamental frequency before and after the plasma spot. This comparison showed the difference of few degrees between these two, much less than the rotation of the 3rd harmonics. We attribute this effect to the $\chi^{(3)}$- and plasma-induced circular birefringence, since this number also corresponds to our estimations of such. This result shows that the 3rd harmonics is influenced by the pump birefringence only insignificantly. Finally, we would like to pay attention to the simplicity of the experiment in comparison to the one utilizing electron direction measurements.

Pump fields in time and frequency representation
------------------------------------------------

In general, we consider one-color field in the form: $$\begin{gathered}
\label{eq:field1}
  E_{x}(t) = E_0e^{-2 \ln 2 \,t^2/\mathcal{T}^2} \cos{\omega_0 t},\\
  E_{y}(t) = -\epsilon E_0e^{-2 \ln 2 \,t^2/\mathcal{T}^2} \sin{\omega_0 t},
  \label{eq:field1a}\end{gathered}$$ where, in our simulations, $\omega_0$ corresponds to the wavelength 0.8 $\mu$m, $E_0=0.07$ au, the pulse duration is $\mathcal{T}=20$ fs and the ellipticity is $\epsilon=0.6$. For the two-color field we have $$\begin{gathered}
\label{eq:field20}
  E_{x}(t) = E_0e^{-2 \ln 2 \,t^2/\mathcal{T}^2} \left(\cos{\omega_0 t} + \cos{2\omega_0 t} \right),\\
  E_{y}(t) = -E_0e^{-2 \ln 2 \,t^2/\mathcal{T}^2} \left(\sin{\omega_0 t}+\sin{2\omega_0 t}\right),
  \label{eq:field2}\end{gathered}$$ with $\omega_0=0.8 \mu$m, $E_0=0.0375$ au, and $\mathcal{T}=20$ fs.

In deriving the equations for the Brunel harmonics response, we used the Jones vector representation for the field harmonics ${\mathrm{\mathbf{E}}}(t) = \sum_n {\mathrm{\mathbf{E}}}_ne^{in\omega}$. The waveshape described in [Eqs. (\[eq:field1\])-(\[eq:field2\])]{}, can be easily translated into ${\mathrm{\mathbf{E}}}_n$. For instance, for an exemplary field $E_x = E_0\cos{(\omega_0 t)}$, $E_y= -E_0\sin{(\omega_0 t)}$ (which corresponds to the fundamental harmonic of the 2-color field) we can easily obtain using the Euler formulas: ${\mathrm{\mathbf{E}}}_1 = E_0/2\cdot(1,-1/i)$, ${\mathrm{\mathbf{E}}}_{-{\mathrm{\mathbf{1}}}} = E_0/2\cdot(1,1/i)$ where the first and second component of the vector are x and y components, respectively. In Results, the unsignificant prefactor $E_0/2$ is removed since only the polarization state is discussed.

TDSE simulations
----------------

Our TDSE equation has the well known form (we use atomic units across the article): $$i \partial_t \psi({\mathrm{\mathbf{r}}},t) = H\psi({\mathrm{\mathbf{r}}},t), \, H=({\mathrm{\mathbf{p}}}+{\mathrm{\mathbf{A}}}(t))^2/2+V(r), 
    \label{eq:1a}
    $$ where $\psi({\mathrm{\mathbf{r}}}, t)$ is the wavefunction of the electron depending on space ${\mathrm{\mathbf{r}}}$ (with $r\equiv |{\mathrm{\mathbf{r}}}|$) and time $t$, $H$ is the electron’s Hamiltonian, with $p$ being the momentum operator, ${\mathrm{\mathbf{A}}}(t)$ the vector field potential of the external field, with the corresponding field strength ${\mathrm{\mathbf{E}}}(t)$. The simulation box size used for the integration was 500 au. Simulations were repeated with the uniform and logarithmic spatial grid, with comparison to ensure the numerical convergence.

The response of the atom, ${\mathrm{\mathbf{E}}}_r$ was calculated as: ${\mathrm{\mathbf{E}}}_r = g\partial_{tt}{\mathrm{\mathbf{p}}} =g\partial_t {\mathrm{\mathbf{j}}}$, where ${\mathrm{\mathbf{p}}}= {\ensuremath{\left\langle\psi\right|}}{\mathrm{\mathbf{r}}}{\ensuremath{\left|\psi\right\rangle}}$ is the atomic polarization, ${\ensuremath{\left|\psi\right\rangle}}$ is the electron wavefunction, ${\mathrm{\mathbf{j}}}=\partial_t{\mathrm{\mathbf{p}}}$ is the current assigned to ${\mathrm{\mathbf{p}}}$, and $g$ is the constant depending on the position of the observer, which we assume $g=1$ for simplicity. In [Fig. \[fig:sy\]]{}, the relevant quantities (polarization direction, ellipticity) were averaged around the frequencies $n\omega_0$ in the range $(0.75n\omega_0,1.25n\omega_0)$ for $n$th harmonic except $n=0$; Polarization direction and ellipticity being weighted with the intensity by averaging. For $n=0$, because of our numerical limitations (huge box in time and space needed, together with the resolution of large number of angular harmonics makes it a very difficult task), we excluded extremely low frequencies $\omega<0.26\omega_0$ from the analysis, thus restricting our predictions to the upper part of the THz frequency range. These frequencies still display a very good agreement with the calculations for the 0th harmonic for $\omega\to 0$.

From TDSE to the delayed Drude model.
-------------------------------------

In the framework of the analytic R-matrix approach [@torlina12; @torlina15] the electron wavefunction can be written as: $\psi({\mathrm{\mathbf{r}}},t) = \sum_{t_s,{\mathrm{\mathbf{p}}}} {\ensuremath{\left\langle{\mathrm{\mathbf{r}}}\right.\left|\, {\mathrm{\mathbf{p}}}\right\rangle}}
  R_{{\mathrm{\mathbf{p}}}} e^{iS({\mathrm{\mathbf{p}}},t, t_s)}$, with $R_{{\mathrm{\mathbf{p}}}}$ encodes the angular structure of the initial state, $t_s$ is the (complex-valued) stationary time (see below), and $S = S_{f} + S_C$, $S_{f}= 1/2\int_{t_s-i\kappa^{-2}}^t({\mathrm{\mathbf{p}}}+{\mathrm{\mathbf{A}}})^2 + I_pt_s$, $S_C=\int_{t_s-i\kappa^{-2}}^tV(\tau)d\tau$, where $\kappa=\sqrt{2I_p}$ and $I_p$ is the ionization potential. The ionization time $t_s$ is inferred via the stationary phase approach from the equation $\partial_{t_s}S=0$. $t_s$ is different from the ionization time $t_0$ in a short range potential defined from $\partial_{t_0}S_f=0$, and hereupon we express $\tau={\textrm{Re}\,}{(t_0-t_s)}$. We note that for hydrogen, the value of $\tau$ can be obtained as a perturbation from the relation $-\tau \partial_{t_s t_s}S_f=\partial_{t_s}S_C$ as $\approx$ 74 attoseconds (for strong enough field). By neglecting the term $\tau\partial_{t_0t_0t_0}S$ in the vicinity of $t_0$ we obtain that $S(t,{\mathrm{\mathbf{p}}},t')\approx S_f(t,{\mathrm{\mathbf{p}}},t')+C(t,{\mathrm{\mathbf{p}}})$ with $C$ being some constant in $t'$. This leads us to the following simplified problem to solve: $\psi({\mathrm{\mathbf{r}}},t) = \sum_{t_0,{\mathrm{\mathbf{p}}}} {\ensuremath{\left\langle{\mathrm{\mathbf{r}}}\right.\left|\, {\mathrm{\mathbf{p}}}\right\rangle}}R_{{\mathrm{\mathbf{p}}}}
  e^{iS_f({\mathrm{\mathbf{p}}},t, t_0-\tau)+iC(t,{\mathrm{\mathbf{p}}})}$. Furthermore, assuming that all electrons are created with ${\mathrm{\mathbf{p}}}=0$ making $C$ an insignificant common phase, allowing to finally obtain the classical [Eq. (\[eq:delayed\_drude\])]{} for ${\mathrm{\mathbf{E}}}_r(t)\equiv \partial_t {\mathrm{\mathbf{j}}}$, with $\rho(t) \approx \int_{-\infty}^t W(t'-\tau)dt'$ is the density of free electrons determined by the ionization rate $W(t)$ delayed by the time $\tau$. In our simple theory for $W(t)$ we used the so called quasi-static formula [@ammosov86].

We note that the delayed Drude model by itself does not fix the ionization time $\tau$, it must be obtained from one other consideration, for instance the one presented above. We note however that this model is also compatible with the Green-function approach [@teeny16; @yakaboylu14] which predicts the small positive delays due to the nonzero values of initial velocity at the exit from the classically-forbidden region. In this case, introducing the delay into the Drude model is in fact an effective shift of the classical trajectory to negative exit times, so that at time $t=0$ such shifted trajectory has already nonzero momentum.

Derivation of the Brunel harmonic polarization in the delayed Drude model
-------------------------------------------------------------------------

Delayed Drude model derived before allows to obtain analytically waveshapes of the Brunel harmonics in dependence on ionization delay and other dynamical properties of the ionization process as the following: To obtain the ionization rate, we assume that before we shifted it in time by $\tau$ or before introducing any further phase shifts, $W_n$ obeys the equation $W_n = W_{-n}$, that is, the corresponding function is symmetric in time in respect to the field extremum at $t=0$. The corresponding ionization rate, if ionization is small enough, obeys then the relation $\rho_n = -\rho_{-n}$. Introducing the time delay $\rho_n\to\rho_ne^{i\omega_0\tau}$ and substituting this expression into the formula for ${\mathrm{\mathbf{E}}}_{r,0}$ above Eq. \[eq:11\], we obtain then ${\mathrm{\mathbf{E}}}_{r,0}=
iE_0(\rho_1+\rho_2)\cdot (\sin{\phi},\cos{\phi})$ where $\phi=-\omega_0\tau$. Neglecting the insignificant prefactor, we obtain immediately Eq. \[eq:11\]. The linear polarization predicted by Eq. \[eq:11\] is easy to explain also invoking the angular momentum conservation. The processes $+2\omega_0-\omega_0-\omega_0$ and $+\omega_0+\omega_0-2\omega_0$ leading to 0th Brunel harmonics both have equal amplitudes and the opposite angular momenta of the low frequency component, giving in sum the angular momentum zero. We remark that the polarizaiton state of 0th harmonics is completely independent on the concrete values of $\rho_n$. That is, it is independent on ionization rate dynamics $W(t)$.

For $n=3$, this is not anymore so; the particular polarization state depends in general on $\rho_2$, $\rho_4$. The expression before [Eq. (\[eq:1\])]{} gives us immediately: ${\mathrm{\mathbf{E}}}_{r,3} = -E_0e^{2i\phi/2}\cdot (1+r/2,i\epsilon (1-r/2))$, where $r=e^{2i\phi}r(4)/r(2)$ which, by removing the common factor irrelevant for the polarization state gives us [Eq. (\[eq:1\])]{}.

[56]{}ifxundefined \[1\][ ifx[\#1]{} ]{}ifnum \[1\][ \#1firstoftwo secondoftwo ]{}ifx \[1\][ \#1firstoftwo secondoftwo ]{}““\#1””@noop \[0\][secondoftwo]{}sanitize@url \[0\][‘\
12‘\$12 ‘&12‘\#12‘12‘\_12‘%12]{}@startlink\[1\]@endlink\[0\]@bib@innerbibempty [****, ()](http://dx.doi.org/10.1038/nphys620) [****,  ()](http://dx.doi.org/10.1038/nature08253) [****,  ()](http://dx.doi.org/10.1038/nphys3436) [****,  ()](\doibase 10.1126/science.1198450) [****,  ()](\doibase
  10.1126/science.1157980) [****, ()](http://dx.doi.org/10.1038/nature10820) [**** (), 10.1126/sciadv.aap8148](\doibase 10.1126/sciadv.aap8148) [****,  ()](\doibase 10.1038/nphys1228) [****,  ()](\doibase 10.1103/PhysRevLett.105.253002) [****, ()](http://josab.osa.org/abstract.cfm?URI=josab-7-4-521) [****, ()](http://stacks.iop.org/1367-2630/13/i=12/a=123029) [****,  ()](\doibase 10.1038/nphoton.2008.153) [****,  ()](http://dx.doi.org/10.1038/nphys982) [****,  ()](\doibase 10.1126/science.1163439) [****,  ()](\doibase 10.1038/nature11025) [****,  ()](http://dx.doi.org/10.1038/nphys2125) [****,  ()](\doibase
  10.1364/OPTICA.1.000343) [****,  ()](\doibase
  10.1126/science.1189401) [****, ()](\doibase 10.1103/PhysRevLett.106.143002) [****,  ()](\doibase 10.1103/PhysRevA.90.012116) [****, ()](\doibase
  https://doi.org/10.1016/j.physrep.2014.09.002) [****,  ()](\doibase 10.1103/PhysRevLett.117.023002) [****,  ()](\doibase
  10.1103/PhysRevLett.116.063003) [****,  ()](\doibase 10.1103/PhysRevLett.119.023201) [****, ()](http://stacks.iop.org/0953-4075/48/i=23/a=234002) [****,  ()](\doibase
  10.1103/PhysRevLett.111.043901) [****, ()](\doibase 10.1103/PhysRevA.87.043811) [****,  ()](\doibase 10.1364/OPTICA.2.000151)

Acknowledgements
================

I.B., A.D. and U.M are thankful to Deutsche Forschungsgemeinschaft (DFG) (projects BA 4156/4-1, MO 850-19/1, MO 850-20/1) for financial support.

Author contributions
====================

I.B. incepted the idea, developed the delayed Drude model and preformed simulations and analytic calculations based on it; A.J.G preformed TDSE simulations; I.B., A.J.G, A.H., F.M., and M.I. analyzed and interpreted the results of the simulations; all authors participated in the article formulation and writing. V.V. preformed the experiments.
