## 1 Introduction  
When two beams of light impinge on an optically nonlinear material simultaneously, under certain conditions, they interact in such a way that the energy is transferred from one beam to the other; this phenomenon is known as two-beam coupling [1]. Optical pump-probe spectroscopy is a technique for studying non-equilibrium electronic dynamics based on such a coupling of pump and probe beams, which finds applications across a range of fields, from physics and chemistry to materials science and biology [2]. This technique involves focusing a short pulse of a strong pump beam and, after a certain delay time, a weak probe beam on a sample to record the transmittance or reflectance of the probe light. By adjusting the delay between the pulses, one can study the transient evolution of the out-of-equilibrium electrons. In particular, different versions of the pump-probe technique have been widely used for studying the electronic response of graphene to fast optical excitations [3, 4, 5, 6, 7, 8, 107, 11]. The nonlinear optical behavior of graphene has been a subject of intense debate ever since the initial theoretical prediction in 2007, suggesting that graphene could exhibit an exceptionally high nonlinear optical response in a broad spectral range [12]. By that time, graphene had already garnered recognition for its remarkable linear optical properties, attributed to its distinctive conical band structure. Demonstrating its nonlinear optical properties would further solidify its standing as an ideal material for optical applications with potential implications for advanced photonics and optoelectronics [13]. Since then, substantial and continuous efforts have been dedicated to the experimental characterization of graphene's nonlinear properties. This includes the observation of higher harmonics [14, 15, 16],the Kerr effect [17; 18], saturable absorption [19; 20], and four wave-mixing [21; 22], along with other nonlinear effects. While there may not be complete agreement on the measured quantities across different articles, the majority of them report high values for graphene's nonlinear properties [23]. Theoretical investigations, particularly concerning the third-order nonlinearity, have taken diverse routes. Early approaches applied the Boltzmann equation to investigate the nonlinear optical response of graphene owing to intraband electronic transitions [12; 24]. To incorporate interband contributions, researchers employed the Bloch equations solved using perturbation theory approximations, which allowed for the calculation of the third-order tensor of conductivity, \(\sigma^{(3)}(\omega_{1},\omega_{2},\omega_{3})\)[25; 26; 27]. Although more intricate, this approach also allows for the calculation of even higher orders in the electric field. Additionally, since perturbation theory enables one to obtain explicit expressions for the conductivity tensor, it facilitates the exploration of all possible effects arising from the third-harmonic generation, difference frequency, and saturable absorption. However, a fundamental problem arises with this approach: perturbation theory's validity hinges on small field intensities. In practical nonlinear optics experiments, the intensities of the involved beams are very high, meaning that the theoretical description of the system goes far beyond the limits of perturbation theory. The electromagnetic fields applied to graphene are able to create populations of carriers that cannot be described assuming equilibrium chemical potentials (\(\mu_{0}\)) and temperature (\(T_{0}\)) for the carriers throughout the calculations, as is the case of the perturbative approaches mentioned above. This seems to be unjustified for typical conditions of nonlinear optical experiments on graphene. In recent experiments, [PERSON] et al. employed high-power THz pulses and showcased graphene's remarkable proficiency in generating higher-order THz harmonics by shining THz radiation onto graphene under normal incidence [28].