# How to Crystallize Glass with a Femtosecond Laser

Ruyue Que, Matthieu Lancry, Maxime Cavillon, Bertrand Poumellec

► To cite this version:

Ruyue Que, Matthieu Lancry, Maxime Cavillon, Bertrand Poumellec. How to Crystallize Glass with a Femtosecond Laser. Crystals, 2024, 14 (7), pp.606. 10.3390/cryst14070606. hal-05466525

HAL Id: hal-05466525

https://hal.science/hal-05466525v1

Submitted on 19 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Article

# How to Crystallize Glass with a Femtosecond Laser

Ruyue Que, Matthieu Lancry, Maxime Cavillon and Bertrand Poumellec *

Institut de Chimie Moléculaire et des Matériaux d'Orsay, Université Paris-Saclay, CNRS, 91405 Orsay, France;
queruyue@hotmail.com (R.Q.); mathieu.lancry@universite-paris-saclay.fr (M.L.);
maxime.cavillon@universite-paris-saclay.fr (M.C.)

* Correspondence: bertrand.poumellec@universite-paris-saclay.fr

**Abstract:** The crystallization of glass through conventional thermal annealing in a furnace is a well-understood process. However, crystallization by femtosecond (fs) laser brings another dimension to this process. The pulsed nature of the irradiation necessitates a reevaluation of the parameters for optimal crystallization and an understanding of the particularities of using fs laser. This includes adjusting the laser pulse energy, the repetition rate, and the writing speed to either initiate nucleation or achieve substantial crystal growth. Therefore, a key challenge of this work is to establish reliable calculations for understanding the link between the size of the crystallized region and an ongoing transition (e.g., solid-to-solid, liquid-to-solid), while accounting for the aforementioned laser parameters. In this context, and based on previous work, a temperature distribution (in space and time) is simulated to model the thermal treatment at any point in the glass. By setting the condition that the temperatures are between the glass transition and melting temperature, the simulated crystallized region size can be compared with experimental observations. For that purpose, knowledge of the beam width at the focal point and of the absorbed beam energy fraction are critical inputs that were extracted from experiments found in the literature. After that, the management of the crystallization process and the width of the crystallization line can be achieved according to pulse energy, e.g., crystallite size, and also the effect of the scanning speed can be understood. Among the main conclusions to highlight, we disclose the laser conditions that determine the extent of the crystallized area and deduce that it is never of interest to increase the pulse energy too much as opposed to the repetition rate for the uniform crystallized line.

**Keywords:** laser induced crystallization in glass, crystallization of oxide glass, femtosecond laser

## 1. Introduction

Three dimensional localized crystallization in glasses brings interest as it introduces the possibility of precipitating crystalline phases that present optical non-linear, ferroelectric, or pyroelectric properties. For example, it leads to the formation of non-center-symmetric crystals in optically suitable glasses [1]. This approach underscores that crystallization, a thermo-kinetics process [2,3], is influenced by temperature (and time), occurring between the glass transition temperature ($T_g$) and the melting temperature ($T_m$) of the given precipitated phase(s), and dictating nucleation and growth rates. While spatially selective crystallization using a femtosecond laser has been extensively demonstrated in various glasses [1], this process appears more complex than using a continuous-wave (CW) laser, arising from the pulsed nature of the laser (up to several MHz), and its high peak intensity triggering non-linear processes. For example, in lithium niobium silicate glasses (LNS) [4], crystallization patterns of the $LiNbO_3$ phase strongly depend on pulse energy or laser polarization.

While some physical aspects of crystallization mechanisms have been understood, we aim to further explain and predict semi-quantitative aspects, particularly the correlation between laser parameters and specific optical or material properties. A goal is to

**Citation:** Que, R.; Lancry, M.; Cavillon, M.; Poumellec, B. How to Crystallize Glass with a Femtosecond Laser. *Crystals* **2024**, *14*, 606. https://doi.org/10.3390/cryst14070606

Academic Editor: Xiaoming Duan

Received: 22 April 2024

Revised: 15 June 2024

Accepted: 22 June 2024

Published: 30 June 2024

<figure><img src="image_2.png" /></figure>

**Copyright:** © 2024 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
minimize the number of experiments by performing calculations that facilitate a selection of relevant parameters to induce crystallization. For this purpose, a simple analytical solution has been developed in a previous work [5] to physically describe the thermal dynamics induced by pulsed laser absorption and the subsequent transformations.

In this work, we first simulate the evolution of the temperature distribution induced by a pulsed fs laser. Once the temperature modulation is known, it can be combined with a time-temperature transformation (TTT) diagram for crystallization. The temperature calculation relies on modeling the energy source as spherical, and requires knowledge of several parameters:

(1) The ratio (labeled $R_τ$) between the pulse period ($τ_{RR} = 1/RR$, with RR being the pulse repetition rate) and the heat diffusion time ($τ_D = w^2/4D_{th}$, w being the light beam waist radius at $1/e$ at the focus, and $D_{th}$ the heat diffusion coefficient or diffusivity). While $R_τ$ depends thus on w, the latter is usually not reliably known due to the complexity of non-linear effects in the light-matter interaction. For this purpose, a compilation of values is extracted from the literature, allowing us to clarify its dependency on the incident pulsed laser energy (please refer to Appendix A).

(2) The fraction of the incident energy effectively absorbed by the material (labeled A). Absorption primarily occurs through non-linear absorption (multiphoton ionization). Plasma reflection or scattering is also expected, reducing the coefficient A. After the absorption of a small part of the pulse energy, the absorption becomes linear due to limitations of the electron plasma density [6] via a combination of multistep absorption, trapping, and electron-electron or electron-phonon collisions. Given the lack of precise information about all these processes, extracted values from the literature (reported in Appendix B) enable one to set A values. As one is not aware of that, A is a function of the incident pulse energy, the pulse duration, and the material composition.

Another problem more specific to crystallization is the location of the crystallization domain in the TTT diagram to the confined condition of the tightly focused fs pulsed laser. Indeed, the specificities of the process may change the limit of the crystallization domain compared to the classical one [7–9]. Due to these uncertainties, the approach is tested using LNS glasses for the laser conditions employed in [4,10], for which these aspects are relatively well known.

In this work, brief recall of the classical crystallization process (Section 2) is provided. Thus temperature calculation from a specific hypothesis (Section 3) is given, describing how crystallization process behaves first within a heat accumulation regime (later defined, see Section 4) and, secondly, without heat accumulation (see Section 5, large temperature oscillations). Finally, the results from the approach are compared with several experimental values according to pulse energy and scanning speed (Section 6) and we draw a conclusion (1) about the optimization of laser parameters, (2) that the width of the crystallization area can be understood, and (3) that the incident light energy fraction may change with the number of pulses.

## 2. Basic: Brief Recall of the Crystallization Process

The appearance of crystals in glass is a two-step process: nucleation and growth [2,9]. Crystallization is generally growth-limited at low temperatures and nucleation-limited at higher ones. Moreover, the domain of crystallization is defined from a TTT diagram [9]. It collects the time (called the incubation time) requested at a given temperature to observe a first seed. This diagram, generally established for isothermal treatment, takes the time origin (t = 0) when the sample is in the condition to crystallize for a given temperature. This temperature is incident between $T_g$ and $T_m$, with crystallization occurring between these temperatures with only one transformation (otherwise, other temperatures have to be considered). A boundary can be plotted in this TTT diagram delimiting the appearance of typ. $10^{-6}$ part crystallized in the glass volume (thick black line in Figure 1). The position of the crystallization domain in the TTT diagram is dependent on the treatment pathway. For example, its position will be shifted if the temperature is continuously decreasing (or
not) from $T_m$ [11]. In that case, the growth is delayed because the nucleation is slowed down. But, for the purpose of this paper, this is ignored and we simply consider the thermal treatment curve experienced by the material at a given spatial position with the coordinates $\alpha_0, \beta_0, \gamma_0$. The time-temperature function is applied to this position ($T(\alpha_0, \beta_0, \gamma_0, t)$). This function of time is simply obtained by the displacement of the temperature distribution function (the variation of the temperature in space) in the static case. This approximation is valid for speeds smaller than the diffusion speed ($2D_{th}/w$, typ. a few cm/s). As recalled above, the thermal treatment should exhibit temperatures in the $T_g$ to $T_m$ range. However, heating a material using a pulsed laser, resulting in an oscillatory temperature in volume impacted by non-linear optical processes, is a complex problem that prevents one from easily drawing conclusions for managing the process. The important aspect of solving this issue is to take into account the effect of these parameters while preserving the main physics associated with the induced transformations. Consequently, in [5] we developed a simplified approach for understanding how $T(\alpha, \beta, \gamma, t)$ changes according to the laser parameters (pulse energy, pulse repetition rate, beam scanning speed, beam width, and material properties).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Crystallization domain in a Time-Temperature Transformation (TTT) diagram (in yellow). Its boundary is black. The glass transition and melting temperatures are lines shown in blue. $T_1$ is an example of isothermal treatment that leads to the incubation time $t_1$ required for crystallization to begin.</figcaption></figure>

### 3. Temperature Distribution and Its Evolution during Laser Writing

For this purpose, a known fraction of the incident energy is assumed to be transformed completely into heat after a period of time that is much shorter than the diffusion time $w^2/D_{th}$, where $w$ is the waist radius at 1/e and $D_{th}$ the thermal diffusivity. For most of the inorganic glasses, the electron-phonon coupling time is around a few ps, and the diffusion time is smaller than µs. This means that the pulse duration is assumed to be smaller than one hundred ns. So, it is largely the case that for the fs laser, the pulse duration is a fraction of a ps. Therefore, the heat diffusion starts after the electron relaxation and occurs only in the lattice phonons if there are no physical or chemical reactions. Consequently, the complicated events occurring at the ps level have no effect on the present calculation; otherwise, if the time averaged temperature, for instance, may change some non-negligible part of the pulse energy, it is not taken into account.

It is also assumed that no endothermic reaction is active in the dark that can use a significant part of the absorbed energy. The band structure of the material is assumed as not changing in the course of the irradiation, as is the non-linear absorption coefficient.

