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

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 T00 and a factor Rτ representing the ratio of the pulse period τp to the diffusion time τd. We show that temperature, oscillating between Tmax and Tmin, 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.


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 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 −11 -10 −14 s) is partially injected into a small focal volume of transparent dielectric solids.This intense laser pulse with high irradiance (>10 13 Wcm −2 ) in the focal region stimulates a series of complex dynamic processes, such as multiphoton ionization, tunneling ionization, inverse bremsstrahlung absorption, and avalanche ionization within an ultrashort time scale [2].Such interactions lead to highdensity 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][7][8][9][10][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][12][13][14][15][16][17][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][20][21][22][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 avalanche 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][15][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/α).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 using Duhamel'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.

CW static cylindrical
Gaussian(r) Beer-Lambert(z) analytical Lax [27] pulsed scanning three axes Gaussian(x,y) Beer-Lambert(z) analytical Sanders [28] pulsed scanning three axes Gaussian(x,y) Beer-Lambert(z) analytical Haba [29] pulsed static spherical Gaussian(r) finite difference Eaton [15] or Zhang [30] pulsed quasi-static cylindrical Gaussian(r) Beer-Lambert(z) analytical one pulse Sakakura [18] CW scanning three axes uniform deposition in parallelepiped volume analytical Miyamoto [31] pulsed static cylindrical Gaussian(r,z) analytical Miyamoto [25] pulsed static cylindrical Gaussian(r) Gaussian(z) numerical Shimizu [33] pulsed scanning cylindrical Gaussian(r) surface absorption analytical Rahaman [34,35] pulsed quasi-static spherical Gaussian(r) analytical this work 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.

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 ∇T where → q is the heat flow (energy per unit area and time).Fourier considers it to be linearly dependent on the temperature gradient.κ 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 •s•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, terms − sink terms.The source term is the laser energy density deposited per unit of time (i.e., absorbed laser power), written symbolically as δρ Q δ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, dt , ρ and C p are the density and specific heat capacity, respectively.
r∂r .Considering a beam moving not too fast, convection can be neglected (i.e., the time derivative of the spatial coordinate), dT dt is thus written as ∂T ∂t .Therefore, we obtain the following equation: 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 •∆T = 0, and Equation (1) becomes: Micromachines 2024, 15, 196 5 of 32 Assuming a Gaussian shape of , 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, w 2 with: After pulse energy deposition, diffusion begins to operate, and Equation (1) becomes: Using the initial and boundary conditions on solutions of Equations ( 2) and ( 4), we obtain: 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 τ p = 1/RR, where RR is the pulse repetition rate: With τ d = w 2 /4D th , r w = r w , Equation (6) reads: We note that the variables involved in Equation (7) are the ratio between the period of the pulses τ p and the diffusion times τ d , while the other laser and material parameters are involved in the amplitude T 00 .Therefore, we introduce the parameter R τ : Therefore, Equation (7) becomes: 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 ≫ 1/R τ .We will show how the temperature changes with the number of pulses according to heat accumulation (hence R τ ), 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.

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 ̸ = 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 .

