Rescaling the complex network of low-temperature plasma chemistry through graph-theoretical analysis

We propose graph-theoretical analysis for extracting inherent information from complex plasma chemistry and devise a systematic way to rescale the network under the following key criteria: (1) maintain the scale-freeness and self-similarity in the network topology and (2) select the primary species considering its topological centrality. Network analysis of reaction sets clarifies that the scale-freeness emerging from a weak preferential mechanism reflects the uniqueness of plasma-induced chemistry. The effect of chemistry rescaling on the dynamics and chemistry of the He + O2 plasma is quantified through numerical simulations. The present chemical compression dramatically reduces the computational load, whereas the concentration profiles of reactive oxygen species (ROS) remain largely unchanged across a broad range of time, space and oxygen admixture fraction. The proposed analytical approach enables us to exploit the full potential of expansive chemical reaction data and would serve as a guideline for creating chemical reaction models.


Introduction
Low-temperature plasmas have been attracting attention because of their great potential in various fields such as surface treatment, material processing and bioplasma applications [1][2][3][4][5]. To understand the underlying operating principles of such plasma systems and optimise their performance in * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. applications, it is important to control plasma chemistry consisting of reactive multispecies. However, as a natural consequence, reaction networks to be treated become more complex when considering more realistic chemistry. Hence, efficient and robust chemical reaction models for studying fundamental chemical processes are worth developing.
Most numerical studies employ chemical reaction models to determine the rate of reaction. The reaction rate, which depends on the concentration of reactants present and the rate constant (rate coefficient), can be estimated numerically by solving differential equations (rate equations) for each reaction. Even if it is limited to fluid-approximation modellings, various scales/sizes of chemical reaction sets (in this paper, the scale/size of a reaction set corresponds to the number of species and the number of reaction processes considered) have been proposed to elucidate the low-temperature plasma chemistry. Zero-dimensional models (volume-averaged models/global models) [6][7][8][9][10][11][12][13][14][15][16][17][18][19], quasi-one-dimensional models (including 0D plug-flow approximation models), [14,[20][21][22][23][24][25][26][27][28][29] and one-dimensional models [30][31][32][33] are useful for simulating detailed chemistry at low computational cost. As done in those studies, one can increase the number of species and reaction processes for analysing the chemistry in lowdimensional modelling. For example, the models for humid air (N 2 + O 2 + H 2 O) plasma [21], Ar with humid air (Ar + N 2 + O 2 + H 2 O) plasma [22,23], He with humid air (He + N 2 + O 2 + CO 2 + H 2 O) plasma [25] and humid air (N 2 + O 2 + CO 2 + H 2 O) plasma contacting water [16] have processed dozens of species and from several hundreds to over one thousand reactions to describe plasmas containing contamination from humid environmental air. In reference [19], several hundreds reactions have been considered for a complex NF 3 + N 2 + O 2 + H 2 plasma. A series of works [13,14,17,29] has developed large reaction models, of which the number of reactions is expanded to over five thousand, to investigate the extensive carbon dioxide/methane chemistry. In parallel with the chemical expansion (the increase in scale/size of the reaction set), multidimensional modelling has been progressing. Two-dimensional models tend to use small reaction sets to balance chemical description with computational costs [34][35][36][37][38][39][40][41][42]. Whereas, some studies have employed large reaction sets consisting of several hundreds reactions even in the 2D framework, for instance, a He + N 2 + O 2 + H 2 O plasma model [43][44][45][46] and a NF 3 + N 2 + O 2 + H 2 plasma model [19]. On the other hand, three-dimensional models covering the core and afterglow regions focus on the macroscopic 3D structure of fluid dynamics and heat transfer phenomena rather than chemistry [47,48]. These multidimensional simulations are computationally intensive, so the scale of the reaction set should be properly examined.
The analysis and reduction of complex chemical mechanisms consisting of hundreds to thousands of species and thousands to tens of thousands of reactions has been extensively studied in the field of combustion science [49]. General methods to gain insight into a dynamical system are sensitivity analysis (SA) and uncertainty analysis (UA) [50][51][52][53][54][55][56][57]. They determine the response of the system to perturbations of individual parameters, so that it allows to identify which parameters contribute most to the solution variability. Another common procedure is an analysis of pathway (a series of reaction sequences) [58][59][60]. The pathway analysis (PWA) uses the list of the chemical reactions of the system, the concentration of all species and the rates of all reactions as input. From the output of the pathway analysis, the important pathways (reaction rates are greater than the user-specified threshold) and less important pathways (small rates) can be selected. A process which eliminates unimportant species and reactions from the detailed mechanism is often called skeletal reduction [61]. Typical skeletal methods include the above-mentioned PWA/SA/UA and a principal component analysis (PCA) [62,63]. In general, a problem in removing chemical species from the detailed reaction model is to distinguish the direct and indirect connections between species. The prediction of the concentration and reaction rate of a certain species strongly depends on the concentration of other species which does not appear in the same reaction equation. One of the techniques to overcome this issue is a directed relation graph method, which can quantitatively evaluate the bond between species and remove unnecessary species from the detailed model [49,64,65]. In this method, a graph theory, on which this paper focuses, is explicitly used.
A number of modelling studies in the low-temperature plasma research include PWA/SA/UA. The following are recent representative works, global models of plasmas with He + H 2 O [12], CO 2 [13,17] [45]. From a practical standpoint for computational studies, it should be useful to develop methods that can extract the primary components in complex chemistry. Based on simulation results obtained with a more comprehensive model, a simplified model can be created. A simple and direct method is to manually select a main species that is denser than the specified threshold under certain operating conditions [8,9,41]. More systematic screenings of a He + O 2 plasma [66,67], a nitrogen plasma [68] and a CO 2 plasma [68] have been performed to model the reasonable-scale reaction sets through the detailed PWA/SA/UA. Furthermore, chemistry reduction has been performed through the PCA for a CO 2 plasma [69] or an intrinsic low-dimensional manifold method (ILDM) [70] for an Ar plasma [71]. Figure 1(a) shows a conceptual illustration of a conventional process of the chemistry analysis/reduction. A number of chemistry analysis/reductions are post-processing or at least parallel processing of system simulations, as long as the output of the simulation is required as an input to determine the importance of the reaction or species. In contrast, as shown in figure 1(b), the present graph-theoretical approach is to examine the reaction model as a pre-process for system simulation. In the first place, the original reaction model is created through the composer's choice (based on experience/intuition). The question is how can we be confident that the reaction set is not too arbitrary or not too far from the inherent nature of plasmainduced chemistry. Therefore, we propose a network analysis that can quantify the topology and relationships of reaction models from a different perspective than most of the previous studies.
If the reaction network is complicated enough to form a weblike structure, we can identify which species play central roles in triggering subsequent reactions through a graph theory [72][73][74]. The graph-theoretical analysis is an analysis of relations using mathematical graphs [75][76][77][78][79] and has been applied to various complex networks in the fields of, for example, geography, economics and sociology, as well as physics, chemistry and biology [80][81][82][83][84][85]. The fundamental features of a complex chemical network appearing in lowtemperature plasmas have been examined in previous works [72][73][74]. The topological properties of methane plasma chemistry have been clarified [72]. It has been suggested that the graph-theory approach is applicable to gaining a rapid and approximate understanding and to the visualization and classification of complex plasma-enhanced chemistry [73]. Furthermore, in an interdisciplinary numerical modelling work, the network of intracellular functions, such as mitochondrial redox homeostasis and energy metabolism, activated by a lowtemperature plasma has been analysed [86]. However, the achieved knowledge has not been applied to the rescaling of plasma-chemistry networks. Besides, it is still elusive how effective the rescaling is with regard to computational plasma research.
The objective of this study is (1) to extract essential information from complex reaction data and devise a way to rescale the reaction network through the graph-theoretical analysis and (2) to quantify the effect of network rescaling on numerical simulations. Here, low-temperature plasmas using He + O 2 and He + humid air are examined. In the first half of this paper, the statistical features of the reaction set as a system and the topological centrality of each species are revealed. Particular emphasis is placed on the identification of similarity through network rescaling, by which a simplified model representing a principal feature can be extracted from the extensive model. In the second half, the numerical simulation on an atmospheric He + O 2 plasma jet is performed to gain confidence in the analytical prediction. The effects of network rescaling on the species concentration and the reaction probability are examined in detail.

Graph-theoretical analysis
In graph theory, the elements of the system and the interactions between these elements are represented as nodes and edges, respectively [75][76][77][78][79]. A directed graph (in which edges have orientations) is an ordered pair G = (V, E) where V is a set of nodes and E is a set of edges. The network structure of n nodes can be represented as an n × n adjacency matrix A.
In plasma chemistry, the nodes correspond to species and the edges (connections among species) represent reaction processes. Table 1 shows the graph-theoretical notation for a typical reaction in He + O 2 plasma as an example, that is,   To analyse a 65-species reaction set, for example, we create a 65 × 65 adjacency matrix to quantify the following topological measures and scale-freeness.
(a) Degree indicates the total number of edges that connect the node of interest with other nodes. As indicated in table 1 and figures 2(a) and (b), for the directed networks, the in-degree of a node is the number of incoming edges (the right-hand side of a reaction formula) and the outdegree is the number of outgoing edges (the left-hand side of a reaction formula). These notations represent the relationship between reactants (also called sources or agents) and products in chemical reactions. (b) The length of a path is the number of edges forming it.
The shortest path length (distance) between two nodes i and j is denoted by L(i, j). As indicated in figure  (c) The closeness centrality of a node is defined as the reciprocal of the average shortest path length. Nodes with high closeness centrality are close to most nodes and are therefore placed in the centre of the network. Therefore, the closeness centrality is a measure of how fast information spreads from a given node to other reachable nodes in the network.
Here, N is the total number of nodes, i and j are the node indicator numbers (i = j) and d i j is the number of edges (distance) from i to j along the shortest path (see figure 2(c)). (d) Betweenness centrality is the ratio of the number of shortest directed paths that pass through a given node to the number of all the directed shortest paths. The betweenness centrality of the node represents the controllability, the informativity or the importance of the role as an intermediate that directly or indirectly bridges the reactants and products.
l indicates a node and j = l = i, s jl is the total number of shortest paths between j and l and s jl (i) is the case of passing through i (see figure 2(c)). (e) The clustering coefficient is used to determine the extent to which the nodes of a network are closely connected with one another (cohesion of the network).
Here, g i is the number of neighbours of i and b i is the number of connected pairs between all neighbours of i. The clustering coefficient is the ratio [the number of edges between the neighbours of i]/[the maximum number of edges that could possibly exist between the neighbours of i].

Scale-freeness.
As shown in figure 3, counting the number of nodes for each degree yields the degree distribution, P(k), which can be defined as the fraction (probability) of nodes in the graph with degree k. The network can be called random when all edges are independent, i.e., the edges cannot communicate with one another. In this case, the degree --

Charged species
He + , He Electrons Electrons distribution forms a normal (Gaussian) distribution because the independent edges spread themselves fairly evenly across the nodes (figure 3(a)). Another common feature of real-world networks is the presence of hubs, i.e., a few nodes are highly connected to other nodes. The presence of hubs will give a degree distribution with a long tail following a power law (figure 3(b)), which is expressed by where k is the number of degrees and γ is the exponent parameter. On the left side of figure 3(b), we can see that most species have a relatively small degree (frequent events). In contrast, on the far right side of the graph, only a few species have a very high degree (rare but important event). The network showing the power-law distribution is formed as a result of a preferential attachment, i.e., a new node is preferentially attached to nodes with high degrees in the generative process [95][96][97][98][99][100].
Since power laws have the same functional form at all scales (the scale invariance of power law functions), the degree distribution remains unchanged when rescaling the independent variable k. Therefore, the network is called 'scale-free' (a lack of characteristic scale). Internet networks, sciencecollaboration/citation networks, bio-chemical/cellular/neural networks show the scale-freeness in real world, for instance [76]. In general, the scale-free network is believed to show robustness against perturbations.

Chemistry analysis
We examine three chemistry xmodels: (1) He + humid air plasma reaction set consisting of 65 species and 1360 reaction processes [25]; (2) He + O 2 plasma reaction set consisting of 20 species and 327 reaction processes (referred to as the 20-species model); and (3) He + O 2 plasma reaction set consisting of 12 species and 104 reaction processes (the 12species model). Table 2 shows a list of species considered in the three models. Figure 1(b) shows a conceptual illustration of the reaction space. The extensive set of chemical reactions has been compiled from the latest available expansive data to develop the present models for He-based plasmas. The scale of chemical networks varies widely; the numbers of species and reactions change from 12 to 65 (a factor of 5.4) and from 104 to 1360 (a factor of 13), respectively.

Chemistry reduction in
He + O 2 plasmas. One important objective of this study is to extract a low-dimensional (small-scale) reaction set containing the primary (principle) component while maintaining the statistical feature of a highdimensional (large-scale) reaction set. We conduct the speciesoriented reduction from the 20-species model to the 12-species model in accordance with the following principles of species selection.
• Centrality indices are taken into account.
• Scale-free topology within the network is maintained.
• The 12-species and 20-species models are implemented in the numerical simulation.

Chemistry extension from He + O 2 plasma to He +
humid air plasma. The He + humid air plasma reaction set used here has been developed to examine the influence of the humid air constitution (nitrogen, water and carbon-related chemical species) on the He + O 2 plasma chemistry [25]. In this paper, comparing the He + humid air plasma reaction model and the He + O 2 plasma reaction model allows us to examine how the essential feature of the network structure is maintained through the species-oriented dimension transform.
Note that our scheme proposed here depends only on the network structure of the reaction set and not on the species concentrations or reaction rates as a result of the system simulation.

Numerical simulation 2.3.1. Numerical domain and calculation condition.
A kHzpower-driven atmospheric-pressure dielectric barrier discharge jet (a kHz jet) is examined as one of the typical example of a low-temperature plasma source [101,102]. Figure 4 shows (a) the schematics of the kHz jet with the electric field streamlines induced and (b) the one-dimensional calculation domain.
The plasma properties vary in the two geometrical components along with the electromagnetic streamlines, the tube longitudinal direction between the powered and grounded electrodes, and the tube radial direction from the centre axis to the electrode covered by the dielectric material. This provides quasi-two-dimensional information of the plasma dynamics and chemistry in the bulk region (parallel to the longitudinal centre axis of the tube) as well as in the region near the wall-electrode boundary (perpendicular to the tube wall). It is assumed that the numerical domain is along the broken red

Governing equations, boundary conditions and numeri-
cal procedures. A fluid flow approximation of Boltzmann's equation is employed to describe the present weakly ionized, collision-dominated plasmas. The nonequilibrium plasma established here is characterized by two temperatures, the electron temperature (T e ) and the heavy-particle gas temperature (T g ). We neglect the momentum and energy conservation equations in the heavy-particle system to focus mainly on the electron system. The governing equations can be reduced to the set of conservation equations of density and energy for electrons with a drift-diffusion approximation [103], and the continuity equation of the mass fraction of each species for a heavy-particle system. The rate equation is combined with the above equations. Electrical properties are determined by solving Poisson's equation.  Conservation of electron number density: where n e is the electron number density, Γ e is the electron flux vector, S e is the electron source term and t denotes time. The flux term is expressed based on the drift-diffusion approximation [103][104][105] as, Γ e = −μ e En e − D e ∇n e , where μ e is the electron mobility, E is the electric field, D e is the electron diffusivity. The electron diffusivity (diffusion coefficient) is defined as D e = (k B T e )/(m e ν m ), where k B is Boltzmann's constant, m e is the electron mass and ν m is the momentum transfer frequency [106]. The source term is defined as S e = j c j R j , where c j is the stoichiometric number for electrons in the reaction j, R j is the reaction rate. The boundary condition on the wall (dielectric surface) for the electron flux is −n · Γ e  2 ν e,th n e − p γ p (Γ p · n), where n represents the outward normal, ν e,th is the electron thermal velocity, γ p is the secondary emission coefficient from the pth positive ion species and Γ p is the ion flux [107].
Conservation of electron energy: where n is the electron energy density, q is the unit (electron) charge, Γ is the flux of the electron energy density, the third term on the left-hand side represents the ohmic heating and S is the source term of the electron energy density. The electron energy density is defined as n = n e¯ , where¯ is the mean electron energy. The flux of the electron energy density is written in a drift-diffusion form as Γ = −n μ E − D ∇n , where μ is the energy mobility and D is the energy diffusivity (diffusion coefficient for electron energy density) [103]. The source term is defined as S e = j c j R j Δ j , where c j is the stoichiometric number, R j is the reaction rate and Δ j is the energy loss/gain due to the reaction j. For the electron energy flux, the surface boundary condition is −n · Γ = 5 6 ν e,th n − p p γ p (Γ p · n), where p is the mean energy of the secondary electrons [107].
Mass balance of species k: where ρ, w k , V k and S k are the density of the mixture, the mass fraction, the multicomponent diffusion velocity and the source term of the kth species, respectively [108][109][110][111]. The mixtureaveraged diffusive flux vector ρW k V k is computed with , M n is the mean molar mass of mixture M n = ( k (w k /M k )) −1 , D k j is the binary diffusion coefficient between species k and j, z k is the charge number and μ m k is the mixture-averaged mobility for the kth species, which is derived using Einstein's relation [112,113], μ m k = (qD m k )/(k B T g ). For the heavy species, ions are lost to the surface and excited neutral particles are quenched owing to surface reactions. As a boundary condition, the heavy particle flux normal to the surface is written as −n · j k = M k R k + M k c k z k μ m k (E · n), where R k is the reaction rate on the dielectric surface, c k is the mass fraction and the second term on the right-hand side is added for ions. The surface reaction scheme is provided as supplementary material.
The electric potential φ is computed using Poisson's equation: where ρ q is the space charge density and r = / 0 , where and 0 are the permittivities of the medium and vacuum, respectively. The surface charge accumulation on the dielectric surfaces is taken into account using the boundary conditions, n · (D 1 − D 2 ) = ρ s and dρ s /dt = n · J ion + n · J e , where D 1 and D 2 are the displacement electric fields on the both sides of the boundary, ρ s is the surface charge density, J ion and J e are the ion and electron current densities on the surface, respectively. The electrical potential is given as φ = 0 on the cathode boundary and φ = V apply on the anode boundary.
The equations are discretized by the Galerkin finite-element method which is one of the most popular methods of weighted residuals [114,115] and the resulting simultaneous linear equations having a sparse coefficient matrix is solved by the parallel direct sparse solver (PARDISO) [116]. The implicit backward differentiation formula time stepping method with order of accuracy varying from one to two is used as a timedependent solver [117]. The simulated portions of plasma and dielectric regions are divided by 185 and 16 elements, respectively. The numerical simulation is conducted using the finite element partial differential equation solver Comsol multiphysics [118].  The effect of the chemical reduction on the degree centrality is examined by comparing the above-described 20-species model and the reduced 12-species model [figures 5(g)-(i)]. As the number of reaction processes decreases from 327 to 104 (about 70% reduction), the absolute value of the degree of each species decreases. Although the reaction set is scaled down, the degree values of O 2 , He, electrons and O( 3 P) are relatively high (351-243). The sum of the degrees of these four species accounts for 65% of the total. The scatter plot shows that the degree of the 12-species model is distributed with a similar tendency to that of the 20-species model; the in-degree broadly changes from 5 to 215, whereas the out-degree is in a narrow range . Note that the species are selected not in order of degree but across all levels of degree to retain the scale-free topology, as will be discussed later in this section.

Results and discussion
The effect of the chemical extension by adding nitrogen, water and carbon species is shown in figures 5(a)-(c). Figures 5(a) and (b) figure 5(a)]. The sum of the degree counts for these six species is 8575, which corresponds to about 47% of the sum total degree of all species (18 098 degrees). Other species such as NO, N( 4 S), O − 2 , O 3 , CO 2 , O − and OH have degree counts from 520 to 320. Figure 5(c) shows that the species in-degree counts are more broadly spread from unity to 10 3 than their out-degree counts. The O 2 , He, N 2 and electrons have high counts of both in-and out-degrees. The comparison among the three kinetic models revealed that the reaction networks have a few key species, such as O 2 , He, O( 3 P) and electrons. They show particularly high centrality in both in-and out-degrees, indicating that they are influential in the network because they serve as both reactants and products. In contrast, some species with low in-degree counts mainly play the role of reactants rather than products.