The 3D Fourier equation, with a spherical energy source at the focal volume, assuming a tightly focused laser beam and temperature-independent physico-chemical properties, is employed. Usually, for 3D direct laser writing, the NA of the focusing objective is above 0.5. This section is mostly drawn from what has been already published in [5]. The approximations taken into account are that the pulse durations considered are such that the thermal diffusion is not yet active, i.e., $\langle\langle 0.1 \ \mu\text{s} \rangle\rangle$ and that the scanning speed ($v$) is much lower than the diffusion speed (typ. few cm/s). Consequently, the shape of the spatial
temperature distribution is not time dependent; the only time dependence is introduced
by the beam scanning taken in the α direction and appears as α + v.t. β and γ are fixed at
the position β₀ and γ₀, with γ direction corresponding to the beam propagation axis. The
scanning speed (v) defines the pulsed number deposited immediately as Nₚ = 2w.RR/v.
The laser pulse energy Eₚ defines the maximum temperature increment T₀₀, introduced by
a single pulse above an initial temperature Tinit: T₀₀ = A:Eₚ / (π²ρCₚw³), with A the fraction of ab-
sorbed energy (previously defined); ρ and Cₚ are the glass density and heat capacity, re-
spectively. Additionally, a normalized radius rᵣ = r/w is defined, with r being a radial dis-
tance from the center of the beam focus.

All calculated temperatures are normalized with respect to T₀₀, as follows: Tₙorm = (T-
Tinit)/T₀₀. Tₙorm is oscillating at periods of 1/RR at any distance from the focus center. This
temperature oscillation is occurrent between a maximum (Tmax) and a minimum (Tmin),
which varies according to the laser irradiation conditions. After a sufficient number of
pulses, both Tmax and Tmin reach a steady state regime. Achieving this steady state regime
is dictated by the ratio Rτ as previously defined and taking the form Rτ = 4Dth/RR w². For
instance, for large Rτ values (i.e., >16 for Tmin < 0.03Tmax) the pulse contributions are inde-
pendent of each other, whereas for smaller Rτ values (<16), heat accumulation from pulse-
to-pulse contribution begins to appear (and thus Tmin overcomes 0.03Tmax).

$N_{ss}$ corresponds to the minimum number of pulses necessary to reach the steady state regime. Its expression is given in Equations (A4) in Appendix D. This calculation can be performed either for $T_{min}$ and $T_{max}$.

The temperature-spatial distributions of Tmax and Tmin (normalized in the form of
Tnorm) for various Rτ values are provided in Figure 2. When Rτ is small (large RR or small
diffusivity), Tmin and Tmax have no large relative difference compared to their average val-
ues because the oscillation amplitude is always limited to exp[−(r<sub>w</sub>)<sup>2</sup>], whereas for Tmax
and Tmin amplitude converges to 2/Rτ at the beam center but converges to π/Rτr<sub>w</sub> erf(r<sub>w</sub>)
elsewhere on heat accumulation. When this situation can be reached, the pulse energy can
be adjusted to compensate for an RR increase. It is also worth noting that, when Rτ de-
creases, the shape of the curve converges to the erf(r<sub>w</sub>)/r<sub>w</sub> curve (see Equation (A6) Appen-
dix D). It is a consequence that the temperature distribution decreases much slower at
large radius, compared to a gaussian one. This last shape is the shape of the energy source
and of the temperature distribution for the case where Rτ is large. Note that in the case
of heat accumulation, the overall maximum temperature (actually at the center of the
beam) is defined by A.E<sub>p</sub>.RR, i.e., roughly the average incident power.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Spatial distribution of *T*<sub>min</sub> (blue dash, by Equation (A3)) and *T*<sub>max</sub> (red, by Equation (A1)) in the normalized form (see text) and as a function of the reduced radius *r*<sub>w</sub> when (a) *R*<sub>*τ*</sub> = 0.1, (b) *R*<sub>*τ*</sub> = 1, and (c) *R*<sub>*τ*</sub> = 10.</figcaption></figure>

For large R<sub>τ</sub> values, corresponding to a small RR and/or large D<sub>th</sub>, the oscillations
are relatively large as pulse contributions are separated and thus T<sub>min</sub> appears to have
small values (See Figure 2c). If one considers negligible where T<sub>min</sub> is smaller than 3% of
T<sub>max</sub>, then R<sub>τ</sub> > 16 (for values smaller than 6% of T<sub>max</sub>, then R<sub>τ</sub> > 7; for R<sub>τ</sub> > 10, T<sub>min</sub> is
smaller than 4.4% of $T_{max}$). This defines the limit of the domain of heat accumulation. The increase in $R_τ$ (Figure 2a–c) clearly shows that the shape of $T_{max}$ is also converging with the increasing $R_τ$ to the shape of the beam energy distribution (a gaussian beam, with an envelope having a spatial distribution following $\exp[-(r_w)^2]$), therefore (and in the limited case) becoming independent of $R_τ$. The maximum absolute temperature at the center is $T_{oo}$ in this case. It is defined by $E_p$ only (independent of $R_τ$ and so especially of RR)

For intermediate $R_τ$ values, exemplified in Figure 2b,c, the temperature oscillations are limited between $T_{max}$ and $T_{min}$. This is shown in particular cases with the energy source volume appearing as a shoe box in [12] for $R_τ = 2$ and 20 or in [13] for $R_τ = 20$.

But, in this paper, it is specified that when $r_w > 2$, the difference between $T_{min}$ and $T_{max}$ vanishes, as seen in [14]. Therefore, the conditions of $r_w$ and $R_τ$ for neglecting the temperature oscillations are deduced and the use of an average temperature $T_{mean}$ might be applicable (Appendix D Equations (A6)–(A8)).

Another important point for this current work is the dependence of the full width at half maximum (FWHM) of the temperature distribution with respect to $R_τ$. Its variation is less than a factor 2 in the relative radius scale ($r_w$). In fact, the largest change is at the base of the distribution.

**4. Case of Low $R_τ$ Values (High RR or Low $D_{th}$), i.e., within Heat Accumulation Regime**

Here, $R_τ$ is such that it is much smaller than 10, e.g., $R_τ = 0.1$, due to the use of a large repetition rate or a weakly thermally diffusing material $T_{min}$ is very close to $T_{max}$ ($R_τ = 10$ is the limit for $T_{min} = 4.4\% T_{max}$). In that case, a pulse contribution has no time to vanish significantly before the next pulse absorption induces another temperature increase. This is the case of heat accumulation, and the behavior is like that of Continuous Wave (CW) lasers.

In this section, we deal with conditions where the pulse number received at the same point ($N_p$ as defined earlier) is above the number of pulses required to reach the steady state ($N_{ss}$). This situation yields two possibilities for which the temperature is higher than $T_g$: (i) the temperature remains lower than $T_m$, and (ii) the temperature overcomes $T_m$ for some distance from the center of the beam focus during the thermal treatment. Each situation is described below.

**(i) Moderate $E_p$ leading to a temperature lower than $T_m**

To describe the temperature evolution and for illustration purposes, the laser parameters are set so that the temperature values fall within an acceptable temperature interval. In this section, we took as an example $T_g = 579$ °C coming from the physico-chemical properties of LNS glass and $T_m = 1257$ °C ($T_m$ is the melting temperature of the formed LiNbO₃ crystals in LNS). Additional data requested for the calculation are provided in the Table A3 in Appendix C. Therefore, low absorption fraction (A) and pulse energy ($E_p$) values are used, as provided in the caption of Figure 3. The beam width at 1/e has been fixed to a "reasonable value" for such a small pulse energy, i.e., $w = 1$ µm. In this condition, similarly to Figure 2a, the oscillations (Equation (A2) in Appendix D) are small enough to be neglected. This can be observed in Figure 3, rendering $T_{mean}$ a valid approximation if requested.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Temperature distribution for $T_{max}$ and $T_{min}$ with respect to the normalized radial distance $r_w$ for $R_τ = \frac{4D_{th}}{RR}w^2$ (called Rtau1 in the figures of this paper) = 0.13, $A = 0.1$, $w = 1.65$ µm, RR = 300 kHz, Ep (called Ep1 in the figures of this paper) = 54 nJ for (a) and 37 nJ for (b); $T_g$ and $T_m$, shifted by room temperature, are also shown. $T_{max}$ and $T_{min}$ are calculated from Equations (A1) and (A3) in Appendix D.</figcaption></figure>

Then, the temperature distribution is transformed into the time domain to obtain thermal treatment curves (Figure 4b,c, T-t dependency), replacing the static normalized radius $r_w$ by the coordinate of the place looked at in the moving material $r_{wd}$, taking the beam center as a reference (see Figure 4a). The relative position in the beam of this point is related to relative time from the following relation: $r_{wd}(t, v, d, \beta) = \sqrt{(v \cdot t/w + d)^2 + \beta^2}$. Then, the time origin of the treatment for crystallization is defined when T starts to overcome $T_g$. This occurs at a distance $d$ from the center of the temperature distribution, i.e., at an abscissa $d$ dependent on $\beta$ (distance from the line axis), i.e., $d^2 = r_{wdg}^0 - \beta^2$. $r_{wdg}^0$; the radius of the isotherm $T_g$ is then obtained by solving $T_{mean}(r_{wdg}^0, R_τ) = T_g/T_{00}$, where $T_{mean}$ is an average temperature defined in Appendix D Equation (A6). $r_{wdg}^0$ is also the maximum width of the crystallized region.

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>

**Figure 4.** Plot of the treatment curves in the time scale from the time when the temperature overcomes $T_g$: (a) scheme for illustrating the various normalized distances used in the calculation of the thermal treatment curve, (b) treatment curve at the center $\beta = 0$ or at $r = 1.2$ w ($\beta = 1.2$), and (c) same as (b) but showing the time in log scale, where $A = 0.1$, $w = 1.65$ µm, $R_\tau = 0.13$, and $E_p = 54$ nJ. v = 100 µm/s. The thermal treatment of the point at coordinates 0,0,0 begins when the temperature overcomes $T_g$, here when the beam center is distant by 1.9 w (third parameter in red), whereas it begins later at 1.5 w for a point in the periphery of the line at coordinates 0,1,2,0.