At the Center
At the center, r w = 0, and thus, from Equation (9): Calling N t = t τ p .The temperature evolutions over the generalized pulse number N t for several R τ are shown in Figure 1.From Figure 1 we observe that: -(0,  ) oscillates between a minimum ( ) and a maximum ( ) in each period between two pulses; - The oscillation amplitude ( ) seems to be the same, whatever  ; - seems to reach a steady state as  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  but larger for small  values.For a larger  , the temporal overlapping of 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 τ ; -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 τ but larger for small R τ values.For a larger R τ , the temporal overlapping of temperature increase contributions from consecutive pulses is weaker, whereas it increases (heat accumulation) when R τ is smaller.

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 − ε), where ε 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)τ p .T min is at the end of the pulse period, just before the (N + 1)th pulse arrival.Using Equation ( 10), we have: (11) T min will be at t = N•τ p − ε, thus not containing the temperature contribution induced by the (N + 1)th pulse, so: and When N ≫ 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 ≫ 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.The Effective Number of Pulses for Reaching the Limit of Tosc( ) When will the temperature in the material reach a stable oscillation?In practice, we can calculate a real number of pulses ( ) to closely reach the oscillation limit.(In notation,  , where sso means steady state of oscillation, 0 means the situation when rw = 0).
Consider when  (0,  ) = (1 − ) •  (0, ∞) , the limit is practically reached, where ε is a small quantity, i.e., a few % (based on the actual situation).Thus, it is: Then, we have: The plot of  is shown in Figure 3, note that the actual number of pulses is the integer part above 1 (as ε should be smaller, bounded by  = ).Some specific parameters are given below for visualizing this value.With ε = 3%, when  = 1 (conceivable combinations of material parameters and laser RR),  = 9.36.So, after 10 pulses, the amplitude of the oscillating temperature reaches 0.97 T00.When  is large, e.g.,  = 10,  = 0.94, so only one pulse rules the oscillation amplitude, and we can understand that pulse contributions are separated.When  is smaller,  increases rapidly, e.g.,  = 0.1 , i.e., 1/ = 10 ,  = 94 .Beyond  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 ( larger than which value) the temperature oscillation is constant since the first pulse.The Effective Number of Pulses for Reaching the Limit of T osc (N 0 sso ) When will the temperature in the material reach a stable oscillation?In practice, we can calculate a real number of pulses (N 0 sso ) to closely reach the oscillation limit.(In notation, N 0 sso , where sso means steady state of oscillation, 0 means the situation when r w = 0).Consider when T osc 0, N 0 sso = (1 − ε)•T osc (0, ∞), the limit is practically reached, where ε is a small quantity, i.e., a few % (based on the actual situation).Thus, it is: Then, we have: The plot of N 0 sso is shown in Figure 3, note that the actual number of pulses is the integer part above 1 (as ε should be smaller, bounded by ε = 1 ).Some specific parameters are given below for visualizing this value.With ε = 3%, when R τ = 1 (conceivable combinations of material parameters and laser RR), N 0 sso = 9.36.So, after 10 pulses, the amplitude of the oscillating temperature reaches 0.97 T 00 .When R τ is large, e.g., R τ = 10, N 0 sso = 0.94, so only one pulse rules the oscillation amplitude, and we can understand that pulse contributions are separated.When R τ is smaller, N o sso increases rapidly, e.g., R τ = 0.1, i.e., 1/R τ = 10, N o sso = 94.Beyond N o sso 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 τ larger than which value) the temperature oscillation is constant since the first pulse.
To conclude on this point, Equations ( 13) and ( 14) provides that the temperature oscillation amplitude is T 00 after N sso 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.To conclude on this point, Equations ( 13) and ( 14) provides that the temperature oscillation amplitude is T00 after Nsso pulses and at an oscillatory steady state.Equation ( 16) provides the effective number of pulses for reaching it.It takes more time when  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.

Tmin and Tmax
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 Tmax and Tmin Tmin: The analytical expression of the minimum temperature is transformed from the sum expression Equation ( 12), with details found in Appendix C. Therefore, we obtain: This expression shows the increase of  according to  and  .It is plotted in Figure 4.

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: This expression shows the increase of T min according to N and R τ .It is plotted in Figure 4.The final limit of  , i.e., when  ≫ 1/ is given below: The same method has been applied to obtain the  limit, and the detail can be also found in Appendix C. We have thus: The final limit of T min , i.e., when N ≫ 1/R τ is given below: 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: When N ≫ 1/R τ , the T max limit is: From these expressions, we see that the difference between T max (0, ∞) and T min (0, ∞) 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 2 R τ .Reintroducing here exceptionally T 00 (Equation ( 3)), we obtain: 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 (κ), 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 0 ssmin , N 0 ssmax ) 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 sso , the first approximation of N ssmin is obtained by solving the following assertion, This cubic equation has three roots, where the physical one is . Therefore, N ssmax : With the same method, we obtain: The N ss for reaching closely the steady state (with ε departure).T osc (0,∞), T min (0,∞), and T max (0,∞) are plotted in Figure 5 according to R τ (with ε = 3%).
From Figure 5, we can see that with 1/R τ 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 τ 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 sso and N ssmax , the value converges to 0 when R τ 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 0 ssmax 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 ε = 0.03 for organic materials and ε = 0.06 for inorganic materials): Nssmax: With the same method, we obtain: The Nss for reaching closely the steady state (with ε departure).Tosc(0,∞), Tmin(0,∞), and Tmax(0,∞) are plotted in Figure 5 according to  (with  = 3%).From Figure 5, we can see that with 1/ increasing, the effective pulse numbers for reaching the steady state increases whether for Tosc, Tmin, or Tmax.For practical use, it is better to define one Nss for calculation.When  is large, since the value of Tmin is almost 0 (no pulse superimposition), it is therefore not meaningful to take it into consideration.For Nsso and Nssmax, the value converges to 0 when  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 Tmax.Therefore,  is the suitable and practical number of pulses needed for reaching the steady state.Some examples for the Nssmax value are shown in Table 2 (with ε = 0.03 for organic materials and ε = 0.06 for inorganic materials):  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 0 ss is (considering the effective number to reach the T max limit): Figure 6 show the plots versus 1 according to three different diffusion times: 0.28 µs (silica, glycine), 1.63 µs for nifedipine, and 4.9 µs for sucrose.
Figure 6.The time in µs to reach the steady state according to 1/ 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.
For small enough values of  , the time reaches the value  / , i.e., 1111 for ε = 3%.Note that for inorganic glass and glycine crystal cited above with  = 0.28 µs, this 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 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 τ .

