Internal Surface Plasmon Excitation as the Root Cause of Laser-Induced Periodic Plasma Structure and Self-Organized Nanograting Formation in the Volume of Transparent Dielectric

Abstract

A computer simulation of the dynamics of an optical discharge produced in the volume of a transparent dielectric (fused silica) by a focused femtosecond laser pulse was carried out taking into account the possibility of developing small-scale ionization-field instability. The presence of small foreign inclusions in the fused silica was taken into account with the model of a nanodispersed heterogeneous medium by using Maxwell Garnett formulas. The results of the calculations made it possible to reveal the previously unknown physical mechanism that determines the periodicity of the ordered plasma-field structure that is formed in each single breakdown pulse and is the root cause of the ordered volume nanograting formation in dielectric material exposed to a series of repeated pulses. Two main points are decisive in this mechanism: (i) the formation of a thin overcritical plasma layer at the breakdown wave front counter-propagated to the incident laser pulse and (ii) the excitation of the “internal surface plasmon” at this front, resulting in a rapid amplification of the corresponding spatial harmonic of random seed perturbations in the plasma and formation of a contrast structure with a period equal to the wavelength of the surface plasmon (0.7 of the wavelength in dielectric).

1. Introduction

The search for a mechanism predetermining the laser-assisted formation of a regular periodic subwavelength structure inside a solid transparent dielectric continues to be one of the unsolved problems in nanophotonics and laser-matter-interaction physics. The nanogratings formed in the volume of dielectric by fs laser pulses are already widely used in various fields of modern nanoengineering and nonlinear optics. Nevertheless, a clear physical model that reveals the nature of the periodicity of the underlying plasma-field structures formed during the optical breakdown of the dielectric and allows the calculation of their most important characteristic — the spatial period, is still missing. As the numerous experiments have shown, the nanograting in the material of dielectric (predominantly, the fused silica) is oriented perpendicular to the laser pulse polarization, has a period less than the wavelength in dielectric, and can be formed only by the multiple repeated breakdown pulses [1,2,3,4,5,6,7,8,9,10]. Though the intervals between successive pulses (~${10}^{-6}-{10}^{-4}$ s) are much greater than the disappearing time of the plasma structure formed in every single pulse, it is evident that only this structure (at high enough contrast of the its inherent profile of the energy deposition density) can be responsible for the periodicity of the dielectric material properties appearing after a series of repeated pulses.

In the last decade, the main attention in studying these (one-pulse-created) ionization structures was paid to numerical simulation of laser breakdown initiated in fused silica by randomly placed small inclusions with a lowered threshold of ionization (nanobubbles or atom structure defects) playing the role of seed ionization centers. In a number of calculations it was shown that multiple plasmoids formed around these centers during optical discharge evolution and showed a tendency to align in more or less ordered rows forming on average a kind of spatial grating with a period about half to a whole wavelength [11,12,13,14]. These results demonstrate the “plasma-field” nature of the considered structure and the “single-pulse” origin of its periodicity; however, being only based on the direct numerical integration of Maxwell’s equations (with the FDTD method) they cannot reveal the concrete physical nature of periodicity and still do not allow to go beyond some common and evident references to wave interference phenomenon. Alternative models of the formation of nanogratings based on the theory of exciton-polariton interaction [15,16] also encounter difficulties, since they require too low electron density in the discharge, and more importantly, the spatial period in these models is actually an external parameter specified from the outside.