Figure 4 displays a time-temperature thermal treatment. Therefore, provided that the crystallization domain is known (TTT diagram, see above), the scanning speed must be adjusted to ensure the thermal curve intersects with the crystallization domain, as in Figure 5, at a speed of $v1$ or $v2$. For a speed higher than $v_{lim}$, the curve does not touch the domain anymore and crystallization is not possible. One can also see that for speed $v1$, crystallization starts from the bottom side of the domain and then crosses it and is thus not limited in growth. However, if $T_{max}$ is not so large ($T_{max} < T_m$), the growth can be controlled. For $v2$, crystallization starts from the growth side and is sensitive to the pre-existence of nucleation seeds.

It is also important to note that $r_{wdg}^0$ is in principle larger than the half-width of the crystallized region. Indeed, $r_{wdg}^0$ is a circle where $t = 0$, but crystallization occurs eventually after some incubation time that could be longer than at the center, due to the decrease in $T_{max}$ when going to the line's edge. This contributes to a narrowing of the line, but this effect is small, as shown in Figure 4 of [15]. In addition, at the speed limit, the width is not vanishing. The disappearance of the crystallization is thus arising from the process described in Figure 5. The width of the crystallization area is weakly sensitive to the scanning speed for a solid-to-solid transformation. This validates the crystallization criterion for $T$ to be between $T_g$ and $T_m$, used in our approach.

<figure><img src="image_8.png" /></figure>

**Figure 5.** Display of the position of the treatment curve versus the crystallization domain in the case of solid glass to crystal transformation ($T_g < T < T_m$).
At the periphery of the line (e.g., for $\beta$ = 1.2 in Figure 4), the temperature is smaller than that at the center, but the crystallization condition remains somewhat similar: the treatment curves cross the crystallization domain almost at the same incubation time. From this viewpoint, the results are rather intuitive. First, a lower scanning speed would lead to a larger crystal growth. Second, there is a minimum laser power (P = Ep.RR) necessary to overcome Tg, with a maximum scanning speed necessary for the thermal treatment curve to cross the crystallization boundary. As the power increases, the speed limit increases until the maximum in space of the average temperature reaches Tm. Beyond this point, the crystallization process is different, and is the topic of the next section.

### (ii) larger *Ep* leading to a temperature larger than *Tm*

A large $E_p$ yields temperatures greater than Tm. In this case, if crystallization occurs, it starts from the liquid phase, with the system potentially crystallizing as the temperature decreases below Tm. Similar to the previous situation, $d(\beta)$ is defined by $d^2 = r_{wdm}^{02} - \beta^2$. $r_{wdm}^0$, which is the radius of the isotherm Tm (i.e., a circle with a radius smaller than the one for $T_g$), and is then obtained by solving the following equation that includes $T_m$: $T_{mean}(r_{wdm}^0, R\tau) = T_m/T_{oo}$. $r_{wdm}^0$ is also the maximum size of the melted region. $\beta < r_{wdm}^0$ is represented by the dashed blue curve in Figure 6. When $r_{wdm}^0 < \beta < r_{wdg}^0$ (green curve), the results of section (i) apply. Then, the scanning speed should be adjusted to ensure the treatment curve reaches the crystallization domain. In this case, the entry into the crystallization domain is necessarily from the top side (nucleation limited). It is sensitive to side stimulation because the side has crystallized before.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** The position of the treatment curve versus the crystallization domain when a part of the trace is melting, i.e., a liquid-to-crystal transformation. The region delimited by the red one is for growth only.</figcaption></figure>

## 5. Case of High $R_\tau$ Values (Low RR and Large $D_{th}$), i.e., Out of Heat Accumulation Regime

Here, $R_\tau$ is such that it is larger than 10, due to the use of a small repetition rate or a strong thermally diffusing material. In that case, a pulse contribution has enough time to vanish significantly before the next pulse absorption arrival. There is no heat accumulation, but in Figure 2c this is true only for when $r_w$ is smaller than a value of approximately 3. For the periphery of the heat-affected region, there is a smoothening of the temperature oscillations due to heat diffusion.

On the contrary, for a smaller reduced radius, $T_{min}$ is close to zero, and the temperature oscillates between $T_{min}$ and $T_{max}$ during each period. Assuming that the pulse number received punctually ($N_p$) is above the pulse number required to reach the stationary state ($N_{ss}$ formula Equation (A4) Appendix D), the temperature oscillations at the steady state can be simulated simply. Note that the oscillation at this stage is due to a series of pulse contributions collected in the function $T_{osc}'(r_w, x, R_\tau)$, where x (between 0 and 1) is the
temporal position in the period and is thus the floating value of *t*.RR. The simulation function of the temperature called *SimpulseT* is *SimpulseT(r<sub>wd</sub>(t, v, d, β), R<sub>τ</sub>, E<sub>p</sub>)/T<sub>00</sub>(E<sub>p</sub>) = T<sub>min</sub>(r<sub>wd</sub>(t, v, d, β), R<sub>τ</sub>) + T<sub>osc</sub>'(r<sub>wd</sub>(t, v, d, β), x(t, RR), R<sub>τ</sub>) (see Equation (A9) Appendix D) with r<sub>wd</sub> defined as in the previous section and related to the time origin for crystallization according to one of the two cases for which the temperature is above or below the melting temperature T<sub>m</sub>.

(i) Moderate *E*<sub>*p*</sub> such that the temperature is not overcoming *T*<sub>*m*</sub>

As it is seen in Figure 7, the punctual thermal treatment is beginning when $T_{max} \cdot T_{00}$ is
overcoming the melting temperature for the first time during the scanning. The relative
position in the beam is translated in relative time from the following relation:
$r_{wd}(t, v, d, \beta) = \sqrt{(v \cdot t/w + d)^2 + \beta^2} \cdot d(\beta) \text{ defined by } d^2 = r_{wdg}^0 - \beta^2 \text{ with } r_{wdg}^0 \text{ the solution of } T_{max}(r_{wdg}^0, R_{\tau}) = T_g/T_{00}$.

However, the existence of oscillations between T_min and T_max is such that there are dead times, defined as the time for which the temperature decreases below T_g. Consequently, the efficiency, taken as the ratio between the active time (during which T_g < T < T_m) and the thermal treatment time, is much smaller than one. This is different from small R_τ values, again for which T_max is close to T_min, and T_mean can be employed.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Thermal treatment curve (T<sub>max</sub>, T<sub>min</sub>, *SimpulseT*) for a point at the center (0,0,0) of the written line. R<sub>τ</sub> = 13, A = 0.15, *w* = 1.65 µm, E<sub>p</sub> = 0.5 µJ, 100 kHz, and *v* = 10 µm/s (*RR* is fixed to 400 Hz just for illustrative purpose). Time is in (a) linear and (b) log scale.</figcaption></figure>

The time for which crystallization is possible (as per the definition of efficiency be-
fore) is computed by solving the following equation: *SimpulseT*(*r*<sub>w</sub>, *x*, *R*<sub>*τ*</sub>) = *T*<sub>*g*</sub>. This gives the
temporal interval *x*1, *x*2, graphically defined in Figure 8a during the pulse period and for
which the temperature is above *T*<sub>*g*</sub> according to *r*<sub>w</sub>, i.e., the normalized distance to the beam
center.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** (a) Positions x1, x2 in the pulse period, between which the temperature is above T<sub>g</sub> (b) the efficiency (definition in the text) as a function of r<sub>w</sub>, x1 red curve, x2 dashed blue curve. The same laser parameters as in Figure 7 are used.</figcaption></figure>
At a given point in the material, the temperature is oscillating between Tmax(r,t) and Tmin(r,t) depending on Rτ and T₀. For contributing to the treatment of crystallization, the temperature must stay between Tg and Tm. Below Tg, the kinetics are too slow and above Tm the material is in the liquid state. When Rτ and T00 is such that T(r,t) is out of this interval, the treatment is inefficient. In this case, just a part of the T oscillation is active in contributing to crystallization treatment; this is defined in Figure 8b. As seen in Figure 8b, the efficiency of the thermal treatment (defined as x1-x2) is a function of the distance of the position away from the beam center (r<sub>w</sub>). This reduction in efficiency means that to reach crystallization around the periphery of the beam, the scanning speed must significantly be decreased. In fact, in those conditions, the thermal treatment time is much smaller than the irradiation time. Furthermore, beyond some departure from the center (e.g., here β > 0.7, i.e., in the direction perpendicular to the scanning direction, with the parameters used in Figure 7), the crystallization is no more possible for kinetics reasons (T < T<sub>g</sub>). r<sub>wdg</sub><sup>0</sup> defines the maximum width of the crystallization region.

(ii) Large *E*<sub>p</sub> such that the temperature is overcoming *T*<sub>m</sub>

In that case, a new treatment begins when a pulse oscillation is decreasing below Tm,
assuming that the material melted before.

The use of the *SimpulseT* function is valid. Using this expression, the position in the
period (x3 in Figure 9a) for which the temperature has overcome *T*<sub>m</sub> and returns below it
is detected. This position is the starting time for thermal treatment that ends at position
x1. But, if the next temperature pulse induces the temperature becoming greater than *T*<sub>m</sub>,
as the material is expected to melt and thus lose its memory, there will be a new time
origin to compute from the next pulse. So, the last pulse defines the beginning of the treat-
ment that goes on afterwards, i.e., when *T*<sub>max</sub> is decreasing below *T*<sub>m</sub> in the course of the
scanning. This is obtained by solving *SimpulseT*(r<sub>w</sub>,x,R<sub>τ</sub>) = *T*<sub>m</sub>. This yields positions x3 and
x4 in the period, as shown in Figure 9b.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Values of *x* for the crossing points according to the distance from the beam center *r*<sub>w</sub>. (a) Definition of the remarkable temporal locations in the period: *x*1 and *x*2, where the temperature equals that of the glass one; *x*3, *x*4, where it equals that of the melting one. (b) Variation of the remarkable temporal locations according to the distance from the beam center. The part of the heat-</figcaption></figure>
affected region where crystallization is possible is also indicated. This computation has been achieved with $E_p = 5$ µJ, $w = 1.65$ µm, $A = 0.2$, $RR = 60$ kHz, and $R_τ = 22$.

