Numerical investigation on the heat deposition characteristics of femtosecond laser pulses undergoing multiple filaments

According to cloud physics theory, there are many physical processes that can affect the formation and growth of snow crystals [27]. Figure 1 briefly summarizes the basic physical effects involved from the filamentation. In this paper, the heat deposition process we consider mainly includes the stage from laser pulse energy deposition to local temperature disturbance (ΔT) and pressure disturbance (Δp) as shown in the yellow area in figure 1, which can be regarded as a very important link between laser physical processes and atmospheric physical processes such as thermal diffusion and turbulent motion.

Zoom In Zoom Out Reset image size

In Fourier space, the evolution of the scalar electric field envelope ${\mathcal E} \left(x,y,z,t\right)$ along the propagation direction z is described by the 3D + 1 unidirectional pulse envelope propagation equation (UPPE) expressed by [28]:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial \tilde{\varepsilon }\left(\omega ,{k}_{\perp },z\right)}{\partial z}=i\left\{{K}_{{\omega }_{0}}\left(\omega ,{k}_{\perp }\right)+\displaystyle \frac{{{\rm{\Delta }}}_{\perp }}{2k\left(\omega \right)}\right\}\tilde{\varepsilon }\\ \,\,\,\,\,\,+i{Q}_{{\omega }_{0}}\left(\omega \right)\left(\displaystyle \frac{\tilde{P}\left(\omega ,{k}_{\perp },z\right)}{{\varepsilon }_{0}}+\displaystyle \frac{i}{\omega }\displaystyle \frac{\tilde{J}\left(\omega ,{k}_{\perp },z\right)}{{\varepsilon }_{0}}\right)\end{array}\end{eqnarray} \tag{ 1 }$$

where $\tilde{ {\mathcal E} }\left(\omega ,{k}_{\perp },z\right)$ denotes the Fourier transform of ${\mathcal E} \left(x,y,z,t\right)$ and ${k}_{\perp }^{2}={k}_{x}^{2}+{k}_{y}^{2}$ is transverse wave vector. The angular frequency of the carrier wave is ω₀ = 2πc/λ₀ (λ₀ = 800 nm). The linear effects of dispersion, diffraction, space-time coupling and normal group-velocity dispersion (GVD) is described by [29]:

$$\begin{eqnarray}&&{K}_{{\omega }_{0}}\left(\omega ,{k}_{x},{k}_{y}\right)\equiv \sqrt{{k}^{2}\left(\omega \right)-{k}_{x}^{2}-{k}_{y}^{2}}\end{eqnarray} \tag{ 2 }$$

The nonlinear dispersion function can be written as:

$$\begin{eqnarray}&&{Q}_{{\omega }_{0}}\left(\omega \right)={\omega }^{2}/2{c}^{2}{K}_{{\omega }_{0}}\left(\omega ,{k}_{x},{k}_{y}\right)\end{eqnarray} \tag{ 3 }$$

In equation (1), the Laplacian operator ${{\rm{\Delta }}}_{\perp }^{2}={\partial }_{x}^{2}+{\partial }_{y}^{2}$ accounts for the beam transverse diffraction. The other two source terms $\tilde{P}\left(\omega ,{k}_{\perp },z\right)$ and $\tilde{J}\left(\omega ,{k}_{\perp },z\right)$ account for the nonlinear polarization and current, respectively. The nonlinear polarization is calculated from the following relation:

$$\begin{eqnarray}&&P=2{\varepsilon }_{0}{n}_{0}{n}_{2}{\left| {\mathcal E} \right|}^{2} {\mathcal E} \end{eqnarray} \tag{ 4 }$$

where n₀ = n(ω₀) and n₂ denote the linear and nonlinear refractive index of the medium, respectively.

The current term in the framework of the Drude model reads:

$$\begin{eqnarray}&&\begin{array}{l}\tilde{J}\left(\omega ,{k}_{\perp },z\right)=c\sigma {\varepsilon }_{0}\left(1+i{\omega }_{0}{\tau }_{c}\right)\rho \tilde{\varepsilon }\\ \,\,\,\,\,+c{n}_{0}{\varepsilon }_{0}\displaystyle \frac{W\left(I\right)K\hslash {\omega }_{0}}{I}\left(1-\displaystyle \frac{\rho }{{\rho }_{at}}\right)\tilde{\varepsilon }\end{array}\end{eqnarray} \tag{ 5 }$$

where σ denotes the cross section for inverse bremsstrahlung and τ_c is the collision time. The intensity dependent ionization rate $W\left(I\right)$ scales as $W\left({\left|\varepsilon \right|}^{2}\right)={\sigma }_{K}{\left|\varepsilon \right|}^{2K}$ and ${\sigma }_{K}={\beta }_{K}/\left(K\hslash {\omega }_{0}{\rho }_{at}\right),$ where $K\equiv \left\langle {U}_{i}/\hslash {\omega }_{0}+1\right\rangle$ denotes the number of photons at the frequency ω₀, U_i is the ionization potential energy, ρ_at is the neutral atom density, and β_K denotes the ionization cross section for MPA with coefficient σ_K corresponding to the number of K-photons. In the laser-induced snow formation experiment, because the distance of laser transmission in the cloud chamber is relatively short, we ignore the cumulative impact of possible existing background water droplets or tiny ice particles on the transmission of fs laser pulses as discussed in [30, 31].

The evolution of electron density ρ is described as:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=\left({\sigma }_{K}{\left|\varepsilon \right|}^{2K}+\displaystyle \frac{\sigma }{U}{\left|\varepsilon \right|}^{2}\rho \right)\left(1-\displaystyle \frac{\rho }{{\rho }_{at}}\right)\end{eqnarray} \tag{ 6 }$$

The two terms on the right-hand side of equation (6) account for multiphoton ionization (MPI) and avalanche ionization (inverse Bremsstrahlung), respectively. A list of parameters and their values used in the simulations is provided in table 1.

In the simulation, equation (1) is solved using full-step (not split-step) Fast Fourier Transformation in transverse space and time dimensions and equation (6) is solved via the fourth-order Runge–Kutta method [32].

A stronger laser filament would induce more heat release through plasma recombination [33]. The instantaneous increase of temperature and pressure has a very important influence on the subsequent thermodynamic processes. Because the pulse transit time (∼ps) and the thermalization time (∼ps) are much shorter than the air flow motion time scales, the process of heating due to filamentation can be considered as an isochoric (constant volume) process [34]. The first law of thermodynamics is used to estimate the temperature variation caused by the heat release, which can be expressed as:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}E={c}_{v}{\rho }_{0}\displaystyle {\int }_{{{\mathbb{R}}}^{3}}\left(T\left(\vec{r}\right)-{T}_{gas}\right){d}^{3}\vec{r}\\ \,\,\,=\,\displaystyle \frac{{c}_{v}}{{k}_{B}}\displaystyle {\int }_{{{\mathbb{R}}}^{3}}\left(p\left(\vec{r}\right)-{p}_{gas}\right){d}^{3}\vec{r}\end{array}\end{eqnarray} \tag{ 7 }$$

where c_v is the gas isochoric heat capacity per molecule. At p = 1013 hPa and T = 300 K, air can be approximately regarded as an ideal diatomic molecular gas and c_v = 5k_B/2, where k_B = 1.38 × 10⁻²³ J · K⁻¹ is Boltzmann constant. ρ₀ = 2.7 × 10²⁷ m⁻³ is the neutral density of gas. In the case of an initial temperature and pressure profile with cylindrical symmetry, equation (7) can be further simplified as:

$$\begin{eqnarray}&&{\rm{\Delta }}E\approx \pi {c}_{v}{\rho }_{0}\displaystyle {\int }_{{\mathbb{R}}}\left({T}_{\max }\left(z\right)-{T}_{gas}\right){r}_{0}{\left(z\right)}^{2}dz\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}E\approx \displaystyle \frac{\pi {c}_{v}}{{k}_{B}}\displaystyle {\int }_{{\mathbb{R}}}\left({p}_{\max }\left(z\right)-{p}_{gas}\right){r}_{0}{\left(z\right)}^{2}dz\end{eqnarray} \tag{ 9 }$$