Time
Behavior out of the Center (r w = r/w ̸ = 0) When r w ̸ = 0, we come back to the expression Equation (10): 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).
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.

T osc , T min , and T max
The Limit of T max and T min (when To calculate the oscillation amplitude T r 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 Nth pulse period.T min is still considered at the end of the Nth one, i.e., when t = N•τ p before the absorption of the (N + 1)th pulse, so:  We observe the following differences according to radius rw = 0, 1, 2: - The amplitude of oscillation is less than 1 (in the unit of T00) for increasing radius; - The maximum temperature during a period is still at the beginning of the pulse deposition for  = 1 with these three   , while at  = 2 the maximum temperature is no more at the beginning.That is because there is time for heat to diffuse from the However, since T max in some situations can be in the middle of the pulse period, we set x m , 0 ≤ x m ≤ 1 to define the place where the T max is.Therefore, T max is at the time t = (N − 1 + x m )τ p : 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 ≫ 1/R τ , x m is shown below, and the results of the cumbersome calculation details can be found in Appendix D: 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 < 3  2 + 2 R τ whatever R τ or R τ small enough (less than 2 r w 2 −1.5 when r w 2 > 1.5).
Situation 2: x m ̸ = 0, i.e., r w > 3 2 + 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 .
Micromachines 2024, 15, 196 14 of 32 (1) For situation 1, when x m = 0, the limit of T osc is described as: The amplitude of the temperature oscillations reaches T 00 •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): 2(1+R τ ) 3/2 is called part 1, and (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 τ is large.
T min is the same as situation 1.For T max , using the trapezoidal rule for approximation as for r w ̸ = 0, we have T max from Equation (26): 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: 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.
It is worth noticing (see Appendix D, Figure A4a,b) that when R τ 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 1 2.44(r w ) 3 for a large R τ (details are shown in Appendix D Figure A4c,d by plotting T osc according to r w and R τ ).Therefore, the range of T osc is given by Equations ( 33) and (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 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 τ 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 τ .
Micromachines 2024, 15,196 17 of 35 It is worth noticing (see Appendix D, Figure A4a,b) that when  increases, Tosc exhibits a small departure from the exact value at around rw = 2, attributed to the existence of a non-zero xm.This departure, if it is generally not negligible, is nevertheless bounded.It is calculated to be .( ) for a large  (details are shown in Appendix D Figure A4c,d by plotting Tosc according to rw and  ).Therefore, the range of Tosc is given by Equations ( 33) and (34): We note that the oscillation amplitude Tosc at situation 1 is  −( ) which is the minimum, while in situation 2, the amplitude is larger due to the influence of xm, with a maximum value of .( ) at the place around rw = 2.By now, with these temperature expressions, we obtain the spatial distribution of the minimum and maximum temperature for a given  at steady state, shown in Figure 8.The temperature is oscillating between these two temperature profiles, and note that at rw = 0, the difference is always 1 regardless of  .
We have now all the information for plotting the T distribution with any Rτ value.From Figure 8, we can see that when  is small (large frequency or small diffusivity), Tmin and Tmax have no large relative difference compared to their average values because the oscillation amplitude is always limited to  −( ) whereas the Tmean amplitude is converging to 2/ (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  , the shape of the curve converges to the erf(rw) curve and decreases much slower than a Gaussian one.
For large values of  (small frequency or large diffusivity), the oscillations are relatively large as the pulses are separated and thus Tmin appears to have small values.It is negligible (<3%) when  > 16 or < 6% for  > 7 ).This limits the domain of heat accumulation.The calculation shows that the shape of Tmax is also converging with increasing  to the shape of the beam energy distribution (Gaussian, here  −( ) )  29)) and T max (red, by Equations ( 27) and ( 30)) according to the relative radius r w when (a) R τ = 0.1, (b) R τ = 1, and (c) R τ = 10.
We have now all the information for plotting the T distribution with any R τ value.From Figure 8, we can see that when R τ 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 τ (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 τ , the shape of the curve converges to the erf (r w ) curve and decreases much slower than a Gaussian one.
For large values of R τ (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 τ > 16 or <6% for R τ > 7).This limits the domain of heat accumulation.The calculation shows that the shape of T max is also converging with increasing R τ to the shape of the beam energy distribution (Gaussian, here exp −(r w ) 2 ) independent of R τ .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 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 ε being a small quantity and , we have < ε.Therefore, the N r sso/ssmin/ssmax solutions are shown below: This expression does not have an analytic root without using the tabulated function W, i.e., the Lambert W function (defined as ωe ω = z, ω = W(z)).In practice, since X 2 ≪ 1, by approximation, it becomes: For N r ssmin and N r ssmax , with the approximation of erf(X The behavior of N r ss 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.
The behavior of  according to  for  = 0 has already been shown in Figure 5, with an overall increase.We have plotted Nsso and Nssmax, according to  , for  = 10 as an example, as shown in Figure 9a, and the plot of the related time for reaching the steady state (using Nssmax and with diffusion time 0.28 µs) in Figure 9b.From Figure 9, we can see that as  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.

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  (, ) = (, ) , and this gives: 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.

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 T(r, N) = 1 τ p pulse period at N T(r, t)dt, and this gives: N.B. due to software problem, the average temperature is sometimes quoted as T and sometimes T mean .They have the same meaning.
The two summations can be permuted as they do not operate on the same variables and are independent.We obtain: This result is obtained without approximation.Then, when N ≫ 1/R τ : We note that here, the steady state temperature at the center will reach: It is the same expression as for T max (0, ∞) or T min (0, ∞) for small R τ values.We can note in Figure 10 that T max (0, ∞) and T min (0, ∞) approach T(0, ∞) when R τ is decreasing.T max (0, ∞) goes to 1 and T min (0, ∞) goes to 0 when R τ is large.From Figure 10, we can also find the heat accumulation bound already defined in Figure 5 (with ε = 0,03).It corresponds to T min (0, ∞) = 0.03 and R τ = 12.On the other hand, when T min (0, ∞) departs from T max (0, ∞) by less than approximately 10%, we can admit that the average T is applicable, i.e., for R τ smaller than 0.17 (purple circle).In this case, we can apply the simple expression Equation (43).  .The plots of T min (0, ∞) (red), T max (0, ∞) (blue dash), and T(0, ∞) (T mean , green dash) according to 1/R τ .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).
The Effective Number of Pulses for Reaching the Limit T(0, ∞) (N 0 ssm ) With the same definition as above, the number of pulses to reach T(0, ∞), i.e., |T(0,∞)−T(0,N)| T(0,∞) < ε, N 0 ssm is obtained: N 0 ssm are compared to N 0 ssmax and N 0 sso in Figure 11.The steady state of the mean temperature is reached as T max .
The Effective Number of Pulses for Reaching the Limit  (0, ∞) ( ) With the same definition as above, the number of pulses to reach  (0, ∞) , i.e.,

| ( , ) ( , )| ( , )
< ,  is obtained: are compared to  and  in Figure 11.The steady state of the mean temperature is reached as Tmax.
Figure 11.The effective number of pulses to reach the limit of Tosc (red), Tmin (blue dash), Tmax (green dash), and Tmean (pink dash) according to Rτ from 0 to 5.
Therefore, the time for reaching the steady state here is not Rτ dependent:  − 1 ≈ = 1111 when  = 0.03.It is the value of the maximum tss for reaching a steady state when  → 0 (Figure 6).With these analytical expressions of temperature at the steady state at the center of the focus ( ,  (0, ∞),  (0, ∞),  (0, ∞)) , and the needed number of pulses ( ,  ,  ) , we have a clear view of how parameter  influences the thermal . 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 τ from 0 to 5.
Therefore, the time for reaching the steady state here is not R τ dependent: τ D It is the value of the maximum t ss for reaching a steady state when R τ → 0 (Figure 6).With these analytical expressions of temperature at the steady state at the center of the focus (T osc , T min (0, ∞), T max (0, ∞), T(0, ∞)), and the needed number of pulses (N 0 sso , N 0 ssmax , N 0 ssm ), we have a clear view of how parameter R τ influences the thermal situation at the focus center.The problem is now to extend these results to any place out of the center.

Temperature out of the Focus Center (T(r, N))
With the definition T(r, N) = 1 τ p pulse period at N T(r, t)dt, we have the average temperature in a period as: This limit when N ≫ 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 er f 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 2 R τ .It is inversely proportional to R τ whatever R τ value.
This limit when  ≫ 1/  is shown in Figure 12 and compared to the Gaussian shape of Tmax when  is large and when  is small.When Tmax(r) is Gaussian for the first case, Tmax(r) has the same shape that Tmean has for the second case.As the  function tends to 1 (already for  > 2), the function tends to be hyperbolic and thus decreases much slower than a Gaussian one (see Figure 12).The amplitude is .It is inversely proportional to  whatever  value.Consistently with the previously calculated Tmax and Tmin, the Tmean curve width is equal to the beam Gaussian for large  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 Tmean (𝑁 )
The effective number of pulses  with the approximation ( •  ) ≈ √  •  as (, ) = √ • < 1, is solved to be: 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.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 r ssm ) The effective number of pulses N r ssm with the approximation erf(X , is solved to be: 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.

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)).
√ π Rτ•r w erf(r w ) = T(r w , ∞) = 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.

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 τ , defined as the ratio of the pulse period τ p to the diffusion time τ 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.27), Figure A2).T osc (r,N) 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).(1) x m = 0, for r w < 3 2 + 2 R τ or R τ small enough (less than 2 r w 2 −1.5 when r w 2 > 1.5).In this situation, x m can be omitted from the expression.
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).We can see from Figure A3a that when  = 1, the two expressions have no obvious differences, i.e., even when  ≠ 0, the approximation by omitting  in the expression is feasible.While when  = 10 (Figure A3b), there are differences appearing around  = 2.The difference distribution is emphasized in Figure A3c.From it, we can see the difference appears around  = 1.6 to 4. The maximum of the difference is around  = 2 .However, it can be not negligible if T00 is large.In addition, we observe that the difference around 2 becomes larger as  increases from 1 to 10. Figure A3d shows the difference according to  at  = 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 xm = 0 is feasible practically, only to have a lower Tmax and smaller Tosc around  = 2, the error will be less than 0.04 T00.
Why is the largest difference at  = 2? That is because when  is small enough, xm We can see from Figure A3a that when R τ = 1, the two expressions have no obvious differences, i.e., even when x m ̸ = 0, the approximation by omitting x m in the expression is feasible.While when R τ = 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 00 is large.In addition, we observe that the difference around 2 becomes larger as R τ increases from 1 to 10. Figure A3d shows the difference according to