As seen in Figure 9, just a part of the pulse is “active” to induce crystallization ($x2-x1 < 1$) and the involvement of just a part of the temperature distribution is also possible. Following the picture in Figure 9b, only the tail of the temperature distribution is used. When $T$ is no longer overcoming $T_m$, and as shown in Figure 10 with $d = 1.92$ for $\beta = 0$, $x3$ and $x4$ do not exist anymore and, until $x1$ and $x2$ do, the temperature is larger than $T_g$.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Temperature evolution over time when pulses have overcome $T_m$, corresponding to a liquid-to-solid transition. The last pulse contribution that does not overcome $T_m$ defines the time origin. The computation is made with the parameters of Figure 9 with a scanning speed of 15 µm/s.</figcaption></figure>

From Figure 9a above, it can be noticed that the efficiency for crystallization is again smaller than 100% and that the writing speed must be decreased to trigger crystallization. Regarding the sides of the line ($\beta \neq 0$), the maximum of $T_{max}$ decreases with $y$ and the process described here will stop when $T_{max}$ will no longer increase above $T_m$ ($r_{wdg}^0 < \beta < r_{wdm}^0$). For these $\beta$, the situation is the same as the one described at the beginning of this section. The transition is solid-solid and starts from the front of the $T$ distribution, i.e., before the line center (in reference to the laboratory clock).

**6. Application to the Case of LNS Glass**

The next section is dedicated to applying the above findings to an experimental case using LNS glasses. The objective of this section is to show how the width of the thermally transformed region is linked to laser parameters, especially to the pulse energy. It is based on a real case for which the beam waist and the fraction of absorbed pulse energy ($A$ in $T_{oo}$ expression, Equation (A10) Appendix D) are crucial information, believed to be much smaller than shown in this paper. After data analysis from different authors, it is pointed out that the dependence on pulse energy of the beam waist $w$ is not very sensitive to the material when comparing silica and the LNS family (Appendix A). On the contrary, the pulse energy dependence of $A$ (Appendix B) is sensitive to the pulse duration and to the material. Therefore, the methodology in this section is to start with the pulse energy, the numerical aperture, and the pulse duration, to show that (1) the beam waist and the fraction of absorbed energy can be defined with enough accuracy for the purpose of this paper, and (2) then, with these quantities and their dependence on the laser parameters established reliably enough, to draw conclusions on the laser-induced crystallization process.

Let us start from a low pulse energy (0.5–1 µJ). The experimental half-width of the crystallized area is measured to be 1.65 µm (using 0.7 µJ in [4]) with the following laser parameters: wavelength 1030 nm, pulse duration 300 fs, $RR = 300$ kHz, $NA = 0.6$, and $v = 5$ µm/s [4]. With these values, solid-solid transformation (crystallization, ref. [16]) is shown to occur. On the other hand, the presence of nanogratings [17] indicates the size of the laser beam as these structures are only created by light. $w = 1.65$ µm is deduced.
However, *A* has to be adjusted. There are different processes that contribute to this fraction. The first one is the reflection of light (from Fresnel, or from the electron plasma mirror that is formed by the electron excitation in the conduction band or also from some light scattering). The second one is that, if a part is absorbed, a part is also transmitted as the absorption is never total in transparent materials with non-linear absorption. Since the modification is a solid-solid transformation, the temperature does not overcome *T*<sub>m</sub> (as shown in the previous sections), and the value of *A* has to be smaller than 0.09 for 0.7 µJ. This value falls within the values found in the literature and presented in Appendix B Figure A2. Once all of these parameters are known, the simulation can be carried out similarly to its presentation in the previous sections.

The results are summarized in Figure 11 and below:

- The possible crystallization would fall within a fraction of a second for v = 5 µm/s provided that the efficiency was 100%, but this is obviously not the case. On the contrary, the maximum efficiency is 22% here (see Figure 11b). The real treatment time is thus (x1-x2) multiplied by the time scale of Figure 11c. On the other hand, as crystallization is actually observed in LNS for v up to 125 µm/s [15] in a 0.5–1 µJ interval, with the following values: 1030 nm, NA = 0.6, τ<sub>RR</sub> = 250 fs, RR = 200 kHz, and w = 1.6 µm (for which the steady state is reached). This leads to a position where the crystallization nose is in the time range of 10 ms, i.e., an order of magnitude lower than the isothermal treatments in a conventional furnace, in a large volume without a stress field [18].

- The calculated crystallized half-width is close to that of the beam waist, with the criterion that the temperature should be above T<sub>g</sub>. It can be concluded that there is almost no line broadening originating from the thermal effect for this pulse energy as it is experimentally observed.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Temperature simulation in LNS glass for *E*<sub>p</sub> = 0.7 µJ, *A* = 0.09, *w* = 1.65 µm, *R*<sub>τ</sub> = 4.4, and scanning speed = 5 µm/s for (a–c): (a) *T*<sub>max</sub> and *T*<sub>min</sub> distribution, (b) the useful part of the T oscillations, (c) the corresponding treatment curve for *v* = 5 µm/s, and (d) the same with *v* = 125 µm/s.</figcaption></figure>

Now, and based on the above discussion, one can wonder how width evolves when
*E*<sub>p</sub> is increased. To answer this aspect, two distinct examples are considered, at low and
higher pulse energies.

**First example:** It is known that solid–solid transformation occurs until 1.0 µJ and that for higher pulse energy, a part of the crystallized region originates from a liquid-to-solid
transformation. For instance, let us consider the pulse energy 1.3 µJ. For this value, the half-width of the crystallized region is 3.73 µm. The laser beam radius is collected from [4] and reported in Figure A1 at 2.6 µm. $R_τ$ thus decreases from 4.4 to 1.8 with the pulse energy increase. The indirect variation of $R_τ$ with the pulse energy was not expected. It pushes the system to be in a heat accumulation regime even if $RR$ is constant. Additionally, $A$ increases itself to 0.15. Using these data, the computation is summarized in Figure 12.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Temperature computation in LNS glass for $Ep = 1.3$ µJ, $w = 2.6$ µm, $A = 0.15$, $R_τ = 1.8$. Laser parameters: $\lambda = 1030$ nm, $\tau_p = 300$ fs, $RR = 300$ kHz, NA = 0.6, and $v = 5$ µm/s. (a) T distribution, (b) efficiency, and (c) simulation on time.</figcaption></figure>

The provided results indicate the following:

- The efficiency is close to 100% due to the effect of heat accumulation that increases $T_{min}$ above $T_g$.

- The calculated crystallized radius is found to be 1.25 $w$. There is thus a small broadening of 25% compared to the beam waist $w$. There is thus a crystallized part that is not submitted to the beam where there is no submicrostructuring induced by the laser light. There is also a melted region smaller than the beam waist (half-width = 0.48 $w$). It is also noted that there is no obvious correlation between the beam width and the melted region.

For completeness, another experimental result available for 1.8 µJ has been computed and reported in Table 1. It shows reasonable agreement with the experimental results.

**Second example:** there is a much higher pulse energy, let us say 4.2 µJ, as used by Veenhuizen et al. [10], with laser parameters as follows: 175 fs, 200 kHz, and a varied scanning speed $v$. The authors of the aforementioned reference observed a crystallization full width of 31.3 µm for 15 µm/s and 9.82 µm for 75 µm/s, i.e., a significant variation. From Figure A1 in Appendix A, $w$ seems to be around 4.5 µm and thus $R_τ = 0.9$; $A$ can be seen to be 0.32 from Figure A2 in Appendix B. The results are shown in Figure 13.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Temperature computation in LNS glass for Ep = 4.2 µJ, w = 4.5 µm, A = 0.32, and Rτ = 0.9. Laser parameters: λ = 1030 nm, τp = 300 fs, RR = 200 kHz, NA = 0.6, and v = 5 µm/s. (a) Temperature distribution, (b) simulation on time for two positions β = 0 (center of the line) and β = 1.7 (periphery).</figcaption></figure>

Figure 13 allows a comparison of the treatment originating from the solid region and from the melted regions. The two curves penetrate the crystallization domain after a similar incubation time. However, in real time, the first one starts in front of the beam at a distance 3.4 w from its center, that is, before the temperature maximum, whereas the second starts 1.7 w after the temperature maximum. Consequently, the thermal treatment curve of the periphery is penetrating the crystallization domain before that of the center. The periphery is thus in a situation to stimulate the nucleation of the melted region. In addition, as both of the treatment curves have almost no oscillations, the efficiency is obviously 100%.

Let us now examine the effect of the scanning speed reported in [10]. A large variation in the width of the crystallization region was found (by a factor of three), whereas the speed varied by a factor of five. The origin of this variation is not primarily caused by the time shift of the treatment curves with the speed, as seen in Figure 5. The effect is quite small. Another possibility could be from the pulse number received punctually and that could be smaller than the quantity required for reaching the steady state. The steady state is reached after $N_{ss} = 7.6 \times 10^4$ pulses (from Equation (A4) Appendix A), whereas the number of pulses received punctually $N_p = 1.2 \times 10^4$ pulses for v = 15 µm/s. However, for v = 75 µm/s, $N_p = 2400$, which is significantly lower than $N_{ss}$ (>30 times less). Consequently, the steady state situation is not reached at those speeds. On the other hand, Figure 2 shows what the width at the beginning of the irradiation is like for a large $R_τ$, i.e. separated pulse contribution and increases in the course of irradiation due to heat accumulation. To compute the width in the transitory stage, it is assumed here that $T(r_{wdg}^0)$ crosses $T_g$ for $r_{wdg}^0 > 2$. The use of $T_{mean}$ is thus a good approximation. Therefore, the equation below has to be solved by deducing $r_{wdg}^0$ according to the number of pulses. It reads as follows (from Equation (A8) Appendix D):

$$ \bar{T}(r_{wdg}^0, R\tau, N) = \frac{\sqrt{\pi}}{R_\tau \cdot r_{wdg}^0} \cdot \left[ \operatorname{erf}(r_{wdg}^0) - \operatorname{erf}\left(\frac{r_{wdg}^0}{\sqrt{1 + N \cdot R\tau}}\right) \right] = \frac{T_g}{T_{00}} \quad \text{with} \quad T_{00} = \frac{A \cdot E_p}{\pi^{3/2} \rho C_p w^3} $$

