Laser-Induced Thermal Stresses in Dense and Porous Silicon Dioxide Films

: The laser-induced thermal stresses in silicon dioxide ﬁlms are calculated using molecular dynamics simulations. The absorption of the laser energy is simulated by the linear temperature growth from room temperature to 1300 K in a time equal to the laser pulse duration. The maximum values of stresses for picosecond pulses are approximately twice as high as for nanosecond pulses. The stresses in highly porous glancing angle deposited ﬁlms are approximately two times lower than in dense ﬁlms. Stress waves caused by picosecond pulses are observed in dense ﬁlms. An increase in the heating temperature to 1700 K leads to an increase in the absolute stress values for picosecond pulses, and a decrease for nanosecond pulses.


Introduction
The investigation of the possible causes of laser-induced damage to optical films is an important task in many applications of high-power lasers [1]. To date, there is still no complete understanding of the effects resulting in laser-induced damage to transparent dielectric coatings. Thermoelastic breakdown, caused by thermal expansion and the associated stresses, is considered as one of such effects [2].
Laser-induced thermal stresses have been investigated in different coatings, both experimentally and theoretically. The propagation of laser-induced elastic waves in fused silica and a single Si crystal was experimentally studied in [3]. It was found that the maximum stress is observed at the interface between the film and the substrate. The interface stress increases with increasing film thickness and decreases with increasing substrate thickness. The surface acoustic waves arising from the laser energy absorption in the thin film-substrate system were investigated both experimentally and theoretically in [4]. Wave-induced stresses were calculated using a finite element model. The dependence of the interfacial stresses on the ratio of the wavelength-to-film thickness was investigated. The correlation between the polarization state of the femtosecond laser radiation and the anisotropy of stresses arising in films during laser writing was experimentally investigated in [5]. Based on the volume expansion model, the maximum stress in fused silica was calculated as 2 GPa [5].
In Ref. [6], the finite element model was developed to simulate laser-induced stress waves in dielectric multilayer coatings. The model allows the study of the dependence of the degradation of material structural integrity, caused by stress waves, on the coating's geometric and material parameters. It was found that closer impedance of the materials of the layers decreases interface stresses by 20%-30%. Also, increasing the thickness of a single layer reduces the maximum interfacial stress. In Ref. [7], the heat conduction and elastic equations were solved in the case of a single Gaussian laser pulse absorbed by the MgF 2 /ZnS double layer. The distributions of the temperature and stresses over the film thickness were calculated. The factors influencing these distributions were identified and discussed. In Ref. [8], the 2D finite element model was used to simulate the laserinduced thermal stresses in different optical thin-film materials. It was found that for all the materials, the induced stress is compressive. The smallest absolute value of stresses was obtained for the MgF 2 film.
Compared with the finite element models, the atomistic simulation methods allow the study of the processes of laser energy absorption at the microscopic level. This is especially important when simulating the fast processes with times in the order of picoseconds or less. The stresses in nickel films, caused by femtosecond laser pulses, were studied using molecular dynamics (MD) simulations in Ref. [9]. It was found that ultrafast heating results in the formation of a 100 nm thick strongly compressed layer below the film surface. This layer remains compressed for less than 50 ps. The laser pulses with low intensity form only a single elastic wave.
The MD method was applied in [10] to study the stresses in growing silicon dioxide films. The stresses were calculated for dense films and for glancing angle deposited (GLAD) films. It was found that the absolute values of stresses in GLAD films are several times less than in dense films. Also, for GLAD films, stress anisotropy in the substrate plane was found.
In this work, the earlier developed MD-based method [11][12][13] was modified and applied to simulate laser-induced thermal stresses in dense and GLAD silicon dioxide films. Laser pulses of different durations were considered. The dependence of the stress values on the thin films structure and on the heating temperature is discussed.

Method
The simulation scheme is shown in Figure 1. The atomistic cluster represents the part of the film and substrate exposed to laser irradiation (bottom of Figure 1). To calculate the laser-induced thermal stresses, the following steps are performed:

•
Simulation of the deposition of SiO 2 thin films as described in the section below. • Equilibration of the deposited structure in the NPT (constant number of particles, pressure and temperature) ensemble. The asymmetric barostat is applied to ensure independent variation in the cluster dimensions during equilibration. This step is necessary to eliminate stresses that appear during thin-film deposition; • Modeling the absorption of laser radiation in the film by increasing the cluster temperature to the value T f . The heating is assumed to be uniform over the film volume since the cluster thickness of 40 nm is insufficient to consider the distribution of the laser radiation intensity along the z-axis inside the film. The calculations of stresses are averaged over the film volume during the film heating; • In the case of the fast laser pulses with durations of τ = 1 ps and τ = 10 ps, the stresses are also calculated during the relaxation time τ r after the end of the film heating. The temperature during this period is kept equal to T f .
The laser-induced thermal stress is studied using the atomistic clusters of the silicon dioxide films obtained in our previous work [10,12,13]. These clusters were deposited using the MD-based step-by-step scheme [11]. Oxygen and silicon atoms are randomly inserted into the upper part of the modeling box in a stoichiometric proportion of 1:2. The initial velocities are directed to the substrate. The values of the velocities correspond to the energies E(Si) = 10 eV and E(Si) = 0.1 eV for high-and low-energy deposition processes, respectively, and E(O) = 0.1 eV in both cases. Then the MD simulation starts in the NVT (constant number of particles, volume and temperature) ensemble. The periodic boundary conditions are applied. The duration of one injection step is equal to 10 ps, and one time step is equal to 0.5 fs. The horizontal substrate dimensions are equal to 10 nm and 30 nm. The substrate thickness in a vertical direction is equal to 6 nm. The thicknesses of the deposited films are about 40 nm. The temperature of the simulation cluster is kept equal to 300 K using the Berendsen thermostat [14]. The energy of the interatomic interactions is calculated using the empirical pairwise DESIL (DEposition of SILica) force field [11]. The deposition is performed at different values of the deposition angle α (Figure 1). The stress tensor σ is calculated as follows [10,15,16]: where V is the modeling box volume, K is the kinetic energy tensor, Ξ is the virial tensor.
where ⊗ is the direct product of two vectors, m i and → υ i are the atomic masses and velocities, → r ij is the radius vector directed from the i-th atom to the j-th atom, → f ij is the force acting from the i-th atom to the j-th atom, and N is the number of atoms.
The simulation is carried out using equipment at the shared research facilities of highperformance computing resources at Lomonosov Moscow State University [18].

