# Usable Analytical Expressions for Temperature Distribution Induced by Ultrafast Laser Pulses in Dielectric Solids

Ruyue Que, Matthieu Lancry, Bertrand Poumellec

► To cite this version:

Ruyue Que, Matthieu Lancry, Bertrand Poumellec. Usable Analytical Expressions for Temperature Distribution Induced by Ultrafast Laser Pulses in Dielectric Solids. Micromachines, 2024, 15 (2), pp.0019. 10.3390/mi15020196. hal-04424108

HAL Id: hal-04424108

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

Submitted on 29 Jan 2024

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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Usable Analytical Expressions for Temperature Distribution Induced by Ultrafast Laser Pulses in Dielectric Solids

Ruyue Que, Matthieu Lancry 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.); matthieu.lancry@universite-paris-saclay.fr (M.L.)

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

**Abstract:** This paper focuses on the critical role of temperature in ultrafast direct laser writing processes, where temperature changes can trigger or exclusively drive certain transformations, such as phase transitions. It is important to consider both the temporal dynamics and spatial temperature distribution for the effective control of material modifications. We present analytical expressions for temperature variations induced by multi-pulse absorption, applicable to pulse durations significantly shorter than nanoseconds within a spherical energy source. The objective is to provide easy-to-use expressions to facilitate engineering tasks. Specifically, the expressions are shown to depend on just two parameters: the initial temperature at the center denoted as $T_{00}$ and a factor $R_{\tau}$ representing the ratio of the pulse period $\tau_{p}$ to the diffusion time $\tau_{d}$. We show that temperature, oscillating between $T_{max}$ and $T_{min}$, reaches a steady state and we calculate the least number of pulses required to reach the steady state. The paper defines the occurrence of heat accumulation precisely and elucidates that a temperature increase does not accompany systematically heat accumulation but depends on a set of laser parameters. It also highlights the temporal differences in temperature at the focus compared to areas outside the focus. Furthermore, the study suggests circumstances under which averaging the temperature over the pulse period can provide an even simpler approach. This work is instrumental in comprehending the diverse temperature effects observed in various experiments and in preparing for experimental setup. It also aids in determining whether temperature plays a role in the processes of direct laser writing. Toward the end of the paper, several application examples are provided.

**Keywords:** temperature distribution; femtosecond pulsed laser; interaction laser-dielectric solid

## 1. Introduction

In the context of an ultrafast laser interacting with solids, temperature plays a special role in the transformation processes. Some of the processes can be thermally activated, others can be temperature driven, such as phase transition but not thermally activated. The objective of this paper is to develop an analytic approximation to predict the behavior of the spatial temperature distribution and the temperature evolution over time according to the key laser parameter combinations and then to deduce their importance. This approach seeks to provide physical insight and semi-quantitative results without relying on heavy and overly detailed finite element calculations. This methodology resonates with the philosophy espoused by Paul Dirac in 1929, as documented in the Proceedings of the Royal Society of London [1]:

> "The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble. It therefore becomes desirable that approximate practical methods of applying

**Citation:** Que, R.; Lancry, M.; Poumellec, A.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

Academic Editor: Antonio Ancona

Received: 1 December 2023

Revised: 16 January 2024

Accepted: 23 January 2024

Published: 27 January 2024

<figure><img src="image_3.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/).
quantum mechanics should be developed, which can lead to an explanation of the main features of complex atomic systems without too much computation".

In the ultrafast laser-matter interaction process, the energy from the laser pulse that has an extremely short pulse duration (10<sup>-11</sup>–10<sup>-14</sup> s) is partially injected into a small focal volume of transparent dielectric solids. This intense laser pulse with high irradiance (>10<sup>13</sup> Wcm<sup>-2</sup>) in the focal region stimulates a series of complex dynamic processes, such as multiphoton ionization, tunneling ionization, inverse bremsstrahlung absorption, and avallanche ionization within an ultrashort time scale [2]. Such interactions lead to high-density electron excitations in condensed matter, creating plasma with high temperatures and pressures. This plasma expands rapidly in the focal zone, resulting in structural modifications as energy relaxes through phonon-electron interactions [3,4].

In the low repetition rate regime, thermal accumulation is usually negligible in the processing. The temperature decreases to the initial degree before the next pulse arrives. The non-linear nature of the optical absorption can confine the formed modifications to the focal volume. These advantages minimize the thermal collateral damage and heat-affected zone [5]. Thus, ultrafast laser direct writing (ULDW) is suggested as a general technique to induce highly localized modifications and optical structures within/near the focus in a wide range of transparent solids [6–11]. In this regime, denoted as non-thermal ULDW, the repetition rate (RR) is usually a few kilohertz, and the fabrication efficiency is also limited by the low pulse RR.

In contrast, when the pulse repetition rate of the ultrafast laser increases, the interval between successive laser pulses is less than the time needed for the absorbed energy to diffuse out of the focal volume and this induces an obvious localized heat accumulation effect [11–18]. In this case, for a given pulse energy, the temperature increases continuously in the focal zone before stabilizing. The final diffusion of the heat into the surrounding material may lead to a material melting beyond the focal volume over a longer time scale. In this regime, denoted as thermal ULDW, the melted modified region is much larger than the focus size. Paralleling to the wide applications of non-thermal ULDW, the localized thermal accumulation has been demonstrated to be important in the ULDW for inducing various phenomena and structures in the transparent solids and improving the performance of the fabricated devices. For example, the thermal accumulation can lead to a higher symmetry of the waveguide cross-section, reducing the propagation loss and enhancing the fabrication efficiency [12,15,16]. Thermal accumulation in the ULDW can also induce unique phenomena, such as elemental redistribution and local crystallization, which are nearly not achievable in the non-thermal ULDW [11,19–23]. In the thermal ULDW regime, the temperature distribution can work as a driving force to redistribute the elements or reorganize the structures in the thermal melted region. The thermal accumulation effect has also been reported to be critical for the formation of periodic nanogratings in some glass systems [24]. Moreover, the thermal accumulation induces a high temperature that can produce thermally excited free electrons, which seeds the avallanche ionization and significantly enhance absorptivity [25]. As a result, more energy can be absorbed, and this further increases thermal accumulation. Until now, thermal accumulation has been established to be an important assistor in many cases to help ULDW to achieve various applications in fundamental science and technological manufacturing [14–16,21,25,26]. Clarifying the principle of thermal ULDW and reviewing its current stage in the applications are highly urgent and significant for guiding future work [11,15,16].

For this aspect of the work, Lax et al. in 1977 [27] published the first paper that described the 3D spatial distribution of the temperature rise induced by the Beer-Lambert absorption of a static Gaussian CW laser beam in cylindrical geometry. Then, Sanders in 1984 [28] described an extension of these calculations for scanning beams and provided analytic expressions. In 1991, Haba et al. [29] described the calculation of a 3D spatial distribution for the Beer-Lambert absorption of a scanning Gaussian pulsed laser in cylindrical geometry. However, even if the expression was quite complete but numerically
solvable, there was no extended discussion on the laser/material parameters. Then, Eaton et al. [15] in 2005 and Zhang et al. [30] in 2007 performed finite difference calculations, for simple pulsed and CW Gaussian beams in spherical geometry, preventing easy material analysis. In 2007, Sakakura et al. [18] solved the Fourier equation in the frame of cylindrical geometry for energy delivered by a Gaussian pulsed fs laser (pulse duration 220 fs, RR 3 Hz, pulse energy < 1 µJ). With such a weak RR, the calculation can be restricted to one pulse as the experimental measurement (a lens effect) was smaller than 1 ms. However, it is not a special case and for material treatment, a large number of pulses are required. That is why Miyamoto et al. [31] in the same year, deduced analytical expressions for scanning uniform pulsed laser energy deposition in a parallelepiped volume of width 2 w corresponding to the scanning CW beam diameter at $1/e$ and length corresponding to the absorption length $(1/\alpha)$. These calculations were used also by Beresna et al. [32] and applied to a particular case, i.e., borosilicate. In 2011, Miyamoto et al. [25] considered a cylindrical source with its full width dependent on z in order to account for the convergence of the beam or the non-linear properties including the self-focusing. In 2012, Shimizu et al. [33] used a static cylindrical Gaussian beam and Gaussian energy deposition in depth for multi-pulsed laser energy deposition but solved the problem numerically. Lastly, in 2019 and 2020, Rahaman et al. [34,35] proposed an analytical solution usingDHamel's theorem and Hankel's transform method, for a transient, two-dimensional thermal model. We summarized the above research in Table 1 below, to compare with our work.

<table><caption>**Table 1.** State of the art of the thermal simulation of laser-matter interactions.</caption><thead><tr><th>Laser Type</th><th>Mode</th><th>Geometry</th><th>Source Shape</th><th>Solving Method</th><th>Refs.</th></tr></thead><tbody><tr><td>CW</td><td>static</td><td>cylindrical</td><td>Gaussian(r)<br>Beer-Lambert(z)</td><td>analytical</td><td>Lax [27]</td></tr><tr><td>pulsed</td><td>scanning</td><td>three axes</td><td>Gaussian(x,y)<br>Beer-Lambert(z)</td><td>analytical</td><td>Sanders [28]</td></tr><tr><td>pulsed</td><td>scanning</td><td>three axes</td><td>Gaussian(x,y)<br>Beer-Lambert(z)</td><td>analytical</td><td>Haba [29]</td></tr><tr><td>pulsed</td><td>static</td><td>spherical</td><td>Gaussian(r)</td><td>finite difference</td><td>Eaton [15] or Zhang [30]</td></tr><tr><td>pulsed</td><td>quasi-static</td><td>cylindrical</td><td>Gaussian(r)<br>Beer-Lambert(z)</td><td>analytical one pulse</td><td>Sakakura [18]</td></tr><tr><td>CW</td><td>scanning</td><td>three axes</td><td>uniform deposition in<br>parallelepiped volume</td><td>analytical</td><td>Miyamoto [31]</td></tr><tr><td>pulsed</td><td>static</td><td>cylindrical</td><td>Gaussian(r,z)</td><td>analytical</td><td>Miyamoto [25]</td></tr><tr><td>pulsed</td><td>static</td><td>cylindrical</td><td>Gaussian(r)<br>Gaussian(z)</td><td>numerical</td><td>Shimizu [33]</td></tr><tr><td>pulsed</td><td>scanning</td><td>cylindrical</td><td>Gaussian(r)<br>surface absorption</td><td>analytical</td><td>Rahaman [34,35]</td></tr><tr><td>pulsed</td><td>quasi-static</td><td>spherical</td><td>Gaussian(r)</td><td>analytical</td><td>this work</td></tr></tbody></table>

In short, the drawback in the available literature is that the authors did not provide simple expressions that allow the reader to easily understand how each parameter of lasers and materials influences the evolution of the temperature distribution, and to control the thermal effect in transparent materials with non-linear optical absorption for which the effect is mainly in volume for a focused beam. However, beyond this step that corresponds to the absorption of a small part of the pulse energy, the absorption becomes linear [36]. This is the reason why we present the analytical approach or link the properties of the materials to the shape of the temperature distribution and use it for explaining the phenomena such as:
- The appearance of several regions in the heat-affected volume including change of the structure of a glass, crystallization, phase separation, thermal erasure while writing providing that energy endo or exo is negligible in front of the laser one;

- The variations in the shape of the interaction volume according to the laser parameters like a change of laser track width, change of laser track morphology.

For this purpose, we restricted ourselves assuming that the physical properties of the material are independent of the temperature, but this does not prevent the possibility of physical deductions. We used the simplest solution of the Fourier equation in spherical geometry, i.e., a Gaussian shape along the perpendicular and longitudinal direction of the beam propagation direction. This applies not only to the sample surface but also to multiphoton absorption by stating the coordinate origin at the geometrical or effective focus. Since the typical application of this model is the thermal accumulation of a high focused beam in a material with non-linear absorption. Namely, the cylindrical symmetry and the Beer-Lambert law along z cannot be considered. We have also considered that the pulse duration (smaller than a few ns) is much smaller than the diffusion time so that the initial temperature distribution is defined by the shape of the absorbed energy source. This is applicable in most cases to femtosecond and nanosecond lasers as the diffusion time is usually of the order of a fraction of µs in inorganic glasses and a few µs in organic materials. In addition, material phase change and non-linear optical effects are not considered in this model except for the presence of coefficient A (see below).

This study was motivated by seeing nowadays that, as the means of simulation are easily accessible, the physical sense is hidden or even lost, which prevents the correct management of the laser parameters according to targeted property modifications.

## 2. Starting Formulation

From a theoretical point of view, the heat deposited at a point by the laser diffuses into the material by following Fourier's law $\vec{q} = -\kappa\vec{\nabla}T$ where $\vec{q}$ is the heat flow (energy per unit area and time). Fourier considers it to be linearly dependent on the temperature gradient. $\kappa$ is the thermal conductivity, in general, a tensor of order 2 which relates the gradient vector of $T$ to the flux. Its dimension is energy ($J/m^2 \cdot s \cdot K$). For isotropic materials, such as glasses, one will suppose that this tensor is reduced to a scalar. To calculate (in principle) the evolution and distribution of $T$, we start with the law of the conservation of energy, $\frac{d\rho_q}{dt} + \vec{\nabla} \cdot \vec{q} = $ source terms – sink terms. The source term is the laser energy density deposited per unit of time (i.e., absorbed laser power), written symbolically as $\frac{\delta\rho_q}{\delta t}$. Its spatial shape defines the symmetry of the problem. For the sake of simplicity for demonstrating physical conclusions, we have assumed spherical symmetry. This means that we do not take into account some changes of focal volume with incident pulse energy due to Kerr self-focusing and electron plasma defocusing. We assume that there were no heat annihilation terms (for example, endothermic chemical transformation), sink terms = 0. Using the definition of specific heat, $\frac{d\rho_q}{dt} = \rho \cdot C_p \cdot \frac{dT}{dt}$, $\rho$ and $C_p$ are the density and specific heat capacity, respectively. $\vec{\nabla} \cdot \vec{q} = \vec{\nabla} \cdot (-\kappa\vec{\nabla}T) = -D_{th} \cdot \Delta T$, with diffusivity $D_{th} = \frac{\kappa}{\rho \cdot C_p}$, $\Delta$ is the Laplace operator written in spherical symmetry, $\Delta = \nabla^2 = \frac{\partial^2}{\partial r^2} + \frac{2\partial}{r\partial r}$. Considering a beam moving not too fast, convection can be neglected (i.e., the time derivative of the spatial coordinate), $\frac{dT}{dt}$ is thus written as $\frac{\partial T}{\partial t}$. Therefore, we obtain the following equation:

$$ \frac{\partial T(r, t)}{\partial t} - D_{th} \cdot \Delta T(r, t) = \frac{1}{\rho C_p} \frac{\delta \rho_q}{\delta t} \quad (1) $$