To provide a sufficiently clear physical model that allows us to identify the specific physical mechanism of the emergence of order from chaos in “multiplasmoid” discharge and calculate (at least qualitatively) the main characteristics of the arising ordered structure and their dependencies on external parameters, it seems natural to use approaches based on the concept of ionization-field instability of a medium exposed to high-intensity electromagnetic radiation. The main tasks in these approaches are: (i) description of the joint evolution of the electric field and plasma density perturbations against the background of some quasi-stationary equilibrium states, (ii) obtaining dispersion equations connecting the growth rate constants (so called increments of instabilities) of these perturbations with characteristics of their spatial structure (spatial periods or wave numbers), and (iii) finding the optimal conditions providing enough fast growth of unstable perturbations in a definite range of the wave numbers at the linear and nonlinear stages of instability (corresponding to small and large perturbations, respectively). As applied to the above nanograting problem this program was realized within the frameworks of simple quasi-one-dimensional models of optical discharge maintained in fused silica (in the absence [17,18] or presence [19] of easily ionizable small inclusions) by a traveling plane wave. These models gave results agreeing with experiments in respect of both grating orientation and observed period range of most quickly growing spatial modulation. In particular, one of the ionization-field instabilities (so-called “plasma-resonance”, one caused by quasi-static mutual amplification of the normal electric field component and plasma density within thin under-critical plasma layers) leads to fast growing spatial modulation of plasma and energy deposition density with a maximum of time growth rate lying within the required wave number range (corresponding to the periods from one-third to two-thirds of the laser wavelength in dielectric) [17,19]. A significant drawback of these models from the point of view of explaining the appearance of a regular periodic structure is, nevertheless, the width of this maximum in the space of wave numbers is too large to allow certain spatial harmonics to stand out in the spectrum of random small (seed) perturbations, which would lead, at the nonlinear stage of instability, to the formation of an ordered periodic structure during one breakdown pulse. To form such a structure, the spectrum of eigenwaves of the considered nonlinear dynamical system should apparently contain, along with the given pump wave, another eigenwave, the excitation of which (if its wave number coincides with the wave number of the spatial spectrum of the initial noise) would lead to the appearance of a corresponding sharp resonant maximum in the spectrum of instability increments. In principle, the role of such an additional eigenwave in an optical discharge plasma could be played by a longitudinal Langmuir wave (bulk plasmon), the resonant excitation of which (due to the nonlocality of the polarization response), as it was shown in the work [20], can be significant in optical breakdown of dense gases. However, in the laser plasma parameters range typical of the experimentally observed formation of nanogratings in a solid dielectric, the nonlocal correction to polarizability is negligible compared to its imaginary part and the longitudinal wave is completely suppressed by absorption due to collisions. Therefore, the bulk plasmon cannot be excited in the experiments, and apparently there is no need to take into account the nonlocality of the plasma polarizability (which ensures the possibility of its excitation) in attempts to elucidate the essence of the physical mechanism of nanograting formation.

The difficulty associated with the absence of an additional eigenwave in the volume of a homogeneous medium does not exist in explaining the mechanism of formation of a periodic structure under conditions of optical breakdown at the dielectric–vacuum interface, where the role of such a wave is played by the surface plasmon excited at the interface, the wavelength of which determines the spatial period of the formed surface nanograting [5]. Such surface plasmon can be excited, however, during breakdown inside a homogeneous medium, if a nonuniformity of the parameter determining the ionization rate of the optical field was created in it previously (for example, as a result of the modification of the medium by previous pulses). This is indicated, in particular, by some results of the numerical calculations described above [11,12], according to which an ordered bulk grating was formed at the junction of two half-spaces with different concentrations of seed ionization centers. The possible role of a surface plasmon excitation on the interface between the unperturbed region of dielectric and the region modified by previous pulses is qualitatively considered also in [21].

In the present work, it is shown that the instability of the discharge with respect to the excitation of the “internal” surface plasmon and the periodic plasma field nanostructure generated by it can arise already at the front of the breakdown wave in a single pulse inside a macroscopically homogeneous medium and does not require the creation of a preliminary inhomogeneity of the ionization rate inside it. To describe this phenomenon, the quasi-one-dimensional electrodynamic models, previously used by us, which virtually ignored the role of the longitudinal (parallel to the wave vector of the pump wave) components of the electric field, which are necessarily present in the surface plasmon, were refined, taking into account the excitation of this component on two-dimensional plasma inhomogeneity. As shown by numerical calculations and estimates based on a relatively simple theoretical model that takes into account the generation of a surface plasmon at the breakdown wave front, an ordered subwavelength nanograting is formed in this case in nanodispersed (multyplasmoid) discharge with an arbitrary noise spectrum of small random perturbations. The work is organized as follows. In the second part, the physical model used is described and the basic equations are formulated that describe the spatiotemporal evolution of the discharge in a nanoporous dielectric, which is considered as a continuous heterogeneous medium that can be described on the basis of Maxwell Garnett formulas. The third part presents the results of a numerical calculation that demonstrates the formation of contrast ordered periodic structures of plasma and energy deposition densities, developing from small initial perturbations randomly distributed in space. In the fourth part, a simple model is proposed that makes it possible to qualitatively explain the occurrence of a periodic structure as a result of the excitation of an unstable surface plasmon. In conclusion, the main results of the work are formulated.

2. The Physical Model. Initial Equations and Approximations

