# Transverse Ultrafast Laser Inscription in Bulk Silicon

M. Chambonneau,1,∗ M. Blothe,1 Q. Li,1 V. Yu. Fedorov,2,3 T. Heuermann,1,4 M.
Gebhardt,1,4 C. Gaida,5 S. Tertelmann,6 F. Sotier,6 J. Limpert,1,4,7 S.
Tzortzakis,2,8,9 S. Nolte1,7 1Institute of Applied Physics, Abbe Center of
Photonics, Friedrich-Schiller-University Jena, Albert-Einstein-Straße 15,
07745 Jena, Germany.
2Science Program, Texas A$\&$M University at Qatar, P.O. Box 23874, Doha,
Qatar.
3P. N. Lebedev Physical Institute of the Russian Academy of Sciences, 53
Leninskiy Prospekt, 119991 Moscow, Russia.
4Helmholtz-Institute Jena, Fröbelstieg 3, 07743 Jena, Germany.
5Active Fiber Systems GmbH, Ernst-Ruska-Ring 17, 07745 Jena, Germany.
6Innolas Photonics GmbH, Justus-von-Liebig-Ring 8, 82152 Krailling, Germany.
7Fraunhofer Institute for Applied Optics and Precision Engineering IOF, Center
of Excellence in Photonics, Albert-Einstein-Straße 7, 07745 Jena, Germany.
8Institute of Electronic Structure and Laser (IESL), Foundation for Research
and Technology—Hellas (FORTH), P.O. Box 1527, GR-71110 Heraklion, Greece.
9Materials Science and Technology Department, University of Crete, 71003
Heraklion, Greece.
<EMAIL_ADDRESS>
###### Abstract

In-volume ultrafast laser direct writing of silicon is generally limited by
strong nonlinear propagation effects preventing the initiation of
modifications. By employing a triple-optimization procedure in the spectral,
temporal and spatial domains, we demonstrate that modifications can be
repeatably produced inside silicon. Our approach relies on irradiation at
$\approx 2$-µm wavelength with temporally-distorted femtosecond pulses. These
pulses are focused in a way that spherical aberrations of different origins
counterbalance, as predicted by point spread function analyses and in good
agreement with nonlinear propagation simulations. We also establish the laws
governing modification growth on a pulse-to-pulse basis, which allows us to
demonstrate transverse inscription inside silicon with various line
morphologies depending on the irradiation conditions. We finally show that the
production of single-pulse repeatable modifications is a necessary condition
for reliable transverse inscription inside silicon.

††preprint: APS/123-QED

##  I Introduction

Ultrafast laser direct writing is a proven technique for in-volume
functionalization of dielectrics in a three-dimensional (3D), fast and
contactless way [1]. Despite its high desirability for applications including
photonics, sensors, photovoltaics, ultrafast microelectronics and quantum
computing, this technique has no equivalent in monolithic silicon [2]. This
lack originates from the formation of filaments during the propagation of
ultrashort laser pulses in this highly nonlinear material [3, 4, 5, 6, 7, 8].
As a consequence of these nonlinear propagation effects in silicon, the energy
deposition is strongly delocalized to the prefocal region at levels below the
breakdown threshold—thus, representing an uttermost challenge for the
production of repeatable modifications which are the basic brick for 3D laser
direct writing. While different approaches have been devised for writing lines
longitudinally (i.e., along the optical axis) in silicon [9, 10, 11, 12, 13,
14, 15], demonstrations of transverse line inscription (i.e., perpendicular to
the optical axis) are seldom, especially in the femtosecond regime. Two
strategies were adopted so far for solving this issue. The first one consists
of using longer pulses (i.e., in the picosecond or nanosecond regime) to
decrease the aforementioned nonlinearities [16, 17, 18, 19, 20, 21, 22]. The
main disadvantage of this approach lies in its thermal nature which limits the
degree of control that can be reached on the inscribed structures. In order to
reduce the heat-affected zone, a second strategy consists of employing
femtosecond pulses with an altered entrance surface, either by depositing an
oxide layer [23], or by employing oil immersion [24]. Nevertheless, in
addition to the impractical aspect of this strategy, surface alteration is
undesirable—if not impossible—in numerous applications.

In this article, we demonstrate contactless femtosecond laser transverse
inscription inside silicon employing a triple-optimization procedure which
combines beneficial effects in the spectral, temporal and spatial domains. The
spectral optimization consists of employing femtosecond laser pulses at
$\approx 2$-µm wavelength, where the peak in the nonlinear refractive index at
this specific wavelength gives rise to enhanced energy deposition [5, 25, 22].
The optimization in the time domain relies on the utilization of temporally-
distorted 750-fs duration pulses. Temporal imperfections were recently shown
to act as a seed inside silicon allowing the material breakdown [26, 27]. We
propose an additional optimization in the spatial domain which consists of
counterbalancing the spherical aberration induced by the focusing lens with
the one provoked by the refractive index mismatch at the air-silicon
interface, as supported by point spread function calculations. This triple-
optimization approach allows us to demonstrate the initiation of repeatable
permanent internal modifications with ultrashort laser pulses, in good
agreement with nonlinear propagation simulations. The laws governing
modification growth after successive irradiations at the same location are
then established. This comprehension on the breakdown behavior leads us to
devise transverse femtosecond laser inscription of continuous lines in the
bulk of silicon. Three distinct line morphologies which depend on the writing
speed and the input pulse energy are exhibited. Ultimately, we demonstrate
that continuous lines can be solely written in the conditions where the
single-pulse modification probability is 100$\%$—thus, shedding light on the
necessity for breakdown reignition on a pulse-to-pulse basis during the
inscription.

##  II Experimental arrangement

The experimental arrangement employed for inducing, detecting and
characterizing permanent modifications inside silicon is depicted in Fig.
1(a). It relies on 750-fs duration pulses [full width at half maximum,
assuming a squared hyperbolic secant (sech2) temporal shape] at 1.95-µm
wavelength emitted at $\Omega_{0}=100$-Hz repetition rate by a Tm-doped fiber
laser source. The input pulse energy $E_{\text{in}}$ is controlled by means of
a half-wave plate and a linear polarizer. This energy is measured in air,
between the focusing optics and the sample (i.e., the Fresnel reflections are
not taken into account), and it may reach $E_{\text{in}}\approx 2.55$ µJ. The
beam size before focusing is adjusted by means of a Galilean telescope so that
the beam diameter is larger than the diameter of the aspherical lens
(Thorlabs, C037TME-D, numerical aperture $\text{NA}=0.85$) employed for
focusing the beam inside silicon. This lens is mounted on a stage allowing its
displacement along the optical axis $z$. It is important to note that, while,
aspheres usually allow eliminating spherical aberrations, the one we have
selected is designed for operation at 9.5-µm wavelength, thus inducing
aberrations at our laser wavelength in addition to those induced by the
refractive mismatch at the air-silicon interface (this will be discussed below
in Sec. III.1). The silicon samples (Siegert Wafer, double-side polished, 1-mm
thickness, $<100>$-oriented, undoped, 200-$\cdot$cm resistivity) can be moved
in the $xy$ plane perpendicular to the optical axis $z$ thanks to a
positioning system (Aerotech, ANT 130).

Figure 1:  (a) Schematic of the arrangement employed for internal breakdown
and transverse line writing in silicon and associated in-situ diagnostics. The
vector $\vec{k}$ indicates the direction of laser propagation. The second-
order autocorrelation traces in (b) and (c) have been obtained with the two
distinct 2-µm laser sources utilized in our study.