where r₀(z) denotes the half width at half maximum (HWHM) of the temperature profile at z [34]. For simplicity, we assume it is approximately equal to radius of the filament (≈150 μm) in this study. Based on the above simplifications, the equations (8) and (9) can be further reduced to:

$$\begin{eqnarray}&&{T}_{\max }\left(z\right)-{T}_{gas}\approx \displaystyle \frac{1}{{c}_{v}{\rho }_{0}\pi {r}_{0}^{2}}\displaystyle \frac{\partial {\rm{\Delta }}E}{\partial z}\left(z\right)\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{p}_{\max }\left(z\right)-{p}_{gas}\approx \displaystyle \frac{{k}_{B}}{\pi {c}_{v}{r}_{0}^{2}}\displaystyle \frac{\partial {\rm{\Delta }}E}{\partial z}\left(z\right)\end{eqnarray} \tag{ 11 }$$

Thus, the maximum local temperature and pressure variation are approximately proportional to $\partial {\rm{\Delta }}E/\partial z,$ defined as the energy deposition rate per unit length along the propagation axis, which can be calculated by [35]:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Delta }}E}{\partial z}=\displaystyle \iiint u(x,y,z,t)dxdydt\end{eqnarray} \tag{ 12 }$$

where $u\left(x,y,z,t\right)$ denotes the density of nonlinear losses. The plasma recombination mainly accounts for heat release [33]. Therefore, according to equation (6), we estimated the density of nonlinear losses with the following equation:

$$\begin{eqnarray}&&u(x,y,z,t)=\left(\sigma \rho I+{\beta }_{K}{I}^{K}\right)\times \left(1-\rho /{\rho }_{at}\right)\end{eqnarray} \tag{ 13 }$$

As above, the maximum local temperature and pressure variation caused by the released heat from the filament can be calculated based on equations (10) to (13).

The initial pulse investigated in this paper is a Gaussian beam and its spatiotemporal envelope of the electric field is given as:

$$\begin{eqnarray} &&{\mathcal E} =\sqrt{2{P}_{in}/\pi {w}_{0}^{2}}\times \exp \left(-\displaystyle \frac{{x}^{2}+{y}^{2}}{{w}_{0}^{2}}-\displaystyle \frac{{t}^{2}}{{t}_{p}^{2}}-i\displaystyle \frac{k\left({x}^{2}+{y}^{2}\right)}{2f}\right)\end{eqnarray} \tag{ 14 }$$

where P_in is the input power, w₀ = 3 mm is the beam width, and t_p = 25 fs is the pulse duration. The energy of the input pulse is ${E}_{in}={P}_{in}{t}_{p}\sqrt{\pi /2}.$ The initial beam passes through a lens with focal length f.

Figure 2 shows variation of on-axis intensity, on-axis electron density, and beam radius along propagation distance z for laser beams with E_in = 4.5 mJ and E_in = 7.1 mJ under external lens of focal length f = 30 cm and f = 50 cm, respectively. As seen from figure 2(a), the on-axis intensity of filament with f = 30 cm lens is higher than that with f = 50 cm lens, which leading to higher electron density in the filament (figure 2(b)). For a given focal length, the on-set distance of filamentation (as shown by the position where the on-axis intensity increases rapidly) is very close for pulses with E_in = 4.5 mJ and E_in = 7.1 mJ. When the input pulse energy increases from 4.5 mJ to 7.1 mJ, the peak on-axis electron density decreases whereas the length of filament tends to increase for pulse with f = 30 cm focal length as roughly seen from the evolution of on-axis intensity. For the pulse with f = 50 cm focal length, the electron density increases whereas the length of filament decreases, which is opposite to the case of f = 30 cm lens. As seen from figure 2(c), the radius of filament increases with the input energy increases. When fixing input energy, the radius tends to be slightly larger under longer external focusing.

Zoom In Zoom Out Reset image size