Our relatively simple (two-dimensional) physical model allows taking into account key electrodynamic factors that determine the possibility of the internal surface plasmon excitation and forming a regular periodic plasma-field structure within a solid dielectric by a single focused laser pulse. The first of these factors is the field amplitude increase in the longitudinal direction inside the beam when approaching the focal plane. It is thanks to this increase that during the ionization process a rather sharp front of the breakdown wave is formed, which can support surface plasmon. The second factor is the relationship between the transverse and longitudinal components of the electric field in the presence of a two-dimensional inhomogeneity of the plasma density. Due to this relationship, the surface plasmon is unstable with respect to arbitrarily weak seed perturbations and imposes its periodicity on the formed structure at the nonlinear stage of instability. Both of these factors are taken into account in our equation:

$$\begin{array}{c}\frac{i}{c^2}\frac{d({\omega }^2\epsilon )}{d\omega }\frac{\partial H_y}{\partial t}\end{array}+\frac{1}{S}\frac{\partial }{\partial z}\left(S\frac{\partial H_y}{\partial z}\right)+\frac{{\partial }^2{H}_y}{\partial x^2}+\frac{{\omega }^2}{c^2}\epsilon H_y-\frac{\nabla \epsilon }{\epsilon }\nabla H_y=0,$$

(1)

for the complex amplitude (slow time envelope) of the magnetic field of a quasi-monochromatic TM wave and the equations for the complex amplitudes of the transverse ($E_x$) and longitudinal ($E_z$) components of the electric field:

$$\frac{\partial E_x}{\partial t}-i\omega E_x=-\frac{1}{\epsilon }\frac{\partial H_y}{\partial z},$$

(2)

$$\frac{\partial E_z}{\partial t}-i\omega E_z=\frac{1}{\epsilon }\frac{\partial H_y}{\partial x}$$

(3)

where all the field components are assumed to be written as real parts of the productions of their complex amplitudes and factor $\exp (-i\omega t)$. In Equations (1)–(3), $\epsilon$ is a dielectric function of the ionized medium. When using Equation (1) for a slow field envelope, it is assumed that the characteristic time of the change in the complex field amplitude and plasma density is large compared to the field period $2\pi /\omega$. The function $S(z)$ in the second term of Equation (1) can be considered here as the effective ray tube cross-section in the near-axis region of the focused wave beam, the field of which is modeled by the above equations. Equation (1) can be obtained from the exact wave equation based on the following approximations: (i) we neglect the second derivative of the field amplitude with respect to time; (ii) we replace the Laplace operator with the mentioned term, which allows us to describe the axial changes in the field amplitude of the focused beam in its axial region under the assumption of a given (converging) ray tube and the field amplitude being constant on the wavefront surface orthogonal to the rays. Similar approximations in the study of the breakdown wave dynamics in a paraxial wave beam were used earlier, for example, in [22]. Modifying the wave equation in this way, we get the opportunity to describe the (important to us) longitudinal changes in amplitude due to the focusing of the beam, at the same time distracting from its large-scale transverse inhomogeneity. Consideration of the latter would be important only in calculating the total transverse size of the emerging structure and is not directly related to the main issue that interests us here, concerning the nature of the small-scale transverse periodicity of the field. In fact, ignoring the large-scale transverse inhomogeneity in the framework of this model, we attribute our study to the structure created in the wide axial region of the long-focus beam. In the calculations below, we take for the ray tube cross-section an expression $S=1+z^2/l_F^2$ corresponding to an unperturbed Gaussian beam with a Rayleigh length $l_F$ and a known axial distribution of the field:

$$E_x^{(0)}=\frac{H_y^{(0)}}{\sqrt{{\epsilon }_s}}=E_F{\left(1+\frac{z^2}{l_F^2}\right)}^{-1/2}\exp \left(i{k}_{s}z-\frac{{(t-z\sqrt{{\epsilon }_s}/c)}^2}{{\tau }^2}\right),k_s=\frac{\omega }{c}\sqrt{{\epsilon }_s},$$

(4)

obtained in this case for a Gaussian pulse (with a duration ${\tau }_p$ and maximum electric field amplitude $E_F$) directly based on the solution of Equations (1)–(3) in the absence of ionization (i.e., in the homogenous dielectric with permittivity ${\epsilon }_s$). Additional formal justifications for the possibility of using the simplified wave Equation (1) are given in Appendix A.

