Preprint becomes especially important when the state of the gas is out of equilibrium, giving rise to non-Boltzmann distributions in the internal energy states of the molecules. For example, this is a common occurrence in the forebody region of spacecraft entering planetary atmospheres in hypersonic environments. The most reliable approach to modeling non-equilibrium flows relies on the solution of the master equation , whereby the population of each internal energy level of molecules is explicitly computed. However, the cost associated with such simulations makes them impractical for computational fluid dynamics (CFD) applications. As such, there is a need for physics-driven data reduction methods that can extract patterns from experimental and simulation data, available only for a few operating conditions, and accurately extrapolate them to new scenarios. Recently, this has led to an interest in using reduced-order modeling techniques and data-driven machine learning to investigate non-equilibrium chemical kinetics . For example, the authors have successfully employed physics-based coarse-graining approaches for constructing reduced-order models (ROMs) . While they have led to significant speed-ups in the simulation of non-equilibrium kinetics and applications to CFD, current ROMs did not bypass the necessity of complete ab initio StS databases. Nevertheless, the calculations of these quantities are computationally burdensome, as they not only depend upon the availability of accurate potential energy surfaces (PESs), but they also require the simulation of billions of collisions via the quasi-classical trajectories (QCT) method . Despite that, few techniques have been developed to reduce the computational cost of QCT calculations for StS kinetic. The inspiring work by Koner et al. was the first attempt to construct and train a machine learning-based surrogate capable of completing a sparse database of reaction rate coefficients for hypersonic applications. However, their regression task was artificially simplified by statistical averaging (i.e., by averaging the StS cross-sections over windows of vibrational and rotational quantum numbers). This solution leads to the loss of fine structural dependencies due to the rotational-vibrational anisotropy, and it finally generates inaccuracies in the microscopic quantities of interest (Fig. ). This brings us to the central questions we aim to address in this work: Is it possible to accurately reconstruct the physics underlying reaction rate coefficients? Can this be done using just a few samples from the StS rate coefficients database? The key to answering these questions lies with (i) identifying the patterns underlying reaction rates, (ii) recognizing the latent physical variables that can be associated with these patterns, and (iii) employing machine learning architectures that allow efficient and automatized extractions of those patterns, together with their representations as functions of the latent physical variables. To this end, we look at reaction rate coefficients from different perspectives. In Fig. , inelastic reaction rate coefficients are separately viewed in the space of the initial, i, and final, j, state indices (left figure), and in the space of internal energy, ε Int , and rotational quantum number, J (right figure). On the one hand, the left plot in Fig. does not highlight any clear pattern, as the i and j indices are only convenience state labeling variables. This fact becomes even clearer when looking at the zoomed-in section of the reaction rate matrix. Characterizing the levels by their indices would imply mere memorization of the matrix and only a partial extraction of the complex underlying physics. Limited generalization capabilities would characterize the resulting surrogate, which would also be fully precluded for generative purposes (i.e., extrapolation to different chemical systems). On the other hand, the right plot in Fig. highlights more interesting patterns. Only a few inelastic reaction rate coefficients (i.e., five of the columns in the figure on the left) are represented on top of the effective diatomic potential, which consists of a centrifugal term added to the rotation-less potential energy . The grey dots in the figure are the different energy states that belong to the diatom, and the oblique strands stretching out from J = 0 to J max represent states having the same vibrational quantum number, v. Consistently, a vertical cut along the z-direction selects all vibrational numbers at a given J. The five distinct ellipsoids in the right plot correspond to five initial states (red dots). The patterns show that the inelastic rate coefficients vary anisotropically in the plane and that the degree of anisotropy is a quantity proper of the initial level . These dependencies have already been observed by Sahai et al. , and their spectral clustering strategy was developed to take these patterns into account. More in detail, for initial states denoted by low vs and low Js, the patterns are stretched out to almost flat lines, with transitions predominantly vibrational invariant. Similarly, in the regions characterized by high vs and low Js, the vibrationally elastic transitions (i.e., reactions that favor changes in J) are still preferred. On the contrary, in the high-J sections of the diatomic potential, the eccentricity of the ellipsoid is significantly lower, showing that the vibrational and rotational transitions are almost equally probable in this region. Importantly, the dependencies highlighted by the right plot in Fig. point to the fact that the physics behind the rate coefficient patterns is strictly connected to the properties of the effective diatomic potential. This paper aims to construct a surrogate capable of reproducing the level-specific ellipsoidal patterns after training with small subsets of a few temperature-dependent reaction matrices. The first approach that might come to mind for accomplishing this task is singular value decomposition (SVD). For example, we could decompose each of the matrices, find the main bases, and represent the level-specific ellipsoids as compositions of such modes. However, compared to classic SVD-based methods, multiple challenges must be tackled. Firstly, the StS reaction data should be represented in the space of the effective diatomic potential properties of the initial and final states rather than being described as a mere sequence of matrices characterized by the convenience indexes i and j. This is necessary, as already explained, to be consistent with the underlying physics and to guarantee generalization capabilities. We highlight that the level-specific properties of the effective diatomic potential are high-dimensional and quantized (i.e., discrete but not arranged on a uniform multidimensional grid). Secondly, the temperature dependencies would also need to be considered (i.e., tensor decomposition). Thirdly, since the ultimate goal is to reduce the computational cost of QCT calculations, we should expect the StS reaction data to be significantly sparse and only a few columns of the reaction matrix in Fig. to be available for training (see Fig. ). To address these requirements, the deep operator network (DeepONet) architecture is employed here to model and predict excitation rate coefficients (i.e., the coefficients for inelastic and exchange reactions). The DeepONet structure was developed initially as a neural operator (i.e., as a general surrogate built upon neural networks to represent the solution of a family of PDEs without requiring numerical integration). While, Preprint in this problem, the reaction rate coefficients are integrated quantities and are not governed by an explicit operator, DeepONet is still an effective emulating architecture. In fact, Venturi and Casey recently showed that DeepONet is interpretable as a continuous extension of SVD that works with multi-dimensional spaces and does not require grid-structured data . The DeepONet architecture comprises two building blocks connected by a dot-product layer. The first block is called the trunk net and outputs the modes of a linear decomposition, while the second block, called the branch net, generates the related coefficients. We stress that the linear combination of modes must accurately reproduce the ellipsoidal kernels (e.g., the patterns in the right plot of Fig. ). Consequently, the modes are represented here in the space of the final state's level properties centered at the initial state's value, and ∆x = x jx i (with x being a generic property) are given as the trunk net's inputs. Instead, the linear combination coefficients are functions of the initial state's properties, and x i are the branch net's inputs. Moreover, this last quantity is added to both trunk and branch inputs to account for the translational temperature dependencies. This paper trains two DeepONets independently for the inelastic and exchange reaction rates due to highly varying magnitudes and different physical governing principles of these two types of interactions. Unlike inelastic collisions, exchange processes are not bound by quantum transition rules, and overall, their rate values are approximately two orders of magnitude smaller than the inelastic ones. Along with energy transfer processes, an NN-based surrogate model is also constructed for the dissociation processes. As discovered in previous work by the authors , the StS dissociation rate coefficients show strong dependencies on the energy-based distance of the state from the centrifugal barrier (see Fig. ). Taking inspiration from this previous work and recognizing that the dissociation rate coefficients need to be mapped back only to initial states, we use a simple feed-forward neural network (FNN) as a surrogate for dissociating collisions. The novel methodology presented in this paper is tested on the kinetics of the O 2 +O collisional system between 1, 500 and 20, 000 K. The choice of the system is based on its relevance to hypersonic flight and the availability of potential energy surfaces for carrying out training and test QCT calculations. The developed methodology reduces by 85% the QCT computational cost for obtaining the StS reaction rate coefficients at the five training temperatures. Such a cost is evaluated by considering the overall number of trajectories simulated for training purposes. Additionally, the SurQCT surrogate can predict any rate coefficient for any temperature between 1,500 and 20,000 K without requiring any additional QCT calculations. This computational speedup is achieved while maintaining an accuracy of within 15% in the QoI, such as quasi-steady state (QSS) dissociation rates and energy relaxation times. This paper is organized as follows. Sec. 1 details on the surrogate models developed for the collision processes and explains the sampling algorithm implemented to select training data points. A discussion of the results is given in Sec. 2. We first assess the quality of the surrogate by comparing the rate coefficients that the model predicts with the reference values from QCT at training and test scenarios. Subsequently, we compare the master equation results produced by rate coefficients obtained from the reference QCT calculations and predicted by SurQCT. Finally, the work is summarized in Sec. 3.