Scale-freeness and self-similarity.
From the statistical perspective, figure 6 shows the degree distributions of (a) He + humid air plasma, (b) He + O 2 plasma, 20-species model and (c) He + O 2 plasma, 12-species model. Horizontal bars added to the dots indicate the range of degree sampling. The profiles of power-law formulae, P(k) ∝ k −2 and P(k) ∝ k −3 , are plotted. Figure 6 indicates that each degree probability profiles can be expressed as (a) P(k) ∝ k −1.1 , (b) P(k) ∝ k −1.2 and (c) P(k) ∝ k −1.0 . The number of degrees k originating from a particular species shows a power-law distribution in all cases, suggesting scale-freeness of the network topology. Similar to the general feature found in real-world networks, the plasmainduced reaction network developed here has a small number of nodes (species) that are highly connected to other nodes (species). As previously discussed, He, O 2 , O( 3 P) and electrons are the most relevant hub species, as determined using the three reaction models. Helium atoms are the feed gas atoms and are the most abundant. Oxygen molecules are seeded for tailoring the target ROS. Ground-state atomic oxygen is the most primary ROS. Electrons are the main triggers for ionization, excitation or dissociation.
Importantly, the exponents of the three formulas, 1.1 ± 0.1, remain unchanged against rescaling. This finding reveals the self-similarity in the topology of the plasma-induced reaction network. The well-known theoretical Barabási-Albert model (BA model) provides a systematic way of generating a scalefree network following a power law of P(k) ∝ k −3 [76]. A principal rule of this model in network growth (one new node joins the other nodes) is preferential attachment, i.e., a new node tends to attach to an existing hub node and the attachment takes place one by one. The strong preferential attachment through the one-by-one rule generates very few huge hubs, which makes the exponent value large. The present exponent value of 1.1 ± 0.1 shown in figure 6 is lower than that of 3 in the BA model. This implies that the generative process of the preference topology in the plasma chemistry is unique. The nodes of the BA model do not experience any 'chemistry' during the growth. When a chemical reaction takes place through the network generation process, such as fragmentation, bonding or polymerization, it renders the preferential growth weak. The marked difference of the plasma-induced chemical network from other networks is that it involves parallel progression of a variety of excitation, ionization, attachment and dissociation processes of reactants as well as products. The rich chemistry induced by high-energy electrons leads to simultaneous attachments of many species onto the existing network configuration rather than the simple one-by-one attachment. Therefore, more species have chances to get edges from the new species, which in turn lessens the preferential attachment feature, i.e., the exponent of the power-law formula decreases.