Three-dimensional in-situ monitoring of the laser-silicon interaction is
performed by means of two customized diagnostics. The first one is a
transmission microscope relying on white light illumination. The light passes
through the sample along the $x$ axis, and images of the modifications in the
$yz$ plane are recorded thanks to an objective lens (Mitutoyo, M Plan Apo NIR-
HR 50$\times$), a tube lens, and an InGaAs camera (Xenics, Bobcat 320). Images
in the $xy$ plane are recorded by a transmission microscope relying on
illumination at 1.03-µm (Spectra-Physics, femtoTrain) along the $-z$ direction
passing through the silicon sample, the focusing asphere, and a tube lens
before being collected by a silicon-based camera. This latter microscope is
utilized for determining the relative position of the focus with respect to
the sample surface by performing surface damage scans at energies near the
breakdown threshold. The focusing depth is obtained by multiplying the
displacement $\Delta z$ of the focusing lens along the optical axis starting
from the entrance surface by the refractive index of silicon at 1.95-µm
($n_{0}\approx 3.45$ [28]). For improved imaging performance, post-mortem
characterizations of the modifications are carried out with a transmission
infrared phase microscope consisting of a customized Mach-Zehnder
interferometer operated at 1.3-µm wavelength whose principle is detailed in
Ref. [29].

Two distinct laser sources have been used in the present study. Both of these
sources exhibit similar spectral and spatial features. However, the temporal
shapes of the pulses emitted by these laser systems are different. The first
laser system (Active Fiber Systems, see details in Ref. [30]) delivers
undistorted bell-shaped pulses as illustrated in Fig. 1(b). In stark contrast,
as shown by the second-order autocorrelation trace in Fig. 1(c), the
femtosecond pulses emitted by the second laser source (Innolas Photonics,
FEMTO*1950-8-T-2500) exhibit a structured temporal profile. This distortion
originates from uncompensated higher-order dispersion resulting in the
existence of satellite peaks. However, one has to emphasize that second-order
autocorrelation does not allow discriminating if the satellite peaks are pre-
or post-pulses.

##  III Results and discussion

###  III.1  Linear propagation calculations

Let us first examine the focusing depth allowing to maximize the energy
deposition inside silicon, i.e., the spatial optimization. To do so, we have
first carried out linear propagation calculations in silicon with our recently
developed in-house vectorial model called “InFocus” [31, 32]. Briefly, this
model gives access to the linear propagation in any medium and relies on point
spread function analysis, allowing the determination of the volumetric field
distribution of light in the focal region. One specificity of this model is
that it accounts not only for the spherical aberration provoked by the
refractive index mismatch, but also for the one induced by the focusing lens
through its associated Zernike polynomials. While this latter spherical
aberration may be neglected for a beam focused with an objective lens, it
becomes significant for thick singlet lenses with short focal length. The
accuracy of our model to predict these two types of spherical aberrations
separately was benchmarked with propagation imaging [31]. While, in principle,
similar measurements could be performed in the present study, these would
require an objective lens with a numerical aperture higher than the one of the
focusing asphere ($\text{NA}=0.85$), which represents a major technical
challenge. In the present study, we thus restrict ourselves to linear
propagation calculations considering the aspherical lens used in the
experiments. As mentioned in Sec. II, because of the difference between the
design wavelength of this asphere (9.5 µm) and the laser wavelength (1.95 µm),
spherical aberration is inevitably induced by the focusing lens.

For simplicity, let us first consider that the propagation is linear, i.e.,
nonlinear effects such as the optical Kerr effect as well as plasma formation
and defocusing are neglected. Under this assumption, the results are
independent of the temporal profile of the pulses—and, thus, valid for the two
employed laser systems [see Figs. 1(b) and (c)]. Figure 2(a) shows the
intensity distributions for various focusing depths in silicon. For short
focusing depths ($\leq 80$ µm), the asymmetry in the intensity distributions
is ascribable to the predominance of lens-induced spherical aberration.
Conversely, for long focusing depths ($\geq 140$ µm), elongated patterns along
the propagation direction indicate the predominance of spherical aberration
provoked by the refractive index mismatch at the air-silicon interface.
Interestingly, for intermediate depths ($\approx 100$ µm), both types of
spherical aberrations compensate each other, leading to localized near-
symmetric intensity distributions. The maximum intensity reached in each
configuration is displayed in Fig. 2(b) and shows a maximum value for 100-µm
focusing depth, where both types of spherical aberrations counterbalance. In
contrast, for focusing depths $\leq 80$ µm and $\geq 120$ µm, the maximum
intensity value reached is significantly lower than at $\approx 100$-µm
focusing depth. This leads us to the conclusion that the spatial optimization
consisting of the compensation of the two types of spherical aberrations is
beneficial for provoking breakdown inside silicon. However, one should
emphasize that this condition on counterbalancing spherical aberrations taken
alone cannot fully explain breakdown experiments as the hypothesis of a linear
propagation is unrealistic given our experimental conditions. Indeed, the
linear and nonlinear refractive indices of silicon at $\lambda\approx 2$-µm
wavelength are $n_{0}=3.45$ [28], and $n_{2}\approx 9.4\times 10^{-18}$
m${}^{2}\cdot$W-1 [33, 34, 35]. The corresponding critical power for self-
focusing $P_{\text{cr}}\approx 1.8962\lambda^{2}/\left(4\pi
n_{0}n_{2}\right)\approx 17.7$ kW is much lower than the experimental peak
powers ($P\geq 1.7$ MW). Therefore, while spherical aberrations exist
experimentally as a linear effect, rigorously, a more sophisticated model
accounting for nonlinear propagation is required for determining the optimal
focusing depth.

Figure 2:  (a) Numerical calculations of the on-axis intensity distribution at
various focusing depths thanks to our linear propagation model [32, 31]. The
vector $\vec{k}$ indicates the direction of laser propagation. The intensity
distributions $I/I_{\text{max}}$ are normalized to the maximum intensity of
the whole set of data. The spatial scale applies to all images. (b) Evolution
of the normalized maximum intensity for each configuration as a function of
the focusing depth.

###  III.2 Nonlinear propagation simulations

In order to consider the nonlinear propagation inside silicon, we performed
numerical simulations based on the unidirectional pulse propagation equation
(UPPE) [36, 37, 6]. As the initial condition at the entrance of the silicon
sample, we took a pulse with the spatial profile calculated by our vectorial
model, which allowed us to include the effect of aberrations into the
simulations. The considered temporal shape corresponds to the autocorrelation
trace in Fig. 1(b). The input pulse energy was chosen to be 1 µJ. Together
with the UPPE, we solved the following equation for the electron density
$\rho$:

$\displaystyle\frac{\partial\rho}{\partial
t}=\sigma_{2}I^{2}(\rho_{\text{nt}}-\rho)+\sigma_{\text{B}}\frac{\mathcal{E}^{2}}{U_{\text{cr}}}\frac{\rho_{\text{nt}}-\rho}{\rho_{\text{nt}}}\rho-\gamma\rho^{3}-\frac{\rho}{\tau_{\text{r}}},$
(1)

where the first term on the right-hand side describes two-photon ionization,
second — avalanche ionization, third — Auger recombination, and fourth —
electron relaxation. Here $\rho_{\text{nt}}=5\times 10^{28}$ m-3 is the
density of neutral atoms,
$\sigma_{2}=\beta_{2}/(2\hbar\omega_{0}\rho_{\text{nt}})$ is the cross-section
of two-photon ionization with the two-photon absorption coefficient
$\beta_{2}=2.97\times 10^{-12}$ m$\cdot$W-1 [33, 34], and $I$ is the pulse
intensity;
$\sigma_{\text{B}}=(e^{2}/m_{e}^{*})\nu_{c}/(\nu_{c}^{2}+\omega_{0}^{2})$ is
the inverse Bremsstrahlung cross section, where $e$, $m_{e}$, and
$m_{e}^{*}=0.5m_{e}$ are the charge, mass, and effective mass of electron,
$\nu_{c}=0.3\times 10^{15}$ s-1 is the collisions frequency [38], $\omega_{0}$
is the central pulse frequency, and $\mathcal{E}$ is the pulse electric field,
while $U_{\text{cr}}=1.5U_{\text{i}}$ with $U_{\text{i}}=1.12$ eV being the
band-gap energy (the coefficient 1.5 accounts for the momentum conservation);
$\gamma=1.1\times 10^{-42}$ m${}^{6}\cdot$s-1 is the Auger recombination
coefficient [39], and finally $\tau_{r}=100$ fs is the electron relaxation
time which we took from Ref. [40] assuming that, due to collisions, the
maximum energy of free electrons cannot be much larger than the band gap.
Ultimately, free-carrier diffusion is neglected as this mechanism is not
expected to play a major role at the time scale of ultrashort pulses [41].