In order to acquire a deeper insight into effects of input pulse energy and external focusing on the temporal dynamics of filaments, we present the evolution of the temporal distribution of the on-axis intensity in figure 3. As shown in these figures, the external focusing strongly affects the distribution of the intensity before and after the focus lens in the clamping region. Clear on-axis pulse splitting effect can be seen at the focus position. Comparing the two figures on the left column of figures 3(a) and (b), the values of their on-axis intensity almost has no difference. However, more obvious pulse splitting (indicated by white dummy arrow) appears after the focusing position for the pulse with input energy E_in = 7.1 mJ. At the post-filamentation region, an intense few-cycle pulse with a duration of ∼50 fs forms in both cases. Comparing the two figures on the right column of figures 3(a) and (b), the pulse with input energy E_in = 7.1 mJ experiences more obvious on-axis pulse splitting than that of pulse with input energy E_in = 4.5 mJ. In particular, at the position z ≈ 0.60 m, the trailing part of the pulse is completely separated from the front of the pulse (indicated by white dotted ellipse). The pulse with input energy E_in = 4.5 mJ forms an intense few-cycle pulse with a duration of ∼40 fs, which lasts longer transmission distance than that of pulse with input energy E_in = 7.1 mJ.

Zoom In Zoom Out Reset image size

The spatial distribution of energy fluence, which is defined as the total energy per unit area at different distances is depicted in figure 4. The energy fluence is calculated by:

$$\begin{eqnarray}&&F\left(x,y,z\right)=\displaystyle {\int }_{-\infty }^{+\infty }{\left|E(x,y,t,z)\right|}^{2}dt\end{eqnarray} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

As shown in figure 4, the highlighted part of the profiles clearly indicates the on-set distance and its distribution of filaments. Only one main defocusing-refocusing cycle is visible in the filament for the four cases. Near the nonlinear focus, clear ripples emitted from the core of the filament are visible. Their peak fluence is all close to 7.0 J · cm⁻². The spatial splitting of the pulse can be clearly observed after the focus position, as indicated by white dummy arrow in left column of figure 4(a). In the right column of figure 4(b), there is also spatial splitting of the pulse occurring before the fluence reaches the maximum value (indicated by white dummy arrow). The different location of multiple filament formation may have important influence on characteristics of heat deposition due to the energy competition between filaments [36]. In addition, due to the formation of multiple filaments under f = 50 cm focal length, the peak fluence slightly decreases.

To appreciate the evolution of the overall structure of the waveform more clearly, we present the transverse profile of fluence at different propagation positions in figure 5. As shown in the figures, the radius of the beam decreases during the self-focusing process (z = 0.2 m in figures 5(a) and (b), z = 0.4 m in figures 5(c) and (d)). With the filament propagation, the radius of the beam gradually becomes larger. The hot spots in the transverse section clearly indicate the formation of a structure made up of multiple filaments. During the filamentation, the intensity and numbers of filaments increase with the growth of initial pulse energy for a given external focal length of lens. When fixing the initial energy, the pulses also tend to form higher intensity and more numbers of filaments under stronger focusing.

Zoom In Zoom Out Reset image size

As the energy deposition of filamentation and subsequent heat release is the primary cause of the air-flow motion [22], which is a key mechanism affecting the efficiency of the snow formation, the influence of the multiple filament formation on these processes was further investigated based on the simulation results. Here, the difference between the input energy and output energy (ΔE = E_in−E_out) is used to define the energy deposition. Meanwhile, the ratio of deposited energy to input energy (E_out/E_in), which is defined as the relative deposited energy, is also used to understand the characteristics of energy deposition.