To study the joint evolution of the field and plasma in the discharge produced by such a pulse, which is our main task here, a numerical solution of Equations (1)–(3) is required together with equations describing the kinetics of the discharge and dielectric function of the ionized medium. In the case of apparently the most practical interest, we should consider the discharge arising in the dielectric (fused silica) as a result of the ionization of multiple foreign inclusions with a reduced ionization threshold that can play the role of primary centers of breakdown. Generally speaking, a fairly complete description of the dynamics of such a discharge is a very complex and time-consuming task that requires taking into account the spatiotemporal evolution of many inhomogeneous plasmoids formed near each of primary ionization centers, describing the processes of their expansion, shape change, and merging (see, for example, [11,12,13,14]). However, the essence of the physical mechanism that determines the periodicity of arising structures, as it seems to us, can be revealed in the framework of the simplified model that we use. We assume that ionization occurs in the volume of many small spheres with a lowered ionization threshold $U=5.2$ eV [12,14,23], the sizes and number of which are fixed. Within this approach, a microscopically inhomogeneous medium can be considered as a continuous heterogeneous one; the fields in the Equations (1)–(3) should be considered as their macroscopic means; and the dielectric function $\epsilon$ should be replaced by its effective value ${\epsilon }_{eff}$, which is determined based on the well-known Maxwell Garnett formulas:

$${\epsilon }_{eff}={\epsilon }_s\frac{1+2f\alpha }{1-f\alpha },\alpha =\frac{{\epsilon }_p-{\epsilon }_s}{\epsilon +2{\epsilon }_s},$$

(5)

$${\epsilon }_p={\epsilon }_s-\frac{N}{N_c}\left(1-i\nu /\omega \right),$$

(6)

where ${\epsilon }_p$ is permittivity of the plasma, $N$ is free electron density, $N_c=m({\omega }^2+{\nu }^2)/(4\pi e^2)$ is its critical value, $\nu$ is the effective frequency of electron collisions, and $f$ is the volume fraction of spherical plasmoids. Figure 1 shows the dependences of the real and imaginary parts of the effective permittivity ${\epsilon }_{eff}$ on the dimensionless plasma density $n=N/({\epsilon }_s{N}_c)$ at the values $\nu /\omega =0.15$ (which seems to us the most realistic for the conditions for producing nanogratings in fused silica based on the results of [23,24,25]) and volume fraction $f=0.3$ (close to that which can be extracted from references [11,12,13,14]).

The change in plasma density is determined by the rate equation

$$\frac{\partial N}{\partial t}=W_{pi}+W_a-\frac{N}{{\tau }_r}+W_n.$$

(7)

The first three terms on the right-hand side of this equation describe multiphoton and avalanche ionization, and recombination. We have used the approximations for dependences of their rates on laser intensity similar to the ones used in the previous studies [23,24,25] of the near-infrared (800 nm) breakdown of the nanoporous fused silica with lowered bandgap $U=5.2$ eV:

$$W_{pi}={\sigma }_3(N_m-N){\left({\left|\eta \right|}^{2}I\right)}^3,W_a=\frac{8\pi \nu e^2}{3\sqrt{{\epsilon }_s}mcU{\omega }^2}\frac{N_m-N}{N_m}{\left|\eta \right|}^{2}IN,$$

(8)

$$\eta ={[1+(1-f)({\epsilon }_p-{\epsilon }_s)/(3{\epsilon }_s)]}^{-1}.$$

Here $I=c\sqrt{{\epsilon }_s}{\left|E\right|}^2/8\pi$ is the laser intensity, $\eta$ is the proportionality coefficient between the field inside the sphere and the averaged macroscopic field, ${\sigma }_3=1.5\times {10}^{-28}$ W⁻³cm⁶s⁻¹, $N_m=2.1\times {10}^{22}{\mathrm{cm}}^{-3}$ is the atom density, and ${\tau }_r=150$ fs [26] is the characteristic recombination time. The last term $W_n$ is a “noise” ionization source introduced into the equation for modeling small seed fluctuations of the plasma density at the initial stage of breakdown in the main region of the discharge ($\left|z\right|<l_F$). It should be noted that the introduction of this source is only a formal technique, allowing us to take into account small initial (seed) fluctuations caused by the statistical spread (even if very small) of the parameters of discrete inclusions in the real discharge. If we try to give this source a more specific and clear meaning and talk about the actual origin of the fluctuations, we can, for example, consider it as a fluctuating component of the ionization rate due to the inevitable presence of small noise fluctuations in the field amplitude in the incident laser pulse.