Figure 5 .
Figure 5.The effective number to reach the limit of Tosc, Tmin, and Tmax according to  from 0.1 to 100 for ε = 0.03.The boundary point is at Rτ = 15.7.

Figure 5 .
Figure 5.The effective number to reach the limit of T osc , T min , and T max according to R τ from 0.1 to 100 for ε = 0.03.The boundary point is at R τ = 15.7.

Figure 7
Figure7shows the temperature evolution based on the above expression over time at two relative distances  = 1 (Figure7a,b) and  = 2 (Figure7c,d).

Figure 9 .
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.

Figure 10
Figure10.The plots of T min (0, ∞) (red), T max (0, ∞) (blue dash), and T(0, ∞) (T mean , green dash) according to 1/R τ .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).

Figure 12 .
Figure12.Plot of normalized T max (r w , R τ ) and T(r w , R τ ) at the steady state for R τ = 0.1 and 100 for T max and 1 for T mean .

Micromachines 2024 ,
15, 196  30 of 35    the one without  .FigureA3a,b displays the spatial distribution of Tmax with the  of the value depending on rw (red) and the one with  = 0 for whatever rw (blue dash).

Figure A3 .
Figure A3.(a,b) Spatial distribution of Tmax with the exact expression ( is a function of  ,  ) and simplified expression (  = 0) .(a)  = 1 , (b)  = 10.(c) The differences between the two expressions of  = 10 (red) and 1 (blue dash).(d) The differences between these two expressions according to  .

Figure A3 .
Figure A3.(a,b) Spatial distribution of T max with the exact expression (x m is a function of R τ , r w ) and simplified expression (x m = 0).(a) R τ = 1, (b) R τ = 10.(c) The differences between the two expressions of R τ = 1 0 (red) and 1 (blue dash).(d) The differences between these two expressions according to R τ .

Table 1 .
State of the art of the thermal simulation of laser-matter interactions.

Table 2 .
Pulse number needed (N ss ) for reaching the steady state in materials, using Equation (23).

Table 3 .
Analytical expressions of final temperature and the effective number for reaching steady state.*: condition of r w , R τ (Equation (

Table A2 .
Thermal physico-chemical data of some materials and processing parameters.