Jean Ebothé
email: jean.ebothe@univ-reims.fr

Magnetic properties dependence of mesoscopic scale thick films on thickness and surface roughness characteristics: Case of nickel samples in Ni/Pt system by

Keywords: Mesoscopic scale thick films, surface topography, thin film's magnetic properties, magnetic microstructures

Thickness (d) and surface roughness (σ) of material thin films are inseparable characteristics engendered by the same deposition process. Both are crucial parameters for magnetic properties (p) of real mesoscopic scale thick ferromagnetic films. When condition (d/σ) <10 2 is validated, sorting out samples of a series is possible only by means of τ = (d/σ) geometrical ratio which enables their coherent comparison. Their magnetic properties (p) are thereby expressed in a single function form p=f(τ). Application of that approach to a series of nanostructured Ni electrodeposits reveals that these samples are consistent with the ratio range 1.23 ≤ τ ≤ 82.00 concurring with the reported description. Their coercivity (Hc) and magnetic domain size (w) behave in agreement with their respective predicted general curve profile when d is substituted by τ, which depicts quite well the similar role of both parameters. The study of these properties evolution using the normalization model indicates a discontinuity in the magnetism of the investigated samples. Bloch magnetic domains (MD)B are associated with mixed domain walls (DWN+DWB) below a critical position (τ0 -1 ) ≈ 0.35, while Néel domain walls (DWN) coexist with mixed magnetic domains (MDB+MDN) beyond that position.

INTRODUCTION

The behaviour of ferromagnetic material thin films submitted under a magnetic field is basically assigned to their chemical composition and structure phase. However, the studies on their actual properties always represent a complicated task either in theory or experiment.

Decreasing the sample size, as currently executed nowadays, enhances the difficulty of entirely matching the results obtained from both different approaches. Once surface adsorption phenomenon and probable internal impurities are controlled, the film global geometry can play a determining role on its behaviour. Thickness and surface topography together represent quite well that geometrical feature. Thickness is the most commonly investigated film characteristic which mainly reflects its core or bulk content configuration. In the frame of nano science and technology, thinner specimens can produce better feed-back usually explored in the fabrication of sensors and many other modern device systems [1][2][3].

Surface topography depicts the material interruption zone, bearing for that reason a typical microstructure and morphology which engenders some modifications under the magnetic field effect compared to bulk. Its impact on magnetic film's properties is nowadays widely evidenced [4][5][6]. By and large, the study of individual role of surface or bulk part on a film physical behaviour is undertaken by approximation, neglecting the little contribution of one of them [8][9][10]. This clearly indicates that their complete separation is impossible, regardless of experimental methods used for the samples deposition and characterization. We have to admit that such an approximation can never be systematic due to various situations linked to elaboration process of real material thin films. In actual fact, usual approximating procedures are very restrictive since one excludes the case of samples having equal impact of both film parts and also the situation when these procedures are impossible. In practice, the easier assumption is sometimes quickly adopted, stating for instance that surface topography aspects have a limited effect on the film behaviour. This explains the little attention paid so far to the mentioned unusual situations. The extensive solicitation of mesoscopic scale thick samples including nano-films and nano-structured specimens in the building of modern device systems prompts to carefully examine all other possible options keeping in mind that surface morphology and bulk interconnect closely as both of them result from the same film growth process. Consequently, only a clear and rational approach to their mutual relationship with film properties can emerge onto consistent interpretations. Theoretical studies on these properties evolution are tremendously improved when they remain as close as possible to a detailed description of the film behaviour, even though that description derive from simple models. As a matter of fact, simplifications or approximations are always associated with strict conditions to be fulfilled in any theory application. Accordingly, reduction of simplifications in analytical treatments makes the theory potentially in better agreement with the film material reality. The calculated property is thereby more adapted to the actual behaviour of the investigated sample.

In the present proposal, the study is focused on the evolution of nano-structured film's magnetic properties, insisting on the particular case for which approximation or simplification in the relationship between bulk and surface contributions is impossible. This is a "non trivial" situation practically ignored in the literature of material thin films. An original approach to that evolution is proposed. Sections 2 and 3 are devoted to a concrete study of Ni films formed using a specific electro-deposition protocol. Application of that approach to these film's magnetic properties is reported, followed by identification of the sample's magnetic microstructures evolved that contribute to their ferromagnetism manifestation.

EXPERIMENTAL DETAILS

The Ni samples investigated are prepared by cathodic voltammetry (C-V) technique described previously [7] using a three-electrode cell system consisted of 1 M Ni(SO4)2 aqueous electrolyte (pH = 3.5) obtained from a FLUKA product. The working electrode or stationary deposition substrate as well as the counter electrode are both made of a 0.5 cm 2 area platinum plates (Goodfellow SARL Company) while a METTLER TOLEDO saturated Calomel one serves as reference for potential measurements. The electrochemical system is completed by a radiometer PGP 201 potentiostat. The film deposition method is based on the drawing of cathodic hemi-cycle current-voltage (I-V) curves within a fixed potential interval [0.00; -1.00] V/SCE, the main variable of the process being the scan rate (r) taken in the range 0.17 ≤ r ≤ 1.67 mV.s -1 usually considered as low values. A regular change of r in that selected range leads to the formation of samples having different thickness values measured from a SRATAGEM (SAMx, France) software associated with X-ray microanalysis, as reported elsewhere [11,12]. The values obtained are corroborated by those provided by the sample cross-section images resulting from scanning electron microscopy (SEM). Global film's magnetic properties are investigated by means of a LAKE SHORE COMPANY vibrating sample magnetometer (VSM). A Nanoscope IIa (VEECO) scanning probe microscope (SPM) set operating in dual-modes serves for examination of local film's characteristics. It is then tuned as atomic force microscope (AFM) for the film topography analyses and as a magnetic force microscope (MFM) for capturing the sample's magnetic phase images from which their magnetic microstructures are measured. The general principle of the SPM experiments in both modes is based on the determining role of the tip-sample surface interaction. Magnetic CoCrcoated tips of approximately 17 nm apex radius are connected to a cantilever of 6 Nm -1 spring coeficient. The system is used with a resonance frequency of about 70 Hz. SPM analyses are performed in the tapping/lift mode which allows the collection of both the topography and magnetic phase images of the same surface area. Further details on the use of each one of the applied modes are provided elsewhere [13].

RESULTS AND DISCUSSION

The impact of deposition substrate on material thin film's properties is well established by now even though it still remains a complex but open research topic due to interference of many other parameters [14][15][16][17][18][19]. We showed elsewhere that gold, Cu and ITO substrates respectively engender typical mechanical and microstructure features to our Ni electrodeposits whilst their fcc unit lattice remains unchanged [13,20,21]. Their measured strain and surface morphology including the grain's shape and size are specific in each case.

One knows how these characteristics affect the film's magnetic properties. Pt substrate used in the present work keeps the same Ni film unit lattice and is then expected to bring its specific impact on the sample's microstructure and magnetic features. In the present section, the study is focused on the particular relationship between thickness (d) and root mean square surface roughness (σ) of these Ni samples [22][23][24]. Only the saturated roughness value precisely defined elsewhere is retained here since it is more representative of the entire sample surface irregularities [25,26]. Implication of that relationship on the films coercivity (Hc) and magnetic domain size (w) is examined. These two film magnetic properties are chosen due to the simplicity of their well known analytical dependence on d and also because they are involved here in investigation on the Ni samples' magnetic microstructures.

However, prior to the results presentation, an overview of our original approach to the film's characteristics seems prerequisite.

Characteristics of real mesoscopic scale thick films in a series

Evolution of a thin film property is commonly investigated from a sample's series made of specimens obtained following a precise protocol. Bringing some change to their deposition conditions is the only way to produce different specimens neatly sorted out by means of a precise discrimination criterion. This objective is reached by acting on the film deposition parameters (x1, x2 ...... xn) including deposition time, t. A change in one of the selected xi ≠ t, keeping fixed all the other parameters, directly affects the growth rate of the film formation.

That gives rise to a series of specimens having their individual (di, σi) couples. A similar result is also obtained by solely changing t value although the origin of modification in (di, σi) couples differs since the film growth rate is maintained here. In both cases, the selected xi simultaneously affects d and σ of the specimens. Accordingly, d1 ≠ d2 ≠ d3 ... ≠ dn and also σ1 ≠ σ2 ≠ σ3 ... ≠ σn. It can be deduced that neither d nor σ is separately controlled by any xi parameter. Besides, if evolution of d can be theoretically predicted like in electrodeposition (i.e. Faraday law), anticipating the one of σ remains very difficult. One concludes that sorting out specimens of a series only by means of their thickness d value is impossible since their surface roughness σ differs. Similarly, taking exclusively roughness criterion σ for the same operation is equally unacceptable due to the variation of their d value. When the condition σ <<< d (with d/σ ≥10 2 ) widely admitted for macroscopic scale thick samples is fulfilled, d is their only one main contributing geometrical characteristic. The sorting out operation is then executed referring to d, as commonly encountered. Under the effect of a constant physical field (E), the dependence of the resulting film property (p) on d is then nearly similar to the one of ideally flat films (d ≠ 0, σ = 0) such as