Results and Discussion
In this work, we investigated the stresses induced by heating to the temperature T f below the softening point of fused silica, 2.0 × 10 3 K [19]. The results of the stresses calculation for low-energy-deposited SiO 2 films are presented in Figure 2. Hereinafter, σ = (σ xx + σ yy )/2, where σ xx and σ yy are the diagonal components of the stress tensor, Equation (1), and the axes directions are shown in Figure 1. The fast oscillation of σ is due to a limited number of particles in the simulation clusters which leads to noticeable thermodynamic fluctuations of the calculated stresses. (see Appendix A).
The temperature of the films increases linearly from an initial temperature of 300 K to T f = 1300 K. We assume that the temperature of 1300 K corresponds to the value of the laser intensity, which is not much lower than the laser-induced damage threshold for SiO 2 films,~4 J/cm 2 [20].
The heating time is denoted as τ and is taken to be 1 ps, 10 ps, and 1 ns, which corresponds to the picosecond and nanosecond laser pulses. The plots in Figure 2 show the stress values averaged over segments of the MD trajectories, with a duration of 0.1 ps.
After the end of heating, the film temperature begins to decrease. To estimate the characteristic time of this decrease, the heat conduction equation is solved for the model task on the condition that there is uniform absorption in the film, with realistic thicknesses of the film and the substrate. It is found that the temperature of the upper edge of the film, deposited on a 1 mm thick substrate, decreases by less than one percent of the maximum value during several microseconds. Increasing the substrate thickness increases the temperature-reducing time. In all the cases, this time is much longer than the times of the MD simulations. For this reason, the temperature is kept equal to T f during the simulation, after the end of the heating.
In all the cases, the sign of the stresses is negative, which corresponds to the compressive stress. In the case of τ = 1 ps, the stress grows for approximately 1-2 ps after the end of heating (see the enlarged plots inside the red rectangle at the top of Figure 2). The stress values change insignificantly when the time of heating increases from τ = 1 ps to τ = 10 ps (compare Figure 2a,b). However, when the time of heating increases to 1 ns, the stress reduces approximately twice (Figure 2c). This could be attributed to the fact that the structure has more time for relaxation under nanosecond heating.
The maximum stress values significantly decrease with an increase in the deposition angle from α = 30 • to α = 60 • . This is due to the fact that deposition of the large angle results in the formation of the glancing angle deposited (GLAD) films with high porosity [21]. The densities of the SiO 2 films for various deposition angles and various deposition energies are shown in Table 1. The density changes insignificantly with an increase in the deposition angle from α = 0 • to α = 30 • , while a further increase in α to 60 • leads to a noticeable decrease in density. The stresses demonstrate a similar dependence on α ( Figure 2).
As mentioned in the Section 1, laser-induced stresses can cause thermoelastic breakdown of the coatings structure. Thus, a decrease in stresses with an increase in the deposition angle can be accompanied by an increase in the laser-induced damage threshold (LIDT). Indeed, in [22] the increase in the LIDT of silicon dioxide films with growth of the deposition angle has been observed experimentally.  The results of the stress calculation for high-energy-deposited silicon dioxide films are shown in Figure 3. The absolute values of the stresses are several times higher than for the low-energy-deposited films ( Figure 2). This can be explained by a decrease in the film porosity with an increase in the deposition energy. Indeed, even for α = 60 • the density of the high-energy-deposited film is higher than the density of the low-energy film deposited at α = 0 • (Table 1). Accordingly, the absolute stress values in all high-energy-deposited films are higher than in low-energy-deposited films. In the GLAD film (α = 60 • ), the stresses are approximately half that of the normally deposited film (α = 0 • ). The formation of film layers, with a high value of compress stress under laser irradiation, was demonstrated earlier in the MD simulations [9].
With fast heating (Figure 3a), periodic changes in stress values are observed in the cases of α = 0 • and α = 30 • . The period of these changes is about 10 ps. We assume that these changes indicate the appearance of stress waves caused by fast heating. Assuming that the wavelength of such waves is close to the characteristic dimension of the simulation area, we find that the wave velocity is about a few km/s. This corresponds to the sound velocity in solids. It can be seen, from the plots in Figures 2 and 3, that the stress waves are formed only in the case of picosecond heating of dense films. The slight difference in the periods at α = 0 • and α = 30 • can be caused by the difference in the film density, since the sound velocity depends on it. As mentioned in the Section 1, the thermally induced elastic waves were observed experimentally [3,4]. These waves can significantly affect the LIDT since their propagation leads to the spalling of thin-surface films [4]. Thus, the LIDT value of GLAD SiO 2 films can be higher due to the absence of stress waves. An increase in the heating temperature to T f = 1700 K leads to an increase in the absolute values of thermal-induced stresses in high-energy-deposited films in the case of fast heating (Figure 4a). Stress waves are also observed in this case, but they decay within 20-30 ps after the end of heating. Apparently, this is due to the beginning of the softening of the structure upon heating to 1700 K. This softening prevents the propagation of elastic waves.
In contrast to fast heating, in the case of τ = 1 ns, an increase in the heat temperature T f from 1300 K (Figure 3c) to 1700 K (Figure 4b) leads to a reduction in the absolute stress value. Moreover, a stress value reaches its maximum when the structure temperature is about 1400 K. A further increase in the temperature to 1700 K leads to a decrease in the stress values. We assume that this effect is also related to structure softening. The calculated values of the stresses are in the interval of the experimental data. The stresses during the laser writing process in silica achieve 2 GPa [5]. The laser-induced compressive stress in silicon dioxide films, deposited by ion beam sputtering, achieve 1.6 GPa in the normal conditions of radiation at 1064 nm wavelength and a pulse duration of 5 ns [23]. The intensity of the laser irradiation in [24] is a laser-induced damage threshold of about 25 J/cm 2 . In [24] it was found that the stresses in SiO 2 /HfO 2 multilayers achieve 5 GPa during the laser-induced damage process.
It is interesting to compare the values of thermal-induced stresses and stresses arising during the growth of silicon dioxide films, calculated for high-energy deposition [10]. The stresses in both cases are compressive. However, the maximum absolute values of thermal-induced stresses are several times higher than the stresses in growing films. In both cases, stresses vary insignificantly with increasing the deposition angle α from 0 • to 30 • . When the deposition angle increases from 30 • to 60 • , the relative reduction in stresses in the growing films is greater than the relative reduction in thermally induced stresses. Finally, the thermally induced stresses at the considered heating temperature are noticeably higher than the stresses in the growing film. On the other hand, the thermal stresses, in contrast to the stresses arising during the film's deposition, affect the film only for a limited time.

