## 2 Methods  
To investigate the propagation dynamics of a femtosecond pulse in a gas medium, we modeled the propagation of a laser pulse with the nonlinear Schrodinger equation (NLSE), which governs the evolution of the electric field envelope \(E\)(x,y,z,t) of the pulse in the reference frame moving at the group velocity \(v_{g}=\partial\omega/\partial k\big{|}_{\omega_{0}}\)[9]:  
\[\frac{\partial E}{\partial z}=\frac{i}{2k_{0}}\Delta_{\perp}E-\frac{ik^{\prime \prime}}{2}\frac{\partial^{2}E}{\partial\tau^{2}}+i\frac{\omega_{0}}{c}n_{2}| E|^{2}E-k_{0}\frac{\omega_{pe}^{2}}{\omega_{0}^{2}}E-\sum_{s=\text{N}_{2},\text{O}_{2},\text{H}_{2}\text{O}}\frac{\beta_{s}^{(K_{s})}}{2}|E|^{2K_{s}-2}E \tag{1}\]  
where \(E\)(x,y,z,t) represents the envelope of the electric field along propagation direction z, and \(\tau\) refers to the retarded time variable, t-z/v\({}_{g}\), according to group velocity as the reference system. \(k_{0}=2\pi\)/\(\lambda_{0}\) and \(\omega_{0}=2\pi\)c/\(\lambda_{0}\) are the central wave number and the angular frequency of the carrier wave corresponding the wavelength \(\lambda_{0}=248\) nm, respectively. The Laplacian operator \(\Delta_{\perp}^{2}=\frac{\partial^{2}}{\partial\tau^{2}}+\frac{\partial^{2}} {\partial\tau^{2}}\) denotes the beam transverse diffraction. The second terms on the right-hand side of Equation (1) account for group velocity dispersion (GVD) with a coefficient of \(k^{\prime\prime}=1.21\times 10^{-28}\) (s\({}^{2}\)/m). The Kerr effect is related to \(n_{2}=8.0\times 10^{-19}\) cm\({}^{2}\)/W [12], which denotes the nonlinear refraction index of air. Plasma defocusing is described by the plasma oscillation frequency, \(\omega_{pe}=\sqrt{q_{e}^{2}\rho/m_{e}\epsilon_{0}}\) (\(q_{\text{e}}\) is the electron charge, \(m_{e}\) is the electron mass, and \(\rho\) is the free electron density). The last term in Equation (1) accounts for multiphoton absorption (MPA). The MPA term accounts for energy absorption due to multiphoton ionization [9]. The coefficient \(\beta_{s}^{(K_{s})}=K_{s}\hbar\omega_{0}\rho_{at,s}\sigma^{K_{s}}\) is related to the multiphoton ionization coefficient, which corresponds to the number of \(K\)-photons needed to ionize the corresponding gas molecule.  
The evolution equation of the electron density \(\rho_{\text{es}}\) is governed by:  
\[\frac{\partial\rho}{\partial t}=\sum_{s}\left(\frac{\beta_{s}^{(K_{s})}}{K_{ s}\hbar\omega_{0}}|E|^{2K_{s}}+\frac{\sigma}{U_{s}}\rho_{\text{es}}|E|^{2} \right)\left(1-\frac{\rho_{\text{es}}}{\rho_{at,s}}\right) \tag{2}\]  
where \(\sigma=5.1\times 10^{-22}\) is the cross section for inverse Bremsstrahlung, \(\hbar=h/2\pi\) is Planck constant, \(U_{s}\) is ionization potential of the gas, and \(\rho_{at,s}\) denotes the initial neutral gas atom density. The coefficients used in Equations (1) and (2) are listed in Table 1.  
The initial pulse investigated in this paper is modeled by a Gaussian beam:  
\[E(x,y,z=0,t)=\sqrt{\frac{2P_{in}}{\pi w_{0}^{2}}}\exp\left(-\frac{x^{2}+y^{2}}{w_{ 0}^{2}}-\frac{t^{2}}{t_{p}^{2}}\right) \tag{3}\]  
where \(P_{in}\), \(w_{0}\), and \(\tau_{p}\) are the input peak power, initial pulse radius, and the pulse duration, respectively.  
A standard split step Crank-Nicholson scheme is used to integrate the linear part of Equation (1) along the propagation axis z and the fourth-order Runge-Kutta method is employed to solve the evolution equation (Equation (2)) [14].