The increase in the width as a function of the pulse number is plotted in Figure 14. There is a broadening with the number of pulses that arises from heat accumulation as $R_τ$ is around 1. $r_{wd}^0 = 3.5$ is found for the speed in consideration, i.e., a decrease from 3.8, the value at the steady state. This does not account for variation by a factor of three of the crystallization width. Even if the beginning of the irradiation is considered ($N_p=1$), the possible variation would only be by a factor of two. On the contrary, the results reported in [15] and used in Figure 11 are in agreement with this explanation (the laser parameters are as follows: $\lambda = 1030$ nm; RR = 200 kHz; NA = 0.6; pulse duration = 250 fs). In that case, $N_{ss}$ is around 50 whereas $N_p > 5280$ for a speed lower than 125 µm/s. The steady state is always reached and the width is weakly varying as experimentally observed (this explanation is given for Figure 5).

The question is thus: are there other reasons to change the crystallization width? The difference between the beginning and the steady state could be the change in absorption
due to the creation of defects in addition to the non-linear one that plays the role of trapping centers of excited electrons produced by previous pulses. It is worth noticing that, only a variation of A from 0.135 to 0.32 may explain that $r_{wdg}^0$ changes from 1.6 to 3.8.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** The half-width variation during the transitory stage (Np < Nss). The maximum width is 3.8 *w* (w = 4.5 µm) for Np > Ns (at steady state), for Ep = 4.2 µJ, Rτ = 0.9.</figcaption></figure>

On the other hand, it is observed in [10] that for a given pulse energy above some speed, the crystallization which occurs is only solid-to-solid transformation. This means that, according to the calculation, $T_{max}$ is no longer overcoming $T_m$. This is in agreement with lower absorption at a larger speed. Furthermore, the observation of the change in this turning point (here called $v_t$) with pulse energy leads one to consider that a particular dose $D_0$ is necessary for the absorption increase and for A to reach the value at the steady state. One would have, therefore, the following relation: $Ep.2w.RR/v_t = D_0$. So, it would be possible to have an increase in temperature disconnected from the heat accumulation but arising from an increase in absorption related to defect creation.

Finally, Table 1 shows that the steady state width of the crystallized region varies with the pulse energy for two reasons: because the beam waist increases with Ep, and the temperature distribution expands relatively with the beam waist, $R_τ$ decreases and the interaction is thus more and more sensitive to heat accumulation, and also because A increases with the pulse energy.

<table><caption>**Table 1.** Summary of experimental and computed width according to the pulse energy.</caption><thead><tr><th>E<sub>p</sub>(μJ)</th><th>Pulse duration τ<sub>p</sub> (fs)</th><th>w(μm)</th><th>A</th><th>R<sub>τ</sub></th><th>r<sub>w</sub>(Crystal) Computed</th><th>r<sub>w</sub>(Crystal) Experimental</th><th>r<sub>w</sub>(Melt) Computed</th><th>r<sub>w</sub>(Melt) Experimental</th><th>Ref.</th></tr></thead><tbody><tr><td>0.5</td><td>300</td><td>1.23</td><td>0.035</td><td>7.9</td><td>1</td><td>uk</td><td>na</td><td>Na</td><td>This work</td></tr><tr><td>0.7</td><td>300</td><td>1.65</td><td>0.09</td><td>4.4</td><td>1.05</td><td>1.00</td><td>na</td><td>Na</td><td>[4]</td></tr><tr><td>1</td><td>300</td><td>2.23</td><td>0.12</td><td>2.4</td><td>1.18</td><td>uk</td><td>na</td><td>Na</td><td>This work</td></tr><tr><td>1.3</td><td>300</td><td>2.6</td><td>0.15</td><td>1.8</td><td>1.47</td><td>1.6</td><td>0.48</td><td>1.1</td><td>This work</td></tr><tr><td>1.8</td><td>300</td><td>3.0</td><td>0.2</td><td>1.3</td><td>2.29</td><td>1.86</td><td>1.02</td><td>Uk</td><td>[4]</td></tr><tr><td>4.2</td><td>175</td><td>4.5</td><td>0.32</td><td>0.9</td><td>3.79</td><td>3.50</td><td>1.71</td><td>Uk</td><td>[10]</td></tr></tbody></table>

N.B.: na means not applicable and uk means unknown.

## 7. Discussion

Considering laser-induced crystallization, an analytical formula can be used to compute the temperature oscillations in time when the steady state is reached. These oscillations are limited between a maximum and minimum temperature, whatever the laser parameters and material properties, and established in the frame of a spherically symmetric
energy source [5]. The broadening effect when increasing the pulse energy or the pulse
density (inversely proportional to the beam scanning speed) was simulated.

Previously, it has been shown how the spatial distribution and specifically its FWHM
is changing during the transition to a steady state depending on the ratio between the
pulse period and the heat diffusion time (defined as Rτ). On the other hand, the tempera-
ture magnitude is defined by the beam waist and the fraction of absorbed incident energy
(A). Thus, when applied to a real case and a reliable simulation is requested, it is necessary
to consider 1) an increase in the beam waist (w) with the pulse energy, and 2) the increase
in the A factor with the pulse energy and the pulse duration. Such constraints that origi-
nate from a non-linear effect in the solid can, for the most part, only currently be extracted
from the experimental results available in the literature (see Appendixes A and B) as
mechanisms in the case of fs laser absorption are not clear enough for a robust simulation.

In such a way, within the validity of the used approximation, the following is de-
duced:

(1) The crystallization management is facilitated in a heat accumulation regime (Rτ << 10). This allows the temperature oscillations to be quite small compared to average temperature values; thus, the laser parameters (i.e., Ep, R) can be easily adjusted to be completely above T<sub>g</sub> and below T<sub>m</sub> during the thermal treatment time that is starting when the temperature is rising above T<sub>g</sub>. Thermal treatment curves are computed. The scanning speed is also straightforward for the thermal curve to efficiently penetrate the crystallization domain. The crystallization process will originate from a solid-solid transition, and the width of the crystallized area is close to the beam waist. For the inorganic oxides, the thermal diffusivity is usually around 10<sup>-6</sup> m<sup>2</sup>/s (within a factor of two) so the diffusion time from a source with a beam waist of 1 µm is around <0.5 µs. The above conditions apply when the repetition rate frequency is much larger than 0.4-0.2 Mhz. However, if the focusing is not so tight, let us say w = 2 µm, the heat diffusion time is four times larger and reaches 1-2 µs. So, the repetition rate has to be much larger than 5-10 MHz. However, in organic materials for which the diffusion time is 10-20 times longer, the heat accumulation regime is more easily reached and the thermal width more easily predictable.