Closeness, betweenness and clustering.
We examine here the importance of the species from a microscopic perspective. Figure 7 shows   In this group, the hydroxyl radical OH, ground-state atomic nitrogen N( 4 S) and NO 2 have relatively high betweenness centrality (higher than 10 −2 ). The species with low betweenness centrality of less than 10 −3 , for instance, the hydrate cluster ions, i.e., O + 2 · H 2 O, NO − 2 · H 2 O and NO − 3 · H 2 O, are located in the outermost regime [figure 7(a)] because these clusters are end-products rather than intermediates. The similar distributions in betweenness-closeness centrality plots in figures 6(b), (e), (h) are also observed in SiH 4 and CH 4 plasma chemistry, as indicated in our previous study [73]. This finding suggests that such trends are in a more generalized framework in molecular plasma chemistry. Figures 7(c), (f ), (i) show the clustering coefficient distribution of the three networks. The clustering coefficient vanishes with k as C clustering ∝ k −0.8±0. 3 . The decline of the clustering coefficient at a higher degree implies the presence of a few primary hubs having a large number of connections. The hub species show broad connectivity and affect the entire net-work. In contrast, species with higher clustering coefficients are found in the lower k regime. Such cluster species strongly connect with neighbours and localise themselves. The powerlaw profile of the clustering distribution is interpreted as an expression of scale-freeness. Figure 8 summarizes the macroscopic relationship among the clustering coefficient, closeness centrality and betweenness centrality of each species. Along with the direction indicated by an α-arrow, the minimum closeness centrality increases from 0.4 to 0.8 with decreasing system size, because the average shortest path length decreases and the network becomes denser upon excluding the species located near the outer edge of the system. When the species betweenness centrality increases from 10 −4 to 10 −2 , the clustering coefficient decreases without apparent change in the closeness centrality [β-arrow]. Species with higher betweenness centrality interact with a wider range of other species and more intermediate reactions rather than communicating tightly with only their neighbours. The species with further higher betweenness centrality (higher than 0.1) also achieve markedly high closeness centrality [γ-arrow]. These central hub species govern the entire kinetic system.   Figure 11 shows the bulk ROS densities calculated with the 20-species model (solid curves) and the 12-species model (dashed curves) as a function of O 2 /(He + O 2 ). The O( 3 P) and O 2 (a 1 D) densities monotonically increase from 10 20 m −3 to 10 21 m −3 with increasing O 2 fraction from 0.01% to 0.1%, where the O 3 density is relatively low. The highest O( 3 P) density of about 10 21 m −3 is obtained under the O 2 fraction of 0.1%-0.5%, at which the O( 3 P) density decreases from its peak. The O 2 fraction giving the peak O 2 (a 1 D) density shifts to 1%-2%. For further O 2 admixture, the O 3 density is comparable to that of O 2 (a 1 D) or becomes more pronounced to 10 22 m −3 at the O 2 fraction of 0.5%. The O 2 admixture dependence of ROS tailoring is common to both reaction models, of which the difference in density is a maximum of a factor of two. Figure 12 shows the reaction rate of the most pronounced production and destruction processes of ROS as a function of O 2 /(He + O 2 ). In the present He + O 2 environment, the O( 3 P) production is mainly driven by electron-induced processes,