Figure 6 presents the evolution of the total deposited energy, absorbed deposited energy rate, estimated local maximum temperature variation, and local maximum pressure variation (see section 2.3) along the propagation distance z. It should be noted that the energy for pulse with E_in = 4.5 mJ is increased by 2.6 mJ for comparison. As shown in figure 6(a), the total deposited energy of pulse with E_in = 7.1 mJ is significantly larger than that of pulse with E_in = 4.5 mJ. When fixing the input pulse energy, the total deposited energy becomes higher for laser pulse with a 30 cm focal length than a 50 cm focal length. The relative deposited energy is 36.76% and 32.75% for pulse with E_in = 4.5 mJ when f = 30 cm and f = 50 cm, respectively. When the input pulse energy increases to 7.1 mJ, the relative deposited energy becomes 43.62% and 40.33% for pulses with f = 30 cm and f = 50 cm. The relative deposited energy increases significantly with input energy. The energy deposition rate is also very different, which can be seen from figure 6(b). The laser pulse with E_in = 7.1 mJ and f = 30 cm can achieve the highest energy deposition rate. Its maximal value is about 95 mJ · m⁻¹, which is much larger than that of laser pulse with E_in = 7.1 mJ and f = 50 cm (58 mJ · m⁻¹), and also larger than laser pulse with E_in = 4.5 mJ and f = 30 cm (62 mJ · m⁻¹).

Zoom In Zoom Out Reset image size

According to the maximum temperature variation range of the transmission path depicted in figure 6(c), the local temperature variation range and the increase amplitude are much higher for the filaments with f = 30 cm. The maximum pressure experiences the similar variation as temperature. This is favorable to the formation of the stronger airflow motion and eventually more mass of snow in the cloud chamber.

In order to adequately understand the influences of input pulse energy on the characteristics of laser filaments, figure 7 shows the simulated dependence of various filamentation parameters including length of filament (L_f), maximum electron density (ρ_max), peak fluence (F_max), total deposited energy (ΔE), maximum energy deposition rate ($\partial {\rm{\Delta }}E/\partial z$), and maximum local temperature variation (ΔT) as a function of input energy. Here, we define the filament length as the length from the onset distance to the value of the peak plasma density less than 10¹⁴ cm⁻³ at the last beam collapse [37]. As shown in figure 7(a), the increase of filament length is limited by the formation of multiple filaments with the growth of the pulse energy. The input energy could also affect the intensity of filament. We see that both the electron density and the peak fluence in the filament increase with the input energy increases (figures 7(b) and (c)). In particular, the increase for cases with f = 30 cm is more obvious when the energy increases from 3.0 mJ to 7.1 mJ. Due to the generation of more electrons with the growth of initial energy, the deposited energy also increases with the input energy increases (figure 7(d)). When fixing the input energy, the maximum electron density of pulses with f = 30 cm is higher than that of pulses with f = 50 cm, which also leads to the higher deposited energy in the former case. We also calculate the rate of absorbed deposited energy by the gaseous medium. The evolution of maximum rate of deposited energy is presented in figure 7(e). It can be seen that the pulses with higher input energy can obtain a larger rate of deposited energy, especially for pulses with f = 30 cm. The larger rate of deposited energy indicates the higher increase of local temperature and pressure (figure 7(f)), which could help to create stronger airflow and more snow in the cloud chamber.

Zoom In Zoom Out Reset image size

Based on the above analysis, we can see that both the input energy and external focusing seriously influence the filaments formation. The numbers and intensity of filaments increases with the increase of input energy. The increase of filament length is limited by the formation of multiple filaments. For pulses with external lens of f = 50 cm, the filament length involves slightly decrease with the input energy increase. For pulses with a given initial energy, the formed filaments under strong external focusing can deposit higher energy with the increase of input energy due to the higher multiphoton ionization and plasma absorption. The deposited energy rate also evolves significantly increase with the input pulse energy under strong external focusing, which indicates the local maximum temperature and pressure variation also increases with the input energy increases. This is beneficial to the increase of airflow velocity and mass of snow with the input energy under f = 30 cm as reported by Sun et al [25]. However, for pulses with external lens of f = 50 cm, the deposited energy rate of the formed filaments experiences less significant increase with the increase of input energy than lens of f = 30 cm. We speculate that this might be related to the formation of multiple filaments before the focus position. As variation of local temperature and pressure can further influence the microphysical processes like turbulence, evaporation, condensation, and freezing as it in nature, we can expect that the differences in airflow velocity and mass of snow for pulses with different input energy may be primarily influenced by the heat deposition of filaments.