Figure 3 summarizes the results of our simulations. Figure 3(a) shows the
dependence of the peak fluence on the propagation distance for different
focusing depths. Here, one can note that the maximum fluence values are almost
independent of the focusing depth. Therefore, based on the fluence data, we
cannot conclude that some focal depths are preferable to others. In turn, Fig.
3(b) shows the dependence of peak intensity $I_{\text{peak}}$ and peak
electron density $\rho_{\text{peak}}$ on the focusing depth. The values of
$I_{\text{peak}}$ and $\rho_{\text{peak}}$ reach their respective maxima at
100-µm focusing depth, in accordance with the predictions made by our
vectorial model (see Fig. 2). However, since the intensity and electron
density reach their peak values at some arbitrary point in space (across the
beam) and time, one cannot fully rely on this data. Moreover, it is known that
criteria based on the peak electron density usually do not provide the most
accurate description of breakdown experiments [42]. Therefore, as a more
reliable breakdown criterion, we decided to calculate the laser energy
delivered to the free-electron subsystem, which then can be transferred to the
crystal lattice [43] (see the details in Appendix A). The peak values of the
delivered energy (estimated by neglecting the surface damage at depths below
40 µm) are displayed in Fig. 3(c) as a function of the focusing depth. It can
be seen that, for 100-µm focusing depth, the peak delivered energy reaches its
maximum, indicating that the breakdown probability for this focusing depth is
the highest.

Figure 3:  Results of numerical simulations of nonlinear propagation in
silicon. (a) On-axis fluence profiles at different focusing depths. (b)
Evolution of the peak intensity (green) and peak electron density (blue) as a
function of the focusing depth. The red region indicates the focusing depths
for which surface damage was experimentally detected. (c) Evolution of the
peak delivered energy as a function of the focusing depth (black circles). The
bars indicate the percentage of the input laser energy spent on ionization
(blue) and acceleration (i.e., heating) of electrons (orange, with the
indicated percent value).

Additionally, since, in our nonlinear propagation model, the laser energy is
only spent on ionization and electron acceleration (heating), by calculating
the total laser energy loss and the total energy delivered to free electrons,
we can separately find the percent of the initial laser energy spent on
ionization and electron heating [see the color bars in Fig. 3(c)]. In
particular, we find that, at 100-µm focusing depth, where the peak delivered
energy reaches its maximum, the percent of laser energy spent on acceleration
of electrons is also the highest and exceeds 90$\%$.

To sum up, both the linear and the nonlinear propagation simulations indicate
an optimal focusing depth near 100 µm, which could thus potentially lead to
continuous transverse line writing inside silicon. In order to verify
experimentally these theoretical predictions, we have first carried out
breakdown experiments in the bulk of silicon using the laser source emitting
temporally-undistorted pulses [see Fig. 1(b)]. These experiments consist of
single-pulse irradiation, and the corresponding binary result (i.e., breakdown
or not) is evaluated by means of the three-dimensional in-situ monitoring
system described in Sec. II. By repeating the experiments on ten different
fresh sites, the bulk modification probability (i.e., the ratio between the
number of breakdown sites and the total tested sites) is evaluated. The
modification probability lies in the 0–30$\%$ range. In other words, no
repeatable modification could be produced in the bulk of silicon. This
observation held even by changing the input pulse energy as well as the
focusing depth. As illustrated by the endeavored transverse line inscription
in Fig. 4 where modifications are randomly produced during writing, transverse
inscription of continuous lines in this stochastic regime is doomed to
failure. Both the size and the depth of the modifications are unrepeatable.
This random behavior for breakdown initiation cannot originate from
fluctuations of the laser power due to the high stability ($\approx 1\%$) of
the employed system. It is thus more likely related to the existence of
precursor defects inside the material [44]. One may note the limited spatial
expansion of the modification along the writing direction (indicated by the
green arrow) in Fig. 4. The modifications can be propagated over $\approx 100$
µm before no more expansion occurs.

Figure 4:  Transmission ($T/T_{0}$) image of transverse line inscription
attempt with pulses corresponding to the autocorrelation trace in Fig. 1(b) at
the writing speed $v=20$ µm$\cdot$s-1. The vector $\vec{k}$ and the green
arrow indicate the direction of laser propagation and the writing direction,
respectively.

Transverse inscription of continuous lines could not be achieved despite the
benefits expected from irradiations at 2-µm wavelength where the deposited
energy is maximized (spectral optimization), and at 100-µm focusing depth
where spherical aberrations of two different origins counterbalance (spatial
optimization). These negative results lead us to add a third optimization in
the temporal domain. In the nanosecond regime, several studies have shown that
the temporal pulse shape strongly influences the initiation [45, 46, 47, 48,
49] and the morphology [50, 51] of laser-produced modifications in fused
silica. Recent results showed that the temporal pulse shape is also a key
parameter in femtosecond-laser-silicon interaction as perfect Gaussian
temporal profile is not suitable for breakdown initiation inside silicon [26].
Hereafter, the temporal optimization used in the present study consists of
employing temporally-distorted pulses associated with the autocorrelation
trace in Fig. 1(c).

###  III.3  Single-pulse breakdown initiation

In order to investigate the possibility offered by our triple-optimization
procedure for inducing internal modifications in silicon, single-pulse
breakdown initiation experiments have been first realized at different
focusing depths and input pulse energies. The corresponding results are shown
in Fig. 5 in terms of bulk breakdown initiation probability as a function of
the focusing depth for various input pulse energies $E_{\text{in}}$. Surface
damage was systematically detected for focusing depths $\leq 40$ µm (not
represented in Fig. 5). No bulk breakdown is observed for all focusing depths
at $E_{\text{in}}=1.5$ µJ, thus defining the modification threshold. The
striking feature in Fig. 5 is the existence of a narrow window (indicated in
pink) for 100–120-µm focusing depth where the breakdown probability shows a
maximum value for $E_{\text{in}}\geq 1.7$ µJ. This result is in excellent
agreement with the theoretical predictions in Secs. III.1 and III.2.
Therefore, this narrow window indicates the focusing depth for which spherical
aberrations induced by the lens and the air-silicon interface counterbalance.
Interestingly, bulk breakdown probability of 100$\%$ is reached for
$E_{\text{in}}\geq 1.9$ µJ. This deterministic character suggests that these
conditions are well-adapted for ultrafast laser direct writing inside silicon,
in opposition to the previous results obtained with temporally-undistorted
pulses. Apart from this narrow window, nonzero modification probabilities
below 100$\%$ are measured, suggesting non-intrinsic breakdown likely
originating from the existence of local precursor defects inside the material
[44]—thus incompatible with the repeatable inscription of functions inside
silicon (see Fig. 4).

Figure 5:  Evolution of the bulk breakdown probability inside silicon with
single laser pulses as a function of the focusing depth for various input
pulse energies. The pink region indicates the depths for which the
modification probability reaches 100$\%$ for an input pulse energy
$E_{\text{in}}\geq 1.9$ µJ.