The equations describing the spatiotemporal evolution of the field and plasma were numerically solved in a rectangular region $0\le x\le L_x$, $L_1\le z\le L_2$. For simplicity, the boundary conditions along the transverse coordinate $x$ were assumed to be periodic with a common period of several wavelengths ($2<L_x/\lambda <15$) in the unperturbed dielectric and significantly exceeding the periods of spatial harmonics of the field and plasma perturbations of interest to us. The boundaries along the longitudinal coordinate $z$ were chosen so that the computational domain includes with sufficient margin the entire region ($\left|z\right|<l_F$) in which the plasma produced can affect the field. In order to get rid of unwanted reflection of waves from the z-boundaries of the computational domain, absorbing layers were placed at these boundaries, providing a sufficiently low reflection coefficient (~10⁻⁵). A noise ionization source generating a random spatial distribution of small initial density fluctuations was specified as a random set of sufficiently large number ($10\times 10=100$) of sinusoidal harmonics:

$$\begin{array}{l}{W}_n=A{\cos }^2\left(\frac{\pi }{2l_F}z\right){\displaystyle \sum _{m=1}^{10}{\displaystyle \sum _{n=1}^{10}\cos \left[\frac{\pi m}{L_x}(x+{\alpha }_{mn}z)+{\varphi }_{mn}\right]}}\mathrm{at}|z|<l_F,t_1<t<t_2,\\ W_n=0\mathrm{at}\left|z\right|\ge l_F,t<t_1,t>t_2,\end{array}$$

(9)

where ${\varphi }_{nm}$ and ${\alpha }_{nm}$ are random variables (with equal probability distributed over the segments $[0,2\pi ]$ and $[-1,+1]$, respectively) that determine the random phases and directions of the wave vectors of spatial harmonics of perturbations. The amplitude $A$ and times of the beginning ($t_1$) and end ($t_2$) of the action of the noise source were set so that by the time when the background (unperturbed) plasma density reached a predetermined small value $N=0.01N_c{\epsilon }_s$, the maximum value of the density fluctuations was 0.01 of this value.

Along with the plasma density and field, we also calculated the spatial distribution of the local energy deposition density in the medium $w(x,z,t)$ by the given time $t$ (see also [27]):

$$w=f{\displaystyle {\int }_{-\infty }^t\left[\nu \frac{N}{N_c}{\left|\eta E\right|}^2+\left(8\pi U+\frac{1}{2}\frac{{\left|\eta E\right|}^2}{N_c}\right)s\frac{\partial N}{\partial t}\right]}dt.$$

(10)

Here, $s=1$ at $\partial N/\partial t>0$ and $s=0$ at $\partial N/\partial t<0$; the first term under integral describes the electron collision losses, the second one takes into account the energy expenditure for the interband transition (overcoming the bandgap $U$) and energy transfer to the newly born free electrons. This function is an important characteristic determining the rate and space profiles of the medium heating and therefore the course of the thermo-mechanical and chemical processes and accumulation effects [14,28] responsible for the volume nanograting formation in the medium by the repeated pulses. The high contrast nanograting can be formed evidently if by the end of the pulse the sufficient energy deposition $w$ is localized in a comparatively narrow regions.

3. Results of Calculation

Numerical simulation was performed at laser pulse parameters typical of the experiments on the formation of nanogratings in the volume of fused silica: the maximum value of the unperturbed intensity at the focus $I_{\max }=c\sqrt{{\epsilon }_s}|E_F{|}^2/(8\pi )$~${10}^{13}-{10}^{14}{\mathrm{W}/\mathrm{cm}}^2$, the dimensionless focal length $k_s{l}_F=35$ (beam convergence angle $\vartheta ={10}^{°}$), and the pulse duration $\tau =100\mathrm{fs}$. Initially, the discharge dynamics were simulated in the absence of a noise ionization source that produces small seed fluctuations in the plasma density $N$ in spherical plasmoids and the effective dielectric function ${\epsilon }_{eff}$ determined by it. In this case, the general scenario for the development of the discharge is characterized by the formation of a one-dimensional breakdown wave, which arises in the focus region and counter-propagates to the incident beam. The nature of the evolution of the longitudinal profile $N_0(z,t)$ of this wave and the maximum density value $N_{0\max }$ achieved at its front depends on the beam intensity $I_{\max }$ and changes significantly when passing through some critical values $I_1\approx 2\times {10}^{13}$ W/cm² and $I_2\approx 7.3\times {10}^{13}$ W/cm². Figure 2 shows the spatiotemporal evolution of dimensionless density profiles $n_0=N_0(z,t)/({\epsilon }_s{N}_c)$ at three intensities. At low intensities $I_{\max }<I_1$ (Figure 2a), i.e., at a relatively low ionization rate, the plasma density in spherical plasmoids does not have time to reach a threshold value $N_{th}$ corresponding to the condition $\mathrm{Re}{\epsilon }_{eff}=0$ during the pulse width. In the region $I_1<I_{\max }<I_2$ (Figure 2b), a thin (compared to wavelength) layer with a suprathreshold density ($N_{0\max }>N_{th}$) is formed, the profile of which remains almost unchanged for quite a long time (actually from the moment the density passes through this threshold until the end of the ionization process).