p ≈ f (d) (1).
In principle, the condition (d/σ ≥10 2 ) as proposed is always expected to be possibly fulfilled, regardless of the film thickness category. In actual fact, only specific growth modes can lead to its achievement with concrete microscopic scale thick samples [27,28]. This achievement becomes very problematic for real mesoscopic scale thick samples, including nanoscopic scale thick ones and nano-structured specimens that exhibit a lower chance to directly satisfy that condition. The literature on the last sample's category seldom discusses on experimental results concurring with the condition d/σ <10 2 , which depicts impossible approximation between σ and d [29][30][31][32][33]. In that precise case, the extremely reduced d measurements combined with unavoidable σ values requires considering the sample's global geometrical feature with both characteristics in term of (d ≠ 0, σ ≠ 0) couple, as mentioned above.

Accordingly, each one of them has to be simultaneously considered for a consistent sample's sorting out operation. This particular situation is practically denied or eluded in the studies of material thin films. A solution can be found for that operation by using the sample non dimensional thickness to surface roughness ratio, simply named by sample geometrical ratio τi, expressed as

τi = (di /σi) = (d /σ )i (2).
In that precise case, τi appears a relevant and well adapted characteristic of each real material thin film. In the series, this ratio enables a safe sample's sorting out operation, keeping or not the initial ordering. Ultimately, τ can serve as the only one coherent variable of p corresponding to the following single function

p = f (τ) (3).
A comparison of individual sample pi values can be safely done here, which permits a study of p evolution in an interval of τ values. As long as σ can never be higher than d, our approach unambiguously demonstrates that only two options correlating d andσ can lead to a sorting out operation in a series of real thin films. The first one, mostly adapted to mesoscopic scale thick samples, uses τ as demonstrated here; the second option, more common to macroscopic and microscopic scale thick samples, uses d. It can thus be deduced that d is an extension of τ in that operation.

By and large, sophisticate developments are praised for interpreting the evolution of thin film properties. We here show that the basic role of film geometrical feature is wrongly neglected most of the time due to some quick approximations. As reported, the impact of τ may markedly affect a sample p value whenever condition d/σ <10 2 is validated, which is not concretely as rare as it seems to be with real mesoscopic scale thick films [29][30][31][32][33]. This impact reflects the material feed-back under a physical field (E) effect and then it directly contributes to the film experimental data. Consequently, the replacement of d by τ in analytical expressions of p, when required, should lead to a better closeness between calculated and measured p values.

Topography and thickness of investigated Ni samples

The visible impact of Pt substrate on Ni film's microstructure feature is provided by the surface topography images of Fig. 1 corresponding to the samples having the two extreme d and τ values, as ordered in Table 1. In both cases, the result obtained quite differs from the needle-like grains or local conglomerates respectively engendered by Cu and ITO substrates shown elsewhere [13]. The grains repartition is here more regular, their size being almost uniform and their global shape almost preserved all over the entire film surface. The reported sample's topography images allow a rapid examination of the grain size evolution. Indeed, one notes in Fig. 1 The results reported in

Study of global and local magnetic characteristics of Ni samples

The significance of the ratio interval 1.23 ≤ τ ≤ 82.00 reported above confirms the fact that the properties of the investigated Ni samples are more correctly studied from the approach expressed by equation (3). In that case, τ affects all their global magnetic properties including the saturation magnetization (Ms), remanence (Mr) and coercivity (Hc). The role of τ is clearly depicted in Fig. 3 from the effect of a parallel magnetization M (H║) of the two samples A confirmation of change is provided in Fig. 4 by the study of film's coercivity and magnetic squareness S = (Mr/Ms). The impact of the film geometrical feature induced by Pt substrate on the behaviour of each magnetic property is here evidenced. One sees in Fig. 4(a) that all the obtained S values are superior to 0.50. This experimental S = f(τ) curve profile can be directly and exclusively compared only with those of similar mesoscopic scale thick samples for which τ characteristic is also established. However, as long as d is considered as an extension of τ for the same material thin film taken from macroscopic till mesoscopic scale thick forms, the evolution of S versus τ or d has the same physical interpretation. One knows that higher squareness values (S > 0.50) of ferromagnetic films indicate that the sample's magnetic reversal is ruled by DW nucleation and motion mechanism [34,35]. Consequently, the S = f(τ) curve profile of Fig. 4(a) has the same meaning. The dependence of film's coercivity Hc on τ in Fig. 4(b) depicts a regular decreasing profile which experimentally corresponds to the exponent law Hc = 85.11 (τ) -0.12 (4).

As it appears, that curve profile concurs quite well with Néel prediction law Hc = A d -n (with A and n, the fitting constants) assigned to ferromagnetic films for which coercivity proceeds from DW motion mechanism [36], as found here. Owing to the close link established above between τ and d, one deduces that experimental Hc = f(τ) curve of Eq. 4 describes the same physical property behaviour .

magnetic domain size (w) using the measurement method proposed by Hsieh et al [39] and also Dumas-Bouchiat et al [37]. The corresponding curve proposed in Fig. 6 clearly shows a dependence obeying the law

w = 147 f (τ 1/2 ) (5)
with K=147, whose profile is quite similar to the one expected from the mentioned Kittel w = f(d) prediction. Indeed, the straight line feature obtained for the large beginning part of this curve depicts a proportionality previously reported in the study of the w = f (d 1/2 ) curve profile of Ni and Co samples [36]. One deduces from the present result that the use of d or τ keeps unchanged the linear dependence of w magnetic property, which confirms the common role played by each one of these two film geometrical characteristics. The deviation from the straight line zone reported at (τ 1/2 ) ≈ 4.60 has been also observed with w = f (d 1/2 ) curves and assigned to some changes in the sample's magnetic patterns [34]. The same phenomenon is then likely to be the origin of the mentioned modification.

Magnetic microstructures evolved in Ni nanostructured films

The behavior of ferromagnetic sample's global and local magnetic characteristics is determined by their magnetic microstructure's types. As observed above from the studies of Ni samples surface morphology and magnetic properties, Pt deposition substrate is expected also to impact the Ni sample's magnetic microstructures. By and large, their types are examined from the film magnetic properties for which thickness has a prominent role. The influence of roughness in that matter got a late attention while the combined contribution of thickness and roughness, as proposed in the present work, still remains poorly investigated.

Regarding our particular mesoscopic scale thick Ni sample's category, the study requires considering the dependence of their magnetic properties on τ, as shown above. This study is here deduced from the well adapted analytical treatment devoted to relative or normalized characteristics of ferromagnetic samples detailed elsewhere [11]. This treatment is based on the behavior of the film's normalized magnetic properties, referring to unit film thickness versus normalized film roughness (σ'= σ/d), which corresponds to reciprocal non dimensional geometrical ratio (τ -1 ) of the present approach. Magnetic microstructure types are then directly correlated to evolution of normalized coercivity (Hc') and magnetic domain size (w').

In that case, the re-written film normalized magnetic domain size (w'=w/d) corresponds to = / (6) where EDW represents the domain wall (DW) surface energy; N, the demagnetizing factor; K1, the in-plane anisotropy constant and Ms, the saturation magnetization. Due to DW nucleation and motion mechanism of the investigated Ni samples, their normalized coercivity (Hc' = Hc/d) is expressed as

′ = + + ! " # $ % & (7)
with δ' the normalized DW thickness, Kv and A being respectively the in-plane volume anisotropy and the exchange constants. The reported experimental result of figure 7 shows that the samples w'= f (τ -1 ) curve profile exhibits two different regimes depending on (τ -1 ) region. A regular increase of w' is observed below (τ0 -1 ) ≈ 0.35, which is followed by a stabilization onset located beyond that value. The study of the samples Hc' versus (τ -1 ) proposed in figure 8 depicts an opposite curve profile since a moderate increasing regime of Hc' value precedes an abrupt enhancement of that property occurred also at (τ0 -1 ) ≈ 0.35 or τ0 ≈ 2.86. One deduces from both resulting curve profiles that different magnetic microstructure types evolved both side hands of (τ -1 )0 value.

In the frame of the analytical treatment considered, evolution of (w') or (Hc') with respect to (τ -1 ) provides a signature of the active magnetic microstructure type. The abrupt increasing trend of w' obtained below (τ0 -1 ) in Fig. 7 is interpreted as the contribution of Bloch magnetic domain (MD)B. The increasing trend of Hc' is rather assigned to the Néel domain walls (DW)N activity. The moderate increasing regime of this property at the same (τ0 -1 ) region reported in Fig. 8, can thus be consistently explained by the simultaneous contributions of (DW)N and Bloch domain walls (DW)B, the last type being gradually replaced by the first one. According to the same treatment, the marked increase of w' versus (τ -1 ) is a signature of active (MD)B while its decrease is associated with Néel magnetic domain (MD)N. Consequently, the moderate increase of w'= f (τ -1 ) curve profile towards a plateau as obtained beyond (τ -1 )0 in

SUMMARIZE AND CONCLUSION

Due to the unavoidable interconnection of thickness d and roughness σ resulting from thin film growth process, the study of nano-structured samples' magnetic properties is problematic when condition (d/σ)<10 2 is fully validated. Defining a new film non dimensional geometrical ratio τ = (d/σ) is prerequisite in that case in order to sort out the specimens of a sample's series. That ratio safely serves as a consistent variable for the study of these properties which can be expressed in a general single function form p = f (τ). That new approach is experimented here for studying coercivity (Hc) and magnetic domain size (w) of Ni electrodeposits for which the condition (d/σ) <10 2 is clearly fulfilled. When τ substitutes d, the results obtained with both properties are in agreement with theoretical predictions, implying that τ is equivalent to d and then well-adapted to the real geometry of investigated samples.

  (b) that image magnification is performed for the smoothest sample (σ = 10 nm) in order to enhance the visibility of the sample surface morphology. This operation appears here necessary to satisfactorily compare it's microstructure with the one of the most disturbed surface (σ = 85 nm) of Fig.1(a). A priori, this fact suggests that the evolution of grain size or surface irregularities is significantly quicker with Pt substrate.

Figure 1 .

 1 Figure 1. -3D-AFM surface topography images of the two Ni samples having the extreme characteristics as shown in Table 1): a) -Sample 1 (with d = 102 nm and σ = 85 nm); b) -Sample 6 (with d = 820 nm and σ = 10 nm).