Since the pulse duration is much less than the diffusion time ($w^2/D_{th}$ with $w$ is the beam waist radius at $1/e$), the latter is at the scale of $10^{-7}$ s and $10^{-6}$ s, the diffusion process
can be considered therefore to be well separated from the deposition process. During the pulse, a deposition of energy density takes place, but the diffusion does not begin, so $D_{th} \cdot \Delta T = 0$, and Equation (1) becomes:

$$ \frac{\partial T(r, t)}{\partial t} = \frac{1}{\rho C_p} \frac{\delta \rho_Q}{\delta t} \qquad (2) $$

Assuming a Gaussian shape of $\frac{\delta \rho_Q}{\delta t}(r, t) = \frac{A \cdot E_p}{\pi^2 w^3} \cdot \exp[-\frac{r^2}{w^2}] \cdot f(t)$, where $w$ is the beam waist radius (at $1/e$), $f(t)$ is the pulse shape (integral of $f(t)$ on the pulse time = 1), $E_p$ is the energy of the pulse, $A$ represents the absorbed fraction of the pulse energy.

Therefore, $T(r, 0) - T_{room} = \frac{1}{\rho C_p} \int_{pulse} \frac{\delta \rho_Q}{\delta t} dt = T_{00} \cdot \exp[-\frac{r^2}{w^2}]$ with:

$$ T_{00} = \frac{A \cdot E_p}{\frac{3}{\pi^2 \rho C_p w^3}} \qquad (3) $$

After pulse energy deposition, diffusion begins to operate, and Equation (1) becomes:

$$ \frac{\partial T(r, t)}{\partial t} - D_{th} \cdot \Delta T(r, t) = 0 \qquad (4) $$

Using the initial and boundary conditions on solutions of Equations (2) and (4), we obtain:

$$ T(r,t) = T_{00} \cdot \frac{w^3}{(w^2 + 4D_{th} \cdot t)^{\frac{3}{2}}} \cdot \exp \left[ - \left( \frac{r^2}{w^2 + 4D_{th} \cdot t} \right) \right] + T_{room} \qquad (5) $$

Equation (5) describes a single-pulse-induced temperature distribution over time. $T_{00}$ is the maximum temperature induced by a laser pulse at the focus center. $T_{room}$ is the ambient temperature, which will be omitted for ease of calculation. The temperature should thus be understood as the temperature increment above the initial sample temperature.

It is important to note that when utilizing the spherical model, the deposited energy volume will consistently yield a higher temperature than reality, as the size along $z$ is usually larger than the waist radius. Given our primary concern lies in assessing the temperature's dependence on various parameters, it is possible to adjust the actual calculated value, which is notably affected by the absorption fraction $A$, to align it more closely with reality.

In the case of the absorption of $N$ pulses, we easily obtain the evolution of the distribution considering the linearity of the differential equation and making up the sum of the solution for one impulsion but shifted in time $\tau_p = 1/RR$, where $RR$ is the pulse repetition rate:

$$ T(r,t) = T_{00} \cdot \sum_{N-1 \text{ integer part}(t \cdot f)}^{N-1} \frac{w^3}{(w^2 + 4D_{th} \cdot t)^{3/2}} \cdot \exp \left[ - \left( \frac{r^2}{w^2 + 4D_{th}(t - n \cdot \tau_p)} \right) \right] \qquad (6) $$

With $\tau_d = w^2/4D_{th}, r_w = \frac{r}{w}$, Equation (6) reads:
$$T(r_w, t) = T_{00} \cdot \sum_{n=0}^{N-1} \frac{1}{\left(1 + \frac{t-n \cdot \tau_p}{\tau_d}\right)^{3/2}} \cdot \exp\left[-\left(\frac{r_w^2}{1 + \frac{t-n \cdot \tau_p}{\tau_d}}\right)\right] \quad (7)$$

We note that the variables involved in Equation (7) are the ratio between the period of the pulses $\tau_p$ and the diffusion times $\tau_d$, while the other laser and material parameters are involved in the amplitude $T_{00}$. Therefore, we introduce the parameter $R_\tau$:

$$R_\tau = \frac{\tau_p}{\tau_d} \qquad (8)$$

Therefore, Equation (7) becomes:

$$\frac{T(r_w, t)}{T_{00}} = \sum_{n=0}^{N-1} \frac{1}{N-1=\text{integer part}\left(\frac{t}{\tau_p}\right) \left[1 + \left(\frac{t}{\tau_p} - n\right) \cdot R_\tau\right]^2} \cdot \exp\left[-\frac{(r_w)^2}{1 + \left(\frac{t}{\tau_p} - n\right) \cdot R_\tau}\right] \quad (9)$$

where N is the number of pulses defined from the time t.

The objective now is to compute the value of the temperature T according to the coordinate $r_w$ when $N \gg 1/R_\tau$. We will show how the temperature changes with the number of pulses according to heat accumulation (hence $R_\tau$), i.e., when T at the end of the period cumulates with the increase induced by the absorption of the next pulse. We will also describe the properties of the temperature on average in the pulse repetition period. We will also show that a steady state can be reached and give the practical number of pulses for that.

### 3. Final Temperatures at Steady State

At first, we separate the temperature problem into two cases: (1) at the center, i.e., $r_w = 0$; (2) for general cases when $r_w \neq 0$ including case (1). Again, for the sake of simplicity, $T_{00}$ will be usually omitted. Therefore, the subsequent temperatures will virtually include $T_{00}$.

#### 3.1. At the Center

At the center, $r_w = 0$, and thus, from Equation (9):

$$T(0, t) = \sum_{n=0}^{N-1} \frac{1}{N-1=\text{integer part}\left(\frac{t}{\tau_p}\right) \left[1 + \left(\frac{t}{\tau_p} - n\right) \cdot R_\tau\right]^2} \quad (10)$$

Calling $N_t = \frac{t}{\tau_p}$. The temperature evolutions over the generalized pulse number $N_t$ for several $R_\tau$ are shown in Figure 1.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Plot of the relative temperature (Equation (10)) at the center $r_w = 0$ according to the generalized pulse number $N_t = \frac{t}{\tau_p}$ with $R_\tau = 0.2, 2, 20$ until several pulse numbers (a) 10, (b) 20, and (c) 100.</figcaption></figure>

From Figure 1 we observe that:

- $T(0, N_t)$ oscillates between a minimum ($T_{min}$) and a maximum ($T_{max}$) in each period between two pulses;

- The oscillation amplitude ($T_{osc}$) seems to be the same, whatever $R_\tau$;

- $T$ seems to reach a steady state as $N_t$ becomes large (already seen in various papers [29,31,32]);

- The number of pulses to reach this 'steady state' appears very small for a large $R_\tau$ but larger for small $R_\tau$ values. For a larger $R_\tau$, the temporal overlapping of temperature increase contributions from consecutive pulses is weaker, whereas it increases (heat accumulation) when $R_\tau$ is smaller.

### 3.1.1. The Oscillation Amplitude $T_{osc}$

We observe the oscillations of temperature on time in Figure 1 on each period. Just after the pulse energy deposition, the temperature experiences a sudden increase and then a slow decrease until the next pulse arrival. It is important to know the amplitude of the temperature oscillations ($T_{osc}$) because when $T_{osc}$ is large, at the beginning of a period, temperature may be high enough for transformation but in a short time, and at the end of the period during a long time duration, the temperature can be low, maybe achieving another transformation. The middle part could therefore often be the most active part.
The Limit of the Temperature Oscillation Amplitude after an Infinite Number of Pulses

The question here is: how do the oscillations evolve in time according to the pulse number N for a given diffusion time? If the period is large ($R_τ$ large), we expect independent pulses and thus the amplitude will be $T_{00}$. However, when the pulse period is small ($R_τ$ small), can we imagine a smaller oscillation? The next calculation provides answers.

For that purpose, we compare the difference between the maximum T and minimum T of the Nth pulse, $T_{max}(0,N) - T_{min}(0,N) = T(0, t_N) - T(0, t_{N+1} - \epsilon)$, where $\epsilon$ is an arbitrary small quantity for ensuring that the number of pulses in the expression (11) is the same. $T_{max}$ is defined just after the deposition of the Nth pulse, so at the beginning of the pulse, $t_N = (N - 1)\tau_p$. $T_{min}$ is at the end of the pulse period, just before the (N + 1)th pulse arrival. Using Equation (10), we have:

$$ 
\begin{aligned}
T_{max}(0,N) &= T(0, t = t_N = (N-1)\tau_p) = \sum_{n=0}^{N-1} \frac{1}{[1 + (N-1-n) \cdot R_τ]^3} \\
&= \sum_{n=0}^{N-1} \frac{1}{[1 + n' \cdot R_τ]^3}
\end{aligned} 
\quad (11) $$

$T_{min}$ will be at $t = N \cdot \tau_p - \epsilon$, thus not containing the temperature contribution induced by the (N + 1)th pulse, so:

$$ 
\begin{aligned}
T_{min}(0,N) &= T(0, t = t_{N+1} - \epsilon = N \cdot \tau_p - \epsilon) = \sum_{n=0}^{N-1} \frac{1}{[1 + (N-n) \cdot R_τ]^3} \\
&= \sum_{n'=1}^{N} \frac{1}{[1 + n' \cdot R_τ]^3}
\end{aligned} 
\quad (12) $$

Therefore,

$$ T_{osc}(0,N) = T_{max}(0,N) - T_{min}(0,N) = 1 - \frac{1}{[1 + N \cdot R_{\tau}]^3} \quad (13) $$

and

$$ T_{osc}(0, \infty) = \lim_{N \to \infty} T_{osc}(0, N) = 1 \quad (14) $$

When $N \gg 1/R_τ$, $T_{osc}$ tends to 1. This means $T_{00}$ in the absolute scale. $T_{osc}$ according to the pulse number is shown in Figure 2. It reaches a maximum value 1, i.e., $T_{00}$, when $N \gg 1/R_τ$. At the beginning of the irradiation, $T_{osc}$ starts with a value smaller than $T_{00}$, where a smaller $R_τ$ leads to a smaller oscillation at the beginning. When the pulse number N increases until some value, $T_{osc}$ reaches $T_{00}$. When $R_τ$ is large, e.g., 10, the amplitude is equal to $T_{00}$ whatever N, as pulse contributions are separated (no overlapping). With the expression of $T_{00}$ (Equation (3)), which is proportional to pulse energy ($E_p$), the temperature oscillation range can be determined.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Temperature oscillation $T_{osc}(0, N)/T_{00}$ (here quoted $T_{osc}(N, R_{\tau})$ dependence according to pulse number and for three values of $R_{\tau} = 0.1, 1$, and 10).</figcaption></figure>

The Effective Number of Pulses for Reaching the Limit of $T_{osc}(N_{ss0}^0)$

When will the temperature in the material reach a stable oscillation? In practice, we can calculate a real number of pulses ($N_{ss0}^0$) to closely reach the oscillation limit. (In notation, $N_{ss0}^0$, where sso means steady state of oscillation, 0 means the situation when $r_w = 0$).

Consider when $T_{osc}(0, N_{ss0}^0) = (1 - \varepsilon) \cdot T_{osc}(0, \infty)$, the limit is practically reached, where $\varepsilon$ is a small quantity, i.e., a few % (based on the actual situation). Thus, it is:

$$ \frac{|T_{ocs}(0, \infty) - T_{osc}(0, N_{ss0}^0)|}{T_{ocs}(0, \infty)} < \varepsilon \qquad (15) $$

Then, we have:

$$ N_{ss0}^0 = \frac{1}{R_{\tau}} \left[ \left( \frac{1}{\varepsilon} \right)^{\frac{2}{3}} - 1 \right] \qquad (16) $$

The plot of $N_{ss0}^0$ is shown in Figure 3, note that the actual number of pulses is the integer part above 1 (as $\varepsilon$ should be smaller, bounded by $\varepsilon = \frac{1}{[1+R_{\tau}]^{\frac{2}{3}}}$. Some specific parameters are given below for visualizing this value. With $\varepsilon = 3\%$, when $R_{\tau} = 1$ (conceivable combinations of material parameters and laser RR), $N_{ss0}^0 = 9.36$. So, after 10 pulses, the amplitude of the oscillating temperature reaches 0.97 $T_{00}$. When $R_{\tau}$ is large, e.g., $R_{\tau} = 10$, $N_{ss0}^0 = 0.94$, so only one pulse rules the oscillation amplitude, and we can understand that pulse contributions are separated. When $R_{\tau}$ is smaller, $N_{ss0}^0$ increases rapidly, e.g., $R_{\tau} = 0.1$, i.e., $1/R_{\tau} = 10$, $N_{ss0}^0 = 94$. Beyond $N_{ss0}^0$ pulses, the oscillation amplitude becomes almost constant. According to the pulse period, we can know the time to reach the constant oscillation amplitude. By comparing with pulse number $N = 1$ (blue dashed line in Figure 3), we can deduce in what condition ($R_{\tau}$ larger than which value) the temperature oscillation is constant since the first pulse.
<figure><img src="image_6.png" /><figcaption>Figure 3. Number of pulses to reach the oscillation amplitude limit $N_{ss0}^0$ according to $R_τ$ from 0.1 to 100, with $ε = 3\%$. For $1/R_τ > 0.1$, weak heat conduction, large RR, and vice versa.</figcaption></figure>

To conclude on this point, Equations (13) and (14) provides that the temperature oscillation amplitude is $T_{00}$ after $N_{ss0}$ pulses and at an oscillatory steady state. Equation (16) provides the effective number of pulses for reaching it. It takes more time when $R_τ$ is very small (slow heat diffusion or high pulse RR) but the requested time remains quite small. In particular, the temperature oscillation amplitude is only relevant to certain laser parameters of a single pulse and material parameters of the energy-to-temperature conversion relationship, independent of RR and diffusion parameters.

$T_{min}$ and $T_{max}$

We demonstrate in Appendix B that the temperature induced by laser pulses will not increase indefinitely but converge to a finite value. This defines a steady state that corresponds to the equilibrium between the energy supplied by the laser and the energy diffusing out of the irradiated voxel.

The Limit of $T_{max}$ and $T_{min}$

$T_{min}$: The analytical expression of the minimum temperature is transformed from the sum expression Equation (12), with details found in Appendix C. Therefore, we obtain:

$$T_{min}(0, N) \approx \frac{1}{2(1 + R_{\tau})^{\frac{3}{2}}} + \frac{1}{2(1 + N \cdot R_{\tau})^{\frac{3}{2}}} + \frac{2}{R_{\tau}} \left[ \frac{1}{\sqrt{1 + R_{\tau}}} - \frac{1}{\sqrt{1 + N \cdot R_{\tau}}} \right] \quad (17)$$

This expression shows the increase of $T_{min}$ according to $N$ and $R_τ$. It is plotted in Figure 4.
<figure><img src="image_7.png" /><figcaption>Figure 4. $T_{min}/T_{00}$ at $r_w = 0$ according to N increasing from 1 to 10,000 when $R_τ = 0.2$, 1, and 5.</figcaption></figure>

The final limit of $T_{min}$, i.e., when $N \gg 1/R_τ$ is given below:

$$T_{min}(0, \infty) = \frac{1}{2(1 + R_{\tau})^{\frac{3}{2}}} + \frac{2}{R_{\tau}} \frac{1}{\sqrt{1 + R_{\tau}}} \approx \frac{2}{R_{\tau}} \frac{1}{\sqrt{1 + R_{\tau}}} \quad (18)$$

The same method has been applied to obtain the $T_{max}$ limit, and the detail can be also found in Appendix C. We have thus:

$$T_{max}(0, N) \approx 1 + \frac{1}{2[1 + R_{\tau}]^{3/2}} + \frac{1}{2[1 + (N - 1)R_{\tau}]^{3/2}} + \frac{2}{R_{\tau}} \left[ \frac{1}{\sqrt{[1 + R_{\tau}]}} - \frac{1}{\sqrt{1 + (N - 1)R_{\tau}}} \right] \quad (19)$$

When $N \gg 1/R_τ$, the $T_{max}$ limit is:

$$T_{max}(0, \infty) = 1 + \frac{1}{2(1 + R_{\tau})^{\frac{3}{2}}} + \frac{2}{R_{\tau}} \frac{1}{\sqrt{1 + R_{\tau}}} \approx 1 + \frac{2}{R_{\tau}} \frac{1}{\sqrt{1 + R_{\tau}}} \quad (20)$$

From these expressions, we see that the difference between $T_{max}(0, \infty)$ and $T_{min}(0, \infty)$ is 1, which is consistent with the oscillation amplitude limitation (Equation (14)).

When $R_τ$ reaches 0 (e.g., by increasing pulse RR or with the material of small thermal conductivity), Equations (18) and (20) are approximately proportional to $\frac{2}{R_τ}$. Reintroducing here exceptionally $T_{00}$ (Equation (3)), we obtain:

$$T_{max}(0, \infty) \sim T_{min}(0, \infty) \sim T_{00} \cdot \frac{2}{R_{\tau}} = \frac{2AE_p}{\pi^{\frac{3}{2}} \rho C_p w^3 R_{\tau}} = \frac{2AE_p \tau_d}{\pi^{\frac{3}{2}} \rho C_p w^3 \tau_p} = \frac{2AE_p f}{\pi^{\frac{3}{2}} D_{th} \rho C_p w} = \frac{2AP}{\pi^{\frac{3}{2}} \kappa w} \quad (21)$$

with P being the average power.

We note that the temperature is now dependent on the incident laser power as is the case for CW lasers, and inversely dependent on the thermal conductivity ($\kappa$), whereas $T_{00}$ was dependent on the incident pulse energy ($E_p$), not on the thermal diffusivity but just on the heat capacity of the material. This is due to large time-overlapping of the pulse contribution when $R_τ$ reaches 0.

Therefore, increasing the pulse RR with constant $E_p$ leads to a faster temperature increase but NOT with constant average power. The same maximal temperature can be achieved with or without heat accumulation. However, $T_{min}$, which is negligible in front of $T_{max}$ for large $R_τ$ values, increases until it almost equals $T_{max}$ for small $R_τ$ values (large RR).
The Effective Number of Pulses for Reaching the Limit of $T_{min}$ and $T_{max}$ ($N_{ssmin}^0, N_{ssmax}^0$)

$N_{ssmin}$: The effective number of pulses to reach the steady state $N_{ss}$ is defined to have temperature reaching $T_{min}$ or $T_{max}$. As the same definition as for $N_{ss0}$, the first approximation of $N_{ssmin}$ is obtained by solving the following assertion, $\frac{|T_{min}(0,N)-T_{min}(0,\infty)|}{T_{min}(0,\infty)} < \varepsilon$, with $\varepsilon$ being a small quantity. Posing $X = \frac{1}{\sqrt{1+N\cdot R_{\tau}}}$, it reads $\frac{X^3}{2} - \frac{2}{R_{\tau}}X + \varepsilon \cdot T_{min}(0, \infty) > 0$. This cubic equation has three roots, where the physical one is $X < \frac{\varepsilon \cdot T_{min}(0, \infty)}{\frac{2}{R_{\tau}}}$. Therefore,

$$N_{ssmin}^0 > \frac{1}{R_{\tau}} \left[ \left( \frac{2}{R_{\tau} \cdot \varepsilon \cdot T_{min}(0, \infty)} \right)^2 - 1 \right] \quad (22)$$

$N_{ssmax}$: With the same method, we obtain:

$$N_{ssmax}^0 > \frac{1}{R_{\tau}} \left[ \left( \frac{2}{R_{\tau} \cdot \varepsilon \cdot T_{max}(0, \infty)} \right)^2 - 1 \right] \quad (23)$$

The $N_{ss}$ for reaching closely the steady state (with $\varepsilon$ departure). $T_{osc}(0, \infty)$, $T_{min}(0, \infty)$, and $T_{max}(0, \infty)$ are plotted in Figure 5 according to $R_{\tau}$ (with $\varepsilon = 3\%$).

<figure><img src="image_8.png" /><figcaption>Figure 5. The effective number to reach the limit of $T_{osc}$, $T_{min}$, and $T_{max}$ according to $R_{\tau}$ from 0.1 to 100 for $\varepsilon = 0.03$. The boundary point is at $R_{\tau} = 15.7$.</figcaption></figure>

From Figure 5, we can see that with $1/R_{\tau}$ increasing, the effective pulse numbers for reaching the steady state increases whether for $T_{osc}$, $T_{min}$, or $T_{max}$. For practical use, it is better to define one $N_{ss}$ for calculation. When $R_{\tau}$ is large, since the value of $T_{min}$ is almost 0 (no pulse superimposition), it is therefore not meaningful to take it into consideration. For $N_{ss0}$ and $N_{ssmax}$, the value converges to 0 when $R_{\tau}$ is large because it can be considered to be already at steady state when the pulses are separated. As observed, the green dash is always higher than the red line when $N > 1$, and the oscillation reaches the steady state faster than $T_{max}$. Therefore, $N_{ssmax}^0$ is the suitable and practical number of pulses needed for reaching the steady state. Some examples for the $N_{ssmax}$ value are shown in Table 2 (with $\varepsilon = 0.03$ for organic materials and $\varepsilon = 0.06$ for inorganic materials):
<table><caption>**Table 2.** Pulse number needed ($N_{ss}$) for reaching the steady state in materials, using Equation (23).</caption><thead><tr><th></th><th>SiO<sub>2</sub></th><th>LNS</th><th>STS</th><th>Glycine</th><th>Zeonex</th><th>Sucrose</th><th>Nifedipine</th></tr></thead><tbody><tr><td>τ<sub>d</sub> (μs)</td><td>0.28</td><td>0.235</td><td>0.04</td><td>0.28</td><td>0.42</td><td>4.9</td><td>1.63</td></tr><tr><td>RR(kHz)</td><td>200</td><td>200</td><td>200</td><td>200</td><td>200</td><td>200</td><td>200</td></tr><tr><td>R<sub>τ</sub></td><td>18</td><td>21</td><td>125</td><td>18</td><td>12</td><td>1</td><td>3</td></tr><tr><td>N<sub>ssmax</sub></td><td>1</td><td>1</td><td>1</td><td>1</td><td>3</td><td>641</td><td>80</td></tr></tbody></table>

We observe that for the inorganic material examples in the table, with $RR = 200$ kHz, there is no heat accumulation, pulse contributions are separated, there is no transient time, and the time variation of $T$ is from one pulse contribution. However, for the organic material examples, except for glycine crystal, with same $RR$, the laser induces heat accumulation. Therefore, not only can we deduce from the known laser and material parameters whether or not there will be heat accumulation, but we can also easily backtrack on how to choose a laser $RR$ that avoids or guarantees heat accumulation in a particular material.

We can thus define the boundary between the two domains by $N_{ssmax}(R_{τ,ε}) = 1$ as shown in Figure 5 by the blue point. When the effective pulse number is equal to 1, the pulse contribution is separated, so it is considered that there is no heat accumulation. This is the heat accumulation definition we propose with a new perspective. $R_{τ}$ varies with the level of sensitivity of the targeted transformation, for $ε = 0.03$, $R_{τ} = 15.7$, for $ε = 0.06$, $R_{τ} = 7$, for instance.

From $N_{ssmax}$ together with the laser pulse $RR$, we know the time needed to reach the steady state. Accordingly, the time for reaching the steady state $t_{ss}^0$ is (considering the effective number to reach the $T_{max}$ limit):

$$t_{ss}^0 = N_{ssmax} \tau_p = \tau_d \left[ \left( \frac{2}{R_{\tau} \cdot \varepsilon \cdot T_{max}(0, \infty)} \right)^2 - 1 \right] \quad (24)$$

Figure 6 show the plots versus $\frac{1}{R_{\tau}} = \frac{RR \cdot w^2}{4 \cdot D_{th}}$ according to three different diffusion times: 0.28 µs (silica, glycine), 1.63 µs for nifedipine, and 4.9 µs for sucrose.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** The time in µs to reach the steady state according to $1/R_{τ}$ value for glycine or silica (red), nifedipine (blue dash), and sucrose (green). The value of the second parameter in $tss$ corresponds to $τ_d$ in Table 2.</figcaption></figure>

For small enough values of $R_{τ}$, the time reaches the value $τ_d/ε^2$, i.e., $1111τ_d$ for $ε = 3\%$. Note that for inorganic glass and glycine crystal cited above with $τ_d = 0.28$ µs, this time is smaller than 1 ms (red profile). However, for sucrose and nifedipine, this time is
5.5 ms and 1.8 ms, respectively. In any case, the important fact is the independency of the transient time with $R_τ$ for small values (see for $R_τ < 1$) and thus it is bounded to quite a small value.

For large enough values of $R_τ$, this time is limited by the period $τ_p$ that increases with $R_τ$.

### 3.2. Time Behavior out of the Center ($r_w = r/w \neq 0$)

When ($r_w \neq 0$, we come back to the expression Equation (10):

$$T(r_w, t) = \sum_{\substack{n=0 \\ N-1=\text{integer part} (t/\tau_p)}}^{N-1} \frac{1}{\left[1 + \left(\frac{t}{\tau_p} - n\right) \cdot R_\tau\right]^{3/2}} \cdot \exp\left[-\frac{(r_w)^2}{1 + \left(\frac{t}{\tau_p} - n\right) \cdot R_\tau}\right]$$

Figure 7 shows the temperature evolution based on the above expression over time at two relative distances $r_w = 1$ (Figure 7a,b) and $r_w = 2$ (Figure 7c,d).

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Plot of the relative temperature ($T/T_{00}$) with $R_τ = 0.2, 2, 20$ at (a,b) $r = w$ from (a) pulse 1 to pulse 20 and (b) pulse 100 to pulse 120 (c,d) $r = 2w$ from (c) pulse 1 to pulse 20 and (d) pulse 100 to pulse 120. Inserts (b,d) zoom of pulse 106–108 of $R_τ = 2$ and 20 at $r = w$ and $r = 2w$, respectively.</figcaption></figure>

We observe the following differences according to radius $r_w = 0, 1, 2$:

- The amplitude of oscillation is less than 1 (in the unit of $T_{00}$) for increasing radius;

- The maximum temperature during a period is still at the beginning of the pulse deposition for $r_w = 1$ with these three $R_τ$, while at $r_w = 2$ the maximum temperature is no more at the beginning. That is because there is time for heat to diffuse from the
center to $r_w$. This renders the following calculation of $T_{max}$ for increasing radius to be more complex.

### 3.2.1. $T_{osc}$, $T_{min}$, and $T_{max}$

The Limit of $T_{max}$ and $T_{min}$ (when $N \gg 1/R_τ$)

To calculate the oscillation amplitude $T_{osc}$, it is the same as the case of $r_w = 0$. In general, we compare the difference between the maximum $T$ and minimum $T$ in the $N$th pulse period. $T_{min}$ is still considered at the end of the $N$th one, i.e., when $t = N \cdot \tau_p$ before the absorption of the $(N+1)$th pulse, so:

$$T_{min}(r_w, N) = \sum_{n=0}^{N-1} \frac{1}{[1 + (N - n) \cdot R_{\tau}]^3} \exp \left[ - \frac{(r_w)^2}{1 + (N - n) \cdot R_{\tau}} \right] \quad (25)$$

However, since $T_{max}$ in some situations can be in the middle of the pulse period, we set $x_m$, $0 \le x_m \le 1$ to define the place where the $T_{max}$ is. Therefore, $T_{max}$ is at the time $t = (N - 1 + x_m)\tau_p$:

$$T_{max}(r_w, N) = \sum_{n=0}^{N-1} \frac{1}{[1 + (N - 1 + x_m - n) \cdot R_{\tau}]^3} \exp \left[ - \frac{(r_w)^2}{1 + (N - 1 + x_m - n) \cdot R_{\tau}} \right] \quad (26)$$

This expression is also a general expression to describe both $T_{max}$ and $T_{min}$, while $T_{min}$ appears at the end of the period, i.e., $x_m = 1$, as well as the case when $r_w = 0$, $T_{max}$ appears at the beginning of the period with $x_m = 0$.

The position of the maximum $x_m$ is solved as a function of $R_τ$, $r_w$, and $N$. When considering the steady state, when $N \gg 1/R_τ$, $x_m$ is shown below, and the results of the cumbersome calculation details can be found in Appendix D:

$$x_m = \frac{\sqrt{R_\tau} \sqrt{9R_\tau + 32(r_w)^2 - 3R_\tau - 8}}{8R_\tau} \quad (27)$$

Based on Equation (27), when considering different values of $r_w$ and $R_τ$, the thermal calculation can be divided into two situations (see Appendix D):

Situation 1: $x_m = 0$, when, $r_w < \sqrt{\frac{3}{2} + \frac{2}{R_τ}}$ whatever $R_τ$ or $R_τ$ small enough (less than $\frac{2}{r_w^2-1.5}$ when $r_w^2 > 1.5$).

Situation 2: $x_m \neq 0$, i.e., $r_w > \sqrt{\frac{3}{2} + \frac{2}{R_τ}}$, in this situation the maximum temperature is in the middle of the period, and the expressions of $T_{max}$ and $T_{osc}$ should contain $x_m$.

(1) For situation 1, when $x_m = 0$, the limit of $T_{osc}$ is described as:

$$T_{osc}(r_w, N) = T_{max}(r_w, N) - T_{min}(r_w, N) = \exp[-(r_w)^2] - \frac{\exp[-\frac{(r_w)^2}{1 + N \cdot R_\tau}]}{(1 + N \cdot R_\tau)^{\frac{3}{2}}} \xrightarrow{N \gg 1/R_\tau} \exp[-(r_w)^2] \quad (28)$$

The amplitude of the temperature oscillations reaches $T_{00} \cdot \exp[-(r_w)^2]$ whatever $R_τ$. It is consistent with our observations in Figure 7, e.g., the amplitudes are 0.368 $T_{00}$ and 0.018 $T_{00}$ at $r_w = 1$ and $r_w = 2$, respectively.

Therefore, for $T_{min}$ and $T_{max}$, using the trapezoidal rule for approximation as for $r_w = 0$, we have $T_{min}$ from Equation (25):
$$T_{min}(r_w, N) \approx \frac{1}{2} \left[ \frac{\exp\left[-\frac{(r_w)^2}{1+R_\tau}\right]}{(1+R_\tau)^{3/2}} + \frac{\exp\left[-\frac{(r_w)^2}{1+N \cdot R_\tau}\right]}{(1+N \cdot R_\tau)^{3/2}} \right] + \frac{\sqrt{\pi}}{R_\tau \cdot r_w} \left[ \operatorname{erf}\left(\frac{r_w}{\sqrt{1+R_\tau}}\right) - \operatorname{erf}\left(\frac{r_w}{\sqrt{1+N \cdot R_\tau}}\right) \right]$$

$$T_{min}(r_w, N) \underset{N \gg 1/R_{\tau}}{\longrightarrow} T_{min}(r_w, \infty) = \frac{\exp\left[-\frac{(r_w)^2}{1+R_{\tau}}\right]}{2(1+R_{\tau})^{3/2}} + \frac{\sqrt{\pi}}{R_{\tau} \cdot r_w} \operatorname{erf}\left(\frac{r_w}{\sqrt{1+R_{\tau}}}\right) \quad (29)$$

$\frac{\exp\left[-\frac{(r_w)^2}{1+R_{\tau}}\right]}{2(1+R_{\tau})^{3/2}}$ is called part 1, and $\frac{\sqrt{\pi}}{R_{\tau} \cdot r_w} \operatorname{erf}\left(\frac{r_w}{\sqrt{1+R_{\tau}}}\right)$ is called part 2 for further use. $T_{max}$ is $T_{min} + T_{osc}$.

(2) For situation 2, even if $x_m$ influences the $T_{max}$ and $T_{osc}$, its influence is bounded. When $x_m = 0$ is used, we calculate the temperature at the beginning of the period and the maximum is thus larger (with a non-zero $x_m$). However, how much larger? In which situations should we care about it? From the analysis, the details are described in Appendix D, we found that the difference appears only around $r_w = 1.6$ to 4 when $R_\tau$ is large.

$T_{min}$ is the same as situation 1. For $T_{max}$, using the trapezoidal rule for approximation as for $r_w \neq 0$, we have $T_{max}$ from Equation (26):

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

$$\underset{N \gg 1/R_{\tau}}{\longrightarrow} \frac{\exp\left[-\frac{(r_w)^2}{1+x_m \cdot R_{\tau}}\right]}{[1+x_m \cdot R_{\tau}]^3} + \frac{\exp\left[-\frac{(r_w)^2}{1+(1+x_m) \cdot R_{\tau}}\right]}{2[1+(1+x_m) \cdot R_{\tau}]^3} + \frac{\sqrt{\pi}}{R_{\tau} \cdot r_w} \left\{ \operatorname{erf}\left[\frac{r_w}{\sqrt{[1+(1+x_m) \cdot R_{\tau}]}}\right] \right\} \quad (30)$$

$\frac{\exp\left[-\frac{(r_w)^2}{1+x_m \cdot R_\tau}\right]}{[1+x_m \cdot R_\tau]^{3/2}}$, $\frac{\exp\left[-\frac{(r_w)^2}{1+(1+x_m) \cdot R_\tau}\right]}{2[1+(1+x_m) \cdot R_\tau]^{3/2}}$, $\frac{\sqrt{\pi}}{R_\tau \cdot r_w} \left\{ \operatorname{erf}\left[\frac{r_w}{\sqrt{[1+(1+x_m) \cdot R_\tau]}}\right] \right\}$ are called part 1, 2, and 3, respectively.

The general expression of $T_{osc}$ (when $N$ tends to infinity or larger than the effective number for reaching the steady state) is given as $T_{max}$ (Equation (30)) minus $T_{min}$ (Equation (29)), and it reads:

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

Part 1 in Equation (29) and part 2 in Equation (30) are smaller than the other parts by a factor 10, so they can be approximately omitted to simplify the expressions in practice.
$$ T_{osc}(R_{\tau}, r_{w}) \approx \frac{\exp\left[-\frac{(r_{w})^{2}}{1+x_{m} \cdot R_{\tau}}\right]}{\left[1+x_{m} \cdot R_{\tau}\right]^{\frac{3}{2}}} + \frac{\sqrt{\pi}}{R_{\tau} \cdot r_{w}} \left\{ \operatorname{erf}\left[\frac{r_{w}}{\sqrt{1+(1+x_{m}) \cdot R_{\tau}}}\right] \right\} - \frac{\sqrt{\pi}}{R_{\tau} \cdot r_{w}} \operatorname{erf}\left(\frac{r_{w}}{\sqrt{1+R_{\tau}}}\right) \quad (32) $$

It is worth noticing (see Appendix D, Figure A4a,b) that when $R_{\tau}$ increases, $T_{osc}$ exhibits a small departure from the exact value at around $r_{w} = 2$, attributed to the existence of a non-zero $x_{m}$. This departure, if it is generally not negligible, is nevertheless bounded. It is calculated to be $\frac{1}{2.44(r_{w})^{3}}$ for a large $R_{\tau}$ (details are shown in Appendix D Figure A4c,d by plotting $T_{osc}$ according to $r_{w}$ and $R_{\tau}$). Therefore, the range of $T_{osc}$ is given by Equations (33) and (34):

$$ T_{osc}(R_{\tau}, r_{w}) \xrightarrow{R_{\tau} \to 0} \exp[-(r_{w})^{2}] \quad (33) $$

$$ T_{osc}(R_{\tau}, r_{w}) \xrightarrow{R_{\tau} \to \infty} \frac{1}{2.44(r_{w})^{3}} \quad (34) $$

We note that the oscillation amplitude $T_{osc}$ at situation 1 is $\exp[-(r_{w})^{2}]$ which is the minimum, while in situation 2, the amplitude is larger due to the influence of $x_{m}$, with a maximum value of $\frac{1}{2.44(r_{w})^{3}}$ at the place around $r_{w} = 2$. By now, with these temperature expressions, we obtain the spatial distribution of the minimum and maximum temperature for a given $R_{\tau}$ at steady state, shown in Figure 8. The temperature is oscillating between these two temperature profiles, and note that at $r_{w} = 0$, the difference is always 1 regardless of $R_{\tau}$.

We have now all the information for plotting the $T$ distribution with any $R_{\tau}$ value.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Spatial distribution of $T_{min}$ (blue dash, by Equation (29)) and $T_{max}$ (red, by Equation (30) and Equation (27)) according to the relative radius $r_{w}$ when (a) $R_{\tau} = 0.1$, (b) $R_{\tau} = 1$, and (c) $R_{\tau} = 10$.</figcaption></figure>

From Figure 8, we can see that when $R_{\tau}$ is small (large frequency or small diffusivity), $T_{min}$ and $T_{max}$ have no large relative difference compared to their average values because the oscillation amplitude is always limited to $\exp[-(r_{w})^{2}]$ whereas the Tmean amplitude is converging to $2/R_{\tau}$ (heat accumulation). This case is interesting if a rather stable temperature is requested. Then, pulse energy can be adjusted for compensating the pulse RR increase. It is also worth noting that, by decreasing $R_{\tau}$, the shape of the curve converges to the $\exp(r_{w})$ curve and decreases much slower than a Gaussian one.

For large values of $R_{\tau}$ (small frequency or large diffusivity), the oscillations are relatively large as the pulses are separated and thus $T_{min}$ appears to have small values. It is negligible (<3%) when $R_{\tau} > 16$ or <6% for $R_{\tau} > 7$). This limits the domain of heat accumulation. The calculation shows that the shape of $T_{max}$ is also converging with increasing $R_{\tau}$ to the shape of the beam energy distribution (Gaussian, here $\exp[-(r_{w})^{2}]$ independent of $R_{\tau}$. This translates that the maximum is almost whatever the radius, at the
beginning of the period. This is not true exactly only around $r = 2w$ where a few % departure from the Gaussian shape of $T_{max}$ is demonstrated.

For $R_τ$ intermediate values, the temperature oscillations are limited between $T_{max}$ and $T_{min}$. This is shown with particular cases with a shoe box in [31] for $R_τ = 2$ and 20 or in [32] for $R_τ = 20$.

However, in this paper, we regard that when $r_w > 2$, the difference between $T_{min}$ and $T_{max}$ is vanishing. We see this in [37].

Therefore, we can deduce in particular, in whatever situation of $r_w$ and $R_τ$, the temperature oscillations can be neglected and the use of an average temperature is applicable.

Other application remarks:

(1) In the intermediate cases around $R_τ = 1$, the center of the heat-affected zone experiences large temperature oscillations whereas the periphery temperature is not oscillating. This may induce differences in the modification structures along the radius. Specifically, the pedestal of the curve, borne by $T_{min}$, increases in width with $R_τ$ as $1.75 + \frac{1}{16}R_τ^{\frac{R_τ+145}{R_τ+20}}$.

(2) For smaller $R_τ$ values, during the transient period (before $N_{ss}$), the width of the temperature distribution starts with the beam waist (Gaussian) and then increases until a size which is defined by $R_τ$. It does not increase indefinitely over time. The order of magnitude is one w per two orders of magnitude on $R_τ$, e.g., the trace width at 1 MHz is twice the one at 10 kHz.

## The Effective Number of Pulses for Reaching the Temperature Limits

Since $x_m$ is not negligible in very limited circumstances, the effective numbers of pulses for reaching the limit of $T_{osc}$, $T_{min}$, and $T_{max}(N_{r'sso/ssmin/ssmax})$ are given in the situation when $x_m = 0$.

With the same definition as we calculated in $r_w = 0$, with $\epsilon$ being a small quantity and $X = \frac{1}{\sqrt{1+N/R_τ}}$, we have $\frac{|T_{osc/max/min}(r_w,N)-T_{osc/max/min}(r_w,\infty)|}{T_{osc/max/min}(r_w,\infty)} < \epsilon$. Therefore, the $N_{r'sso/ssmin/ssmax}$ solutions are shown below:

$$N_{sso}^r > \frac{-3 \cdot W \left[ -\frac{1}{3} e^{-\frac{1}{3}r_w^2} \cdot r_w^2 \cdot \epsilon^{\frac{2}{3}} \right] - 2r_w^2}{3 \cdot R_τ \cdot W \left[ \left( -\frac{1}{3} e^{-\frac{1}{3}r_w^2} \cdot r_w^2 \cdot \epsilon^{\frac{2}{3}} \right) \right]} = \frac{1}{R_τ} \left[ -\frac{2r_w^2}{3 \cdot W \left[ \left( -\frac{1}{3} e^{-\frac{1}{3}r_w^2} \cdot r_w^2 \cdot \epsilon^{\frac{2}{3}} \right) \right]} - 1 \right] \quad (35)$$

This expression does not have an analytic root without using the tabulated function $W$, i.e., the Lambert W function (defined as $\omega e^{\omega} = z$, $\omega = W(z)$). In practice, since $X^2 \ll 1$, by approximation, it becomes:

$$N_{sso}^r > \frac{1}{R_τ} \left[ \left( \frac{1}{\epsilon \cdot \exp[-(r_w)^2]} \right)^{\frac{2}{3}} - 1 \right] \quad (36)$$

For $N_{ssmin}^r$ and $N_{ssmax}^r$, with the approximation of $erf(X \cdot r_w) \approx \frac{2}{\sqrt{\pi}} X \cdot r_w$,

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

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

The behavior of $N_{ss}^r$ according to $R_τ$ for $r_w = 0$ has already been shown in Figure 5, with an overall increase. We have plotted $N_{sso}$ and $N_{ssmax}$, according to $r_w$, for $R_τ = 10$ as
an example, as shown in Figure 9a, and the plot of the related time for reaching the steady state (using $N_{ssmax}$ and with diffusion time 0.28 µs) in Figure 9b.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** (a–c) Effective number of $N_{ssmax}$ (red), $N_{ssmin}$ (blue dash), and $N_{sso}$ (green) for reaching a steady state according to $r_w$ from 0 to 5 when (a) $R_τ = 10$; (b) $R_τ = 1$; (c) $R_τ = 0.1$. (d) The time (in µs) for reaching the steady state according to $r_w$ from 0 to 5 when $R_τ = 0.1, 1$, and 10.</figcaption></figure>

From Figure 9, we can see that as $r_w$ increases, it takes more pulses (three orders of magnitude more) and this corresponds to a longer time to reach the steady state. Therefore, in reality, even though at the exact center of the beam, the temperature is stable, the periphery is still evolving. In particular, in the case of a moving beam, the maximum speed of scanning is limited by the change at the focus periphery increasing from zero.

4. The Mean Temperature in the Period between Two Pulses

For many transformations induced by laser irradiation (fictive temperature, crystallization, erasure of previously induced structures, stress relaxation, and so on), the large temperatures occurring within a pulse period are so brief that the system has no time to significantly respond. On the contrary, for smaller temperatures occurring at the end of the period, the system may have time to respond if the temperatures are not too small (this is the case for overlapping pulse contribution, i.e., heat accumulation). Therefore, the system responds efficiently predominantly for intermediate temperatures in the main part of the period. On the other hand, when $R_τ$ is small (large pulse RR versus diffusion time), temperature oscillations are relatively small whatever the radius, or for large radius values whatever $R_τ$ values (see Figure 8), the temperature oscillation can be neglected. For these, the use of an average temperature is relevant. In any case, the average value can be a guide for following the temperature distribution in space and its evolution. Hence, this section is devoted to simple expressions of average temperature values in the function of material and laser parameters.

We define the averaging by $\bar{T}(r, N) = \frac{1}{\tau_p} \int_{\text{pulse period at N}} T(r, t) dt$, and this gives:
$$ \bar{T}(r, N) = \frac{1}{\tau_p} \int_{\text{pulse period at } N} T(r, t) dt \qquad (39) $$

N.B. due to software problem, the average temperature is sometimes quoted as $\bar{T}$ and sometimes $T_{mean}$. They have the same meaning.

### 4.1. Temperature at the Center ($\bar{T}(0, N)$)

$$ \bar{T}(0, N) = \frac{1}{\tau_p} \int_{\frac{t}{\tau_p} = N-1}^{\frac{t}{\tau_p} = N} \sum_{n=0}^{N-1} \frac{1}{\left[1 + \left(\frac{t}{\tau_p} - n\right) \cdot R_{\tau}\right]^{\frac{3}{2}}} dt \qquad (40) $$

The two summations can be permuted as they do not operate on the same variables and are independent. We obtain:

$$ \bar{T}(0, N) = \frac{1}{\tau_p} \sum_{n=0}^{N-1} \int_{\frac{t}{\tau_p} = N-1}^{\frac{t}{\tau_p} = N} \frac{1}{\left[1 + \left(\frac{t}{\tau_p} - n\right) \cdot R_{\tau}\right]^{3/2}} dt = \frac{1}{R_{\tau}} \sum_{n=0}^{N-1} \left[ - \frac{2}{\left[1 + \left(\frac{t}{\tau_p} - n\right) \cdot R_{\tau}\right]^{1/2}} \right]_{\frac{t}{\tau_p} = N-1}^{\frac{t}{\tau_p} = N} = \frac{2}{R_{\tau}} \left(1 - \frac{1}{\sqrt{1 + N \cdot R_{\tau}}}\right) \qquad (41) $$

This result is obtained without approximation. Then, when $N \gg 1/R_{\tau}$:

$$ \bar{T}(0, \infty) = \lim_{N \gg 1/R_{\tau}} \bar{T}(0, N) = \frac{2}{R_{\tau}} \qquad (42) $$

We note that here, the steady state temperature at the center will reach:

$$ \frac{2 \cdot T_{00}}{R_{\tau}} = \frac{2A \cdot P}{\pi^{3/2} \kappa W} \qquad (43) $$

It is the same expression as for $T_{max}(0, \infty)$ or $T_{min}(0, \infty)$ for small $R_{\tau}$ values. We can note in Figure 10 that $T_{max}(0, \infty)$ and $T_{min}(0, \infty)$ approach $\bar{T}(0, \infty)$ when $R_{\tau}$ is decreasing. $T_{max}(0, \infty)$ goes to 1 and $T_{min}(0, \infty)$ goes to 0 when $R_{\tau}$ is large. From Figure 10, we can also find the heat accumulation bound already defined in Figure 5 (with $\epsilon = 0,03$). It corresponds to $T_{min}(0, \infty) = 0,03$ and $R_{\tau} = 12$. On the other hand, when $T_{min}(0, \infty)$ depends from $T_{max}(0, \infty)$ by less than approximately 10%, we can admit that the average T is applicable, i.e., for $R_{\tau}$ smaller than 0.17 (purple circle). In this case, we can apply the simple expression Equation (43).
<figure><img src="image_13.png" /><figcaption>Figure 10. The plots of $T_{min}(0, \infty)$ (red), $T_{max}(0, \infty)$ (blue dash), and $\bar{T}(0, \infty)$ ($T_{mean}$, green dash) according to $1/R_{\tau}$. The defined boundary points of heat accumulation and negligible oscillation in the system are marked by a blue circle and purple circle, respectively (for definitions, see text).</figcaption></figure>

The Effective Number of Pulses for Reaching the Limit $\bar{T}(0, \infty)$ ($N_{ssm}^0$)

With the same definition as above, the number of pulses to reach $\bar{T}(0, \infty)$, i.e., $\frac{|\bar{T}(0, \infty)-\bar{T}(0,N)|}{\bar{T}(0, \infty)} < \varepsilon$, $N_{ssm}^0$ is obtained:

$$N_{ssm}^0 > \frac{1}{R_{\tau}} \left[ \left( \frac{1}{\varepsilon} \right)^2 - 1 \right] \quad (44)$$

$N_{ssm}^0$ are compared to $N_{ssmax}^0$ and $N_{sso}^0$ in Figure 11. The steady state of the mean temperature is reached as $T_{max}$.

<figure><img src="image_14.png" /><figcaption>Figure 11. The effective number of pulses to reach the limit of $T_{osc}$ (red), $T_{min}$ (blue dash), $T_{max}$ (green dash), and $T_{mean}$ (pink dash) according to $R_{\tau}$ from 0 to 5.</figcaption></figure>

Therefore, the time for reaching the steady state here is not $R_{\tau}$ dependent: $\tau_D \left[ \left( \frac{1}{\varepsilon} \right)^2 - 1 \right] \approx \frac{\tau_D}{\varepsilon^2} = 1111\tau_D$ when $\varepsilon = 0.03$. It is the value of the maximum $t_{ss}$ for reaching a steady state when $R_{\tau} \to 0$ (Figure 6).

With these analytical expressions of temperature at the steady state at the center of the focus ($T_{osc}$, $T_{min}(0, \infty)$, $T_{max}(0, \infty)$, $\bar{T}(0, \infty)$), and the needed number of pulses ($N_{sso}^0$, $N_{ssmax}^0$, $N_{ssm}^0$), we have a clear view of how parameter $R_{\tau}$ influences the thermal
situation at the focus center. The problem is now to extend these results to any place out of the center.

## 4.2. Temperature out of the Focus Center ($T(r, N)$)

With the definition $\bar{T}(r, N) = \frac{1}{\tau_p} \int_{\text{pulse period at N}} T(r, t) dt$, we have the average temperature in a period as:

$$ \bar{T}(r_w, N) = \frac{1}{\tau_p} \int_{\tau_p}^{\tau} \sum_{n=0}^{N-1} \frac{1}{\left[1 + \left(\frac{t}{\tau_p} - n\right) \cdot R\tau\right]^2} \exp\left[-\frac{(r_w)^2}{1 + \left(\frac{t}{\tau_p} - n\right) \cdot R\tau}\right] \cdot dt = \frac{\sqrt{\pi}}{R\tau \cdot r_w} \cdot \left[\text{erf}(r_w) - \text{erf}\left(\frac{r_w}{\sqrt{1 + N \cdot R\tau}}\right)\right] \underset{N \gg R\tau}{\longrightarrow} \frac{\sqrt{\pi}}{R\tau \cdot r_w} \text{erf}(r_w) \quad (45) \\ \text{So } \bar{T}(r_w, \infty) = \frac{\sqrt{\pi}}{R\tau \cdot r_w} \text{erf}(r_w) $$

This limit when $N \gg 1/R_τ$ is shown in Figure 12 and compared to the Gaussian shape of $T_{max}$ when $R_τ$ is large and when $R_τ$ is small. When $T_{max}(r)$ is Gaussian for the first case, $T_{max}(r)$ has the same shape that $T_{mean}$ has for the second case. As the erf function tends to 1 (already for $r_w > 2$), the function tends to be hyperbolic and thus decreases much slower than a Gaussian one (see Figure 12). The amplitude is $\frac{2}{R_τ}$. It is inversely proportional to $R_τ$ whatever $R_τ$ value.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Plot of normalized $T_{max}(r_w, R_τ)$ and $\bar{T}(r_w, R_τ)$ at the steady state for $R_τ = 0.1$ and 100 for $T_{max}$ and 1 for $T_{mean}$.</figcaption></figure>

Consistently with the previously calculated $T_{max}$ and $T_{min}$, the $T_{mean}$ curve width is equal to the beam Gaussian for large $R_τ$ values and to the curve limit given in Equation 45 and shown in Figure 8a which is wider.

The Effective Number of Pulses for Reaching the Limit of $T_{mean}$ ($N_{ssm}^r$)

The effective number of pulses $N_{ssm}^r$ with the approximation $\text{erf}(X \cdot r_w) \approx \frac{2}{\sqrt{\pi}} X \cdot r_w$ as $X(N, R_τ) = \frac{1}{\sqrt{1+N \cdot R_τ}} < 1$, is solved to be:

$$ N_{ssm}^r(r_w) > \frac{1}{R_τ} \left[ \left( \frac{2}{R_τ \cdot \varepsilon \cdot \bar{T}(r_w, R_τ)} \right)^2 - 1 \right] = \frac{1}{R_τ} \left[ \left( \frac{2 \cdot r_w}{\varepsilon \cdot \sqrt{\pi} \cdot \text{erf}(r_w)} \right)^2 - 1 \right] \quad (46) $$

From the expression above, we see that the periphery of the distribution is stabilized later than the center as we noticed already in the previous section.
## 5. Application Examples

To demonstrate the practical significance of the aforementioned calculations, we are discussing below several problems where we can apply these equations to analyze the temperature effects.

**Laser-induced crystallization.** It is known that for crystallizing a glass, it is necessary to control temperature and time in order to penetrate the crystallization domain [38]. A method for reaching it with a pulsed laser is described in [39]. It is shown that crystallization with a single pulse is possible from the solid state if the beam scanning speed is sufficiently low according to the nucleation time and crystallization growth rate. For a larger scanning speed, it is necessary to increase the pulse energy or the RR to reach the crystallization domain. This is for the formation of nanocrystals that are orientable with laser polarization. The decrease in the speed leads to the growth of the nanocrystals. In turn, crystallization is still possible if the speed is increased but the pulse energy should be increased. In such a way, the temperature overcomes the melting one during a time long enough in the pulse period and the material is melted in such a way that crystallization does not progress more after each pulse. From the calculations in this paper, the best method appears to be a high RR with moderate pulse energy in order to maintain T (control of $T_{max}$ and $T_{min}$) around the crystallization temperature.

**Erasure process during laser writing.** In the case of pure silica, there is a first regime called type I for which the refractive index increases for pure silica glasses [40,41]. It is partly produced by a change in fictive temperature [42,43]. For that, the time for the temperature to decrease until a given value of temperature has to be larger than the relaxation time of the glass. This time is roughly defined by the cooling time which is itself defined by the moving spatial curve for one pulse [44]. Pulse energy can be adjusted consequently.

In the case of the materials in which type II transformation is achievable giving rise to a large birefringence based on self-organized nanograting (NG) and nanopores, there is a pulse-energy-RR-scanning-speed-related domain [45,46]. This domain is limited for large pulse energies depending on RR. In addition, for large RR values, the retardance decreases until it is no longer possible to write NGs. One of the hypotheses is the following: NG is based on the existence of nanopores distributed in a self-organized NG [9]. Recent work [47] shows that the thermal stability of such an object is defined by the viscosity that itself depends on T. Therefore, as T increases when pulse energy and RR are increased, they have to be limited to avoid an in-pulse erasure after creation during the pulse.

**Concurrent processes.** In organic materials, according to pulse energy and RR for the same mean power, two different processes are observed whereas we could believe that modification is just dependent on mean power (i.e., dose) [48]. One is the destruction of the material at low RR and high $E_p$, while the other is the creation of luminophores at high RR and low $E_p$. We explain this by looking at the amplitude of the T oscillations; we can say that for the first case, the oscillations are large whereas for the other case, they are small. In the first case, the temperature is overcoming the decomposition temperature of the material but not in the second case.

**Size of the crystallized trace.** In [49], we find a size of the heat-affected zone that is much larger for the glass called Silica-SrTiO$_3$ than for Silica-LiNbO$_3$ for comparable laser parameters. $R_τ$ for the first is 84 and for the second is 11. The first remark is that there is almost no heat accumulation (separated pulse contribution to the temperature) especially for the STS glass. For this glass, it is even possible to use the formula for one contribution. Then, the size is defined by the lowest maximal temperature according to the radius (Equation (30)) at least larger than $T_g$. Note that due to the expected size of the trace, the use of $T_{mean}$ is possible for intermediate $R_τ$ values (Equation (45)).

$$ \frac{\sqrt{\pi}}{R\tau \cdot r_w} erf(r_w) = \bar{T}(r_w, \infty) = \frac{T_g}{T_{00}} $$
**Crown effect.** In [49], we also find an example of a crystallized shell (the center is not crystallized). Similarly, as above, there is a highest maximal temperature that is equal to $T_{melting}$. Above this temperature, the viscosity decreases strongly and other processes may appear (see [49]) to be blocking crystallization on cooling.

Speed effect on the laser trace width. In glasses for mid IR [50], the energy threshold for the appearance of a sudden spatial broadening depends on the scanning speed, so we can deduce that it is not related to temperature as the writing speed is not involved in the thermal diffusion for a speed lower than a few m/s.

There are also some remarks that we can deduce from the calculation:

- If a process is actually independent of RR, it does not depend on temperature (see [46]);

- The number of pulses received by the material locally depends on the scanning speed. As the number of pulses for reaching a steady state is different at the center than at the periphery, it is possible that the appearance of the trace on the edge depends on the scanning speed during the transient stage;

- However, the transient stage is not dependent on the pulse energy.

## 6. Conclusions

This study derives analytical expressions for the temperature distribution at the steady state induced by ultrafast multi-pulses within a spherical geometry, based on the laser and constant material parameters. These expressions depend only on two parameters: the initial temperature at the center (denoted as $T_{00}$) and a quantity $R_{\tau}$, defined as the ratio of the pulse period $\tau_p$ to the diffusion time $\tau_d$. We recall that temperature oscillates between $T_{max}$ and $T_{min}$, eventually reaching a steady state, and we calculate the minimum number of pulses required to attain this state. For ease of use, a geometry with spherical symmetry was chosen for the energy deposition in order to lead to simple temperature expressions compiled in Table 3. The Table 4 is a further simplification usable in most of the cases. This approach also facilitates a clear and precise definition of the onset of heat accumulation.

<table><caption>**Table 3.** Analytical expressions of final temperature and the effective number for reaching steady state. *: condition of $r_w$, $R_\tau$ (Equation (27), Figure A2).</caption><tbody><tr><td colspan="3" style="text-align:center; font-weight:bold;">r<sub>w</sub> = 0</td></tr><tr><td></td><td style="text-align:center; font-weight:bold;">T(0, N)</td><td style="text-align:center; font-weight:bold;">T(0, N → ∞)</td><td style="text-align:center; font-weight:bold;">N<sup>0</sup><sub>ss</sub>(ε, R<sub>τ</sub>)</td></tr><tr><td style="text-align:left; font-weight:bold;">T<sub>osc</sub>(0, N)</td><td style="text-align:center;">1 - X<sup>3</sup></td><td style="text-align:center;">1</td><td style="text-align:center;">1/R<sub>τ</sub> [((1/ε)<sup>2/3</sup> - 1)]</td></tr><tr><td style="text-align:left; font-weight:bold;">T<sub>min</sub>(0, N)</td><td style="text-align:center;">1/2 [X1<sup>3</sup> + X<sup>3</sup>] + 2/R<sub>τ</sub> [X1 - X]</td><td style="text-align:center;">1/2 X1<sup>3</sup> + 2/R<sub>τ</sub> X1</td><td style="text-align:center;">1/R<sub>τ</sub> [((2/R<sub>τ</sub> · ε · T<sub>min</sub>(∞))<sup>2</sup> - 1)]</td></tr><tr><td style="text-align:left; font-weight:bold;">T<sub>max</sub>(0, N)</td><td style="text-align:center;">1 + T<sub>min</sub>(N)</td><td style="text-align:center;">1 + T<sub>min</sub>(∞)</td><td style="text-align:center;">1/R<sub>τ</sub> [((2/R<sub>τ</sub> · ε · T<sub>max</sub>(∞))<sup>2</sup> - 1)]</td></tr><tr><td style="text-align:left; font-weight:bold;">T&#773;(0, N)</td><td style="text-align:center;">2/R<sub>τ</sub> (1 - X)</td><td style="text-align:center;">2/R<sub>τ</sub></td><td style="text-align:center;">1/R<sub>τ</sub> [((1/ε)<sup>2</sup> - 1)]</td></tr><tr><td colspan="4" style="text-align:center; font-weight:bold;">r<sub>w</sub> ≠ 0 *</td></tr><tr><td></td><td style="text-align:center; font-weight:bold;">T(r, N)</td><td style="text-align:center; font-weight:bold;">T(r, N → ∞)</td><td style="text-align:center; font-weight:bold;">N<sup>r</sup><sub>ss</sub>(ε, R<sub>τ</sub>)</td></tr><tr><td style="text-align:left; font-weight:bold;">T<sub>osc</sub>(r, N)</td><td style="text-align:center;">exp[-(r<sub>w</sub>)<sup>2</sup>] - X<sup>3</sup>exp[-(X · r<sub>w</sub>)<sup>2</sup>]</td><td style="text-align:center;">exp[-(r<sub>w</sub>)<sup>2</sup>]</td><td style="text-align:center;">1/R<sub>τ</sub> [((1/ε · exp[-(r<sub>w</sub>)<sup>2</sup>])<sup>2/3</sup> - 1)]</td></tr><tr><td style="text-align:left; font-weight:bold;">T<sub>min</sub>(r, N)</td><td style="text-align:center;">1/2 [X1<sup>3</sup>exp[-(X1 · r<sub>w</sub>)<sup>2</sup>]<br/>+ X<sup>3</sup>exp[-(X · r<sub>w</sub>)<sup>2</sup>]]<br/>+ √π/R<sub>τ</sub> · r<sub>w</sub> [erf(X1 · r<sub>w</sub>) - erf(X · r<sub>w</sub>)]</td><td style="text-align:center;">1/2 X1<sup>3</sup>exp[-(X1 · r<sub>w</sub>)<sup>2</sup>]<br/>+ √π/R<sub>τ</sub> · r<sub>w</sub> erf(X1 · r<sub>w</sub>)</td><td style="text-align:center;">1/R<sub>τ</sub> [((2/R<sub>τ</sub> · ε · T<sub>min</sub>(∞))<sup>2</sup> - 1)]</td></tr></tbody></table>
<table><caption>**Table 4.** Practical approximated analytical expressions of final temperature and the effective number for reaching steady state. With $X(N, R_τ) = \frac{1}{\sqrt{1+N·R_τ}}$ and $X1(R_τ) = \frac{1}{\sqrt{1+R_τ}}$. *: condition of $r_w, R_τ$ (Equation (27), Figure A2).</caption><tr><td>T<sub>max</sub>(r, N)</td><td>exp[-(r<sub>w</sub>)<sup>2</sup>] + T<sub>min</sub>(N)</td><td>exp[-(r<sub>w</sub>)<sup>2</sup>]+T<sub>min</sub>(∞)</td><td>1/R<sub>τ</sub> [((2/(R<sub>τ</sub>·ε·T<sub>max</sub>(∞)))<sup>2</sup> - 1)]</td></tr><tr><td>T&#773;(r, N)</td><td>√π/R<sub>τ</sub>·r<sub>w</sub> [erf(r<sub>w</sub>) - erf(X·r<sub>w</sub>)]</td><td>√π/R<sub>τ</sub>·r<sub>w</sub> erf(r<sub>w</sub>)</td><td>1/R<sub>τ</sub>τ [(2·r<sub>w</sub>/(ε√π·erf(r<sub>w</sub>)))<sup>2</sup> - 1]</td></tr></table>

<table><thead><tr><th></th><th colspan="2">r<sub>w</sub> = 0</th><th colspan="2">r<sub>w</sub> ≠ 0 *</th></tr><tr><th></th><th>T(0, N → ∞)</th><th>N<sup>0</sup><sub>ss</sub>(ε, R<sub>τ</sub>)</th><th>T(r, N → ∞)</th><th>N<sup>r</sup><sub>ss</sub>(ε, R<sub>τ</sub>)</th></tr></thead><tbody><tr><td>T<sub>osc</sub>(r, N)</td><td>1</td><td>1/R<sub>τ</sub> [(1/ε)<sup>2/3</sup> - 1]</td><td>exp[-(r<sub>w</sub>)<sup>2</sup>]*</td><td>1/R<sub>τ</sub> [(1/ε·exp[-(r<sub>w</sub>)<sup>2</sup>])<sup>2/3</sup> - 1]</td></tr><tr><td>T<sub>min</sub>(r, N)</td><td>2/R<sub>τ</sub> X1</td><td>1/R<sub>τ</sub> [(2/(R<sub>τ</sub>·ε·T<sub>min</sub>(∞)))<sup>2</sup> - 1]</td><td>√π/R<sub>τ</sub>·r<sub>w</sub> erf(X1·r<sub>w</sub>)</td><td>1/R<sub>τ</sub> [(2/(R<sub>τ</sub>·ε·T<sub>min</sub>(∞)))<sup>2</sup> - 1]</td></tr><tr><td>T<sub>max</sub>(r, N)</td><td>1+T<sub>min</sub>(∞)</td><td>1/R<sub>τ</sub> [(2/(R<sub>τ</sub>·ε·T<sub>max</sub>(∞)))<sup>2</sup> - 1]</td><td>exp[-(r<sub>w</sub>)<sup>2</sup>]+T<sub>min</sub>(∞)</td><td>1/R<sub>τ</sub> [(2/(R<sub>τ</sub>·ε·T<sub>max</sub>(∞)))<sup>2</sup> - 1]</td></tr><tr><td>T&#773;(r, N)</td><td>2/R<sub>τ</sub></td><td>1/R<sub>τ</sub> [(1/ε)<sup>2</sup> - 1]</td><td>√π/R<sub>τ</sub>·r<sub>w</sub> erf(r<sub>w</sub>)</td><td>1/R<sub>τ</sub> [(2·r<sub>w</sub>/(ε√π·erf(r<sub>w</sub>)))<sup>2</sup> - 1]</td></tr></tbody></table>

We analyzed the distribution of the temperature oscillations relative to the radius from the center and the parameter $R_τ$. Oscillations are large at the center for large $R_τ$ values but decrease strongly for a large radius $r_w > 2$, i.e., for the periphery where the light intensity decreases almost by a factor of 10. On the contrary, oscillations are minimal everywhere for small $R_τ$ values (i.e., high frequency or low thermal diffusivity). In such conditions, the average of the temperature from the last period can be used, yielding even simpler expressions. Additionally, we found that the periphery of the focus reaches the steady state later than the center. By examining the pulse number required for the steady state according to the radius, we can better control transformations in these regions and understand the variations from the center.

This work aids in understanding how temperature variations influence different experimental observations, mentioned at the end. It can also be helpful to detect if temperature is acting on the processes of direct laser writing.

Future work includes refining this approach by considering the asymmetry of fs focus, making differences between transversal radius and depth to deduce how the trace changes over time. Another interesting point is the T dependence of the physical-chemical parameters, but this cannot be investigated without finite element calculation. Our approach allows us to choose the most representative parameters for applying such a calculation. In addition, the asymmetry of the focal volume, along and perpendicular to the propagation axis, is a refinement of the present calculations that we could be interested in to include variations with the pulse energy (like Kerr focusing, plasma density defocusing, plasma mirror).

**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, R.Q.; writing—original draft, B.P. and R.Q.; writing—review and editing, R.Q., B.P. 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. R.Q. acknowledges the China Scholarship Council (CSC) for the funding No. 201808440317 of her PhD fellowship.
**Data Availability Statement:** Data is contained within the article.

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

## Appendix A

<table><caption>**Table A1.** Glossary of laser and material parameters.</caption><thead><tr><th>Parameters</th><th>Definitions</th><th>Units</th></tr></thead><tbody><tr><td>A</td><td>Fraction of reflected light by the plasma</td><td>none</td></tr><tr><td>α</td><td>Light absorption</td><td>μm<sup>-1</sup></td></tr><tr><td>ε</td><td>A small quantity of computational needs</td><td>none</td></tr><tr><td>τ<sub>p</sub></td><td>Period of the pulses</td><td>μs</td></tr><tr><td>τ<sub>D</sub></td><td>Heat diffusion time τ<sub>D</sub> = <sup>w<sup>2</sup></sup>&frasl;<sub>4D<sub>T</sub></sub></td><td>μs</td></tr><tr><td>R<sub>τ</sub></td><td>τ<sub>p</sub>/τ<sub>D</sub></td><td>none</td></tr><tr><td>w</td><td>Beam waist radius (at 1/e)</td><td>μm</td></tr><tr><td>D<sub>T</sub></td><td>Thermal diffusivity D<sub>T</sub> = <sup>κ</sup>&frasl;<sub>ρ·c<sub>p</sub></sub></td><td>m<sup>2</sup>/s</td></tr><tr><td>κ</td><td>Thermal conductivity</td><td>W/(m·K)</td></tr><tr><td>E<sub>p</sub></td><td>Pulse energy</td><td>J</td></tr><tr><td>f</td><td>Pulse repetition rate</td><td>MHz</td></tr><tr><td>ρ</td><td>Density</td><td>kg/m<sup>3</sup></td></tr><tr><td>C<sub>p</sub></td><td>Specific heat capacity</td><td>J/(kg·K)</td></tr></tbody></table>

<table><caption>**Table A2.** Thermal physico-chemical data of some materials and processing parameters.</caption><thead><tr><th></th><th>ρ (kg/m<sup>3</sup>)</th><th>C<sub>p</sub> (J/(kg·K))</th><th>κ (W/(m·K))</th><th>D<sub>T</sub> (m<sup>2</sup>/s)</th><th>τ<sub>D</sub> (μs)</th><th>Melting Point (K)</th></tr></thead><tbody><tr><td>STS glass</td><td>3887</td><td>410</td><td>10.1</td><td>6.34 × 10<sup>-6</sup></td><td>0.039</td><td>1585</td></tr><tr><td>LNS glass</td><td>3830</td><td>650</td><td>2.65</td><td>1.06 × 10<sup>-6</sup></td><td>0.235</td><td>1530</td></tr><tr><td>SiO<sub>2</sub>(glass)</td><td>2200</td><td>703</td><td>1.38</td><td>8.92 × 10<sup>-7</sup></td><td>0.28</td><td>1983</td></tr><tr><td>Borosilicate (Schott D263) [40]</td><td>2510</td><td>820</td><td>0.96</td><td>4.66 × 10<sup>-7</sup></td><td>0.534</td><td>1324</td></tr><tr><td>Glycine</td><td>1160.7</td><td>1266</td><td>1.3 [51]</td><td>8.85 × 10<sup>-7</sup></td><td>0.283</td><td>506 (decomp.)</td></tr><tr><td>Zeonex</td><td>1010</td><td>1000</td><td>0.045</td><td>4.445 × 10<sup>-8</sup></td><td>5.624</td><td>553</td></tr><tr><td>Nifedipine</td><td>1300</td><td>1000</td><td>0.2</td><td>1.54 × 10<sup>-7</sup></td><td>1.63</td><td>446</td></tr><tr><td>Sucrose</td><td>1587</td><td>1243.1</td><td>0.1</td><td>5.07 × 10<sup>-8</sup></td><td>4.93</td><td>458 (decomp.)</td></tr></tbody></table>

## Appendix B. The Existence of a Bound for $T_{max}$ and $T_{min}$

Does $T_{max}$ theoretically have a limitation?

Since the temperature is the sum of the temperature contribution induced by each pulse, as shown in Equation (10), we need to know the convergence of this sum when the time or pulse number $N$ increases to infinity. A condition necessary but not sufficient is that the increase between the $T_{max}$ just after $N + 1$ pulses and just after $N$ pulses, i.e., $t_N = (N - 1) \tau_p$ and $t_{N+1} = (N) \tau_p$, is tending to 0. It is:

$$ 
\begin{aligned}
T_{max}(0, t_{N+1}) - T_{max}(0, t_N) &= \sum_{n=0}^{N} \frac{1}{[1 + (N - n) \cdot R_{\tau}]^{\frac{3}{2}}} - \sum_{n=0}^{N-1} \frac{1}{[1 + (N - 1 - n) \cdot R_{\tau}]^{\frac{3}{2}}} \\
&= \frac{1}{[1 + N \cdot R_{\tau}]^{\frac{3}{2}}}
\end{aligned} 
\quad (A1) $$

When $R_{\tau}$ is fixed, Equation (A1) goes to 0 when $N$ tends to infinity. It is the same if we consider the difference of $T_{min}$.
On the other hand, the sum expression of Equations (11) and (12) can be proved to be convergent by using the p-series test, since the p value in them (p = 3/2) is larger than 1 [52].

## Appendix C. Expression Approximation

To calculate the limit of the temperature, we start with the evolution of the minimum temperature on each period as long as N increases to infinity (so, with Equation (12), the limit of $T_{min}(0, t_{N+1})$ with N ≫ 1/R$_{τ}$).

Here, we introduce the approximation derived from the trapezoidal rule for the calculation of the integral, which is:

$$ \sum_{n=0}^{N-1} f(n) \approx \frac{f(0) + f(N-1)}{2} + \int_{0}^{N-1} f(n) dn \quad (A2) $$

Here,

$$ f(n) = \frac{1}{[1 + (N - n) \cdot R_{\tau}]^{\frac{3}{2}}} $$

with an error smaller than $(\frac{N-1}{12N^3}) f''(n)n$.

Therefore, with this approximation, we obtain:

$$ T_{min}(0, N) \approx \frac{1}{2(1 + R_{\tau})^{\frac{3}{2}}} + \frac{1}{2(1 + N \cdot R_{\tau})^{\frac{3}{2}}} + \frac{2}{R_{\tau}} \left[ \frac{1}{\sqrt{1 + R_{\tau}}} - \frac{1}{\sqrt{1 + N \cdot R_{\tau}}} \right] \quad (A3) $$

When N goes to infinity, $T_{min}(0, N)$ goes to:

$$ T_{min}(0, \infty) = \frac{1}{2(1 + R_{\tau})^{\frac{3}{2}}} + \frac{2}{R_{\tau}} \frac{1}{\sqrt{1 + R_{\tau}}} \approx \frac{2}{R_{\tau}} \frac{1}{\sqrt{1 + R_{\tau}}} \quad (A4) $$

$T_{max}$: The same method has been applied to obtain the $T_{max}$ limit but with the last term from the sum extracted out (due to the nature of the trapezoidal rule, there is a sharp increment in the temperature of the last pulse, which should not be averaged out):

$$ \sum_{n=0}^{N-1} f(n) \approx f(N-1) + \frac{f(0) + f(N-2)}{2} + \int_{0}^{N-2} f(n) dn \quad (A5) $$

$$ f(n) = \frac{1}{[1 + (N - 1 - n) \cdot R_{\tau}]^{\frac{3}{2}}} $$

We have thus:

$$ T_{max}(0, N) \approx 1 + \frac{1}{2[1 + R_{\tau}]^{3/2}} + \frac{1}{2[1 + (N - 1)R_{\tau}]^{3/2}} + \frac{2}{R_{\tau}} \left[ \frac{1}{\sqrt{[1 + R_{\tau}]} - \sqrt{[1 + (N - 1)R_{\tau}]}} \right] \quad (A6) $$

$T_{max}$ can also be obtained by the $T_{min}$ Equation (A4) + $T_{osc}$ (Equation (13)), thus:

$$ T_{max}(0, N) = 1 + \frac{1}{2[1 + R_{\tau}]^{3/2}} - \frac{1}{2[1 + N \cdot R_{\tau}]^{3/2}} + \frac{2}{R_{\tau}} \left[ \frac{1}{\sqrt{[1 + R_{\tau}]} - \sqrt{[1 + N \cdot R_{\tau}]}} \right] \quad (A7) $$
Both Equation (A6) or Equation (A7), when $N$ tends to infinity:

$$T_{max}(0, \infty) = 1 + \frac{1}{2(1 + R_{\tau})^{\frac{3}{2}}} + \frac{2}{R_{\tau}} \frac{1}{\sqrt{1 + R_{\tau}}} \approx 1 + \frac{2}{R_{\tau}} \frac{1}{\sqrt{1 + R_{\tau}}} \quad (A8)$$

The errors: For $T_{min}$, the error is smaller than $(N-1)^3 / 12N^3$ $f''(n)$ with $f''(n) = \frac{15 \cdot R_{\tau}^2}{4[1+(N-n) \cdot R_{\tau}]^{7/2}}$ with $0 \le n \le N - 1$ [53,54] and has to be compared with $T_{min}(0, \infty) = \frac{1}{2[1+R_{\tau}]^{3/2}} + \frac{2}{R_{\tau}} \frac{1}{\sqrt{1+R_{\tau}}}$. The error is smaller than 1.1% for any value of $R_{\tau}$, as shown in Figure A1a.

For $T_{max}$, we have an error smaller than $(N-2)^3 / 12(N-1)^3$ $f''(n)$ with $f''(n) = \frac{15 \cdot R_{\tau}^2}{4[1+(N-1-n) \cdot R_{\tau}]^{7/2}}$ with $0 \le n \le N - 2$. The error has to be compared with $\frac{T_{max}(0, \infty)}{T_0} = 1 + \frac{2}{2[1+R_{\tau}]^{3/2}} + \frac{2}{R_{\tau}} \frac{1}{\sqrt{1+R_{\tau}}}$. The error is smaller than 1.7% for any value of $R_{\tau}$, as shown in Figure A1b.

<figure><img src="image_16.png" /><figcaption>**Figure A1.** The error in using the trapezoidal rule for the approximation of computing. (a) $T_{min}$ and (b) $T_{max}$ versus $R_{\tau}$.</figcaption></figure>

## Appendix D. About Situation 2 in Section 3.2.1: When the $T_{max}$ is in the Middle of a Pulse Period

To calculate the oscillation amplitude $T_{osc}$ for $r \ne 0$, same as the case of $r_w = 0$, in general, we compare the difference between the maximum $T$ and minimum $T$ of the $N$th pulse period. $T_{min}$ is considered still at the end of the $N$th pulse, i.e., when $t = N \cdot \tau_p$ without the $(N+1)$th pulse, so:

$$T_{min}(r_w, N) = \sum_{n=0}^{N-1} \frac{1}{[1 + (N - n) \cdot R_{\tau}]^{\frac{3}{2}}} \exp \left[ - \frac{(r_w)^2}{1 + (N - n) \cdot R_{\tau}} \right] \quad (A9)$$

while $T_{max}$ in some situations can be in the middle of the pulse period, we thus set $x_m$, $0 \le x_m \le 1$ to define the place where the $T_{max}$ is. Therefore, $T_{max}$ is at $t = (N - 1 + x_m)\tau_p$:

$$T_{max}(r_w, N) = \sum_{n=0}^{N-1} \frac{1}{[1 + (N - 1 + x_m - n) \cdot R_{\tau}]^{\frac{3}{2}}} \exp \left[ - \frac{(r_w)^2}{1 + (N - 1 + x_m - n) \cdot R_{\tau}} \right] \quad (A10)$$

This expression is the general expression to describe both $T_{max}$ and $T_{min}$, while $T_{min}$ appearing at the end of the period means $x_m = 1$, as well as the case when $r_w = 0$, $T_{max}$ appears at the beginning of the period with $x_m = 0$.
To find the place of the maximum in the period, we search for the root $x_m$ of the derivative of Equation (A10). We obtain the root $x_m$ as a function of $R_τ$ and $r_w$ and N, but to know only the final situation, and also for mathematical solvability, $x_m$ is described as Equation (A11) as $N \gg 1/R_τ$. The $x_m$ values according to $R_τ$, $r_w$ are plotted in Figure A2.

$$x_m = \frac{\sqrt{R_τ\sqrt{9R_τ + 32(r_w)^2} - 3R_τ - 8}}{8R_τ} \quad (A11)$$

<figure><img src="image_17.png" /><figcaption>**Figure A2.** 2D map of the solution $x_m$ according to $R_τ$ and the reduced radius ($r/w$). Lines correspond to iso $x_m$.</figcaption></figure>

In Figure A2, $x_m = 0$ is located in the center of this map (at $r_w = 0$), the boundary is between dark green and green, below it $x_m$ is negative and it should be considered to be 0 (this is the maximum T location), i.e., the blue and purple area is also $x_m = 0$. In this parameter region ($r_w, R_τ$), the maximum of temperature is at the beginning of the period. From green to yellow then orange, the $x_m$ value increases to 0.5 then to 1, i.e., the maximum of temperature is in the middle of the period or even at the end. The boundary between yellow-green and orange is the boundary of $x_m = 1$, thus with the parameters in the orange region and red region, $x_m = 1$. When in the parameter region of ($r_w, R_τ$), the $x_m > 1$ is also considered to be 1. The thermal calculation can be divided into two situations:

(1) $x_m = 0$, for $r_w < \sqrt{\frac{3}{2} + \frac{2}{R_τ}}$ or $R_τ$ small enough (less than $\frac{2}{r_w^2-1.5}$) when $r_w^2 > 1.5$). In this situation, $x_m$ can be omitted from the expression.

(2) $x_m \neq 0$, for $r_w > \sqrt{\frac{3}{2} + \frac{2}{R_τ}}$ or $R_τ > \frac{2}{r_w^2-1.5}$ when $r_w^2 > 1.5$, $x_m$ will be appearing in the temperature expressions.

To further analyze the condition in situation 2 that $x_m$ can be set to 0, since $x_m$ only influences the value of $T_{max}$ and $T_{osc}$, we analyze $T_{max}$ by Equation (A10) with $x_m$ in it and
the one without $x_m$. Figure A3a,b displays the spatial distribution of $T_{max}$ with the $x_m$ of the value depending on $r_w$ (red) and the one with $x_m = 0$ for whatever $r_w$ (blue dash).

<figure><img src="image_18.png" /><figcaption>**Figure A3.** (a,b) Spatial distribution of $T_{max}$ with the exact expression ($x_m$ is a function of $R_\tau, r_w$) and simplified expression ($x_m = 0$). (a) $R_\tau = 1$, (b) $R_\tau = 10$. (c) The differences between the two expressions of $R_\tau = 10$ (red) and 1 (blue dash). (d) The differences between these two expressions according to $R_\tau$.</figcaption></figure>

We can see from Figure A3a that when $R_\tau = 1$, the two expressions have no obvious differences, i.e., even when $x_m \neq 0$, the approximation by omitting $x_m$ in the expression is feasible. While when $R_\tau = 10$ (Figure A3b), there are differences appearing around $r_w = 2$. The difference distribution is emphasized in Figure A3c. From it, we can see the difference appears around $r_w = 1.6$ to 4. The maximum of the difference is around $r_w = 2$. However, it can be not negligible if $T_{oo}$ is large. In addition, we observe that the difference around 2 becomes larger as $R_\tau$ increases from 1 to 10. Figure A3d shows the difference according to $R_\tau$ at $r_w = 2$ and 3. From Figure A3d, the differences seem to have a limit which is less than 0.04. Therefore, using the expressions of situation 1 to approximately simulate situation 2, i.e., always $x_m = 0$ is feasible practically, only to have a lower $T_{max}$ and smaller $T_{osc}$ around $r_w = 2$, the error will be less than 0.04 $T_{oo}$.

Why is the largest difference at $r_w = 2$? That is because when $r_w$ is small enough, $x_m$ is 0; while when $r_w$ is very large, the temperature is quickly decreasing, not mentioning
the case when $R_τ$ is large. Therefore, $r_w = 2$ is the middle value that is not small enough to cross the $x_m = 0$ boundary, and not so far away to have a low temperature.

For calculating the analytical expression of the amplitude of temperature, $T_{osc}$ at $r > 0$, due to a non-zero $x_m$ value, the parts in the summation of $T_{max}$ at this point can no longer be eliminated with the parts of $T_{min}$, except in situation 1. We give the analytical expression for situation 1 first, and the general expression will be presented after having given the expressions for $T_{min}$ and $T_{max}$.

When $x_m = 0$, the limit of $T_{osc}$ is described as:

$$T_{osc}(r, N) = T_{max}(r, N) - T_{min}(r, N) = \exp[-(r_w)^2] - \frac{\exp[-\frac{(r_w)^2}{1 + N \cdot R_τ}]}{(1 + N \cdot R_τ)^{\frac{3}{2}}} \underset{N \gg 1/R_τ}{\longrightarrow} \exp[-(r_w)^2] \quad (A12)$$

The amplitude of the oscillations appears to be $T_{00} \cdot \exp[-(r_w)^2]$.

For $T_{min}$ and $T_{max}$, using the trapezoidal rule for approximation as for $r_w \neq 0$, we have $T_{min}$ from Equation (A9):

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

$$\underset{N \gg 1/R_τ}{\longrightarrow} \frac{\exp\left[-\frac{(r_w)^2}{1 + R_τ}\right]}{2(1 + R_τ)^{3/2}} + \frac{\sqrt{\pi}}{R_τ \cdot r_w} \erf\left(\frac{r_w}{\sqrt{1 + R_τ}}\right) \quad (A14)$$

We call $\frac{\exp\left[-\frac{(r_w)^2}{1 + R_τ}\right]}{2(1 + R_τ)^{3/2}}$ part 1, and $\frac{\sqrt{\pi}}{R_τ \cdot r_w} \erf\left(\frac{r_w}{\sqrt{1 + R_τ}}\right)$ part 2.

Furthermore, $T_{max}$ from Equation (A10):

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

$$\underset{N \gg 1/R_τ}{\longrightarrow} \frac{\exp\left[-\frac{(r_w)^2}{1 + x_m \cdot R_τ}\right]}{[1 + x_m \cdot R_τ]^{3/2}} + \frac{\exp\left[-\frac{(r_w)^2}{1 + (1 + x_m) \cdot R_τ}\right]}{2[1 + (1 + x_m) \cdot R_τ]^{3/2}} + \frac{\sqrt{\pi}}{R_τ \cdot r_w} \left\{ \erf\left[\frac{r_w}{\sqrt{1 + (1 + x_m) \cdot R_τ}}\right] \right\} \quad (A15)$$

$\frac{\exp\left[-\frac{(r_w)^2}{1 + x_m \cdot R_τ}\right]}{[1 + x_m \cdot R_τ]^{3/2}}$, $\frac{\exp\left[-\frac{(r_w)^2}{1 + (1 + x_m) \cdot R_τ}\right]}{2[1 + (1 + x_m) \cdot R_τ]^{3/2}}$, $\frac{\sqrt{\pi}}{R_τ \cdot r_w} \left\{ \erf\left[\frac{r_w}{\sqrt{1 + (1 + x_m) \cdot R_τ}}\right] \right\}$ are called part 1, 2, and 3, respectively.

The general expression of $T_{osc}$ (when $N$ tends to infinity or larger than the effective number) is given as Equations (A14) and (A15), and it reads:
$$ T_{osc}(R_{\tau}, r_{w}) = \frac{\exp\left[-\frac{(r_{w})^{2}}{1+x_{m} \cdot R_{\tau}}\right]}{\left[1+x_{m} \cdot R_{\tau}\right]^{\frac{3}{2}}} + \frac{\exp\left[-\frac{(r_{w})^{2}}{1+(1+x_{m}) \cdot R_{\tau}}\right]}{2\left[1+(1+x_{m}) \cdot R_{\tau}\right]^{\frac{3}{2}}} + \frac{\sqrt{\pi}}{R_{\tau} \cdot r_{w}} \left\{ \operatorname{erf}\left[\frac{r_{w}}{\sqrt{\left[1+(1+x_{m}) \cdot R_{\tau}\right]}}\right] \right\} \\ - \frac{\exp\left[-\frac{(r_{w})^{2}}{1+R_{\tau}}\right]}{2(1+R_{\tau})^{\frac{3}{2}}} - \frac{\sqrt{\pi}}{R_{\tau} \cdot r_{w}} \operatorname{erf}\left(\frac{r_{w}}{\sqrt{1+R_{\tau}}}\right) \quad (A16) $$

Part 1 in Equation (A14) and part 2 in Equation (A15) (i.e., $\frac{\exp\left[-\frac{(r_w)^2}{1+(1+x_m) \cdot R_\tau}\right]}{2(1+(1+x_m) \cdot R_\tau)^{3/2}}$ and $\frac{\exp\left[-\frac{(r_w)^2}{1+(1+x_m) \cdot R_\tau}\right]}{2(1+(1+x_m) \cdot R_\tau)^{3/2}}$) are smaller than the other parts by a factor of 10, so they can be approximately omitted to simplify the expressions in practice. Then, we obtain:

$$ T_{osc}(R_{\tau}, r_{w}) \approx \frac{\exp\left[-\frac{(r_{w})^{2}}{1+x_{m} \cdot R_{\tau}}\right]}{\left[1+x_{m} \cdot R_{\tau}\right]^{\frac{3}{2}}} + \frac{\sqrt{\pi}}{R_{\tau} \cdot r_{w}} \left\{ \operatorname{erf}\left[\frac{r_{w}}{\sqrt{\left[1+(1+x_{m}) \cdot R_{\tau}\right]}}\right] \right\} - \frac{\sqrt{\pi}}{R_{\tau} \cdot r_{w}} \operatorname{erf}\left(\frac{r_{w}}{\sqrt{1+R_{\tau}}}\right) \quad (A17) $$

The plots in Figure A4a are $T_{osc}$ according to $r_w$ at $R_\tau = 0.1, 1$, and $10$, accompanied with $\exp[-(r_w)^2]$ (the result of situation 1) for comparison, while the distribution difference between the exact and the approximation is shown in Figure A4b. It is the same as the difference for $T_{max}$ as in Figure A3c. It is also plotted according to $R_\tau$, with $r_w = 2$ and 3, accompanied with $\exp[-(2)^2]$ and $\exp[-(3)^2]$.

<figure><img src="image_19.png" /><figcaption>Figure A4. (a) Spatial distribution of $T_{osc}$ with the real expression ($x_m$ is a function of $R_\tau, r_w$) when $R_\tau = 10, 1$, and $0.1$, and the reduced expression ($x_m = 0$) in the pink dash. (b) The differences between the two expressions according to $r_w$. (c) $T_{osc}$ with the real expression ($x_m$ is a function of ($R_\tau, r_w$), and the reduced expression ($x_m = 0$) according to $R_\tau$ at $r_w = 1, 2, 3$. (d) A zoom of $R_\tau$ from 0 to 10. (e) $T_{osc}$ distribution when $R_\tau = 1000$ (red), comparing with Tosc in situation 1 (blue dash) and the maximum $T_{osc}$ value when $R_\tau \to \infty$ according to $r_w$ (green).</figcaption></figure>

Similarly to $T_{max}$, we note that from Figure A4a, when $R_\tau$ increases, attributed to the non-zero value of $x_m$, compared to situation 1 ($x_m = 0$), the oscillation amplitude has only
a small increase around $r_w = 2$ which is sometimes not negligible. The $R_τ$ dependence is shown in Figure A4c,d, and we note that in Figure A4c, the oscillation amplitude at $r_w = 2$ (green dash) and $r_w = 1$ (cyan) increases when $R_τ$ increases, until it reaches a limit. These limits are calculated to be $\frac{\exp[-\frac{3}{2}]}{[\frac{2}{3}(r_w)^2]^{3/2}}$ when $R_τ \to \infty$ (from Equation (A11), $x_m \cdot R_τ \xrightarrow{R_τ \to \infty} \frac{2}{3}(r_w)^2 - 1$). Figure A4e shows this limit value according to $r_w$ around $r_w = 2$.

Therefore, the range of $T_{osc}$ is given by Equations (A18) and (A19):

$$T_{osc}(R_{\tau}, r_{w}) \xrightarrow{\substack{x_m = 0 \text{ and in particular when } R_{\tau} \to 0 \text{ whatever } r_w \\ \xrightarrow{R_{\tau} \to \infty} \frac{2}{3}(r_w)^2 - 1}} \exp[-(r_w)^2] \quad (A18)$$

$$T_{osc}(R_{\tau}, r_{w}) \xrightarrow{R_{\tau} \to \infty \text{ and } r_{w} \text{ around } 2} \frac{1}{2.44(r_{w})^{3}} \quad (A19)$$

We note that the oscillation amplitude $T_{osc}$ at situation 1 is $\exp[-(r_w)^2]$ which is the minimum, while in situation 2, the amplitude is larger due to the influence of $x_m$, with a maximum value of $\frac{1}{2.44(r_w)^3}$ at a place around $r_w = 2$. Therefore, in practice and for accuracy, consistent with Figure A4e, when $r_w < 1.2$ or $R_τ$ is small, we can apply Equation (A18). Furthermore, when $r_w > 1.2$ or $R_τ$ is large, Equation (A19) can be applied for $T_{osc}$. Even though, by using situation 1 for all the situations, the error will be less than the limits shown in Figure A4c.

## References

1. Dirac, P.A.M. Quantum mechanics of many-electron systems. Proc. R. Soc. London Ser. A Contain. Pap. Math. Phys. Character 1929, 123, 714–733.

2. 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.

3. Bäuerle, D. Laser Processing and Chemistry; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013.

4. Seuthe, T.; Mermillod-Blondin, A.; Grehn, M.; Bonse, J.; Wondraczek, L.; Eberstein, M. Structural relaxation phenomena in silicate glasses modified by irradiation with femtosecond laser pulses. Sci. Rep. 2017, 7, 43815.

5. Chichkov, B.N.; Momma, C.; Nolte, S.; Von Alvensleben, F.; Tünnermann, A. Femtosecond, picosecond and nanosecond laser ablation of solids. Appl. Phys. A 1996, 63, 109–115.

6. Tan, D.; Sharafudeen, K.N.; Yue, Y.; Qiu, J. Femtosecond laser induced phenomena in transparent solid materials: Fundamentals and applications. Prog. Mater. Sci. 2016, 76, 154–228.

7. Gattass, R.R.; Mazur, E. Femtosecond laser micromachining in transparent materials. Nat. Photonics 2008, 2, 219–225.

8. Davis, K.M.; Miura, K.; Sugimoto, N.; Hirao, K. Writing waveguides in glass with a femtosecond laser. Opt. Lett. 1996, 21, 1729–1731.

9. Shimotsuma, Y.; Kazansky, P.G.; Qiu, J.; Hirao, K. Self-organized nanogratings in glass irradiated by ultrashort light pulses. Phys. Rev. Lett. 2003, 91, 247405.

10. Glezer, E.N.; Mazur, E. Ultrafast-laser driven micro-explosions in transparent materials. Appl. Phys. Lett. 1997, 71, 882–884.

11. Fernandez, T.; Sakakura, M.; Eaton, S.; Sotillo, B.; Siegel, J.; Solis, J.; Shimotsuma, Y.; Miura, K. Bespoke photonic devices using ultrafast laser driven ion migration in glasses. Prog. Mater. Sci. 2018, 94, 68–113.

12. Schaffer, C.B.; Brodeur, A.; García, J.F.; Mazur, E. Micromachining bulk glass by use of femtosecond laser pulses with nanojoule energy. Opt. Lett. 2001, 26, 93.

13. Minoshima, K.; Kowalevicz, A.M.; Hartl, I.; Ippen, E.P.; Fujimoto, J.G. Photonic device fabrication in glass by use of nonlinear materials processing with a femtosecond laser oscillator. Opt. Lett. 2001, 26, 1516–1518.

14. Schaffer, C.; García, J.; Mazur, E. Bulk heating of transparent materials using a high-repetition-rate femtosecond laser. Appl. Phys. A Mater. Sci. Process. 2003, 76, 351–354.

15. Eaton, S.M.; Zhang, H.; Herman, P.R.; Yoshino, F.; Shah, L.; Bovatsek, J.; Arai, A.Y. Heat accumulation effects in femtosecond laser-written waveguides with variable repetition rate. Opt. Express 2005, 13, 4708–4716.

16. Eaton, S.M.; Zhang, H.; Ng, M.L.; Li, J.; Chen, W.-J.; Ho, S.; Herman, P.R. Transition from thermal diffusion to heat accumulation in high repetition rate femtosecond laser writing of buried optical waveguides. Opt. Express 2008, 16, 9443–9458.

17. Bérubé, J.-P.; Bernier, M.; Vallée, R. Femtosecond laser-induced refractive index modifications in fluoride glass. Opt. Mater. Express 2013, 3, 598–611.
18. Sakakura, M.; Terazima, M.; Shimotsuma, Y.; Miura, K.; Hirao, K. Heating and rapid cooling of bulk glass after photoexcitation by a focused femtosecond laser pulse. *Opt. Express* **2007**, *15*, 16800–16807.

19. Liu, Y.; Zhu, B.; Wang, L.; Qiu, J.; Dai, Y.; Ma, H. Femtosecond laser induced coordination transformation and migration of ions in sodium borate glasses. *Appl. Phys. Lett.* **2008**, *92*, 121113.

20. Liu, Y.; Shimizu, M.; Zhu, B.; Dai, Y.; Qian, B.; Qiu, J.; Shimotsuma, Y.; Miura, K.; Hirao, K. Micromodification of element distribution in glass using femtosecond laser irradiation. *Opt. Lett.* **2009**, *34*, 136–138.

21. Shimizu, M.; Sakakura, M.; Kanehira, S.; Nishi, M.; Shimotsuma, Y.; Hirao, K.; Miura, K. Formation mechanism of element distribution in glass under femtosecond laser irradiation. *Opt. Lett.* **2011**, *36*, 2161–2163.

22. Zhu, B.; Dai, Y.; Ma, H.; Zhang, S.; Lin, G.; Qiu, J. Femtosecond laser induced space-selective precipitation of nonlinear optical crystals in rare-earth-doped glasses. *Opt. Express* **2007**, *15*, 6069–6074.

23. Dai, Y.; Zhu, B.; Qiu, J.; Ma, H.; Lu, B.; Cao, S.; Yu, B. Direct writing three-dimensional Ba2TiSi2O8 crystalline pattern in glass with ultrashort pulse laser. *Appl. Phys. Lett.* **2007**, *90*, 181109.

24. Zhang, B.; Tan, D.; Liu, X.; Tong, L.; Kazansky, P.G.; Qiu, J. Self-organized periodic crystallization in unconventional glass created by an ultrafast laser for optical attenuation in the broadband near-infrared region. *Adv. Opt. Mater.* **2019**, *7*, 1900593.

25. Miyamoto, I.; Cvecek, K.; Schmidt, M. Evaluation of nonlinear absorptivity in internal modification of bulk glass by ultrashort laser pulses. *Opt. Express* **2011**, *19*, 10714–10727.

26. Miyamoto, I.; Okamoto, Y.; Tanabe, R.; Ito, Y. Characterization of plasma in microwelding of glass using ultrashort laser pulse at high pulse repetition rates. *Phys. Procedia* **2014**, *56*, 973–982.

27. Lax, M. Temperature rise induced by a laser beam. *J. Appl. Phys.* **1977**, *48*, 3919–3924.

28. Sanders, D.J. Temperature distributions produced by scanning Gaussian laser beams. *Appl. Opt.* **1984**, *23*, 30–35.

29. Haba, B.; Hussey, B.W.; Gupta, A. Temperature distribution during heating using a high repetition rate pulsed laser. *J. Appl. Phys.* **1991**, *69*, 2871–2876.

30. Zhang, H.; Eaton, S.M.; Li, J.; Herman, P.R. Heat accumulation during high repetition rate ultrafast laser interaction: Waveguide writing in borosilicate glass. *J. Phys. Conf. Ser.* **2007**, *59*, 682–686.

31. 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.

32. Beresna, M.; Gertus, T.; Tomašiūnas, R.; Misawa, H.; Juodkazis, S. Three-dimensional modeling of the heat-affected zone in laser machining applications. *Laser Chem.* **2008**, *2008*, 976205.

33. Shimizu, M.; Sakakura, M.; Ohnishi, M.; Yamaji, M.; Shimotsuma, Y.; Hirao, K.; Miura, K. Three-dimensional temperature distribution and modification mechanism in glass during ultrafast laser irradiation at high repetition rates. *Opt. Express* **2012**, *20*, 934–940.

34. Rahaman, A.; Kar, A.; Yu, X. Thermal effects of ultrafast laser interaction with polypropylene. *Opt. Express* **2019**, *27*, 5764–5783.

35. Rahaman, A.; Du, X.; Zhou, B.; Cheng, H.; Kar, A.; Yu, X. Absorption and temperature distribution during ultrafast laser microcutting of polymeric materials. *J. Laser Appl.* **2020**, *32*, 022044.

36. 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.

37. 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.

38. Neuville, D.R. *From Glass to Crystal. Nucleation, Growth and Phase Separation: From Research to Applications*; chemistry/materials; EDP Sciences: 2017; 640p.

39. 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.

40. 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.

41. Sudrie, L. Propagation Non-Linéaire des Impulsions Laser Femtosecondes Dans la Silice. Université de Paris Sud XI Orsay. 2002. Available online: http://www.ensta.fr/ilm/Archives/Theses pdf/L_Sudrie (access on 1 January 2012).

42. Bressel, L.; de Ligny, D.; Sonneville, C.; Martinez, V.; Mizeikis, V.; Buividas, R.; Juodkazis, S. Femtosecond laser induced density changes in GeO₂ and SiO₂ glasses: Fictive temperature effect [Invited]. *Opt. Mater. Express* **2011**, *1*, 605–613.

43. Lancry, M.; Poumellec, B.; Chahid-Erraji, A.; Beresna, M.; Kazansky, P.G. Dependence of the femtosecond laser refractive index change thresholds on the chemical composition of doped-silica glasses. *Opt. Mater. Express* **2011**, *1*, 711–723.

44. Guerfa-Tamezait, D. Modifications d'une Couche Mince sur un Substrat, Induites par L'irradiation Mono-Impulsionnelle UV ns. Ph.D. Thesis, Université Paris-Saclay, Bures-sur-Yvette, France, 2022.

45. Yao, H.; Xie, Q.; Cavillon, M.; Neuville, D.R.; Pugliese, D.; Janner, D.; Dai, Y.; Poumellec, B.; Lancry, M. Volume nanogratings inscribed by ultrafast IR laser in alumino-borosilicate glasses. *Opt. Express* **2023**, *31*, 15449–15460.

46. 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.

47. Cavillon, M.; Wang, Y.; Poumellec, B.; Brisset, F.; Lancry, M. Erasure of nanopores in silicate glasses induced by femtosecond laser irradiation in the Type II regime. *Appl. Phys. A* **2020**, *126*, 876.
48. Que, R.; Houel-Renault, L.; Temagoult, M.; Herrero, C.; Lancry, M.; Poumellec, B. Space-selective creation of photonics functions in a new organic material: Femtosecond laser direct writing in Zeonex glass of refractive index change and photoluminescence. *Opt. Mater.* **2022**, *133*, 112651.

49. He, X.; Liu, Q.; Lancry, M.; Brisset, F.; Poumellec, B. Space-Selective Control of Functional Crystals by Femtosecond Laser: A Comparison between SrO-TiO₂-SiO₂ and Li₂O-Nb₂O₅-SiO₂ Glasses. *Crystals* **2020**, *10*, 979.

50. Delullier, P.; Calvez, L.; Druart, G.; De La Barrière, F.; Humbert, C.; Poumellec, B.; Lancry, M. Photosensitivity of Infrared Glasses under Femtosecond Laser Direct Writing for mid-IR Applications. *Appl. Sci.* **2022**, *12*, 8813.

51. Simmons, J.A. Thermal Conductivity of Glycine. *Nature* **1967**, *216*, 1302.

52. Wang, X. Convergence-Divergence of p-Series. *Coll. Math. J.* **2002**, *33*, 314–316.

53. Weideman, J.A.C. Numerical integration of periodic functions: A few examples. *Am. Math. Mon.* **2002**, *109*, 21–36.

54. Atkinson, K.E. *An Introduction to Numerical Analysis*; John Wiley & Sons: Hoboken, NJ, USA, 2008.

**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.