Besides the statistical aspects related to ultrafast laser-induced breakdown
inside silicon, we have investigated the associated material changes. To do
so, the morphology of single-pulse modifications has been characterized with
different methods. First, in-situ transverse observations reveal that the
modifications consist of a few micro-bubbles along the optical axis $z$, as
shown in Fig. 6(a). The number of individual bubbles varies from one to four
in the tested input energy range, and their typical size is $\approx 4$ µm.
Post-mortem phase measurements have also been carried out in the $xy$ plane,
as shown in Fig. 6(b). The phase shift—and, thus, the refractive index
change—associated with the laser-produced modifications is positive. Two main
types of material change may explain such an observation. First, the positive
refractive index change may originate from local laser-induced strain inside
the material, analogously to breakdown produced in the nanosecond and
picosecond regimes [52, 19, 12, 13]. Second, the modifications may consist of
an allotropic change of silicon, from monocrystalline to polycrystalline [53],
or other polymorphs [54].

Figure 6:  Morphology of a single-pulse modification produced at
$E_{\text{in}}=2.1$-µJ. (a) In-situ transmission microscopy image of the
modification in the $yz$ plane. (b) Post-mortem phase image of the
modification in the $xy$ plane. The phase shift $\Delta\varphi$ is mapped out
thanks to a four-step procedure described in Ref. [29]. The vector $\vec{k}$
indicates the direction of laser propagation.

###  III.4  Multi-pulse modification growth

Establishing the laws governing modification growth (i.e., the increase in the
modification size on a pulse-to-pulse basis) in silicon is a prerequisite for
envisioning multi-pulse functionalization processes. While, previously, we
have determined the conditions for which spherical aberrations induced by the
lens and the air-silicon interface counterbalance, leading to repeatable
breakdown (see Fig. 5), experiments on the evolution of the modification size
after successive irradiations have been carried out. Such a growth sequence is
illustrated by transmission microscopy images in the $xz$ plane after
different numbers of applied pulses in Fig. 7(a). The first pulse focused at
120-µm depth initiates a modification similar to the one in Fig. 6. This
barely detectable single-shot modification expands spatially after several
irradiations to form a strongly absorbing damage site. The modification
preferentially grows toward the entrance surface of the sample, and only
marginally after the focus. This behavior is ascribable to enhanced absorption
in the modified zone. One may note the formation of a tail growing rapidly
toward the entrance surface for a number of pulses $N\geq 30$, as indicated by
white arrows in Fig. 7(a). While deeper investigations out of the scope of the
present study are necessary for establishing the physical origin of this tail,
it must be emphasized that the expansion of this tail does not provoke any
detectable surface damage. Quantitative analyses of the evolution of the
modification lengths $L_{x}$ and $L_{z}$ indicated in Fig. 7(a) along the $x$
and the $z$ axis, respectively, are shown in Fig. 7(b) as a function of the
number of pulses $N$. For all $N$ values, $L_{z}>L_{x}$, meaning that the
modifications are more elongated along the optical axis than in the transverse
direction. Both $L_{x}$ and $L_{z}$ scale logarithmically with the number of
pulses $N$. Remarkably, this evolution for bulk modifications is slower than
for surface modifications in fused silica, which usually exhibit linear and
exponential growth on the entrance and the exit surface, respectively [55, 56,
57].

Figure 7:  Laser-induced modification growth sequence inside silicon for
$E_{\text{in}}=2.1$ µJ. (a) Transmission images of a modification for
different numbers of irradiations at the same location (120-µm below the
entrance surface of the sample). The transmission $T/T_{0}$ is obtained by
dividing a transmission optical image after $N$ pulses at the same location by
a reference image. The spatial scale applies to all images. The vector
$\vec{k}$ indicates the direction of laser propagation. The dotted yellow
lines indicate the sample entrance surface. (b) Evolution of the modification
length $L_{x}$ and $L_{z}$ indicated in (a) along the $x$ and the $z$ axis,
respectively, as a function of the number of applied pulses.

Additional growth experiments have been carried out for different input pulse
energies. The corresponding results are displayed in Fig. 8. The evolution of
the modification lengths $L_{x,y}$ and $L_{z}$ are shown in Figs. 8(a) and
8(b), respectively, as a function of the number of pulses $N$. Similar to the
trends in Fig. 7(b), each data set is fairly well fitted by a logarithmic
function as

$L_{i}=A_{i}\text{ln}(N)+B_{i},$ (2)

where $i=x$, $y$ or $z$, and $A_{i}$ and $B_{i}$ are the growth coefficient
and the single-pulse modification length along $i$, respectively. Given that
the characterizations in Fig. 6 have revealed symmetric modifications in the
$xy$ plane, one could intuitively expect a conservation of this symmetry in
the multi-pulse configuration—and, thus, $A_{x}=A_{y}$. As shown in Fig. 8(a),
this is not the case, and the modifications are more elongated along $x$ than
along $y$. Consecutive irradiations thus produce elliptical modifications
whose major and minor axes are parallel and perpendicular to the laser
polarization, respectively, as illustrated by the inset in Fig. 8(a) where the
modification is 60$\%$ smaller along $y$ than along $x$. Interestingly,
elliptical modifications were also observed in silicon for multi-pulse
picosecond irradiations in the bulk [58, 26, 25, 59], as well as for
nanosecond irradiations on the exit surface [60]. However, the physical
explanation for the origin of this elliptical morphology remains open and
requires a dedicated study.

Figure 8:  Laws of modification growth inside silicon. (a) Evolution of the
transverse modification dimensions $L_{x}$ (circles) and $L_{y}$ (triangles)
along the $x$ and the $y$ axis, respectively, as a function of the number of
applied pulses $N$ for different input pulse energies. The inset is an optical
microscopy observation of a modification produced after 1,000 irradiations at
$E_{\text{in}}=2.4$ µJ. The vectors $\vec{k}$ and $\vec{E}$ indicate the
direction of laser propagation and the polarisation, respectively. (b)
Evolution of the modification dimension $L_{z}$ along the optical axis $z$ for
the same input pulse energies as in (a). In both (a) and (b), the curves with
indicated equations are logarithmic fits of the experimental data according to
Eq. (2). The values of the coefficients $A_{i}$ and $B_{i}$ ($i=x,y$ or $z$)
are reported in Table 1.

The relatively small dispersion in the $L_{x,y,z}$ values in Fig. 8 for
different input pulse energies $E_{\text{in}}$ indicate that the laws of
modification growth described by Eq. (2) do not depend much on this parameter
in the tested input pulse energy range ($E_{\text{in}}=1.7$–2.4 µJ). This is
in agreement with the saturation of the laser intensity during the propagation
of ultrashort pulses in silicon [4, 5]. A statistical analysis of $A_{i}$ and
$B_{i}$ coefficients ($i=x$, $y$ or $z$) appearing in Eq. (2) is recapitulated
in Table 1. Consistent with Figs. 7 and 8, $A_{z}>A_{x,y}$ and
$B_{z}>B_{x,y}$, showing again that the modifications grow preferentially
along $z$.

Table 1:  Statistical analysis of the $A_{i}$ and $B_{i}$ coefficients ($i=x$, $y$ or $z$) appearing in the laws of modification growth described in Eq. (2). The data have been obtained by analyzing 288 modifications initiated at 120-µm focusing depth by 1, 10, 100 and 1,000 laser pulses with $E_{\text{in}}=1.7$–2.4 µJ. |  $A_{x}$ (µm) |  $B_{x}$ (µm) |  $A_{y}$ (µm) |  $B_{y}$ (µm) |  $A_{z}$ (µm) |  $B_{z}$ (µm)
---|---|---|---|---|---|---
Average | 2.34 | 1.81 | 1.31 | 2.03 | 7.18 | 8.91
Standard deviation | 0.12 | 0.07 | 0.08 | 0.13 | 0.31 | 1.43
Minimum | 2.12 | 1.71 | 1.14 | 1.86 | 6.82 | 6.49
Maximum | 2.48 | 1.90 | 1.38 | 2.25 | 7.66 | 10.54