Figure 2 .

 2 Figure 2. -Dependence of sample surface roughness (σ) and geometrical ratio (τ) on the investigated Ni film thickness value (d).

numbered 1 (Figure 3 .

 13 Figure 3. -Normalized magnetization (M/Ms) versus applied parallel magnetic field (H║): hysteresis loops of the two Ni samples having the extreme geometrical ratio values (τ = 1.23 and 82.00).

Figure 4 . 5 .Fig. 1 .

 451 Figure 4. -Dependence of the global magnetic properties of Ni samples on their geometrical ratio: a) -Evolution of magnetic squareness (S=Mr/Ms); b) -Evolution of film's coercivity (Hc).

Figure 5 .

 5 Figure 5. -Local SPM analyses of Ni samples related to the two extreme geometrical ratio values: (I) -(a) topography image and profile; (b) magnetic phase image and profile for the sample labeled τ = 1.23. (II) -(c) topography image and profile; (d) magnetic phase image and profile for the sample labeled τ = 82.00.

Figure 6 .

 6 Figure 6. -Dependence of the investigated Ni samples magnetic domain size (w) on their geometrical ratio (τ).

(

  Film geometrical ratio) 1/2 , √(τ)

Figure 7 .

 7 Figure 7. -Dependence of the normalized magnetic domain size (w') of the investigated Ni samples on their reciprocal geometrical ratio (τ -1 ).

Fig. 7

 7 Fig. 7 can be assigned to the presence of both magnetic domain types. The reported stabilization onset of w' in the same figure hence corresponds to an attenuation of (MD)B activity which is gradually replaced by an enhanced role of (MD)N type. The marked increase

Figure 8 .

 8 Figure 8. -Dependence of the normalized coercivity (Hc') of the investigated Ni samples on their reciprocal geometrical ratio (τ -1 ).

  Updated normalized treatment applied to these Ni films shows a discontinuity in the magnetism of these samples at τ0 ≈ 2.86. Below (τ0 -1 ), Bloch magnetic domains (MD)B coexist with mixed domain walls (DWN + DWB) while Néel domain walls (DW)N are associated with mixed magnetic domains (MDB + MDN) beyond that critical value. This preliminary contribution unambiguously reveals that τ concurs quite well with the study of 0 0.09 0.24 0.38 0.61 0.93 Norm. coercivity, H c ' (Oe/nm) Reciprocal film geometrical ratio, (τ -1 )

Table 1 -

 1 Table1have some common features with those observed earlier with the other mentioned substrates, which are mainly assigned to the film deposition technique involved. Although the thickness limits always differ with the substrate used, the values obtained here are located in the interval 102 ≤ d ≤ 820 nm, which is typical of mesoscopic scale thick samples (d <1 μm). Their measured surface roughness values are ranged in the region 10 ≤ σ ≤ 85 nm. As depicted in Fig.2, thinner samples have the higher roughness values and vice versa in agreement with previous experimental observations[13]. We already explained that evolution by the determining role of scan rate in the dynamics of the film deposition process, which confirms the relevancy of (C-V) technique used for elaboration of smoother thick samples. The behavior of sample geometrical ratio τ versus d reported in the same figure provides a new feature of the film characteristic's modification. Although τ Thickness, roughness and ordering of the investigated Ni samples sorted out by their geometrical ratio characteristic (τ).

	d (nm)	102	127	160	320	530	820
	σ (nm)	85	72	56	43	25	10
	τ = (d/σ)	1.23	1.76	2.85	7.44	21.22	82.00
	Samples
	order	1	2	3	4	5	6

regularly follows the increase of d, its evolution exhibits two different regimes in the investigated d values. The first one corresponds to a moderate increase of that ratio below d ≈ 530 nm while the second step depicts a tremendous enhancement of its value beyond that sample thickness. As shown in Table

1

, the thinnest sample leads to τ = 1.23, which evidences the fact that σ cannot be neglected compared to d. The ratio τ = 82.00 of the thicker sample implies that approximation σ <<< d is not yet acceptable although that condition is about to be fulfilled, knowing that τ = 100 represents our accepted minimal limit of negligibility. One deduces from the obtained ratio evolution that the individual role of d and σ is never identical for all the Ni samples of the series. The range of the reported τ values implies that σ can never be neglected compared to d for all our investigated samples. This result clearly justifies the use of our approach of sub-section 3.1 and leads to neatly sort out the Ni samples of the series as proposed in Table

1

referring to τ, their only one common characteristic.

ferromagnetic samples for which approximation between d and σ is impossible. An appropriate use of τ in theoretical analyses is expected leading to more accurate results in better closeness with experimental data. In the light of the present approach, examination of different ferromagnetic material thin films remains the main concern of next investigation steps aims at confirming the real practical interest of the present proposal.


domain: spi chim phys

text:
Ihar Babushkin

Álvaro Jiménez Galán

José Ricardo Cardoso De Andrade

Anton Husakou

Felipe Morales

Martin Kretschmar

Tamas Nagy

Virgilijus Vaičaitis

Liping Shi

David Zuber

A J Galan

J R C Andrade

L Bergé

S Skupin

I A Nikolaeva

N A Panov

D E Shipilo

O G Kosareva

A N Pfeiffer

A Demircan

M J J Vrakking

U Morgner

M Ivanov

All-optical attoclock for imaging tunnelling wavepackets

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Recent experiments on measuring time-delays during tunnelling of cold atoms through an optically created potential barrier are reinvigorating the controversial debate regarding possible time-delays during light-induced tunnelling of an electron from an atom. Compelling theoretical and experimental arguments have been put forward to advocate opposite views, confirming or refuting the existence of finite tunnelling time delays. Yet, such a delay, whether present or not, is but a single quantity characterizing the tunnelling wavepacket; the underlying dynamics are richer.

Here we propose to augment photo-electron detection in laser-induced tunnelling with detection of light emitted by the tunnelling electron -the so-called Brunel radiation. Using a combination of single-color and two-color driving fields, we identify the all-optical signatures of the re-shaping of the tunnelling wavepacket as it emerges from the tunnelling barrier and moves away from the core. This reshaping includes not only an effective time-delay but also time-reversal asymmetry of the ionization process, which we describe theoretically and observe experimentally. We show how both delay and reshaping are mapped on the polarization properties of the Brunel radiation, with different harmonics behaving as different hands of a clock moving at different speeds. The all-optical detection paves the way to time-resolving optical tunnelling in condensed matter systems, e.g. tunnelling across bandgaps in solids, on the attosecond time-scale.

INTRODUCTION

Tunnelling of an electron through the potential barrier created by an oscillating electric field and the binding potential of the core is a key resource in attosecond science [START_REF] Corkum | Attosecond science[END_REF] and is at the heart of high harmonic generation and high harmonic spectroscopy [START_REF] Baker | Probing proton dynamics in molecules on an attosecond time scale[END_REF][START_REF] Smirnova | High harmonic interferometry of multi-electron dynamics in molecules[END_REF][START_REF] Bruner | Multidimensional high harmonic spectroscopy[END_REF][START_REF] Torres | Revealing molecular structure and dynamics through high-order harmonic generation driven by mid-ir fields[END_REF][START_REF] Pedatzur | Attosecond tunnelling interferometry[END_REF]. Generally, high harmonic generation is associated with radiative recombination following the return of the laser-driven electron to the parent ion. Yet, even when the electron does not return to the core, harmonic radiation is still emitted. Often referred to as the "Brunel radiation" or "Brunel harmonics" [START_REF] Brunel | Harmonic generation due to plasma effects in a gas undergoing multiphoton ionization in the high-intensity limit[END_REF][START_REF] Babushkin | Terahertz and higher-order Brunel harmonics: from tunnel to multiphoton ionization regime in tailored fields[END_REF][START_REF] Balciunas | Optical and THz signatures of sub-cycle tunneling dynamics[END_REF], it is associated with bursts of current triggered by laser-induced tunnelling and is ubiquitous in atoms, molecules, and solids [START_REF] Lanin | Mapping the electron band structure by intraband high-harmonic generation in solids[END_REF][START_REF] Silva | Topological strong-field physics on sub-lasercycle timescale[END_REF][START_REF] Jürgens | Origin of strong-field-induced low-order harmonic generation in amorphous quartz[END_REF][START_REF] Jürgens | Characterization of laser-induced ionization dynamics in solid dielectrics[END_REF]. Sub-cycle tunnelling bursts occur every time the oscillating electric field goes through a maximum. Both odd and even harmonics emerge when the incident laser field contains both fundamental radiation and its second harmonic, including the so-called '0-th' order Brunel harmonic corresponding to terahertz (THz) emission, see e.g. [START_REF] Babushkin | Terahertz and higher-order Brunel harmonics: from tunnel to multiphoton ionization regime in tailored fields[END_REF][START_REF] Kim | Coherent control of terahertz supercontinuum generation in ultrafast laser-gas interactions[END_REF][START_REF] Babushkin | Tailoring terahertz radiation by controlling tunnel photoionization events in gases[END_REF]. To produce THz radiation, the waveshape must be at least asymmetric in respect to the change E → -E. The origin of this radiation can be traced already in a classical picture: According to Maxwell's equations, a change of the electric current dJ/dt produces radiation. A small volume emits, in a small time interval dt, Brunel radiation δE Br ∝ dJ produced by the change of free electron current dJ ∝ ρ(t)E(t)dt due to acceleration of electrons with density ρ(t) in the driving field E(t) in the situation when the ρ(t) is time-dependent [START_REF] Jürgens | Characterization of laser-induced ionization dynamics in solid dielectrics[END_REF][START_REF] Geissler | Light propagation in field-ionizing media: Extreme nonlinear optics[END_REF]. Additional terms may arise due to, for instance, nonzero birth velocity of electrons (see Methods), but the above expression is rather good approximation for the lowest harmonics [START_REF] Babushkin | Generation of terahertz radiation from ionizing two-color laser pulses in ar filled metallic hollow waveguides[END_REF][START_REF] Zhang | Continuum electron giving birth to terahertz emission[END_REF].

Here we show that different Brunel harmonics generated in elliptically polarized single-color and two-color laser fields provide a detailed picture of light-induced tunnelling. Our approach to imaging the ionization dynamics is distinctly different from the currently used attoclock approaches based on detecting photo-electrons [START_REF] Eckle | Attosecond angular streaking[END_REF][START_REF] Eckle | Attosecond ionization and tunneling delay time measurements in helium[END_REF][START_REF] Huismans | Timeresolved holography with photoelectrons[END_REF][START_REF] Shafir | Resolving the time when an electron exits a tunnelling barrier[END_REF][START_REF] Kienberger | Atomic transient recorder[END_REF][START_REF] Pfeiffer | Attoclock reveals natural coordinates of the laser-induced tunnelling current flow in atoms[END_REF][START_REF] Landsman | Ultrafast resolution of tunneling delay time[END_REF][START_REF] Yakaboylu | Wigner time delay for tunneling ionization via the electron propagator[END_REF][START_REF] Torlina | Interpreting attoclock measurements of tunnelling times[END_REF][START_REF] Ni | Tunneling ionization time resolved by backpropagation[END_REF][START_REF] Camus | Experimental evidence for quantum tunneling time[END_REF][START_REF] Sainadh | Attosecond angular streaking and tunnelling time in atomic hydrogen[END_REF][START_REF] Eicke | Attosecond-scale streaking methods for strong-field ionization by tailored fields[END_REF]. It allows us to introduce a complementary, alloptical measurement protocol, which should enable one to extend measurements of tunnelling dynamics into the bulk of solids, where photo-electron spectroscopy is not readily available.

In the strong field regime, the conventional photoelectron attoclock [START_REF] Eckle | Attosecond angular streaking[END_REF][START_REF] Eckle | Attosecond ionization and tunneling delay time measurements in helium[END_REF][START_REF] Ivanov | Strong-field ionization of he by elliptically polarized light in attoclock configuration[END_REF] measures the deflection in the photo-electron spectrum generated by a few-cycle, nearly circularly polarized pulse, with respect to the pulse polarization ellipse. Interpretation of these observations uses mapping between the instant at which the electron "escapes the tunnelling barrier" and its final velocity at the detector (see Fig. 1a). This interpretation, and the possible time-delays associated with the ionization, are a subject of active discussions [START_REF] Landsman | Ultrafast resolution of tunneling delay time[END_REF][START_REF] Yakaboylu | Wigner time delay for tunneling ionization via the electron propagator[END_REF][START_REF] Torlina | Interpreting attoclock measurements of tunnelling times[END_REF][START_REF] Ni | Tunneling ionization time resolved by backpropagation[END_REF][START_REF] Camus | Experimental evidence for quantum tunneling time[END_REF][START_REF] Landsman | Attosecond science and the tunnelling time problem[END_REF][START_REF] Teeny | Ionization time and exit momentum in strong-field tunnel ionization[END_REF][START_REF] Kaushal | Spin-orbit larmor clock for ionization times in one-photon and strong-field regimes[END_REF][START_REF] Eicke | Attoclock with counter-rotating bicircular laser fields[END_REF], just as the general problem of tunnelling times, see e.g. the latest beautiful experiment on this subject [START_REF] Ramos | Measurement of the time spent by a tunnelling atom within the barrier region[END_REF]. The difficulty stems from the need to deconvolve ionization and post-ionization dynamics in the presence of the electron-core interaction, which deflects the outgoing electron and thus distorts the mapping between the tunnelling event and the measurement at a faraway detector. Moreover, possible under-the-barrier delays and the post-barrier electron velocity make the one-to-one mapping from the state at detector and emission time impossible [START_REF] Camus | Experimental evidence for quantum tunneling time[END_REF][START_REF] Teeny | Ionization time and exit momentum in strong-field tunnel ionization[END_REF]. Optical techniques, relying on measuring recombination-based high harmonic emission [START_REF] Bruner | Multidimensional high harmonic spectroscopy[END_REF][START_REF] Pedatzur | Attosecond tunnelling interferometry[END_REF] instead of electron detection, face an additional challenge: the need to decouple the ionization and photorecombination steps.

Therefore, it is advantageous to place the "observation point" as close to the tunnel exit as possible. Brunel radiation is an ideal candidate: it is produced by the current generated by the tunneling electron right near the barrier exit. We find that already the lowest, 0thorder Brunel harmonic, encodes the attosecond-scale delays and the wavepacket deflection by the core potential. The information is contained in the rotation of the polarization ellipse of this harmonic. We prove this analytically and numerically by comparing the proposed new protocol against the standard attoclock setup.

Moreover, we find theoretically and confirm experimentally that higher-order Brunel harmonics encode the temporal profile of the ionization burst, providing several connected "clock hands" rotating at different speed. Using the third-order Brunel harmonic generated in helium, we are able to reconstruct the dynamical reshaping of the electronic wavepacket during ionization. Interestingly, we find asymmetry in the ionization dynamics relative to the instantaneous peak of the driving electric field. We attribute this effect to pre-and post-ionization electron dynamics.

PHYSICAL PRINCIPLE AND THEORETICAL ANALYSIS

We first validate the key idea behind the all-optical attoclock: the vectorial properties of the emitted light, such as the rotation of the polarization ellipse, are determined by the vectorial properties of the current generated by the tunnelling electron and therefore reflect the tunnelling dynamics.

We consider two field arrangements. In the first one, intense circularly polarized infrared (IR) pump with frequency ω 0 is combined with its co-rotating second harmonic 2ω 0 , see Fig. 1b. This generates a total electric field with one well-pronounced maximum per cycle. The orientation of this maximum, here set along the x-axis, is controlled by the relative phase between the two colors and provides the reference direction for the optical attoclock. In the second arrangement, such reference direction is provided by the major axis of the single-color elliptically polarized driving field.

All-optical measurement of strong-field ionization delays

We begin with the first arrangement. The nonlinear response contains even and odd harmonics nω 0 , including the "zeroth-order" harmonic, where the signal is dominated by the Brunel radiation [START_REF] Meng | Enhancement of terahertz radiation by using circularly polarized two-color laser fields[END_REF][START_REF] Tulsky | Boosting terahertz-radiation power with twocolor circularly polarized midinfrared laser pulses[END_REF][START_REF] Tailliez | Terahertz pulse generation by two-color laser fields with circular polarization[END_REF]. Let a classical free electron be injected by strong-field ionization into the atomic continuum with some (not necessarily zero) velocity and be accelerated in the laser field and the potential of the core. The electric field associated with the Brunel radiation, emitted by this electron, reads E Br (ω) ∝ dp(t) dt e iωt dt, with p being the electron momentum. For the 0-th harmonic ω → 0, and E Br (0) ∝ p(t → ∞). Thus, the polarization direction of the 0-th Brunel harmonic tracks the final velocity of the electron after the interaction, i.e. performs the same measurement as the photo-electron attoclock, but by alloptical means. Note that this conclusion is independent of assumptions regarding the ionization process. While exactly zero frequency cannot be observed in the experiment, measurements in the low-frequency (THz) range are feasible and discussed below.

Our conclusion also holds within the fully quantummechanical analysis, see Methods. In this case, E Br (0) tracks the center of mass of the electronic wavepacket. Specifically,

φ[E Br (0)] = φ τ = φP (φ)dφ, (1)
where φ[E Br (0)] denotes the azimuthal angle formed by the vector E Br (0) in the xy-plane, orthogonal to the light propagation direction z, measured clockwise from the positive y direction; P (φ) = |ψ(t → FIG. 1. Photoelectron (a) and optical (b-d) attoclock. In (a), the angle-resolved photo-electron spectrum generated by the driving field (red line) is registered at a detector, revealing attosecond delays and deflections of electronic wavepacket shifting its maximum by φτ , interpreted as an effective delay τ = φτ /ω0. In (b), this delay is probed, instead of electrons, by the 0th-order Brunel radiation in a two-color field (magenta line). Such harmonic generated in gas or solid (blue box) can be selected by a low-pass filter (brown box). It is nearly linearly polarized (brown line) and rotated, having the same effective τ as in (a). The inset shows schematically the harmonic response spectrum, including the pump (blue) and harmonic of interest (red). (c) Polarization state of the 0th harmonic for the equivalent intensity 0c|Emax| 2 /2 = 150 TW/cm 2 , obtained from TDSE simulations with Yukawa potential, Coulomb potential, and using the simple-man classical Drude model. (d) Effective delay τ as a function of intensity obtained from the polarization rotation of the 0th harmonic as determined by TDSE simulations and compared to the results of the photo-electron attoclock: the center of mass position of the electronic wavepacket (TDSE-CM), as well as the simulations of [START_REF] Torlina | Interpreting attoclock measurements of tunnelling times[END_REF] (where the plane z = 0 instead of center of masses was analyzed) and the experiment of [START_REF] Sainadh | Attosecond angular streaking and tunnelling time in atomic hydrogen[END_REF].

∞, p, φ, θ)| 2 p 2 sin θdpdθ is the photo-electron distribution integrated over the momentum and the polar angle, and ψ(t, p, φ, θ) is the electron wavefunction in prepresentation. In the photo-electron attoclock, the angle φ τ is mapped onto the effective ionization time-delay τ as φ τ = ω 0 τ , where ω 0 is the carrier defining the "standard attoclock mapping" M, M(φ τ ) = φ τ /ω 0 = τ . For our two-color field configuration, Fig. 1b, the mapping is slightly different (see Methods): τ = M(φ τ ) φ τ /ω 0 , where ω 0 = 4/3ω 0 . This mapping should not be confused with the physical interpretation of τ , which requires the analysis of the deflection of the tunnelling electron in the core potential, possible energy gain during optical tunnelling and possibly non-zero initial electron velocity when exiting the barrier [START_REF] Landsman | Ultrafast resolution of tunneling delay time[END_REF][START_REF] Yakaboylu | Wigner time delay for tunneling ionization via the electron propagator[END_REF][START_REF] Torlina | Interpreting attoclock measurements of tunnelling times[END_REF][START_REF] Ni | Tunneling ionization time resolved by backpropagation[END_REF][START_REF] Camus | Experimental evidence for quantum tunneling time[END_REF][START_REF] Landsman | Attosecond science and the tunnelling time problem[END_REF][START_REF] Teeny | Ionization time and exit momentum in strong-field tunnel ionization[END_REF][START_REF] Kaushal | Spin-orbit larmor clock for ionization times in one-photon and strong-field regimes[END_REF][START_REF] Eicke | Attoclock with counter-rotating bicircular laser fields[END_REF].

Putting together Eq. ( 1) and the definition of M for the two-color field configuration, the polarization vector of the 0-th Brunel harmonic is predicted to be

E Br (0) ∼ (sin(ω 0 τ ), cos(ω 0 τ )), ( 2
)
where τ is the same effective ionization delay as obtained from the photo-electron spectrum. We now verify this conclusion by using ab-initio timedependent Schrödinger equation (TDSE) simulations to compute the radiated field. We start with the short-range Yukawa potential, which allows us to remove the effects of the long-range Coulomb tail of the core. In this case, the TDSE simulations show that the THz radiation is indeed nearly linearly polarized and that the major axis of its polarization ellipse is orthogonal to the major axis of the Lissajous figure of the driving two-color field, Fig. 1b,c. Thus, we confirm that, for a short-range potential supporting a single bound state, both the photo-electron attoclock and the optical attoclock measure zero ionization delays. This result also coincides with the result of classical simple-man Drude model, where Eq. ( 3) was used with ionization rate ρ(t) calculated using Yudin-Ivanov formula [START_REF] Yudin | Nonadiabatic tunnel ionization: Looking inside a laser cycle[END_REF]. The finite ellipticity seen in the TDSE results is due to the finite spectral width of the response.

For the Hydrogen atom, TDSE simulations show rotation of the polarization ellipse of the 0-th harmonic, Fig. 1c. Figure 1d compares the time-delay extracted from this rotation with the numerical results of [START_REF] Torlina | Interpreting attoclock measurements of tunnelling times[END_REF] and experimental results of [START_REF] Sainadh | Attosecond angular streaking and tunnelling time in atomic hydrogen[END_REF] for the photo-electron attoclock, as well as the positions of the centers of mass of the electronic wavepackets. The data points in Fig. 1d show the average polarization as discussed in Methods. This average is weighted with the spectral energy density, and coincides with the rotation angles at ω → 0 (see Supplementary). Again, we confirm that the two-color all-optical attoclock demonstrates essentially the same effective time delay τ as the single-color photo-electron one, especially if we compare the rotation angles with the center of mass of the electronic wavepacket, as our theory suggests. The broader polarization ellipse for the Coulomb potential is due to the broader spectrum of the THz response and the fact that higher frequencies are more elliptically polarized (see Supplementary).

Extending optical attoclock using higher-order Brunel harmonics

The attoclock mapping φ τ → τ establishes the relation between only two numbers, φ τ and τ . We now extend the optical attoclock beyond the measurement of a single parameter to access the sub-cycle ionization yield. To this end, we consider the second field arrangement with a single-color, elliptically polarized field, see Fig. 2a, and focus on the rotation of the polarization ellipse of the higher-order Brunel harmonics.

Consider again the simple model of a classical electron extracted with zero velocity. The Brunel radiation, given by the electrons born with ionization rate W (t) and producing free current J(t) is described by a simple-man (so-called Drude) model [START_REF] Babushkin | Generation of terahertz radiation from ionizing two-color laser pulses in ar filled metallic hollow waveguides[END_REF][START_REF] Kim | Generation of coherent terahertz radiation in ultrafast laser-gas interactions[END_REF]:

E Br (t) ∝ dJ/dt ∝ E(t)ρ(t) ∝ E(t) t -∞ W (t )dt , (3)
where ρ(t) is the free electron density.

On the other hand, we can also use the photo-electron spectrum and the attoclock mapping between the polar detection angle φ of electrons, measured at infinity, to define an effective ionization time t and ionization rate W e (t),

W e (t) = P (φ)dφ/dt, φ = M -1 (t). ( 4
)
Note that such definitions have nothing to do with the interpretation of the ionization delays in either atto-clock measurement.

As shown in the Methods section using the R-matrix approach [START_REF] Torlina | Interpreting attoclock measurements of tunnelling times[END_REF][START_REF] Torlina | Time-dependent analytical rmatrix approach for strong-field dynamics. i. one-electron systems[END_REF], the simple-man equation Eq. ( 4) remains fully valid also beyond the classical approximation, at least for few lowest order harmonics, if we assume W (t) = W e (t). Thus, as long as the measured harmonics are associated with the Brunel radiation, observing E Br (t) gives all-optical access to P (φ), imaging the tunneling wavepacket by all-optical means. Below we confirm this conclusion by comparing W (t) obtained optically and directly from the photo-electron distribution.

To reconstruct W (t) from the Brunel radiation, we rewrite Eq. (3) in the frequency domain:

E Br,n ∝ i ω 0 j E j W n-j n -j , (5)
where W k and E j are the Fourier components of the ionization rate and of the driving field, respectively, and E Br,n are n-th Brunel harmonics. Already the polarization state of the 3rd harmonic gives rich information about the ionization dynamics. For the two co-rotating circular fields used above, the polarization of the 3rd harmonic is very close to circular, and thus no information can be extracted from its polarization direction since the latter is undefined. For the single-color elliptical pump (with ellipticity = 1), see Fig. 2a, Eq (5) yields

E Br,3 ∝ (1 + r/2, i (1 -r/2)) , r = W 4 /W 2 . (6)
Eq. [START_REF] Pedatzur | Attosecond tunnelling interferometry[END_REF] shows that the polarization state of the 3rd harmonic depends critically on the 2nd and 4th harmonics of the ionization burst W (t). Most importantly, it allows us to determine whether the ionization burst is centered at the field maximum and symmetric around this center, or delayed and/or asymmetric. In the first case, the ratio r = W 4 /W 2 is real-valued and the rotation of the polarization of the third harmonic is zero.

FIG. 2.

Extracting photoionization information from higher order Brunel harmonics for the Hydrogen atom. In (a), a single-color elliptically polarized pump (red line) produces the 3rd Brunel harmonic (magenta line), with the polarization state (ellipticity and polarization direction) encoding the sub-cycle dynamics of the ionization process. Inset shows schematically the harmonic response spectrum, including the pump (blue) and harmonic of interest (red). (b) Polarization of the 3rd harmonic for the single-color driving pulse (|Emax| = 36.3 GV/m, = 0.6, I = 280 TW/cm 2 ) calculated using TDSE with Coulomb and Yukawa potentials as well as simple-man Drude model (denotations as in Fig. 1c). Inset shows the attoclock mapping M and its inverse. (c) Optical reconstruction of the ionization dynamics (Wopt(t), red solid line) compared to the reconstruction from the photoelectron spectrum (We(t), red dot-dashed line) for the Coulomb potential. The estimation of error is given in Methods. Optical reconstruction for the Yukawa potential (WY (t), blue line) is also presented. Black dot-dashed line shows the electron wavepacket distribution P (φ). The attoclock delay ω0τ given by the the position of the center of mass (CM) of We(t) ≈ Wopt(t) is also indicated (red vertical dotted line). (d) The positions of the maxima of the electronic spectra P (φ) (black markers), their CM(blue markers) and the effective delays τ = M(φ) reconstructed from the photo-electron spectra (red markers) for different FWHM pulse durations and ellipticities of the driving field. Horisontal red line shows the attoclock delay extracted optically using the two-color configuration (Fig. 1) for the corresponding peak intensity. Asymmetry of the wavepacket reveals itself from the different positions of the maximum and CM of the photo-electron spectra (see the dotted arrows tracing these for an exemplary case of = 0.5).

In the second case, the ratio r is complex-valued and the polarization is rotated. Thus, extracting the phase of r from the rotation of E Br,3 and comparing it with the ionization delay allows one to identify the presence of asymmetry in the ionization burst W (t) with respect to its maximum. Additional information is, of course, provided by the higher-order harmonics and the phase of r = W 6 /W 4 , r = W 8 /W 6 etc. In the experiment, we focus on E Br,3 .

To evaluate this approach, we first make TDSE simulations for the Hydrogen atom in our single-color field configuration (Fig. 2a) and observe the polarization of the 3rd harmonic, see Fig. 2b. As expected, for the Yukawa potential the rotation of E Br,3 is noticeably smaller than for the Hydrogen atom and is different from the one expected from the time-delay measured by the 0-th harmonic of the two-color optical attoclock. Together, the information extracted from both measurements is already sufficient to reconstruct W (t) defined by Eq. ( 4) up to 4th harmonic,

W opt (t) = W 0 + W 2 cos(2ω 0 t + δ 2 ) + W 4 cos(4ω 0 t + δ 4 ), (7)
where, we changed the notations as 2|W j | → W j , arg(W j ) → δ j with respect to Eq. ( 5). The subscript in W opt (t) stresses the optical route for the reconstruction. The coefficients W 2 , W 4 , δ 2 , δ 4 are reconstructed uniquely if we assume that the center of mass of W opt (t) appears at t = τ (see details in Methods).

The resulting W opt (t) is compared in Fig. 2c with W e (t) extracted directly from the photo-electron spectra P (φ) (red solid and red dot-dashed curves, respectively). We find excellent agreement between both, proving the viability of the all-optical method.

Figure 2c shows a clear asymmetry in W (t) and P (φ) with respect to time reversal, i.e., the field maximum, in the case of the Coulomb potential. This asymmetry disappears completely for the Yukawa potential, proving that the Coulomb tail not only introduces the time-delay, but also distorts the outgoing wavepacket. While in retrospect such conclusion appears transparent, hindsight is always 20/20: to the best of our knowledge, this is the first time when such asymmetry has been both reported and quantified. The presence of both non-trivial ionization phase and the distortion of the outgoing electronic wavepacket is important not only for the Brunel, but also for the recollision-based harmonics in molecules, where it will manifest as channel-dependent ionization phase associated with different ionic states between ionization and recombination, and thus different effective core potentials. Such channel-dependent ionization phase has indeed been confirmed experimentally in [START_REF] Bruner | Multidimensional high harmonic spectroscopy of polyatomic molecules: detecting sub-cycle laser-driven hole dynamics upon ionization in strong mid-ir laser fields[END_REF]. Interestingly, our simulations (see Fig. 2d) also show that the asymmetry becomes small for nearly circularly polarized and few-cycle pulses, suggesting that asymmetry could be also related to intermediate excitations building up inside the potential well, expected in the regime of intermediate Keldysh parameters [START_REF] Ivanov | Anatomy of strong field ionization[END_REF]. This interpretation is further confirmed by the pronounced Freeman resonances in the rotation angle of E Br,3 as a function of laser intensity; Besides, for high input intensities, saturation effects start to play a role (see Supplementary information for the case of He atom).

We note that the ionization delay τ extracted from the photo-electron spectra is nearly independent of the ellipticity and pulse duration and coincides with that extracted optically, using the two-color configuration, see Fig. 2d. We see thus that the ionization delay is, as we can judge, waveform-independent, in contrast to the deflection of the electronic cloud.

IMAGING IONIZATION DYNAMICS OF HELIUM

We now move to the experimental results. The predicted rotation of the polarization ellipse of the nonlinear response is confirmed by the experimental measurements with setup depicted in Fig. 3a. Here we report the rotation angle and the ellipticity of the third harmonic generated by a strong pulse with ellipticity = 0.61 focused in helium. The experimental approach is described in detail in the Methods section. The polarization of the pump was measured and found rotate by a negligible angle of less than 0.5 degrees upon the propagation, thus excluding the influence of plasma-induced birefringence. We note nevertheless, that the Brunel mechanism of 3rd harmonic generation competes with the χ (3) -based one. Since those two mechanisms give contributions with relative phase π/2, this results in additional polarization rotation (see Supplementary Material). The correction induced by the χ (3) nonlinearity is shown as a blue shaded region above the theoretical curves in Fig. 3. Its impact decreases with intensity and becomes almost negligible above I ≈ 1500 TW/cm 2 .

The experimentally measured intensity-dependent parameters of the polarization ellipse (rotation angle and ellipticity) are shown in Fig. 3b, compared with the results of TDSE simulations for the helium atom. The details regarding the error bars are described in the Methods. Experiment and simulation are in a very good agreement. The polarization of the 3rd harmonic is rotated by nearly 5 degrees as the intensity increases from 1000 to 2000 TW/cm 2 , while the ellipticity moderately increases from around 0.56 to 0.62. Combining this information with the theoretical method developed above, we are able to reconstruct the subcycle ionization dynamics, taking the strong-field ionization delay τ of He from the experiments performed in [START_REF] Eckle | Attosecond ionization and tunneling delay time measurements in helium[END_REF].

The coefficients W 2 , W 4 , δ 2 , δ 4 , as defined in Eq. ( 7) and reconstructed from our measurement, are shown in Fig. 3(c,d) as a function of the pump intensity. We find that W 2 and W 4 remain nearly constant as intensity increases, whereas the angles δ 2 and δ 4 experience significant variation. The fact that W i barely change shows that ionization remains "equally fast" for the range of intensities in the experiment. On the other hand, the growth of |δ i | indicates the changing asymmetry of the effective sub-cycle ionization rate W (t) with respect to its maximum. In Fig. 3(e,f), W (t) is visualized for two values of the pump intensity, showing that the asymmetry increases with increasing intensity. While the asymmetry remains fairly small even for large intensities, its presence in the measurements is clearly seen from Fig. 3(c,d): the measured values of |δ i | deviate substantially from zero, well beyond the experimental error bars. Thus, even small asymmetries lead to large deviations of the rotation angle, which is well detected by the all-optical method.

CONCLUSIONS AND OUTLOOK

In conclusion, we have established a firm quantitative link between measuring the photo-electron spectra in strong-field ionization and measuring the Brunel radiation generated by electrons on their way to the continuum, revealing the reshaping of the electronic wavepacket during laser-induced tunnelling. Such reshaping is mapped on the effective ionization delays and the asymmetry of the instantaneous ionization rate with respect to the peak of the field, as it is imaged by the Brunel harmonics. According to this, Brunel harmonics in THz and UV contain signatures of attosecondand sub-angstrom-scale electron dynamics, deeply subwavelength for these frequencies. Remarkably, we observe that the ionization delay is largely independent on the driving field waveshape, in contrast to the spatial reshaping of the electronic wavepacket. We attribute the origin of the ionization asymmetry to the dynamics of the electron wavepacket during and after tunnelling and, for high intensities, saturation effects. Our experiment in helium shows that the all-optical method is very sensitive even to small asymmetries of such kind. This demonstrates promising capabilities of the proposed new protocol for imaging tunnelling, which allows one to explore attosecond-scale wavepacket reshaping in systems where photo-electron detection is not readily available. One of such systems is the bulk solids, where the detection of light is much simpler than the detection of electrons and where attosecond-scale delays were recently reported [START_REF] Hofmann | Noninstantaneous polarization dynamics in dielectric media[END_REF][START_REF] Sommer | Attosecond nonlinear polarization and light-matter energy transfer in solids[END_REF][START_REF] Hassan | Optical attosecond pulses and tracking the nonlinear response of bound electrons[END_REF]. We foresee that Brunel harmonics of yet higher order will allow resolving electron dynamics even closer to the core.

METHODS

Experimental setup

The laser used in the experiments is a Ti:Sapphire-based regenerative amplifier (Spectra Physics Spitfire) delivering pulses with 3 mJ of energy and the spectrum supports a 43 fs of pulse duration at the driving wavelength of 800 nm. The beam was delivered to the experimental setup via an evacuated hollow-core fibre with 70% transmission. This ensured a clean beam profile with an M 2 < 1.1 and a 5.5 mm beam diameter. The delivered energy was selected using a pair of thin-film polarizers in combination with a broadband half-wave plate. To further attenuate the pulses and ensure the polarization purity of the electric field a fused silica wedge was used in Brewster's angle and the s-polarized pulse was used in the experiment. This allowed a tuning of the pulse energy between 40 and 225 µJ with a high linear polarization contrast. An achromatic quarterwave plate covering the whole range of the laser spectrum was used to convert the linear polarization into elliptical. The polarization was analyzed using a Glan-Laser calcite polarizer. The measured ellipse had an intensity ratio of 0.372 between the minor and major axes. The light was then guided into a vacuum chamber containing a parabolic gold mirror with a focal distance of 100 mm. After the focus the beam was recollimated using an Aluminium curved mirror and guided out of the chamber using further Aluminium coated mirrors. A dichroic beam-splitter was used to separate the thirdharmonic from most of the fundamental. The IR light rejected by the dichroic beam splitter was directed onto a power meter in order to analyze power fluctuations during the measurements. The reflected light was used in the polarization characterization portion of the experiment. The polarization state of the third harmonic was analyzed using an α-BBO Glan-Laser polarizer (ThorLabs GLB10-UV) mounted on a home-built electronic rotation stage. Detection was performed after a dielectric HR mirror for 266 nm light allowing further separation of the UV component from the IR. The 800 nm transmitted through the 266 nm HR mirror was detected on a Silicon photodiode while the UV was measured using a Silicon Carbide (SiC) photodiode. The SiC detector is not sensitive to visible or infrared radiation, which allowed a measurement where the background is not influenced by any co-propagating IR light. All three measurement devices were connected to analogical boxcar integrators and the signals were then digitized using a 16-bit ADC connected to an arduino nano. This setup allowed to characterize the polarization state of the UV elliptical light, and to compare the direction of its long axis to that of the exciting IR field. The chamber was evacuated and subsequently filled with He at a pressure of 1.3 bar. No other nonlinearities from the optical components were observed.

Pump fields in time and frequency representation, and the attoclock mapping M

We consider driving pulses propagating in z-direction, with electric field in OXY plane: E = {Ex, Ey} and vector potential A{Ax, Ay}, A(t) = -t -∞ E(t )dt in the form:

Ax(t) = E 0 ω 0 e -2 ln 2 t 2 /T 2 cos ω 0 t, (8)
Ay(t) = - E 0 ω 0 e -2 ln 2 t 2 /T 2 sin ω 0 t, (9)
where, in our simulations ω 0 corresponds to the carrier wavelength of 0.8 µm, E 0 is the field amplitude. The pulse duration T and the ellipticity were chosen to be 20 fs and 0.61 for helium (Fig. 3). For Fig. 2c, T and were 10 fs and 0.6, correspondingly. For the two-color field we have

Ax(t) = - E 0 ω 0 e -2 ln 2 t 2 /T 2 sin ω 0 t + 1 2 sin 2ω 0 t , (10)
Ay(t) = - E 0 ω 0 e -2 ln 2 t 2 /T 2 cos ω 0 t + 1 2 cos 2ω 0 t , (11)
with the same ω 0 , and T = 20 fs.

The attoclock mapping M -1 which maps time t to the direction of the vector potential A(t), that is, φ = M -1 (t) ≡ arctan (Ay(t)/Ax(t)), with φ defined such that it respects the angle definition after Eq. ( 1), is then determined as follows: For the pulse given by Eqs. ( 8)- [START_REF] Balciunas | Optical and THz signatures of sub-cycle tunneling dynamics[END_REF] we have (neglecting slow envelope effects):

t = M(φ) ≡ 1 ω 0 arctan ( tan(φ)). (12)
For small φ can be approximated as M(φ) ≈ φ/ω 0 , ω 0 = ω 0 / . For nearly circularly polarized pulses ( ≈ 1) this expression approaches the "traditional" mapping τ = M(φ) = φ/ω 0 . Both M and M -1 are shown in the inset to Fig. 2b. For Eqs. ( 10)- [START_REF] Silva | Topological strong-field physics on sub-lasercycle timescale[END_REF], the mapping M(φ) is difficult to obtain analytically, but its inverse M -1 (t) is given by:

φ = M -1 (t) ≡ arctan 2 cos (ω 0 t) + cos (2ω 0 t) 2 sin (ω 0 t) + sin (2ω 0 t) , (13)
which is rather good approximated by M -1 (τ ) ≈ 4ω 0 τ /3. That is, τ = M(φ) ≈ φ/ω 0 , ω 0 = 4ω 0 /3.

Jones vector formalism

In deriving the equations for the Brunel harmonics response, we used the Jones vector representation for the field harmonics E(t) = n Ene -inω 0 t . The waveshape described in Eqs. ( 8)- [START_REF] Silva | Topological strong-field physics on sub-lasercycle timescale[END_REF], can be easily translated into En. For instance, a singlecolor field with carrier (fundamental) frequency ω 0 defined by Ex = E 0 cos (ω 0 t), Ey = -E 0 sin (ω 0 t) (which corresponds to the fundamental frequency component of the 2-color field) we can easily obtain using the Jones formalism as:

E 1 = E 0 /2 • (1, -1/i), E -1 = E 0 /2 • (1, 1/i)
where the first and second component of the vector are its x-and y-components, respectively. In the main text of the manuscript, the insignificant pre-factor E 0 /2 is removed since only the polarization state matters.

TDSE simulations

The TDSE equation has the well-known form (we use atomic units here and in the following subsection):

i∂tψ(r, t) = Hψ(r, t), H = (p + A(t)) 2 /2 + V (r), (14)
where ψ(r, t) is the wavefunction of the electron depending on space r (with r ≡ |r|) and time t, H is the electron's Hamiltonian, with p being the momentum operator, A(t) the vector potential of the external (laser) field with the electric field strength E(t) = -∂tA.

The simulation box size used for the integration in the two color case (Fig. 1) was 500 au. In this case, to ensure the numerical convergence, simulations were repeated with both uniform and logarithmic spatial grids. For the single color case, the photo-electron spectrum (PES) in Fig. 2 was extracted by means of the iSURF method [START_REF] Morales | isurf: a family of infinite-time surface flux methods[END_REF] using projection to Volkov states. The box size was taken to be up to 300 au in this case. The convergence of PES was also tested by projecting to the Coulomb states, obtained in a simulation with a box of 3000 au.

Reconstruction of harmonic response from the TDSE simulations

The field Er resulting from the atomic response was calculated as Er = g∂ttp = g∂tj, where p = ψ| r |ψ is the atomic polarization, |ψ is the electron wavefunction, j = ∂tp is the current assigned to p, and g is a constant depending on the position of the observer, which we assume equal to unity for simplicity. Note that this response contains not only Brunel but also all other nonlinearity mechanisms. For the results shown in Fig. 3, the relevant quantities (polarization direction, ellipticity) were averaged in the range (0.75nω 0 , 1.25nω 0 ) for the nth harmonic except for n = 0. For n = 0 (Fig. 1), the average was made in the range (0, 0.25ω 0 ), that is from zero to approx. 100 THz. In all cases, the average was performed using weighting with the intensity distribution. Examples of the spectral dependencies for 0th harmonic are shown in the Supplementary (Fig. 8).

The Drude-like model from R-matrix theory

In this section we derive the Drude-like model Eq. ( 5) using the framework of the analytic R-matrix approach [START_REF] Torlina | Interpreting attoclock measurements of tunnelling times[END_REF][START_REF] Torlina | Time-dependent analytical rmatrix approach for strong-field dynamics. i. one-electron systems[END_REF]. In this approach, the electron wavefunction is written as: ψ(r, t) = ts,p r | p Rpe iS(p,t,ts) where Rp encodes the angular structure of the initial state, ts is a (complex-valued) stationary time (see below), and S = S f + S C , S f = 1/2 t ts-iκ -2 (p + A(t )) 2 dt -Ip(ts -t), S C = t ts-iκ -2 V (r(p, τ, ts))dτ , where κ = 2Ip and Ip is the ionization potential of the atom. The ionization time ts is inferred via the stationary phase approach from the equation ∂t s S = 0. ts is different from the ionization time t 0 in a short range potential defined from ∂t 0 S f = 0, and hereupon we express τ = Re (ts -t 0 ). According to [START_REF] Torlina | Interpreting attoclock measurements of tunnelling times[END_REF], the value of τ can be obtained as a perturbation from the relation -τ ∂t sts S f = ∂t s S C . The approach of [START_REF] Torlina | Interpreting attoclock measurements of tunnelling times[END_REF] can be also seen in the following way: by neglecting second and higher order derivatives of S C as well as third and higher order derivatives of S f , in the vicinity of t 0 we obtain S(p, t, ts) ≈ S f (p, t, t 0 + τ ) + C(t, p). That is, by changing variables t + τ → t we transform (approximately) S → S f + const(t 0 ). Here, we extend this approach, by noting that there always exists a (in general complex) function τ (t, p, t 0 ), that makes the mapping mentioned above exact:

S(p, t, ts) = S f (p, t, t 0 + τ (t, p, t 0 )) + C(t, p) (15)
for arbitrary t 0 . That is, introducing an effective delay τ should effectively map the Coulomb-driven dynamics after the exit of the barrier to the dynamics without the Coulomb tail. Formally speaking, this mapping is general enough to connect arbitrary functions S and S f which are analytic everywhere except at isolated points. This leads us to the following expression for the free wavefunction

ψ(r, t) = t 0 ,p r | p R p e iS f (p,t,t 0 +τ (t,p,t 0 )) , (16)
with R p = iC(t,p) . The physical meaning of τ requires that τ (t, p, t 0 ) → 0 as t → t 0 , since if we consider electrons directly after the tunnel exit, they have not seen the Coulomb tail yet.

On the other hand, as the modulus of T = t -t 0 becomes large, τ (t, p, t-T ) should become independent of T (meaning that at large times after the electrons were born they do not see the Coulomb tail anymore). That is, τ (t, p, t -T ) approaches some limiting function τ (t, p) at large |T |. In this situation, the expression for free electron density ρ(t) can be most straightforwardly and transparently derived assuming a constant C(t, p) ≈= C and only time dependent τ (t, p) ≈ τ (t) at least for those p which contribute mostly to Eq. ( 16) (see below the discussion how these assumptions can be relaxed). In this case we obtain:

ρ(t) = C 0 ρ f (t -τ (t)), (17)
where C 0 = e -2Im C , and ρ f (t) is the free electron density in a short range potential (that is, assuming C = 0, τ = 0 in Eq. ( 16)). Neglecting ∂tτ , we also obtain that the ionization rate W (t) ∝ dρ/dt which obeys the same rule:

W (t) = C 0 W f (t -τ (t)),
and W f (t) is the corresponding ionization rate for the short-range potential. Note that, independently from the simple derivation above, Eq. ( 17) provides a general transformation to map the two electron densities ρ(t) and ρ f (t), except for pathological cases, if ρ(t → -∞) = ρ f (t → -∞) = 0 holds. In this case, τ (t) can be calculated directly from ρ(t) and ρ f (t). That means, Eq. ( 17) is valid also beyond the simplifications used in the derivation above.

In that more general case, the expression relating the delay τ in Eq. ( 17) and in Eq. ( 15) will be more involved.

To complete the derivation of our Drude-like model, we take into account that since the mapping Eq. ( 15) has been applied, we have to consider electrons propagating free from the action of the Coulomb tail. Neglecting quantum interference effects in the electronic wavepacket and the term responsible for the photon absorption [START_REF] Jürgens | Origin of strong-field-induced low-order harmonic generation in amorphous quartz[END_REF][START_REF] Jürgens | Characterization of laser-induced ionization dynamics in solid dielectrics[END_REF][START_REF] Geissler | Light propagation in field-ionizing media: Extreme nonlinear optics[END_REF], we can write the current J(t) formed by free electrons at large enough times t as a sum of "elementary currents" from electrons born at particular times:

J(t) = t -∞ δJ(t, t )dt , with δJ(t, t )dt = -n 0 C 0 W f (t -τ (t ))(A(t) -A(t ) + v 0 (t ))
)dt represent the current from electrons born between t and t + dt with initial velocity v 0 and atomic density n 0 . Differentiating this expression and taking into account Eq. ( 17), we obtain:

dJ dt = n 0 (ρ(t)E(t) -W (t)v 0 (t)) . (18)
Using the relation dρ ≈ W (t)dt, this expression can be reformulated as:

dJ ∝ ρ(t)E(t)dt -v 0 (t)dρ, (19)
which more explicitly underlines that, whereas the latter term is responsible for the birth of new electrons with nonzero velocity, the former one describes acceleration of already born electrons. The harmonics which are typically referred to as "Brunel harmonics" appear due to this former term in the situation when ρ(t) is time-dependent, since if ρ is constant, the corresponding emission contains no harmonics of the driving field. According to typical analytical theories of tunneling (for instance [START_REF] Yudin | Nonadiabatic tunnel ionization: Looking inside a laser cycle[END_REF]) v 0 ≈ 0, thus yielding the Drude-like model Eq.( 3). We also checked that for the low frequencies considered in the present article and the values of v 0 , obtained from the Wigner trajectory approach [START_REF] Yakaboylu | Wigner time delay for tunneling ionization via the electron propagator[END_REF][START_REF] Camus | Experimental evidence for quantum tunneling time[END_REF], the term containing v 0 is also negligible. A figure showing relative influence of these parameters can be found in the Supplementary. The expression Eq. ( 19) with v 0 = 0 is the one shown in the Introduction.

From the expression for J above, one can also see that the average propagation direction at large times is given by J∞ = -C 0 W f (t -τ (t))A(t)dt. That is, the number of electrons dn, born between times t and t + dt flying in the direction φ determined by A(t) is given by C 0 W f (t -τ (t))dt. At the same time by the definition of P , the same number dn is given by P (φ)dφ. Comparing these two expressions, we obtain:

C 0 W f (t -τ (t))dt/dφ = W (t)dt/dφ = P (φ), (20)
which gives us Eq. ( 4), taking into account the attoclock mapping t(φ) = M(φ). In the other words: if we define a quantity We(t) according to Eq. ( 4), Eq. ( 20) shows us that W (t) defined above is identical to We(t). That is, the effective ionization rate We(t) obtained by mapping the angular electron spectrum into time with M is equivalent to the ionization rate, we observe from the Brunel harmonics, as given by the Drude-like model Eq. ( 18) or Eq. ( 3). As we see from the derivation above for low harmonics there is a one-to-one correspondence between wavefunction shape and the ionization rate, reflecting the effective attoclock delay τ . We see also that the delay imaged by Brunel harmonics is the same as the attoclock delay extracted from the electron distribution. This attoclock delay is also in a close connection with the Wigner delay -both are different (but closely related) measures of the group delay of the electronic wavepacket at large times. As we have seen here, if we consider τ as a function of time rather than a single number, this allows to track not only delay but also asymmetry in ionization.

As a short outlook concluding this section, for higher Brunel harmonics we can neglect neither v 0 in Eq. ( 18), nor T -dependence in τ (t, p, t -T ), thus making harmonics dependent on the electron dynamics even closer to the core. This consideration is however beyond the scope of this paper.

Details of reconstruction of W (t) and estimation of errors for hydrogen and helium

For the reconstruction of the W i in Eq. ( 7), Fig. 2, Fig. 3 we use the following procedure: for known delay τ (which we take from [START_REF] Eckle | Attosecond ionization and tunneling delay time measurements in helium[END_REF]) and W 0 we numerically find the values of W 2 , W 4 , δ 2 , δ 4 such that the measured value of r is achieved, and the center of mass of the reconstructed distribution is located at τ . This procedure involves numerical minimization, and gives us a unique solution for the above parameters. The value of W 0 was estimated from the ionization rate at the minimum strength of our elliptical pump field. The particular value of W 0 does not influences significantly the reconstructed dynamics. In particular, despite the analytical formula gives the symmetric ionization rate, this does not influences the asymmetry reconstructed using our procedure. We checked this with different ionization models (tunnel, ADK, Yudin-Ivanov) and obtained very similar results for the rest of W i .

Although the observation of the 3rd harmonic provides only limited information, it is already sufficient to reconstruct the ionization 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 potential and |r| = 0.27 for the Coulomb potential. It is reasonable to expect, that the next term in the expansion of the ionization rate is of the order of |r| 2 which gives us 7% uncertainty from the unknown contribution W 6 for the Coulomb potential and 2% for the Yukawa one. This determined error bars along y-axis shown in Fig. 2c. Besides, the uncertainty in δ 6 can lead to shifts along the time axis, which was taken into account by adding error bars corresponding to deviations along that axis with the amplitude equal to |r| 2 . In Fig. 3 the reconstruction of ionization dynamics in helium was influenced by the impact of bound-bound transitions (described in the text), errors to the experimental reconstruction were estimated using analytic signal approach discussed in Supplementary Information. In Fig. 3, the pulse energy measured in the experiment can not be directly converted to the peak intensity because of the strong influence of defocusing. We used rescaling delivered by computation of nonlinear unidirectional pulse propagation equation for 140 µJ, giving the peak intensity of 1000 TW/cm 2 and the waist radius 5.8 µm. Still, we note that nonlinear propagation effects in the pump introduce certain uncertainty. To take this into account, we introduced intensity error bars of 10%. The details of the numerical simulations can be found in the Supplementary Information. This energy-to-intensity rescaling was used for the whole intensity range.

FIG. 3 .

 3 FIG. 3. Experimental reconstruction of the subcycle ionization dynamics in He, compared to the theoretical simulations. (a) Experimental setup (see Methods for more details): A 800-nm, 43-fs long elliptically-polarized pump pulse (P) is focused into a plasma spot (S) where 3rd harmonic (H) is generated. The polarization components are carefully separated and detected. (b) Experimentally measured rotation of the long axis of the 3rd harmonic with respect to the long axis of the pump (solid blue curve, left axis) and ellipticity of the 3rd harmonic (solid red curve), in comparison with the results of TDSE calculations (dashed blue and red lines). The shaded region shows uncertainty introduced by influence of bound-bound transitions. Further errors are discussed in Methods. (c),(d) Reconstruction of the coefficients Wi, δi describing ionization dynamics as given by Eq. (7) (with denotations as in (b)). (e),(f) Ionization dynamics W (t) reconstructed for an exemplary intensities I = 1030 TW/cm 2 (e) and I = 1820 TW/cm 2 (f) according to experiment (thick red solid lie) and theory (blue dashed line). Blue dotted line indicates the theoretical curve with bound-bound-transitions correction included. Thin red solid line indicates the time inverted t → -t version of W (t) and yellow-dashed region between red lines thus indicates the asymmetry of W (t) (in (e) the asymmetry is almost invisible).

ACKNOWLEDGMENTS I.B., A.D. and U.M are thankful to Deutsche Forschungsgemeinschaft (DFG) (projects BA 4156/4-2, MO 850-19/2) as well as Cluster of Excellence PhoenixD (EXC 2122, Project ID 390833453) for financial support. O.G.K., I.A.N., N.A.P.,D.E.S. thank the support from Russian Science Foundation, grant Nr. 21-49-00023 and National Natural Science Foundation of China 12061131010. SS acknowledges support by the Qatar National Research Fund (Grant NPRP 12S-0205-190047) and HPC resources from GENCI (Grant # A0080507594).

AUTHOR CONTRIBUTIONS

I.B. suggested the idea, developed the quantum Drude-like model, performed simulations and analytics based on this model, performed TDSE simulations; A.J.G performed TDSE simulations for the two-color case; S.S.,O.G.K.,I.A.N., N.A.P., and D.E.S. performed simulations of unidirectional propagation equations; J.R.C.A, M.K., D.Z., L.S. and V.V., performed the experiment; I.B., A.J.G, A.H., F.M., L.B., S.S., U.M., A.D., T.N., M.V. and M.I. analyzed and interpreted the results of the simulations and experiments; all authors participated in the article formulation and writing.