If the pulse energy is increased such that the temperature is becoming larger than the
melting temperature around the center of the written line, the process changes; the start-
ing time of the incubation is no longer when the temperature becomes larger than T<sub>g</sub> but
when T decreases below T<sub>m</sub>, as the process is a liquid-to-solid transformation. The crystal-
lization domain is penetrated from the growth side, which renders the process sensitive
to external stimulation, in particular to crystallization induction from the periphery where
the temperature does not overcome Tm. When departing from the center, the temperature
maximum of the treatment curve decreases below T<sub>m</sub>; a solid-to-solid transformation is
occurring again. There are two regions in the crystallized line. It is shown that the solid-
to-solid transformation region crystallizes sooner than the melted one and thus is a source
of nucleation for the melted one. The size of the crystallized region becomes larger than
the beam waist in that case. The texture orientation of the line (if any) is thus defined by
the peripheral region outside of the beam and so defined mainly by the temperature gra-
dient orientation. The large broadening with the pulse energy is defined in this case
mainly by the base of the temperature distribution (which follows a hyperbolic shape
(T ~ 1/r<sub>w</sub>), although the beam waist and absorption fraction also give rise to a broadening
contribution. We also pointed out an unexpected result: the stimulation of the heat accu-
mulation regime on a pulse energy increase and not only on RR.

(2) On the other hand, when heat accumulation is small or negligible (Rτ > 10, small RR, large diffusion time), the above processes are the same according to the pulse energy, but the temperature oscillates significantly. As a consequence, just a small part of the oscillation is efficient (a few or ten %). This contributes to a proportional decrease in the scanning speed in addition to a repetition rate that could here be smaller.

However, the weak efficiency of the growth rate and a moderate pulse energy allows
for the penetration of the crystallization on the growth-limited side of the crystallization
domain, and may allow for the production of independent nanocrystals that have been revealed to be orientable with laser polarization in the case of LNS glasses [19].

(3) An absorption increase may be observed in the course of the irradiation and produce a large increase in temperature, changing the crystallization process (from solid/solid to liquid/solid). This effect could be dose dependent ($D_0 = E_p \cdot 2w \cdot RR/v_t$) and leads to a transition speed $v_t$ varying with the pulse energy. We suggest that this phenomenon arises from the formation of defects that are more easily ionized.

## 8. Conclusions

Using analytical formulae computed from very simplified approximations (spherical approximation), it is possible to simulate crystallization processes by irradiation with a fs laser. However, to obtain compatible results, the beam waist and the light absorption fraction have to be approximately known apart from the physico-chemical data of the glass. After that, the management of the crystallization process and the width of the crystallization line can be obtained according to the pulse energy, e.g., crystallite size, and also the effect of the scanning speed can be understood. For example, in particular that of LNS glass, the crystallization process may change with the scanning speed.

Lastly, the application of the temperature calculation can be extended to other processes depending on the temperature: the distribution of the fictive temperature, the viscoplastic effect (stress distribution), any chemical reaction that is thermally activated, chemical migration, the maximum energy for not erasing a type X or type II modification, the creation of luminescent defects, or the decomposition temperature of organic material. This management shall be even more accurate when the fs laser is applied to organic matter in order to avoid complete material decomposition.

**Author Contributions:** Conceptualization, B.P.; funding acquisition, M.L.; investigation, R.Q. and B.P.; methodology, B.P.; project administration, B.P. and M.L.; resources, M.L.; supervision, B.P.; validation, R.Q.; visualization, M.C.; writing—original draft preparation, B.P.; writing—review and editing, R.Q., B.P., M.C., and M.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Agence Nationale de la Recherche (ANR), FLAG-IR Project, award number ANR-18-CE08-0004-01, and REFRACTEMP project, award number ANR-22-CE08-0001-01.

**Data Availability Statement:** The original contributions presented in the study are included in the article material; further inquiries can be directed to the corresponding authors.

**Acknowledgments:** R.Q. acknowledges the China Scholarship Council (CSC) for the funding (No. 201808440317) of her PhD fellowship.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## Appendix A. The Beam Width from the Literature and from Our Experiments

For this, the width of modification lines achieved by a non-thermal effect is used. The most obvious one is a line width of type II nanogratings as their structure is produced only by the laser light, whatever the material’s chemical composition.

Considering first the beam focused with a length with NA and the effect of self-focusing, the formula for the effective width at the focus in the material is given by Marburger et al. [20], Ashcom et al. [21], and Schaffer et al. [22], but their formula below does not take into account the defocusing effect due to electron excitation. The defocusing effect is dependent on the excited electron density during the pulse; thus, a too-short lifetime or a weak pulse energy in the excited state of the electron may lead to the absence of this effect. The effective beam waist radius is given below:

$$w_{eff} = \sqrt{\frac{\lambda^2/0.9\pi}{\frac{2\pi NA^2}{1-NA^2} + \frac{I\lambda^2}{P_{cr}}} \quad \text{with} \quad P_{cr} = \frac{3.77\lambda^2}{8\pi n_0 n_2} \quad n_2 \text{ being the non-linear refractive index (3.2} \cdot 10^{-16} \text{ cm}^2/\text{W [23,24])}, n_0 \text{ the relative one, I the light intensity.}$$
Increasing the pulse energy, the intensity increases and thus the width decreases due to non-linear effects. It is the contrary of what it is observed in the case of the laser leading to electron excitations [25].

In Lei et al. [26], the width increases with the energy above 0.2 µJ, although the electron density is of the order of 0.01 excited electrons.nm⁻³. Even for 0.1 excited electrons.nm⁻³, it is said that there is a weak defocusing effect and the half-width is 0.4 µm. However, without the defocusing effect, the width should decrease with a pulse energy increase.

These results and others are collected in Table A1 below.

<table><caption>**Table A1.** Some half-widths according to pulse energy, numerical aperture, and materials.</caption><thead><tr><th>Wavelength (nm)</th><th>Pulse Energy (µJ)</th><th>NA</th><th>Half Width (µm)</th><th>Material</th><th>Ref.</th></tr></thead><tbody><tr><td>1030</td><td>0.1</td><td>0.6</td><td>0.31</td><td>silica</td><td>Ashcom et al. [21]</td></tr><tr><td>1030</td><td>0.5</td><td>0.6</td><td>0.24</td><td>silica</td><td>idem</td></tr><tr><td>1030</td><td>1</td><td>0.6</td><td>0.1</td><td>silica</td><td>idem</td></tr><tr><td>1030</td><td>0.19</td><td>0.3</td><td>0.4</td><td>silica</td><td>Lei et al. Figure S2 of [27]</td></tr><tr><td>1030</td><td>0.55</td><td>0.3</td><td>1.25</td><td>silica</td><td>idem</td></tr><tr><td>1025</td><td>0.7</td><td>Not known</td><td>1.</td><td></td><td>Beresna et al. [28] (in air) Figure 4</td></tr><tr><td>1030</td><td>1</td><td>0.6</td><td>2</td><td>silica</td><td>Qiong et al. [29]</td></tr><tr><td>1030</td><td>0.5</td><td>0.6</td><td>1.3</td><td>G1?</td><td>idem</td></tr><tr><td>1030</td><td></td><td>0.6</td><td>0.7</td><td>B33</td><td>idem</td></tr><tr><td>1030</td><td></td><td>0.6</td><td>1.65</td><td>ULE</td><td>idem</td></tr><tr><td>1030</td><td>3</td><td>0.6</td><td>2.1</td><td>ULE</td><td>idem</td></tr><tr><td>1030</td><td>1</td><td>0.6</td><td>1.25</td><td>BGG-ZnO</td><td>idem</td></tr><tr><td>1030</td><td>2</td><td>0.6</td><td>2.25</td><td>BGG-ZnO</td><td>idem</td></tr><tr><td>1030</td><td>0.7</td><td>0.6</td><td>1.65</td><td>LNS</td><td>Cao et al. [4]</td></tr><tr><td>1030</td><td>1.3</td><td>0.6</td><td>2.6</td><td>LNS</td><td>idem</td></tr><tr><td>1030</td><td>0.5</td><td>0.6</td><td>1.6</td><td>LNSB</td><td>Muzi et al. [15]</td></tr><tr><td>1030</td><td>0.2</td><td>0.6</td><td>0.5</td><td>Suprasil</td><td>Our work</td></tr><tr><td>1030</td><td>5</td><td>0.6</td><td>5</td><td>Suprasil</td><td>Our work</td></tr><tr><td>1030</td><td>0.15</td><td>0.6</td><td>0.45</td><td>Suprasil</td><td>Our work</td></tr><tr><td>1030</td><td>2</td><td>0.6</td><td>1.75</td><td>Suprasil</td><td>Our work</td></tr><tr><td>800</td><td>0.18</td><td>1.35</td><td>0.4</td><td>glasses</td><td>Gamaly et al. [30] Figure 3</td></tr><tr><td>1030</td><td>0.8</td><td>0.6</td><td>2</td><td>Zn BGG</td><td>Yao et al. [31] Figure 2</td></tr><tr><td>1030</td><td>0.19</td><td>0.3</td><td>0.4</td><td>silica</td><td>Figure S2 of [27]</td></tr><tr><td>1030</td><td>0.34</td><td>0.23</td><td>0.65</td><td>silica</td><td>idem</td></tr><tr><td>1030</td><td>0.55</td><td>0.16</td><td>2.2</td><td>silica</td><td>idem</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>**Figure A1.** The half beam width according to the pulse energy from Table A1 for various compounds, $\lambda$ = NIR, various NA, often at the steady state. The dots are experimental measurements, the line is just a guide for the eye.</figcaption></figure>

## Appendix B. Fraction of Absorbed Beam Energy

The incident beam intensity from the laser is a parameter that is dependent on the material properties as well as the beam propagation, and especially on the focusing strength. Indeed, for working in volume, materials do not absorb linearly but non-linearly with an exponent of the intensity light that depends on the materials and on the pulse duration because it depends on the excited electron density that occurs in the few 10 s fs [6]. The fraction of incident light absorbed is thus dependent on the intensity at the focus and around (before and after). In addition, the excited electron density may behave like a mirror and reflect a part of the incident intensity [32]. Moreover, for the largest incident power, the excited electron density may have a defocusing effect that decreases the intensity and thus the non-linear absorption[25]. All of these processes are intricate and complex and lead to difficult prediction without heavy computations. Needing an estimate of this quantity for our computation, various results are collected in the Table A2 below.

<table><caption>**Table A2.** Fraction of absorbed energy according to pulse energy, pulse duration, numerical aperture of the focusing lens, and materials.</caption><thead><tr><th>Wavelength (nm)</th><th>Pulse Energy (μJ)</th><th>Pulse Duration (fs)</th><th>NA</th><th>Fraction</th><th>Material</th><th>Ref.</th></tr></thead><tbody><tr><td>800</td><td>0.1–0.7</td><td>50</td><td>0.5</td><td>0.2</td><td>silica</td><td>Couairon et al. [33]</td></tr><tr><td>800</td><td>1</td><td>90</td><td>Not mentioned</td><td>0.35</td><td>silica</td><td>Wu et al. [34]</td></tr><tr><td>800</td><td>2</td><td>90</td><td>Not mentioned</td><td>0.45</td><td>silica</td><td>idem</td></tr><tr><td>800</td><td>3</td><td>90</td><td>Not mentioned</td><td>0.54</td><td>silica</td><td>idem</td></tr><tr><td>800</td><td>0.06–0.12</td><td>100</td><td>Not mentioned</td><td>&lt;0.12</td><td>metallic surface</td><td>Meng et al. [35] metal surface</td></tr><tr><td>800</td><td>0.25</td><td>100</td><td>Not mentioned</td><td>0.25</td><td>silica</td><td>Zhang et al. [36]</td></tr><tr><td>800</td><td>1/7.4</td><td>100</td><td>Not mentioned</td><td>&lt;0.1</td><td>silica</td><td>idem</td></tr><tr><td>1026</td><td>4.2</td><td>175</td><td></td><td>0.58</td><td>silica</td><td>Computed from Veen-huisen et al. [10], see text</td></tr><tr><td>1030</td><td>0.6</td><td>190</td><td>0.16</td><td>0.25</td><td>silica</td><td>Wang et al. [37]</td></tr><tr><td>1030</td><td>1</td><td>190</td><td>0.16</td><td>0.36</td><td>silica</td><td>idem</td></tr><tr><td>1030</td><td>0.2</td><td>300</td><td>0.3</td><td>0.04 to 0.06</td><td>silica from 1 to 30 pulses</td><td>Wang et al. [37]</td></tr></tbody></table>
<table><tr><td>1030</td><td>0.6</td><td>300</td><td>0.16</td><td>0.05</td><td>silica</td><td>idem</td></tr><tr><td>1030</td><td>0.69</td><td>300</td><td>0.16</td><td>&lt;0.15</td><td>silica</td><td>idem</td></tr><tr><td>1030</td><td>1</td><td>300</td><td>0.16</td><td>0.25</td><td>silica</td><td>idem]</td></tr><tr><td>1030</td><td>0.5</td><td>300</td><td>0.6</td><td>0.035</td><td>LNS</td><td>This work</td></tr><tr><td>1030</td><td>0.7</td><td>300</td><td>0.6</td><td>0.09</td><td>LNS</td><td>This work</td></tr><tr><td>1030</td><td>1.3</td><td>300</td><td>0.6</td><td>0.16</td><td>LNS</td><td>This work</td></tr><tr><td>1030</td><td>1.8</td><td>300</td><td>0.6</td><td>0.2</td><td>LNS</td><td>This work</td></tr><tr><td>1026</td><td>4.2</td><td>175</td><td>0.6</td><td>0.32</td><td>LNS</td><td>This work</td></tr><tr><td>1030</td><td>0.7</td><td>500</td><td>0.16</td><td>0.15</td><td>silica</td><td>Wang et al. [37]</td></tr><tr><td>1030</td><td>1</td><td>500</td><td>0.16</td><td>0.19</td><td>silica</td><td>idem</td></tr><tr><td>1030</td><td>1</td><td>700</td><td>0.16</td><td>0.1</td><td>silica</td><td>idem</td></tr></table>

<figure><img src="image_19.png" /><figcaption>**Figure A2.** Plot of the absorbed energy fraction versus the pulse energy according to pulse duration (fs). The colors of the points make reference to Table A2: the blue curve is for silica and the red one is for LNS glasses. However, lines are just guides for the eyes.</figcaption></figure>

The absorbed energy fraction is very sensitive to the pulse duration for <3.7 microJ but is less so for higher values (Zhang et al. [36]).

## Appendix C

<table><caption>**Table A3.** Physico-chemical and laser parameters.</caption><thead><tr><th></th><th>LNS Glass: 33Li<sub>2</sub>O-<br>33Nb<sub>2</sub>O<sub>5</sub>- 34SiO<sub>2</sub> (mol %)</th><th>Ref.</th></tr></thead><tbody><tr><td>Specific density (ρ, kg/m<sup>3</sup>)</td><td>3830</td><td>Interpolated between SiO<sub>2</sub> and NbLiO<sub>3</sub>.</td></tr><tr><td>Heat capacity (C<sub>p</sub>, J/kg.K)</td><td>787</td><td>idem</td></tr><tr><td>K = D<sub>th</sub>.p.C<sub>p</sub> (J/K.m.s)</td><td>2.65</td><td>Computed from diffusivity</td></tr><tr><td>D<sub>th</sub> (m<sup>2</sup>/s)</td><td>9.1 10<sup>-7</sup></td><td>[38,39]</td></tr><tr><td>T<sub>g</sub> (°C)</td><td>579</td><td>[40,41]</td></tr><tr><td>T<sub>m</sub> (°C)</td><td>1257</td><td>[40,41]</td></tr></tbody></table>

## Appendix D

After reaching the steady state (the number of pulses for that is given in Table A3 of [5]), the expressions of the temperature maximum, minimum and oscillations are as follows:
$$ T_{max}(r_w, R_τ) = \frac{\exp\left[-\frac{(r_w)^2}{1+x_m \cdot R_τ}\right]}{\left[1+x_m \cdot R_τ\right]^3} + \frac{\exp\left[-\frac{(r_w)^2}{1+(1+x_m) \cdot R_τ}\right]}{2\left[1+(1+x_m) \cdot R_τ\right]^3} + \frac{\sqrt{\pi}}{R_τ \cdot r_w} \left\{ \operatorname{erf}\left[\frac{r_w}{\sqrt{1+(1+x_m) \cdot R_τ}}\right] \right\} \quad (A1) $$

Where $x_m = \frac{\sqrt{R_τ} \sqrt{9R_τ + 32(r_w)^2 - 3R_τ - 8}}{8R_τ}$ and where $r_w > \sqrt{\frac{3}{2} + \frac{2}{R_τ}}$ ($x_m$ limited to 1) or $x_m = 0$, when $r_w < \sqrt{\frac{3}{2} + \frac{2}{R_τ}}$ for whatever value of $R_τ$ or $R_τ$ less than $\frac{2}{r_w^2 - 1.5}$ when $r_w^2 > 1.5$).

