## 2 Modeling Electromagnetic Fields in Nonlinear Media  
The Equations (1)-(3) model the laser beam propagation by means of the electromagnetic fields \(\vec{E}(\vec{r},t)\) and \(\vec{H}(\vec{r},t)\). Throughout the light propagation path, fields generate and couple with the density of free charges as free electrons \(n_{fe^{-}}(\vec{r},t)\), which in turn interact with the local field. In the given form, nonlinearities are determined by the electric displacement field \(\vec{D}(\vec{r},t)\) in the Ampere-Maxwell Equation (1). Meanwhile, the evolution of the free-electron density or medium conductivity is described by Equation (3), where \(G(\vec{r},t)\) is the carrier generation rate and \(R(\vec{r},t)\) is the carrier recombination rate. \[\
abla\wedge\vec{H}(\vec{r},t) = \frac{\partial\vec{D}(\vec{r},t)}{\partial t} \tag{1}\] \[\
abla\wedge\vec{E}(\vec{r},t) = -\mu\frac{\partial\vec{H}(\vec{r},t)}{\partial t}\] (2) \[\frac{\partial n_{fe^{-}}(\vec{r},t)}{\partial t} = G(\vec{r},t)-R(\vec{r},t) \tag{3}\]The electric displacement field is composed by the electric field contribution plus the polarization vectors \(\vec{D}(\vec{r},t)=\varepsilon_{0}\varepsilon_{\infty}\vec{E}(\vec{r},t)+\sum \vec{P}(\vec{r},t)\)[2]. The magnitude \(\varepsilon_{\infty}\) is the non-unity high-frequency relative permittivity of the media. On the right hand side of the displacement electric field, the sum of polarization vectors in the model that we are considering \(\sum\vec{P}(\vec{r},t)\) is the contribution of four terms, the metal polarization \(\vec{P}_{Metal}(\vec{r},t)\), the Kerr effect polarization \(\vec{P}_{Kerr}(\vec{r},t)\), the Raman polarization \(\vec{P}_{Raman}(\vec{r},t)\) and the free charge polarization vector \(\vec{P}_{fe^{-}}(\vec{r},t)\) such as \(\sum\vec{P}(\vec{r},t)=\vec{P}_{Metal}(\vec{r},t)+\vec{P}_{Kerr}(\vec{r},t)+ \vec{P}_{Raman}(\vec{r},t)+\vec{P}_{fe^{-}}(\vec{r},t)\) with the detail of each term below:  
\[\vec{P}_{Metal}(\vec{r},t) = \frac{\varepsilon_{0}a}{j\omega+b}\vec{E}(\vec{r},t)\ \ \ \ \ a,b\in\mathbb{C} \tag{4}\] \[\vec{P}_{Kerr}(\vec{r},t) = \alpha e_{0}\chi_{0}^{(3)}\vec{E}(\vec{r},t)\int_{-\infty}^{t} \delta(t-\tau)|\vec{E}(\vec{r},\tau)|^{2}d\tau\] (5) \[\vec{P}_{Raman}(\vec{r},t) = \Psi\vec{E}(\vec{r},t)\int_{-\infty}^{t}e^{-\frac{(t-\tau)}{ \tau_{2}}}\sin(\frac{t-\tau}{\tau_{1}})U(t-\tau)|\vec{E}(\vec{r},\tau)|^{2}d\tau\] (6) \[\vec{P}_{fe^{-}}(\vec{r},t) = \int_{-\infty}^{t}\sigma_{fe^{-}}(\vec{r},t)\vec{E}(\vec{r},t)d\tau \tag{7}\]  
where \(\Psi=\varepsilon_{0}\chi^{(3)}(1-\alpha)\frac{\tau_{1}^{2}+\tau_{2}^{2}}{\tau _{1}\tau_{2}\tau_{2}\tau_{2}\tau_{1}\tau_{2}\tau_{2}\tau_{1}\tau_{2}\tau_{2} \tau_{1}\tau_{2}\tau_{2}\tau_{1}\tau_{2}\tau_{2}\tau_{1}\tau_{2}\tau_{2}\tau_{ 1}\tau_{2}\tau_{2}\tau_{1}\tau_{2}\tau_{2}\tau_{1}\tau_{2}\tau_{2}\tau_{1}\) is a coefficient gathering parameters depending on the considered medium. Expressions Equations (5) and (6) providing the Kerr effect polarization and Raman polarization, respectively, are taken from [9]. Here, \(\alpha\) is a real-valued constant in the range \(\alpha\in[0,1]\) that parameterizes the relative strengths of Kerr and Raman contributions, \(\delta(t)\) refers to the Dirac delta function that models the Kerr nonresonant virtual electronic transitions, which is medium dependent and in the order of few femtoseconds for fused silica, \(U(t)\) is the Heaviside unit step function, \(e^{\frac{-t}{\alpha}}\sin(\frac{t}{\tau_{1}})\) models the impulse response of a single Lorentzian relaxation centered on the optical phonon frequency \(1/\tau_{1}\) and having a bandwidth of \(1/\tau_{2}\). Finally \(\chi_{0}^{(3)}\) is the strength of the third-order nonlinear electric susceptibility. The upper limit of the integrations in Equations (5)-(7) extends only up to \(t\) because the response function must be zero for \(\tau>t\) to ensure causality [3]. Moreover, free-charge polarization is introduced by using a free-charge density of current, which will be defined later.