Numerical simulation
and electronically excited heavy-particle-induced processes, The dominant O( 3 P) destruction processes are a three-body reaction including background He and O 2 and a collision with O 3 in dissociative excitation of O 2 to produce O 2 (a 1 D).
The major O 2 (a 1 D) production processes are reaction (19) and through direct electron impact excitation of O 2 , The main destruction process of O 2 (a 1 D) is quenching through collisions with O( 3 P) and O 3 [reaction (17)].
Of the two-body reaction rate constants used in the numerical simulations, those with values higher than 10 −15 m 3 s −1 are shown in decreasing order in figure 13. The reaction constants contained in the two models are broadly distributed from 10 −12 to 10 −15 m 3 s −1 in a similar manner. After comparing the rate constants included in the two models, it is perhaps not surprising that those species concentrations [figures 9, 10 and 11] and reaction rates [ figure 12] agree well across a wide variety of species, time, space and oxygen admixture fraction. It is worth emphasizing that the similarity in reaction probability is retained through the proposed network rescaling by the species selection across all levels of topological centrality, not by a simple choice from the top.
Under the typical numerical conditions set herein, the calculation time taken for a simulation has been markedly shortened, for instance, from 9 h using the 20-species model to 90 min with the 12-species model, corresponding to an approximately 85% reduction in running time. Consequently, the computational costs can be significantly reduced while maintaining accuracy in plasma-induced chemistry by means of the present complex network analysis.