Conclusions
The laser-induced stresses in dense and porous silicon dioxide films are calculated using MD simulations. The absorption of the laser energy is simulated by the linear temperature growth from room temperature to a final temperature of 1300 K and 1700 K in a time equal to the laser pulse duration. The stresses in high-energy-deposited films are approximately two times higher than in low-energy-deposited films. This difference is caused by the low porosity of the low-energy-deposited films. At a given deposition energy, the stresses in the GLAD films are less than the stresses in normally deposited films. This result agrees with the experimental data.
The picosecond laser pulses produce stress waves in the high-energy-deposited dense films. These waves are not observed in low-energy-deposited films or in high-energy GLAD films. Stress waves can significantly affect LIDT as they lead to spalling of thinsurface films.
For high-energy-deposited films, an increase in the heating temperature to 1700 K leads to an increase in the absolute stress values in the case of picosecond pulses, and to a decrease in these values in the case of nanosecond pulses. We assume that the latter effect is caused by the softening of the films structure, since the temperature of 1700 K is close to the experimentally known softening temperature of fused silica, 2 × 10 3 K.

Appendix A
The characteristic time of the temperature decrease can be calculated using the heat conduction equation, as follows: where a 2 = χ/(сρ), с, χ and ρ are the specific heat capacity, heat conductivity and mass density, respectively, I(z) is the function proportional to the absorbed energy. The direction of the z-axis is shown in Figure A1. The film and substrate thickness are denoted as d and L, respectively ( Figure A1). The following boundary and initial conditions are applied: It is assumed in Equation (A3) that the temperature of the lower edge of the substrate is kept equal to the initial temperature before heating. At the initial moment, the film starts to absorb energy according to the uniform distribution function ( Figure A1). The solution of Equation (A1) is given by the following function: where: At the end of absorption, the temperature distribution T 1 (z, t) is found as the solution of the heat equation with the following conditions: where τ is the laser pulse duration, T(z, τ) is given by Equation (A5). The solution of Equation (A7) is given by the following function: The relative reduction in the temperature at the upper edge of the film (z = 0) can be characterized by the following function: The material parameters are с = 10 3 J/(kg·К), ρ = 2.2 × 10 3 kg/m 3 , and χ = 2.0 W/(m·К) [25].
The graph of the function R(t) is shown on the left side of Figure A2. The substrate thickness is L = 1 mm, the film thickness is d = 1 mkm, which corresponds to the technological dimensions. As can be seen from the plot, the temperature of the upper edge of the film changes insignificantly during the first several tens of nanoseconds. This means that in the MD simulation with nanoseconds trajectories, the temperature after the end of heating is practically constant. If the substrate thickness exceeds 1 mm, then the temperature decreases slower than in the considered case.
The dependence of the temperature on the vertical coordinate is shown on the right in Figure A2. As can be seen, 10 ns after the end of heating, the temperature inside the film layer with z < 0.5 mkm is practically constant. So, the temperature inside the 40 nm thick layer can also be considered as constant during the simulation time of 1 ns or less. Figure A2. Dependence of R(t), Equation (A12), versus time t after the end of film heating, and the dependence of the film temperature T, Equation (A11), along the vertical coordinate z (see Figure 1) at moment t = 10 ns after the end of film heating.