###  III.5  Transverse inscription

By constantly moving the sample perpendicularly to the laser beam at various
writing speeds $v$, transverse inscription has been endeavored inside silicon
at 120-µm focusing depth and $E_{\text{in}}=2.5$ µJ. As shown by the overview
in Fig. 9(a), our triple-optimization procedure allows the first demonstration
of contactless transverse inscription of lines in bulk silicon with
femtosecond laser pulses. Finer characterizations of the line morphology in
different imaging planes are displayed in Figs. 9(b) and 9(c), revealing three
distinct morphologies. The first one appearing for $v=400$ µm$\cdot$s-1 is the
discrete morphology, consisting of a track of individual single-pulse
modifications. This morphology associated with high writing speeds is in
excellent agreement with the prediction of the phenomenological model detailed
in Appendix B, which allows the determination of the threshold minimum speed
value for discrete line inscription $v_{\text{min}}\approx 262$ µm$\cdot$s-1.
In stark contrast, for $v\leq v_{\text{min}}$ the lines are continuous. In
particular, for $40\leq v\leq 200$ µm$\cdot$s-1, the line morphology is
steady, i.e., with a nearly constant track width. Finally, for low speeds
($v\leq 20$ µm$\cdot$s-1), bulgy clusters are exhibited, and the line
morphology is erratic (i.e., the track width is not constant). This complex
morphology likely originates from cumulative effects as suggested in Fig. 9(c)
for $v=20$ µm$\cdot$s-1. In this case, the high number of pulses per point
($N\geq 10$) associated with the slow inscription process gives rise to
modification growth, provoking shielding and scattering effects at the focus
during the subsequent irradiations. As a consequence, the energy deposition
during the subsequent irradiations is drastically reduced so that no new
modification is initiated at the focus until the growth-induced shielding and
scattering effects stop. Finally, the evolution of the line width in the $xy$
plane is displayed in Fig. 9(d) as a function of the writing speed. The line
width is the same for $v\geq 100$ µm$\cdot$s-1 (i.e., $N\leq 2$), and
corresponds to the dimensions of single-pulse modifications (see Fig. 6). This
demonstrates the absence of cumulative effects for high writing speeds. In
contrast, for $v<100$ µm$\cdot$s-1, the line width increases logarithmically
as the writing speed decreases, thus defining the threshold speed for the
aforementioned cumulative effects.

Figure 9:  Transverse ultrafast-laser-inscribed lines inside silicon. (a)
Overview of 4-mm long laser-written lines at 120-µm focusing depth. (b) Post-
mortem transmission microscopy images in the $yz$ plane of lines inscribed at
$E_{\text{in}}=2.5$ µJ for different writing speeds $v$. (c) In-situ
transmission microscopy images in the $xy$ plane of lines inscribed at
$E_{\text{in}}=2.5$ µJ for different writing speeds $v$. In each sub-figure,
the spatial scale applies to all images, the vector $\vec{k}$ indicates the
direction of laser propagation, and the green arrow shows the writing
direction. (d) Evolution of the line width in the $xy$ plane as a function of
$v$ for $E_{\text{in}}=2.5$ µJ. The different colors correspond to the
indicated morphology domains.

Ultimately, we have studied the energy dependence of the three different
morphologies exhibited in Fig. 9 (i.e., discrete, steady and erratic) by
characterizing transversely-inscribed lines for various input pulse energies
and writing speeds. The corresponding results are displayed in Fig. 10. As
discussed previously, high input pulse energies lead to discrete and erratic
morphology for excessively high and low speeds, respectively. Interestingly,
for intermediate speeds, the steady morphology solely exists for
$E_{\text{in}}>1.8$ µJ, corresponding to a single-pulse breakdown probability
of 100$\%$ (see the red curve in Fig. 10). Conversely, the line morphology is
discrete for $v=40$–100 µm$\cdot$s-1 and $E_{\text{in}}\approx 1.7$ µJ, which
corresponds to a nonzero modification probability strictly lower than $100\%$
(see Fig. 5). Given that modification growth occurs at $E_{\text{in}}\approx
1.7$ µJ (see Fig. 8), one can conclude that continuous line inscription is not
related to the reignition of modifications enabling their transverse
expansion, but to intrinsic material breakdown. In other words, a single-pulse
modification probability of $100\%$ is a necessary condition for the
transverse inscription of continuous lines (either steady or erratic). This
conclusion is all the more supported by the observations in Fig. 4.

Figure 10:  Morphology domains of transversely-inscribed lines as a function
of the input pulse energy and the writing speed. The focusing depth is 120 µm.
The black points indicate the conditions for which line inscription has been
endeavored, and the color zones correspond to the different morphology
domains. The red curve corresponds to a sigmoid fit of the single-pulse
breakdown probability data in Fig. 5.

##  IV Summary

We have demonstrated the possibility to contactlessly inscribe continuous
transverse lines in bulk silicon with ultrafast laser pulses. While this
achievement is usually unattainable due to nonlinear propagation effects, we
have proposed a triple-optimization procedure for circumventing these
limitations. This approach consists of using femtosecond pulses (i) at the
beneficial wavelength of $\approx 2$ µm, (ii) exhibiting temporal distortions,
and (iii) focused at a depth for which spherical aberrations of different
origins counterbalance. The concomitance of these spectral, temporal and
spatial optimizations was shown to be indispensable for producing repeatable
breakdown—and, thus, transverse inscription—inside silicon.

Our methodology first relies on the identification of the conditions that
allow obtaining $100\%$ breakdown probability with single-pulse irradiations.
Deterministic breakdown was observed for input pulse energies
$E_{\text{in}}\geq 1.9$ µJ and focusing depths $\approx 100$–120 µm below the
entrance surface. This narrow focusing depth window is in good agreement with
point spread function analyses as well as nonlinear propagation simulations.
In addition, we have studied the spatial expansion of the modifications on a
pulse-to-pulse basis. By applying different numbers of pulses for various
input pulse energies, we have determined the laws governing modification
growth. In all spatial directions, the growth behavior was found to scale
logarithmically with the number of pulses applied, with no strong dependence
on the input pulse energy. Finally, lines have been transversely inscribed in
the bulk of silicon. Depending on the writing speed and the input pulse
energy, the transversely-written lines may show an erratic, steady or discrete
morphology. The discrete morphology appears either for input pulse energies
leading to single-pulse breakdown probabilities $<100\%$, or for too high
writing speeds.

The demonstration of transversely-inscribed lines paves the way to contactless
3D laser direct writing inside silicon. While we have concentrated our efforts
on transverse line inscription at a single focusing depth (120 µm), we
anticipate that phase-control elements (e.g., spatial light modulator) can be
employed for precompensating spherical aberrations at arbitrary depths—thus,
allowing ultrafast laser inscription in all three spatial directions. In the
future, a broad range of in-chip applications can be addressed, including
optical functionalization, wafer dicing, and microfluidics.

###### Acknowledgements.

This research has been supported by the Bundesministerium für Bildung und
Forschung (BMBF) through the NUCLEUS project, grant No. 03IHS107A, as well as
the glass2met project, grant No. 13N15290. It has also been supported by the
National Priorities Research Program grant No. NPRP11S-1128-170042 from the
Qatar National Research Fund (member of The Qatar Foundation). The authors
gratefully acknowledge R. Herda and A. Zach (Toptica Photonics) for technical
support.

##  Appendix A Calculations of the delivered energy

In order to calculate the laser energy delivered to the free-electron
subsystem, we first consider the work $W$ done by the electromagnetic field
force $\vec{F}$ acting on the electrons of total charge $q$ during time
interval $dt$ [61]:

$\displaystyle
W=\vec{F}\cdot\vec{dl}=q(\vec{\mathcal{E}}+\vec{v}\times\vec{\mathcal{B}})\cdot(\vec{v}\,dt)=q(\vec{\mathcal{E}}\cdot\vec{v})dt.$
(3)