A further increase in intensity $I_{\max }>I_2$ (Figure 2c) leads again to an expansion of the density profile and a decrease in its maximum because of the formation of an extended but not sufficiently ionized region in the front (far extending toward the beam) part of the discharge. The absorption of the incident wave in this region, due to its large extent, leads to a strong decrease in the field and the actual stop of the breakdown process at a relatively low level of plasma density already in far approaches to the focal region.

The main characteristics of the suprathreshold plasma layer formed in the discharge as functions of the intensity are illustrated in Figure 3, which shows the intensity dependences of the dimensionless maximum density,$n_{0\max }=N_{0\max }/({\epsilon }_s{N}_c)$, maximum suprathreshold layer thickness $L_{\max }$, (Figure 3a), and its lifetime $T_0$, (Figure 3b). As can be seen from the figure, in the entire indicated range of intensities, the formed layer is characterized by a relatively small thickness (shorter than the wavelength) and a sufficiently long lifetime (~60–80 fs), which contributes to the excitation of an internal surface plasmon on this layer and to the formation of a periodic structure (see below).

The spatiotemporal evolution of the discharge in the presence of small seed noise (Figure 4), which makes it possible for the instability we are interested in, is illustrated in Figure 5, Figure 6 and Figure 7 for the intensity value $I_{\max }=4\times {10}^{13}$ W/cm² lying in the intermediate region of greatest interest $(I_1,I_2)$, where the peak density at the front of the unperturbed breakdown wave goes over a threshold value $N_{th}$. It is after such a transition (Figure 5) that transverse density modulation (with a period approximately equal to $0.67\lambda$) arises and rapidly increases on this front, leading to the formation of a contrasting periodic structure. We note that randomly placed plasmoids are not shown in our figures. Their coordinates do not appear in our consideration; their sizes are much smaller than the characteristic scale of the shown density distributions. The fact that there is a close correlation between the effects of the formation of a thin layer of above-threshold density (with $N_{0\max }>N_{th}$, i.e., $\mathrm{Re}{\epsilon }_{eff}<0$) and its strong transverse modulation, together with the results of a qualitative analysis based on the simple model presented below, allows us to conclude that the excitation of the internal surface plasmon is crucial.

The appearance of this plasmon as a new natural wave in the system under study, as shown in Section 4, leads to a resonant increase in the instability increment of the spatial harmonic of perturbations synchronous with it (whose wave vector coincides with the surface plasmon wave vector). The released harmonic relatively quickly leads to a strong corrugation of the breakdown wave front with a period equal to plasmon wave length; due to the strong modulation of the field amplitude (and hence the ionization rate) along the front, the discharge decays into separate layers with suprathreshold density (Figure 6 and Figure 7). At the final stage of the process, these layers (in fact the same as in the case of an initially purely harmonic perturbation with a given period [18,27]) lead to the formation of a contrast grating structure characterized by a strong nonuniformity in the transverse direction distribution of the energy deposition density (Figure 8). Within the parameter range used in calculations, the spatial period of this structure $\mathsf{\Lambda }$ turned out to be close to $0.67{\lambda }_s\approx 350$ nm, which fits into the framework determined by the experimental data. It is also important to note that, as the simulation showed, this quantity (in the interval of parameters where an ordered structure appeared) did not depend on the concrete realization of the random initial distribution of small perturbations and the sizes of the calculation region. However, there are some factors that are not accounted for in our model but can affect the plasmon wavelength and, therefore, the period of the structure. Probably the non-sphericity of the small plasmoids formed under real conditions in the nanodispersed discharge belongs to them in the first place. This factor can be analyzed within the framework of a more advanced model (based on the same approaches as the model we proposed), in which the discharge in fused silica is described as anisotropic nanodispersed media, whose permittivity tensor components are calculated on the basis of extended Maxwell Garnett theory. Also, the plasmon wavelength can be influenced by the curvature of the breakdown wave front, which arises due to the smooth transverse inhomogeneity of the laser beam field and by the inhomogeneity of the ionized medium due to its heating by a series of multiple laser pulses. Taking these factors into account, however, significantly complicates the model and goes beyond the scope of our work, the purpose of which is to clearly demonstrate, based on the simple model, the physical mechanism determining the periodicity of the nanogratings.