$$ T_{osc}(r_w, R_τ) \approx \frac{\exp\left[-\frac{(r_w)^2}{1+x_m \cdot R_τ}\right]}{\left[1+x_m \cdot R_τ\right]^3} + \frac{\sqrt{\pi}}{R_τ \cdot r_w} \left\{ \operatorname{erf}\left[\frac{r_w}{\sqrt{1+(1+x_m) \cdot R_τ}}\right] \right\} - \frac{\sqrt{\pi}}{R_τ \cdot r_w} \operatorname{erf}\left(\frac{r_w}{\sqrt{1+R_τ}}\right) \quad (A2) $$

$$ T_{min}(r_w, R_τ) = \frac{\exp\left[-\frac{(r_w)^2}{1+(1+x_m) \cdot R_τ}\right]}{2\left[1+(1+x_m) \cdot R_τ\right]^3} + \frac{\sqrt{\pi}}{R_τ \cdot r_w} \left\{ \operatorname{erf}\left[\frac{r_w}{\sqrt{1+R_τ}}\right] \right\} \quad (A3) $$

$$ N_{ssmax}^r > \frac{1}{R_τ} \left[ \left( \frac{2}{R_τ \cdot \varepsilon \cdot T_{max}(r_w, \infty)} \right)^2 - 1 \right] \quad N_{ssmin}^r > \frac{1}{R_τ} \left[ \left( \frac{2}{R_τ \cdot \varepsilon \cdot T_{min}(r_w, \infty)} \right)^2 - 1 \right] \quad (A4) $$

$$ T_{max}(r_w, R_τ, N) \approx \frac{\exp\left[-\frac{(r_w)^2}{1+x_m \cdot R_τ}\right]}{\left[1+x_m \cdot R_τ\right]^3} + \frac{\exp\left[-\frac{(r_w)^2}{1+(1+x_m) \cdot R_τ}\right]}{2\left[1+(1+x_m) \cdot R_τ\right]^3} + \frac{\exp\left[-\frac{(r_w)^2}{1+(N-1+x_m) \cdot R_τ}\right]}{2\left[1+(N-1+x_m) \cdot R_τ\right]^3} + \frac{\sqrt{\pi}}{R_τ \cdot r_w} \left\{ \operatorname{erf}\left[\frac{r_w}{\sqrt{1+(1+x_m) \cdot R_τ}}\right] - \operatorname{erf}\left[\frac{r_w}{\sqrt{1+(N-1+x_m) \cdot R_τ}}\right] \right\} \quad (A5) $$

When the mean T can be used, $\bar{T}(r, N) = \frac{1}{\tau_{RR}} \int_{N}^{\tau_{RR}} T(r, t) dt$; the formula is even simpler to use after reaching the steady state. $\tau_{RR}$ is the period = 1/RR

$$ \bar{T}(r_w, R_τ, \infty) = \frac{\sqrt{\pi}}{R_τ \cdot r_w} \operatorname{erf}(r_w) \quad (A6) $$

$$ N_{ssm}^r(r_w) = \frac{1}{R_τ} \left[ \left( \frac{2 \cdot r_w}{\sqrt{\pi \cdot \varepsilon \cdot \operatorname{erf}(r_w)}} \right)^2 - 1 \right] \quad (A7) $$

$$ \bar{T}(r_w, R_τ, N) = \frac{\sqrt{\pi}}{R_τ \cdot r_w} \cdot \left[ \operatorname{erf}(r_w) - \operatorname{erf}\left(\frac{r_w}{\sqrt{1+N \cdot R_τ}}\right) \right] \quad (A8) $$

The formula built for the simulation of the temperature oscillations at the steady state:

$$ \text{SimpulseT}(r_w, x, R_τ, E_p) = [T_{min}(r_w, R_τ) + T_{osc}'(r_w, x, R_τ)] \cdot T_{00}(E_p) \quad (A9) $$

where $T_{osc}'(r_w, x, R_τ)$ is just the function $T_{osc}(r_w, R_τ)$ where $x_m$ is replaced by x.

The maximum temperature induced by one pulse is $T_{00} = \frac{A \cdot E_p}{\frac{3}{\pi^2} \rho c_p w^3}$ with A the fraction of absorbed energy, (A10) $\rho$ and $c_p$ the glass density and heat capacity, $w$ the beam waist at 1/e, and $E_p$ the pulse energy.

Note that the time is introduced in the formula by $r_{wd}$ which is a "moving" $r_w$ (as the beam is scanning) and in x by rendering it periodic, i.e., in taking only the floating value of $t.RR$.

## References

1. Komatsu, T.; Honma, T. Laser patterning and characterization of optical active crystals in glasses. *J. Asian Ceram. Soc.* **2013**, *1*, 9–16. https://doi.org/10.1016/j.jasper.2013.02.006.

2. Neuville, D.; Cormier, L.; Caurant, D.; Montagne, L. *From Glass to Crystal. Nucleation, Growth and Phase Separation: From Research to Applications*; EDP Sciences: Paris, France, 2017; p. 640.

3. Teng, Y.; Zhou, J.; Sharafudeen, K.; Zhou, S.; Miura, K.; Qiu, J. Space-selective crystallization of glass induced by femtosecond laser irradiation. *J. Non-Cryst. Solids* **2014**, 383, 91–96. https://doi.org/10.1016/j.jnoncrystol.2013.04.015.

4. Cao, J.; Lancry, M.; Brisset, F.; Mazerolles, L.; Saint-Martin, R.; Poumellec, B. Femtosecond Laser-Induced Crystallization in Glasses: Growth Dynamics for Orientable Nanostructure and Nanocrystallization. *Cryst. Growth Des.* **2019**, *19*, 2189–2205. https://doi.org/10.1021/acs.cgd.8b01802.

5. Que, R.; Lancry, M.; Poumellec, B. Usable Analytical Expressions for Temperature Distribution Induced by Ultrafast Laser Pulses in Dielectric Solids. *Micromachines* **2024**, *15*, 196. https://doi.org/10.3390/mi15020196.
6. Lancry, M.; Groothoff, N.; Poumellec, B.; Guizard, S.; Fedorov, N.; Canning, J. Time-resolved plasma measurements in Ge-doped silica exposed to infrared femtosecond laser. *Phys. Rev. B* **2011**, *84*, 245103. https://doi.org/10.1103/PhysRevB.84.245103.

7. Uhlmann, D.R. A kinetic treatment of glass formation. *J. Non-Cryst. Solids* **1972**, *7*, 337–348. https://doi.org/10.1016/0022-3093(72)90269-4.

8. Weinberg, M.C.; Birnie, D.P.; Shneidman, V.A. Crystallization kinetics and the JMAK equation. *J. Non-Cryst. Solids* **1997**, *219*, 89–99. https://doi.org/10.1016/S0022-3093(97)00261-5.

9. Deubener, J.; Allix, M.; Davis, M.J.; Duran, A.; Höhe, T.; Honma, T.; Komatsu, T.; Krüger, S.; Mitra, I.; Müller, R.; et al. Updated definition of glass-ceramics. *J. Non-Cryst. Solids* **2018**, *501*, 3–10. https://doi.org/10.1016/j.jnoncrysol.2018.01.033.

10. Veenhuizen, K.; McAnany, S.; Nolan, D.; Aitken, B.; Dierolf, V.; Jain, H. Fabrication of graded index single crystal in glass. *Sci. Rep.* **2017**, *7*, 44327. https://doi.org/10.1038/srep44327.

11. Yinnon, H.; Uhlmann, D.R. Glass Science and Technology; Uhlmann, D.R., Kreidl, N.J., Eds.; Academic Press: New York, NY, USA, 1983; Volume 1, p. 465.

12. Miyamoto, I.; Horn, A.; Gottmann, J.; Wortmann, D.; Yoshino, F. Fusion Welding of Glass Using Femtosecond Laser Pulses with High-repetition Rates. *J. Laser Micro Nanoeng.* **2007**, *2*, 57–63.

13. Beresna, M.; Gertus, T.; Tomašiūnas, R.; Misawa, H.; Juodkazis, S. Three-Dimensional Modeling of the Heat-Affected Zone in Laser Machining. *Appl. Laser Chem.* **2008**, *2008*, 976205.

14. Miyamoto, I.; Horn, A.; Gottmann, J. Local Melting of Glass Material and Its Application to Direct Fusion Welding by Ps-laser Pulses. *J. Laser Micro Nanoeng.* **2007**, *2*, 7–14.

15. Muzi, E.; Cavillon, M.; Lancry, M.; Brisset, F.; Que, R.; Pugliese, D.; Janner, D.; Poumellec, B. Towards a Rationalization of Ultrafast Laser-Induced Crystallization in Lithium Niobium Borosilicate Glasses: The Key Role of The Scanning Speed. *Crystals* **2021**, *11*, 290. https://doi.org/10.3390/cryst11030290.

16. Cost246 Reliability of fiber Bragg gratings 1997

17. Bricchi, E.; Klappauf, B.G.; Kazansky, P.G. Form birefringence and negative index change created by femtosecond direct writing in transparent materials *Opt. Lett.* **2004**, *29*, 119–121. 10.1364/OL.29.000119

18. Vigouroux, H. Etude de Vitrocéramiques Optiques pour le Doublement de Fréquence. Doctoral Dissertation, University of Bordeaux, Bordeaux, France, 2012.

19. Cao, J.; Matthieu, L.; François, B.; Bertrand, P. Orientable Nonlinear Optical Crystals and Periodic Nanostructure by Femtosecond Laser Irradiation. *Laser Optoelectron. Prog.* **2022**, *59*, 1516001. https://doi.org/10.3788/lop202259.1516001.

20. Marburger, J.H. Self-focusing: Theory. *Prog. Quantum Electron.* **1975**, *4* Pt 1, 35–110. https://doi.org/10.1016/0079-6727(75)90003-8.

21. Ashcom, J.B.; Gattass, R.R.; Schaffer, C.B.; Mazur, E. Numerical aperture dependence of damage and supercontinuum generation from femtosecond laser pulses in bulk fused silica. *J. Opt. Soc. Am. B* **2006**, *23*, 2317–2322. https://doi.org/10.1364/JOSAB.23.002317.

22. Schaffer, C.B.; Brodeur, A.; Mazur, E. Laser-induced breakdown and damage in bulk transparent materials induced by tightly focused femtosecond laser pulses. *Meas. Sci. Technol.* **2001**, *12*, 1784. https://doi.org/10.1088/0957-0233/12/11/305.

23. Kato, T.; Suetsugu, Y.; Nishimura, M. Estimation of non-linear refractive index in various silica-based glasses for optical fibers. *Opt. Lett.* **1995**, *20*, 2279–2281.

24. Le Boiteux, S.; Segonds, P.; Canioni, L.; Sarger, L.; Cardinal, T.; Duchesne, C.; Fargin, E.; Le Flem, G. Nonlinear optical properties for TiO₂ containing phosphate, borophosphate, and silicate glasses. *J. Appl. Phys.* **1997**, *81*, 1481–1487.

25. Couairon, A.; Mysyrowicz, A. Femtosecond filamentation in transparent media. *Phys. Rep.* **2007**, *441*, 47–189. https://doi.org/10.1016/j.physrep.2006.12.005.

26. Lei, Y. Laser Modification of Materials at Micro- and Nanoscale for Photonics and Information Technology. Doctoral Dissertation, University of Southampton, Southampton, UK, 2023.

27. Wang, H.; Lei, Y.; Shayeganrad, G.; Svirko, Y.; Kazansky, P.G. Increasing Efficiency of Ultrafast Laser Writing via Nonlocality of Light-Matter Interaction. *Laser Photonics Rev.* **2024**, 2301143. https://doi.org/10.1002/lpor.202301143.

28. Beresna, M.; Kazansky, P.G.; Svirko, Y.; Barkauskas, M.; Danielius, R. High average power second harmonic generation in air. *Appl. Phys. Lett.* **2009**, *95*, 121502. https://doi.org/10.1063/1.3232235.

29. Xie, Q.; Cavillon, M.; Pugliese, D.; Janner, D.; Poumellec, B.; Lancry, M. On the Formation of Nanogratings in Commercial Oxide Glasses by Femtosecond Laser Direct Writing. *Nanomaterials* **2022**, *12*, 2986. https://doi.org/10.3390/nano12172986.

30. Gamaly, E.G.; Juodkazis, S.; Nishimura, K.; Misawa, H.; Luther-Davies, B.; Hallo, L.; Nicolai, P.; Tikhonchuk, V.T. Laser-matter interaction in the bulk of a transparent solid: Confined microexplosion and void formation. *Phys. Rev. B* **2006**, *73*, 214101.

31. Zaiter, R.; Lancry, M.; Fargues, A.; Adamietz, F.; Dussauze, M.; Rodriguez, V.; Poumellec, B.; Cardinal, T. Optical and structural characterization of femtosecond laser written micro-structures in germanate glass. *Sci. Rep.* **2023**, *13*, 11050. https://doi.org/10.1038/s41598-023-35730-3.

32. Garcia-Lechuga, M.; Haahr-Lillevang, L.; Siegel, J.; Balling, P.; Guizard, S.; Solis, J. Simultaneous time-space resolved reflectivity and interferometric measurements of dielectrics excited with femtosecond laser pulses. *Phys. Rev. B* **2017**, *95*, 214114. https://doi.org/10.1103/PhysRevB.95.214114.

33. Couairon, A.; Sudrie, L.; Franco, M.; Prade, B.; Mysyrowicz, A. Filamentation and damage in fused silica induced by tightly focused femtosecond laser pulses. *Phys. Rev. B* **2005**, *71*, 125435. https://doi.org/10.1103/PhysRevB.71.125435.

34. Wu, A.Q.; Chowdhury, I.H.; Xu, X. Femtosecond laser absorption in fused silica: Numerical and experimental investigation. *Phys. Rev. B* **2005**, *72*, 085128. https://doi.org/10.1103/PhysRevB.72.085128.
35. Meng, Y.; Gong, A.; Chen, Z.; Wang, Q.; Guo, J.; Li, Z.; Li, J. Atomistic-Continuum Study of an Ultrafast Melting Process Controlled by a Femtosecond Laser-Pulse Train. *Materials* **2023**, *17*, 185. https://doi.org/10.3390/ma17010185

36. Zhang, S.; Menoni, C.; Gruzdev, V.; Chowdhury, E. Ultrafast Laser Material Damage Simulation—A New Look at an Old Problem. *Nanomaterials* **2022**, *12*, 1259. https://doi.org/10.3390/nano12081259

37. Wang, H.; Lei, Y.; Wang, L.; Sakakura, M.; Yu, Y.; Shayeganrad, G.; Kazansky, P.G. 100-Layer Error-Free 5D Optical Data Storage by Ultrafast Laser Nanostructuring in Glass. *Laser Photonics Rev.* **2022**, *16*, 2100563. https://doi.org/10.1002/lpor.202100563

38. Kingery, W.D. Thermal Conductivity: XII, Temperature Dependence of Conductivity for Single-Phase Ceramics. *J. Am. Ceram. Soc.* **1955**, *38*, 251–255. https://doi.org/10.1111/j.1151-2916.1955.tb14940.x

39. Morgan, R.A.; Kang, K.I.; Hsu, C.C.; Koliopoulos, C.L.; Peyghambarian, N. Measurement of the thermal diffusivity of nonlinear anisotropic crystals using optical interferometry. *Appl. Opt.* **1987**, *26*, 5266–5271. https://doi.org/10.1364/AO.26.005266

40. Sigaev, V.N.; Golubev, N.V.; Stefanovich, S.Y.; Komatsu, T.; Benino, Y.; Pernice, P.; Aronne, A.; Fanelli, E.; Champagnon, B.; Califano, V.; et al. Second-order optical non-linearity initiated in Li₂O-Nb₂O₅-SiO₂ and Li₂O-ZnO-Nb₂O₅-SiO₂ glasses by formation of polar and centrosymmetric nanostructures. *J. Non-Cryst. Solids* **2008**, *354*, 873–881. https://doi.org/10.1016/j.jnoncrysol.2007.08.033

41. Shimada, M.; Honma, T.; Komatsu, T. Laser patterning of oriented LiNbO₃ crystal particle arrays in NiO-doped lithium niobium silicate glasses. *Int. J. Appl. Glass Sci.* **2018**, *9*, 518–529. https://doi.org/10.1111/ijag.12390

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.