where $\vec{\mathcal{E}}$ and $\vec{\mathcal{B}}$ are the electric and
magnetic fields, respectively, and $\vec{v}$ is the electron velocity. For
continuous charge distribution $q=\rho\,dV$, with electron density $\rho$ and
volume element $dV$, we obtain

$\displaystyle W$
$\displaystyle=(\rho\,dV)(\vec{\mathcal{E}}\cdot\vec{v})dt=\vec{\mathcal{E}}\cdot(\rho\vec{v})\,dVdt$
$\displaystyle=(\vec{\mathcal{E}}\cdot\vec{J})\,dVdt,$ (4)

where $\vec{J}=\rho\vec{v}$ is the current. Thus, the work done by the
electromagnetic field force per unit time, per unit volume, that is, the
energy delivered to the free electrons per unit time, per unit volume, is
given by

$\displaystyle\frac{W}{dVdt}=\vec{\mathcal{E}}\cdot\vec{J}.$ (5)

In turn, we are interested in the energy delivered per unit volume, which we
calculate as

$\displaystyle\frac{W}{dV}=\int\vec{\mathcal{E}}\cdot\vec{J}\,dt.$ (6)

In our simulations we use this delivered energy, measured in units of
J$\cdot$m-3, to estimate the breakdown probability inside silicon [43].

Figure 11:  Simulated fluence and delivered energy distributions for different
focusing depths indicated. The green shading marks the region below 40 µm
corresponding to surface damage in the experiments.

In order to provide an intuition, Fig. 11 shows the fluence and delivered
energy distributions, calculated by Eq. (6), for different focusing depths.
One can see that, compared to the fluence, the delivered energy distribution
is much more localized in space (note the logarithmic scale for the delivered
energy distributions in Fig. 11). Following the experiments, we are interested
only in bulk breakdown in silicon and, therefore, exclude from our analysis
the energy delivered at depths below 40 µm (marked by the green shading in
Fig. 11).

##  Appendix B Phenomenological model for transverse inscription

Let us examine theoretically the conditions allowing transverse inscription of
continuous lines inside silicon. To do so, we develop a phenomenological model
relying on the laws of modification growth established in Sec. III.4 and
geometrical considerations. Intuitively, high writing speeds should lead to a
non-continuous (i.e, a discrete) line morphology consisting of individual
single-pulse modifications next to each other. The discrete morphology thus
solely exists when there is no interplay between individual modifications. In
other words, the modification produced by one pulse does not provoke any
shielding nor scattering effect during the subsequent irradiation.

Without limiting the generality of the model, let us consider a sample
movement—and, thus, an inscription—along the $y$-axis, as schematically
depicted in Fig. 12 where the laser irradiation (in red) is subsequent to a
laser-produced modification (in black). The two consecutive irradiations are
spatially separated by a distance $\Delta y=v/\Omega_{0}$, where $v$ is the
writing speed, and $\Omega_{0}=100$ Hz is the laser repetition rate. In this
figure, we introduce the angle $\varphi$ between the focal plane and the point
of the modification with the closest distance to the optical axis $z$. The
condition for discrete line inscription is thus

$\varphi\leq\frac{\pi}{2}-\theta,$ (7)

where $\theta=\text{arcsin}(\text{NA}/n_{0})$ is the maximum half-angle of the
cone of light, NA is the numerical aperture of the focusing lens, and $n_{0}$
is the refractive index of silicon at the laser wavelength. The angle
$\varphi$ reads

$\varphi=\text{arctan}\left(\frac{L_{z}}{\Delta y-L_{y}/2}\right).$ (8)

Here, $L_{z}$ and $L_{y}$ are the maximum lengths of the modification along
$z$ and $y$, respectively. In the framework of our study, as established in
Sec. III.4 for modification growth at the same position, these lengths both
scale logarithmically with the number of pulses as

$L_{y,z}=A_{y,z}\text{ln}(N)+B_{y,z},$ (9)

where $A_{y,z}$ and $B_{y,z}$ are the growth coefficient and the length of a
single-pulse modification along $y$ and $z$, respectively, and $N\geq 1$ is
the number of pulses applied during the growth sequence. In the case of
transverse laser inscription where the focus position is continuously moved
perpendicular to the laser, the number of pulses per point can be defined as
$N=2w_{0}\Omega_{0}/v$, where $w_{0}$ is the beam waist. Recalling that
$\text{tan}(\pi/2-x)=1/\text{tan}(x)$, the combination of Eqs. (7) and (8)
yields

$\frac{L_{z}}{v/\Omega_{0}-L_{y}/2}\leq\frac{1}{\text{tan}\left(\text{arcsin}(\text{NA}/n_{0})\right)}.$
(10)

Figure 12:  Geometrical considerations associated with the phenomenological
model for transverse inscription. The modification (in black) is produced by a
pulse preceding the current irradiation (in red).

Finally, a condition on the writing speed $v$ can be expressed by combining
Eqs. (9) and (10). This condition takes the form

$\alpha v+\beta\text{ln}(v)\geq\gamma$ (11)

where