Conclusions
We explored the complex chemical network in lowtemperature plasmas using the graph-theoretical analysis to extract its inherent information and devise a systematic way for its rescaling.
Network analysis from a statistical perspective has clarified that the scale-freeness emerging from a weak preferential mechanism in the network generative process reflected the richness in the plasma-induced chemistry sustained by the variety of excitation, ionization, attachment and dissociation processes. The identification of the species topological centrality revealed that the system was governed and mediated by a few hub species, such as electrons, the main discharge driver; He, O 2 , N 2 and H 2 O, the background species acting as both reactants and products; and O( 3 P), the highly reactive primary excited species, in He + humid air plasmas. Thus, the analytical way of rescaling the chemical model maintaining the common scale-freeness and self-similarity in network topology was devised. Furthermore, we quantified the effect of chemistry rescaling on the dynamics and chemistry of the He + O 2 plasma through one-dimensional numerical simulations. Implementing a simplified reaction model may not guarantee exact reproducibility, but it provided the plasma characteristics similar to those obtained with the extended model. The concentration profiles across a broad range of species, the temporal and spatial evolutions of ROS and their O 2 admixture dependence remained largely unchanged, because the similarity of the reaction probability was maintained through the primary components selection. The present chemistry compression dramatically reduced the computational time.
The proposal here, i.e., 'take centrality as an index of species importance', is an idea depending only on the network structure, not directly on the physico chemical properties of species. Therefore, the centrality indices of each species might not necessarily match the importance expected from a variety of viewpoints of discharge sciences or possible applications. In addition, this work could not cover all possible species or reaction processes. Nevertheless, we believe that the present range of network rescaling is wide enough to clarify the inherent nature of plasma-induced chemistry. The proposed analytical approach provides valuable insights into complex systems and enables us to exploit the full potential of expansive chemical reaction data, hence it serves as a guideline for creating chemical reaction models. Combining this analytical approach with existing state-of-the-art rescaling techniques such as PWA/SA/UA/PCA/ILDM is expected to lead to more efficient and robust chemical systems.