Outside the intensity interval $(I_1,I_2)$ (i.e., at $N_{0\max }<N_{th}$) it can lead, during the evolution of the discharge, to the formation of small-scale high-density regions, but the regular periodic structure of interest to us does not form here.

4. Instability of a Thin Plasma Layer with Respect to the Excitation of a Surface Plasmon

As a structure simulating an unperturbed plasma density profile at the moment of formation of a narrow density peak in the breakdown wave (Figure 2b), we consider a stationary discharge supported in the form of a thin uniform plasma layer by a given external field $x_0{E}_0\exp (-i\omega t)$ parallel to its boundaries $z=\pm l$ (see Figure 9a); the stationary plasma density in the layer is $N_0$. Its dielectric constant ${\epsilon }_0$, in view of the relatively small value of the parameter $\nu /\omega =0.15$, used in our calculations, is assumed to be purely real in the simplified model under consideration. As is known, various types of surface waves can propagate along the plasma layer. Among them, the main interest for us is the “antisymmetric” surface plasmon, in which the field component $E_x~\exp \{-i(\omega t-hx)\}$ is an odd function of the coordinate $z$ perpendicular to the layer. As follows from the dispersion equation:

$$D(h)\equiv \tanh {\kappa }_{0}l+{\epsilon }_0\frac{{\kappa }_e}{{\kappa }_0}=0,{\kappa }_0=\sqrt{h^2-k_0^2{\epsilon }_0},{\kappa }_e=\sqrt{h^2-k_0^2},k_0=\frac{\omega }{c},$$

(11)

(see, for example, [29]), for a small layer thickness ($k_{0}l<<1$) this plasmon can exist even if the stationary plasma density $N_0$ slightly exceeds a threshold value $N_{th}$ (i.e., for $-{\epsilon }_0<<1$). Consequently, as the plasma density in the layer increases from zero, this plasmon begins to propagate first (see Figure 9b), and at the threshold of occurrence (${\epsilon }_0\approx -2k_{0}l$) it has a wavelength ${\lambda }_0=2\pi /h_0=\lambda /\sqrt{2}$, regardless of the specific values of the layer parameters.

We now turn to an analysis of the instability of the discharge in the layer with respect to perturbations of the plasma density and the field amplitude of the form $N_1,E_1~\exp (\gamma t)\cos (qx)$ and write down the equations relating these quantities in the plasma in the linear approximation [17,19]

$$\begin{array}{c}\Delta E_{1x}+k_0^2{\epsilon }_0{E}_{1x}-E_0{k}_0^2\frac{N_1}{N_c}-E_0\frac{{\partial }^2}{\partial x^2}\frac{N_1}{N_c}=0,\\ \frac{\partial N_1}{\partial t}=\beta \frac{E_1}{E_0}-\alpha N_1, \alpha =N_c\frac{\partial W(E_0,N_0)}{\partial N_0}, \beta =\frac{E_0}{N_c}\frac{\partial W(E_0,N_0)}{\partial N_0},\end{array}$$

(12)

where $W$ is the total ionization rate. For a given dependence of perturbations on $x$, Equation (12) (taking into account the boundary conditions for the continuity of the tangential components of the electric and magnetic fields at the layer boundaries) leads to the following problem of finding the eigenfunctions $E_{1x}(z)$ and eigenvalues $p_i$:

$$\begin{array}{c}\frac{{\partial }^2}{\partial z^2}{E}_{1x}-p_i^2{E}_{1x}=0,\\ {\left(\frac{\partial E_{1x}}{\partial z}\pm \frac{p_0^2}{\epsilon p_e}{E}_{1x}\right)|}_{z=\pm l}=0,\\ p_0^2=q^2-k_0^2\epsilon ,p_e^2=q^2-k_0^2,\\ p_i^2=p_0^2+\left(1-\frac{q^2}{{\epsilon }_0}\right)\frac{\beta }{\gamma -\alpha }.\end{array}$$

(13)

The indicated problem has two types of solutions (symmetric and antisymmetric in the $z$ coordinate), the eigenvalues of which $q_i$ are found as the roots of the equation