$\begin{split}\alpha=&2/\Omega_{0},\\\
\beta=&A_{y}+2A_{z}\text{tan}\left(\text{arcsin}(\text{NA}/n_{0})\right),\\\
\gamma=&A_{y}\text{ln}(2w_{0}\Omega_{0})+B_{y}\\\
&+2(A_{z}\text{ln}(2w_{0}\Omega_{0})+B_{z})\text{tan}(\text{arcsin}(\text{NA}/n_{0}).\\\
\end{split}$ (12)

The threshold minimum speed $v_{\text{min}}$ for discrete line inscription
appears in Eq. (11) when the left- and the right-hand sides are equal. In this
case, the analytical solution is:

$v_{\text{min}}=\frac{\beta}{\alpha}W\left(\frac{\alpha}{\beta}\text{exp}\left({\frac{\gamma}{\beta}}\right)\right).$
(13)

Here, the notation $W$ denotes the Lambert $W$-function [62]. Inserting in Eq.
(12) the set of parameters determined experimentally in Sec. III.4, we find a
threshold minimum speed value for discrete line inscription
$v_{\text{min}}\approx 262$ µm$\cdot$s-1.

## References

  * Beresna  _et al._ [2014] M. Beresna, M. Gecevičius, and P. G. Kazansky, Ultrafast laser direct writing and nanostructuring in transparent materials, Advances in Optics and Photonics 6, 293 (2014), 1303.6810 .
  * Chambonneau  _et al._ [2021a] M. Chambonneau, D. Grojo, O. Tokel, F. Ö. Ilday, S. Tzortzakis, and S. Nolte, In-volume laser direct writing of silicon – Challenges and opportunities (2021a), unpublished.
  * Kononenko  _et al._ [2012] V. V. Kononenko, V. V. Konov, and E. M. Dianov, Delocalization of femtosecond radiation in silicon, Optics Letters 37, 3369 (2012).
  * Kononenko  _et al._ [2016] V. V. Kononenko, E. V. Zavedeev, and V. M. Gololobov, The effect of light-induced plasma on propagation of intense fs laser radiation in c-Si, Applied Physics A 122, 293 (2016).
  * Zavedeev  _et al._ [2016] E. V. Zavedeev, V. V. Kononenko, and V. I. Konov, Delocalization of femtosecond laser radiation in crystalline Si in the mid-IR range, Laser Physics 26, 016101 (2016).
  * Chanal  _et al._ [2017] M. Chanal, V. Y. Fedorov, M. Chambonneau, R. Clady, S. Tzortzakis, and D. Grojo, Crossing the threshold of ultrafast laser writing in bulk silicon, Nature Communications 8, 773 (2017).
  * Mareev  _et al._ [2020] E. I. Mareev, K. V. Lvov, B. V. Rumiantsev, E. A. Migal, I. D. Novikov, S. Y. Stremoukhov, and F. V. Potemkin, Effect of pulse duration on the energy delivery under nonlinear propagation of tightly focused Cr:forsterite laser radiation in bulk silicon, Laser Physics Letters 17, 015402 (2020).
  * Chambonneau  _et al._ [2021b] M. Chambonneau, Q. Li, V. Y. Fedorov, M. Blothe, K. Schaarschmidt, M. Lorenz, S. Tzortzakis, and S. Nolte, Taming Ultrafast Laser Filaments for Optimized Semiconductor–Metal Welding, Laser & Photonics Reviews 15, 2000433 (2021b).
  * Chambonneau  _et al._ [2016] M. Chambonneau, Q. Li, M. Chanal, N. Sanner, and D. Grojo, Writing waveguides inside monolithic crystalline silicon with nanosecond laser pulses, Optics Letters 41, 4875 (2016).
  * Pavlov  _et al._ [2017] I. Pavlov, O. Tokel, S. Pavlova, V. Kadan, G. Makey, A. Turnali, Ö. Yavuz, and F. Ö. Ilday, Femtosecond laser written waveguides deep inside silicon, Optics Letters 42, 3028 (2017).
  * Matthäus  _et al._ [2018] G. Matthäus, H. Kämmer, K. A. Lammers, C. Vetter, W. Watanabe, and S. Nolte, Inscription of silicon waveguides using picosecond pulses, Optics Express 26, 24089 (2018).
  * Chambonneau  _et al._ [2019a] M. Chambonneau, X. Wang, X. Yu, Q. Li, D. Chaudanson, S. Lei, and D. Grojo, Positive- and negative-tone structuring of crystalline silicon by laser-assisted chemical etching, Optics Letters 44, 1619 (2019a).
  * Kämmer  _et al._ [2019] H. Kämmer, G. Matthäus, K. A. Lammers, C. Vetter, M. Chambonneau, and S. Nolte, Origin of Waveguiding in Ultrashort Pulse Structured Silicon, Laser & Photonics Reviews 13, 1800268 (2019).
  * Turnali  _et al._ [2019] A. Turnali, M. Han, and O. Tokel, Laser-written depressed-cladding waveguides deep inside bulk silicon, Journal of the Optical Society of America B 36, 966 (2019).
  * Alberucci  _et al._ [2020] A. Alberucci, N. Alasgarzade, M. Chambonneau, M. Blothe, H. Kämmer, G. Matthäus, C. P. Jisha, and S. Nolte, In-Depth Optical Characterization of Femtosecond-Written Waveguides in Silicon, Physical Review Applied 14, 024078 (2020).
  * Ohmura  _et al._ [2007] E. Ohmura, F. Fukuyo, K. Fukumitsu, and H. Morita, Modified-layer formation mechanism into silicon with permeable nanosecond laser, International Journal of Computational Materials Science and Surface Engineering 1, 677 (2007).
  * Verburg  _et al._ [2014] P. C. Verburg, G. R. B. E. Römer, and A. J. Huis in ’t Veld, Two-photon–induced internal modification of silicon by erbium-doped fiber laser, Optics Express 22, 21958 (2014).
  * Tokel  _et al._ [2017] O. Tokel, A. Turnalı, G. Makey, P. Elahi, T. Çolakoğlu, E. Ergeçen, Ö. Yavuz, R. Hübner, M. Zolfaghari Borra, I. Pavlov, A. Bek, R. Turan, D. K. Kesim, S. Tozburun, S. Ilday, and F. Ö. Ilday, In-chip microstructures and photonic devices fabricated by nonlinear laser lithography deep inside silicon, Nature Photonics 11, 639 (2017).
  * Chambonneau  _et al._ [2018a] M. Chambonneau, D. Richter, S. Nolte, and D. Grojo, Inscribing diffraction gratings in bulk silicon with nanosecond laser pulses, Optics Letters 43, 6069 (2018a).
  * Wang  _et al._ [2019] X. Wang, X. Yu, H. Shi, X. Tian, M. Chambonneau, D. Grojo, B. DePaola, M. Berg, and S. Lei, Characterization and control of laser induced modification inside silicon, Journal of Laser Applications 31, 022601 (2019).
  * Wang  _et al._ [2020a] X. Wang, X. Yu, M. Berg, B. DePaola, H. Shi, P. Chen, L. Xue, X. Chang, and S. Lei, Nanosecond laser writing of straight and curved waveguides in silicon with shaped beams, Journal of Laser Applications 32, 022002 (2020a).
  * Richter  _et al._ [2020] R. A. Richter, N. Tolstik, S. Rigaud, P. D. Valle, A. Erbe, P. Ebbinghaus, I. Astrauskas, V. Kalashnikov, E. Sorokin, and I. T. Sorokina, Sub-surface modifications in silicon with ultra-short pulsed lasers above 2 $\backslash$textmu m, Journal of the Optical Society of America B 37, 2543 (2020).
  * Nejadmalayeri  _et al._ [2005] A. H. Nejadmalayeri, P. R. Herman, J. Burghoff, M. Will, S. Nolte, and A. Tünnermann, Inscription of optical waveguides in crystalline silicon by mid-infrared femtosecond laser pulses, Optics Letters 30, 964 (2005).
  * Sreenivas  _et al._ [2012] V. V. P. Sreenivas, M. Bülters, and R. B. Bergmann, Microsized subsurface modification of mono-crystalline silicon via non-linear absorption, Journal of the European Optical Society: Rapid Publications 7, 12035 (2012).
  * Chambonneau  _et al._ [2019b] M. Chambonneau, L. Lavoute, D. Gaponov, V. Fedorov, A. Hideur, S. Février, S. Tzortzakis, O. Utéza, and D. Grojo, Competing Nonlinear Delocalization of Light for Laser Inscription Inside Silicon with a 2-$\backslash$textmu m Picosecond Laser, Physical Review Applied 12, 024009 (2019b).
  * Wang  _et al._ [2020b] A. Wang, A. Das, and D. Grojo, Temporal-contrast imperfections as drivers for ultrafast laser modifications in bulk silicon, Physical Review Research 2, 033023 (2020b).
  * Wang  _et al._ [2020c] A. Wang, A. Das, and D. Grojo, Ultrafast Laser Writing Deep inside Silicon with THz-Repetition-Rate Trains of Pulses, Research 2020, 1 (2020c).
  * Li [1980] H. H. Li, Refractive index of silicon and germanium and its wavelength and temperature derivatives, Journal of Physical and Chemical Reference Data 9, 561 (1980).
  * Li  _et al._ [2016] Q. Li, M. Chambonneau, M. Chanal, and D. Grojo, Quantitative-phase microscopy of nanosecond laser-induced micro-modifications inside silicon, Applied Optics 55, 9577 (2016).
  * Baudisch  _et al._ [2018] M. Baudisch, M. Beutler, M. Gebhardt, C. Gaida, F. Stutzki, S. Hadrich, R. Herda, K. Zawilski, P. Schunemann, A. Zach, J. Limpert, and I. Rimke, 2.3-12 $\mu$m tunable, sub-10 optical cycle, ZnGeP2-based OPA directly pumped by a Tm:fiber laser at 1.96 $\mu$m and 100 kHz, in  _Laser Congress 2018 (ASSL)_, Vol. Part F121- (OSA, Washington, D.C., 2018) p. AW4A.3.
  * Li  _et al._ [2021] Q. Li, M. Chambonneau, M. Blothe, H. Gross, and S. Nolte, A flexible, fast and benchmarked vectorial model for focused laser beams (2021), accepted for publication in Applied Optics.
  * InF [2021] _InFocus_ (2021), https://github.com/QF06/InFocus.
  * Bristow  _et al._ [2007] A. D. Bristow, N. Rotenberg, and H. M. van Driel, Two-photon absorption and Kerr coefficients of silicon for 850–2200nm, Applied Physics Letters 90, 191104 (2007).
  * Lin  _et al._ [2007] Q. Lin, J. Zhang, G. Piredda, R. W. Boyd, P. M. Fauchet, and G. P. Agrawal, Dispersion of silicon nonlinearities in the near infrared region, Applied Physics Letters 91, 021111 (2007).
  * Wang  _et al._ [2013] T. Wang, N. Venkatram, J. Gosciniak, Y. Cui, G. Qian, W. Ji, and D. T. H. Tan, Multi-photon absorption and third-order nonlinearity in silicon at mid-infrared wavelengths, Optics Express 21, 32192 (2013).
  * Kolesik and Moloney [2004] M. Kolesik and J. V. Moloney, Nonlinear optical pulse propagation simulation: From maxwell’s to unidirectional equations, Physical Review E 70, 036604 (2004).
  * Fedorov  _et al._ [2016] V. Y. Fedorov, M. Chanal, D. Grojo, and S. Tzortzakis, Accessing Extreme Spatiotemporal Localization of High-Power Laser Radiation through Transformation Optics and Scalar Wave Equations, Physical Review Letters 117, 043902 (2016).
  * Mouskeftaras  _et al._ [2014] A. Mouskeftaras, A. V. Rode, R. Clady, M. Sentis, O. Utéza, and D. Grojo, Self-limited underdense microplasmas in bulk silicon induced by ultrashort laser pulses, Applied Physics Letters 105, 191103 (2014).
  * Jonsson  _et al._ [1997] P. Jonsson, H. Bleichner, M. Isberg, and E. Nordlander, The ambipolar Auger coefficient: Measured temperature dependence in electron irradiated and highly injected n-type silicon, Journal of Applied Physics 81, 2256 (1997).
  * Tanimura  _et al._ [2019] H. Tanimura, J. Kanasaki, K. Tanimura, J. Sjakste, and N. Vast, Ultrafast relaxation dynamics of highly excited hot electrons in silicon, Physical Review B 100, 035201 (2019).
  * Mouskeftaras  _et al._ [2016] A. Mouskeftaras, M. Chanal, M. Chambonneau, R. Clady, O. Utéza, and D. Grojo, Direct measurement of ambipolar diffusion in bulk silicon by ultrafast infrared imaging of laser-induced microplasmas, Applied Physics Letters 108, 041107 (2016).
  * Chimier  _et al._ [2011] B. Chimier, O. Utéza, N. Sanner, M. Sentis, T. Itina, P. Lassonde, F. Légaré, F. Vidal, and J. C. Kieffer, Damage and ablation thresholds of fused-silica in femtosecond regime, Physical Review B 84, 094104 (2011).
  * Zhokhov and Zheltikov [2018] P. A. Zhokhov and A. M. Zheltikov, Optical breakdown of solids by few-cycle laser pulses, Scientific Reports 8, 1824 (2018).
  * Bloembergen [1974] N. Bloembergen, Laser-induced electric breakdown in solids, IEEE Journal of Quantum Electronics 10, 375 (1974).
  * Glebov  _et al._ [1984] L. B. Glebov, O. M. Efimov, G. T. Petrovskiĭ, and P. N. Rogovtsev, Influence of the mode composition of laser radiation on the optical breakdown of silicate glasses, Soviet Journal of Quantum Electronics 14, 226 (1984).
  * Smith and Do [2008] A. V. Smith and B. T. Do, Bulk and surface laser damage of silica by picosecond and nanosecond pulses at 1064 nm, Applied Optics 47, 4812 (2008).
  * Batavičiūtė  _et al._ [2013] G. Batavičiūtė, E. Pupka, V. Pyragaitė, L. Smalakys, and A. Melninkaitis, Effect of longitudinal laser mode beating in damage probability measurements (2013) p. 88851M.
  * Diaz  _et al._ [2014] R. Diaz, M. Chambonneau, R. Courchinoux, P. Grua, J. Luce, J.-L. Rullier, J.-Y. Natoli, and L. Lamaignère, Influence of longitudinal mode beating on laser-induced damage in fused silica, Optics Letters 39, 674 (2014).
  * Lamaignère  _et al._ [2017] L. Lamaignère, R. Diaz, M. Chambonneau, P. Grua, J.-Y. Natoli, and J.-L. Rullier, Correlation between laser-induced damage densities of fused silica and average incubation fluences at 1064 nm in the nanosecond regime, Journal of Applied Physics 121, 045306 (2017).
  * Chambonneau  _et al._ [2014] M. Chambonneau, R. Diaz, P. Grua, J.-L. Rullier, G. Duchateau, J.-Y. Natoli, and L. Lamaignère, Origin of the damage ring pattern in fused silica induced by multiple longitudinal modes laser pulses, Applied Physics Letters 104, 021121 (2014).
  * Chambonneau  _et al._ [2018b] M. Chambonneau, J.-L. Rullier, P. Grua, and L. Lamaignère, Wavelength dependence of the mechanisms governing the formation of nanosecond laser-induced damage in fused silica, Optics Express 26, 21819 (2018b).
  * Iwata  _et al._ [2018] H. Iwata, D. Kawaguchi, and H. Saka, Crystal structures of high-pressure phases formed in Si by laser irradiation, Microscopy 67, 30 (2018).
  * Verburg  _et al._ [2015] P. C. Verburg, L. A. Smillie, G. R. B. E. Römer, B. Haberl, J. E. Bradby, J. S. Williams, and A. J. Huis in ’t Veld, Crystal structure of laser-induced subsurface modifications in Si, Applied Physics A 120, 683 (2015).
  * Rapp  _et al._ [2015] L. Rapp, B. Haberl, C. Pickard, J. Bradby, E. Gamaly, J. Williams, and A. Rode, Experimental evidence of new tetragonal polymorphs of silicon formed through ultrafast laser-induced confined microexplosion, Nature Communications 6, 7555 (2015).
  * Norton  _et al._ [2006] M. A. Norton, E. E. Donohue, M. D. Feit, R. P. Hackel, W. G. Hollingsworth, A. M. Rubenchik, and M. L. Spaeth, Growth of laser damage on the input surface of SiO2 at 351 nm, in  _Proceedings of SPIE_, Vol. 6403 (2006) p. 64030L.
  * Sozet  _et al._ [2016] M. Sozet, J. Neauport, E. Lavastre, N. Roquin, L. Gallais, and L. Lamaignère, Laser damage growth with picosecond pulses, Optics Letters 41, 2342 (2016).
  * Chambonneau and Lamaignère [2018] M. Chambonneau and L. Lamaignère, Multi-wavelength growth of nanosecond laser-induced surface damage on fused silica gratings, Scientific Reports 8, 891 (2018).
  * Mori  _et al._ [2015] M. Mori, Y. Shimotsuma, T. Sei, M. Sakakura, K. Miura, and H. Udono, Tailoring thermoelectric properties of nanostructured crystal silicon fabricated by infrared femtosecond laser direct writing, Physica Status Solidi (a) 212, 715 (2015).
  * Das  _et al._ [2020] A. Das, A. Wang, O. Uteza, and D. Grojo, Pulse-duration dependence of laser-induced modifications inside silicon, Optics Express 28, 26623 (2020).
  * Luong  _et al._ [2020] K. P. Luong, R. Tanabe-Yamagishi, N. Yamada, and Y. Ito, Micromachining of the back surface and interior of Si using infrared femtosecond laser pulses, Journal of Laser Applications 32, 012017 (2020).
  * Griffiths [2013] D. J. Griffiths,  _Introduction to Electrodynamics, Fourth Edition, Chap. 8, pp. 357_ (Pearson Education Inc., 2013).
  * Corless  _et al._ [1996] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, On the Lambert W function, Advances in Computational Mathematics 5, 329 (1996).