## 1 Introduction  
The finite difference time-domain (FDTD) method enables the computation and modeling of light propagation and scattering processes in linear and nonlinear dispersive materials possessing wavelength-dependent and intensity-dependent properties. This ranks the method among the most universal and powerful numerical tools for optics, electrodynamics, antennas, and waveguides theory. In this context, the general vector auxiliary differential equation (GVADE) method, which has been reported in [1] and improved in [2], is one of the most popular methods to deal with nonlinear media. However, the method suffers from some drawbacks that limit its speed and stability due to the semi-implicit updating relation for the current polarization densities. In summary, the method is not appropriate or recommended for a broad range of problems, especially those that address the designs with complex 3D topology, where iterative solutions of coupled nonlinear equations can be unstable and resource-demanding [3]. Some direct or explicit methods take into consideration the Kerr effect in the Born-Oppenheimer approximation [4]. The Kerr effect arises from the nonlinear third-order electric susceptibility. In [3], the authors take into consideration the Kerr effect in an explicit form. Other direct methods [5] require excessive computational burdens, especially for the consideration of complex geometries. When metal-dielectric interfaces are treated, especially at nanometric scales, plasmonic effects appear [6, 7] and generate high concentrations or high local densities of the electromagnetic field [8]. Since dielectrics experience nonlinear optical processes, the treatment of these regions in the interfaces is particularly challenging. In this manuscript, we develop an algorithm that permits solving a nonlinear dielectric medium in which [PERSON] and [PERSON]'s effects [9] occur simultaneously. The Raman effect arises from the scattering caused by the conversion of photon energies into vibrational energy of molecules. The scattered photons have less energy than the incident photons corresponding, therefore lowering frequency fields [2]. We also assume that this takes place in an inhomogeneous dielectric environment that can potentially include metallic heterogeneities with specific optical properties. Along with the laser propagation, provided that the intensity is sufficient to ensure nonlinear absorption, the medium can also turn into a metallic state as a free electron population is generated during multiphoton absorption (MPA) or tunneling and avalanche ionization processes. It is well known that the electron transition from the valence band (VB) to the conduction band (CB) is followed by various processes, including heating in the CB, impact ionization, and various kinds of electron collisions. The purpose of this work consists of computing all these events in a full explicit form. The carrier generation rate \(G(\vec{r},t)\) has been computed owing to the BVkP approach implying Bloch-Volkov states [10]. The algorithms presented in this manuscript are especially suitable for simulations dealing with confinements and guiding high-energy pulsed beams. In summary, in this work, we coupled Maxwell equations and the electric charge carriers conservation equation in which the quantum ionization and recombination are governed by BVkP self-consistently. In the next sections, by following and enriching some pioneers' works [11, 12, 13], we detail a fully explicit, stable algorithm based on GVADE FDTD, able to solve the curl Maxwell's equation in a frequency-dependent way for a nonlinear medium. In order to illustrate the combination of all these nonlinear effects occurring simultaneously in a simulation, we have modeled a Bessel beam traveling through dielectric fused silica, where some metal inclusions are modeled in a similar form as in Ref. [14]. In Section 2, we introduce the theoretical model, which describes the nonlinear effects.