$$\tanh p_{i}l=-{\left(\frac{\epsilon p_e{p}_i}{p_0^2}\right)}^{\pm 1},$$

(14)

(plus and minus correspond to symmetric and antisymmetric mods), determining (based on the last of the Equation (13)) time increments of perturbations as a function of their wave numbers $q$ (or periods $\mathsf{\Lambda }=2\pi /q$); in particular for antisymmetric perturbation

$$\gamma (q)=\alpha +\frac{l^3{\left(q^2-k_0^2\epsilon \right)}^{3/2}}{3{\epsilon }_0}\frac{\beta }{D(q)}.$$

(15)

An analysis of Equation (14) shows that, at low concentrations (${\epsilon }_0>0$), the increments of both types of perturbations are finite, however, with increasing concentration in the layer, the possibility of excitation of a surface plasmon appears, as a result of which the instability increment $\gamma (q,{\epsilon }_0)$ becomes infinitely large for an antisymmetric perturbation with a wave number equal to the wave number of the surface plasmon $q=h_0$. Under realistic conditions (at a finite frequency of collisions) this infinity is obviously absent, however, as estimates made on the basis of Equation (15) show, the perturbation increment with a period equal to the plasmon wavelength $\gamma (h_0,{\epsilon }_0)~\beta \omega /\nu$ significantly exceeds the increments of all other perturbations, this harmonic quickly goes into the nonlinear stage and after its standing out determines the further evolution of the discharge as a whole. Thus, the excitation of a surface plasmon is apparently the main factor ensuring the separation of a single harmonic of the perturbation and leads to the appearance of an ordered periodic structure from initially chaotic perturbations.

5. Conclusions

We have suggested the physical model that explains the laser-pulse-induced nanograting formation within the volume of a solid dielectric as a result of ionization-field instability of the optical discharge in respect to excitation of the “internal surface plasmon” generated at the narrow leading edge of the breakdown wave counter-propagated to the incident pulse inside the dielectric. The instability we consider, like other ionization-field instabilities, is due to the effects of mutual amplification of small perturbations in the field amplitude and the plasma density. The mechanism of its occurrence can be briefly explained as follows. At the leading edge of the breakdown wave counter-propagating with the incident beam, a thin plasma layer is formed in which the effective permittivity is negative. A surface wave (surface plasmon) with a wavelength shorter than the incident wavelength can propagate along such a layer. This plasmon is not excited by the incident electromagnetic wave as long as the layer remains uniform in the longitudinal direction. However, with the appearance of an infinitesimal periodic perturbation of the plasma density along the layer with a spatial period equal to the wavelength of the plasmon itself, it is effectively excited (in the form of a standing surface wave). The maxima of the amplitude of the total electric field (and, consequently, the maxima of the ionization rate) arising from the interference between the fields of the plasmon and the incident wave, that excited it, are located in the same places where the maxima of the density perturbation are located. Due to an increase in the ionization rate in these places, the density maxima will continue to grow, and thus, a positive feedback loop arises, indicating the appearance of instability. The role of the initial (seed) perturbation is played by the corresponding spectral component, which is always present in the continuous spectrum of weak random fluctuations of parameters (in our case, plasma density), which are always present in any real system. We have carried out numerical simulation based on the equation system including the rate equation for the plasma density and wave equation describing the spatiotemporal evolution of the average field in nanoporous fused silica that we have considered as a heterogeneous medium with effective dielectric permittivity calculated using the Maxwell Garnett formulas. We have found that a thin plasma layer with negative (and close to zero) effective permittivity was formed in the head part of the breakdown wave counter-propagated to the incident laser pulse. Surface plasmon have excited at this layer resulting in strong modulation of plasma density along it (in the direction of laser polarization), strong corrugation of the breakdown wave front and formation of contrast periodic structure that can be considered probably as a root cause of subwavelength nanograting formation in the dielectric material. The spatial period of the formed periodical nanostructure (coinciding according to our qualitative model with the wave length of the internal surface plasmon) did not depend on the structure of random small initial perturbations playing the role of a small seed which is necessary to any instability. Within the parameter range used in our calculations this period was about 0.7 of the laser wavelength in unperturbed dielectric which fits within the framework of experimentally measured values. The experimentally observed dependences of the nanograting period on laser pulses parameters could be explained probably only based on a more complicated model (in particular, out of Maxwell Garnett approximations and with thermo-chemical processes between pulses taken into account).