Causal emergent principles and relations for mechanical properties of covalent and ionic crystals

Material chemistry concerns composition–structure–property triangle relationships. As shown in Fig. 1, three kinds of structure–property, composition–structure, and composition–property analytic relationships occur naturally. Among them, the structure–property relationship has been intensively and extensively investigated with experimental and theoretical methods. In addition to heavy experimental syntheses and characterizations, lattice constants, the volume of the unit cell (V₀), bond length (d), bulk modulus (B₀), phonon spectra, bandgap, and other physical properties are able to be obtained via first-principle calculations from the atomic number, mass, and size of constituent elements as input,^1,2 which requires more complex software and significant computation power. Besides ab initio total-energy approaches, empirical models take advantage of their applicability to a broad class of materials and concise and useful trends for material designing.^3–8 As exemplified by models No. 1–6 in Table I, some analytic relations between mechanical properties and structural parameters were proposed for covalent and ionic crystals. It is seen from Table I that empirical models connect mechanical properties with some lattice structural features and some elemental features, which exhibit various dependencies on the dimension of length or scale of bond length and should be clarified for these discrepancies. Nonetheless, quantitative structure–property relations remain little known for material chemistry due to some structural features not being digitized and most material data under operation conditions being missing. Furthermore, it is widely accepted that targeted performances of material are regulated in practice predominantly by changing compositions and secondly by controlling microstructures like refining crystallite sizes and so on. The composition–structure and composition–property relations set up in nature are a substantial cornerstone to regulate those desired properties for engineering applications as well as to discover new materials.^9–11 Since the features of composition, including chemical formulas and chemical bonds, are not yet digitized, both analytic composition–structure and composition–property relations wait for unveiling, while experimental observations and theoretical calculations just get a simple mapping relationship for materials.

Besides the error-and-trial and screening methods, material design in a predictive way is also developing to speed up the development and deployment of new materials.^11–15 Physical prediction models for designing materials require at least three key elements:¹⁴ (1) physical and chemical principles, (2) causal analytic relations on changes, and (3) desired targets for every application. Schrödinger equations accompanying the symmetry principle describe the dynamic behaviors of atoms and atom close packing systems successfully, while the thermodynamic equation of Gibbs free energy describes phase formability, structural phase transformation, and chemical reaction’s direction and extent, which set up a substantial physical and chemical foundation for materials. The desired targets are determined by engineering application or technological development trends while discovering analytic relations between targeted performance and chemical composition/microstructure (object system), processing parameters (manufacturing system), and operation conditions (engineering system), which thus play a key role in establishing physical prediction models. Wherein, digital feature descriptors are substantial not only for the interpretability of physical models but also for the accuracy and efficiency of predictions. Currently, AI technology represented by machine learning is being introduced into the material community to undertake tasks such as selecting features, discovering relations, and solving equations for accelerating the discovery of materials.^15–18 In recent years, machine learning has also been tried to find the hidden relations behind material data (AI for materials)^16,17 or to find solutions to partial differential equations (AI for science).¹⁸ Different from traditional mathematical models, machine learning provides a data-centric modeling method that constructs models by learning representations or rules from data and works with statistical information rather than structural rules about how systems work.¹³ Most of the impressive progress obtained by deep learning was built on big data for small tasks, while AI for materials encountered the challenge of small data for multi-objective tasks.¹⁹ Machine learning is thus inefficient at tackling the dynamic and emergent behaviors of complex and large systems like materials. For the utmost goal, causal relationships are most substantial in both nature and data science to understand how a system works, why it works in this way, and what happens when it is disturbed externally.¹³ Therefore, when facing the reality of materials being small data, creating knowledge and data-synergized models and algorithms is alternatively promising and unavoidable to bridge the gap between data-driven models and mechanistic models, and developing causally intelligent computations could balance the efficiency of statistical modeling with the conciseness of physical models.^{13,14,20–23}

Material research with traditional paradigms of experimental observation, theoretical modeling, and computational simulation provides data, models, and theories, while the data-science paradigm provides new insights and methods for data analysis and data generation.²⁴ To unveil analytic relations hidden behind material data, digitizing chemical composition as feature descriptors plays a vital role in understanding the performance of those known materials as well as in developing new materials. Reductionism usually treats the features of chemical elements and lattice structure as descriptors, which inevitably result in high-dimensionally analytic relations, while emergentism says that the features of atoms (ions) and electrons cannot directly describe the collective behaviors and emergent properties of the atom close packing system.^25,26 Therefore, it is necessary to construct digital ensemble descriptors from the features of chemical elements for crystalline materials, which naturally reduces the dimension for the next regression operation.^14,20–23 In accordance with system theory and the quantum mechanical Schrödinger equation, as shown in Fig. 2, the inputs for designing crystalline materials are the feature data of chemical elements contained in the Periodic Table, element combination, and concentration, of which the choices of element and concentration are the content of inverse designing for material; the outputs are both structural/thermodynamic state parameters and physical/chemical/biological response parameters, of which their comprehensive combination varies with respect to engineering application and are the content of forward prediction. In contrast to correlative relations among the state and response parameters, causal analytic relations of the state/response parameters with respect to digital ensemble descriptors are much more substantial, on the basis of which an intelligent prediction model is physical and feasible with respect to interpretability, generalizeability, and accuracy. Through unveiling causal emergence and regulation principles behind every target property, physical prediction models and intelligent generation algorithms will be established consequently by combining ensemble descriptors and causal analytic relations. Only then is AI for material science physically interpretable: (1) forward prediction models to discover those governing rules, and (2) inverse design models capturing the physical and chemical world of materials accurately.^13,18

It is known from condensed matter physics that macroscopic physical properties are emergent phenomena of the atom close packing system, so they are not able to be described directly using the features of the atom such as atomic number, weight, radius, charge, spin, and so on.²⁶ In crystalline solids, chemical composition and atom geometric arrangement in a unit cell are the features of matter dimension. Herein, two digital ensemble descriptors are creatively constructed for covalent and ionic solids with the features of elements including atomic weight, covalent radius, effective nuclear charge, the number of atoms in a unit cell as^14,20–23 $μ× r 1 r 2 / a B$ (μ is the reduced mass of primitive cell,^20,21 r_i covalent radius of constituent atoms,²⁷ and a_B Bohr radius), and $Z 1 * Z 2 * / Z H *$ (⁠$Z i *$ and $Z H *$ are effective nuclear charge of constituent atoms and hydrogen,²⁸ respectively), in which a_B and $Z H *$ are included here just for dimensionless convenience. Adopting geometric averaged covalent radii and effective nuclear charge is due to two sublattices in these diamond-, zinc-blende-, and rocksalt-type crystallographic structures. As illustrated in Fig. 2, digitizing chemical composition and structural trends as ensemble descriptors makes it possible to establish causal quantitative analytic relations with structural and thermodynamic state/physical response parameters.

So far, it has been seen that, through adopting the data science paradigm and hackling the connotation and denotation of material data, the relationships behind material data are discriminated into two kinds: causality and correlation. With digital ensemble descriptors constructed by the author, causal analytic relations between chemical composition and state/response parameters could be found via small-data mining, and thereafter, causal emergence and regulation principles for forward prediction and inverse design of crystalline materials could be unveiled. In this section, some covalent and ionic crystals were exemplified to illustrate this architecture of intelligent computation for materials, and a physical prediction model for superhard materials was finally illustrated on the basis of causal analytic relations.

Within the experimental observation paradigm, symmetry was treated as a passive result of the interaction of the atom close packing system. On the contrary, symmetry plays an active role in dictating the interaction of atom close packing systems within the theoretical modeling paradigm.^29,30 As illustrated in Figs. 1 and 2, choosing chemical elements and concentration is active in designing materials, while the stable configuration of atom geometric arrangement is a resulted state with a minimum free energy of the atom close packing system dictated by lattice symmetry under a certain thermodynamic boundary condition. In consequence, lattice structure, electronic structure, thermodynamic state, and physical properties are inherently embodying information of symmetry and correlative each other, which vary with changing constituent elements and thermodynamic conditions and reflect various aspects of state/response to external stimuli. Figure 3 illustrates such correlations among bulk moduli (B₀), hardness (H_v), the volume of the unit cell (V₀), and bond length (d) for diamond-, zinc-blende-, and rocksalt-structured covalent and ionic solids with various ionicities in a range of 0–1.00.³¹ Through a least-square regression, a linear relation between B₀ ∝ d⁻³ ∝ 1/V₀ and H_v ∝ d⁻³ ∝ 1/V₀ was obtained, and those proportional coefficients are summarized in Table II for various symmetric crystals.

On bulk modulus (B₀), the proportion coefficient (⁠$C B d fit$⁠) as a function of d⁻³ was obtained as 1824, 1655, 1009, and 635 for diamond-structured solids of group IV, zinc blende solids of group III-V, II-VI, and rocksalt solids of group I-VII, respectively. As illustrated by Fig. 3, this classification is similar to Cohen’s observation (model No. 1 in Table I) of the proportional coefficient determined by the combination of crystallographic structure and constituent elements.^1,2 However, in contrast to Cohen’s empirical model with different scales for tetrahedral- and octahedral-coordinated solids, a universal scale law of B₀ = (1912 − 421λ)/d³ occurs to depict the quantitative bulk modulus of covalent and ionic solids with various ionicities in a range of 0–1.00, of which the proportion coefficient has a simple dimension of energy and λ = 0, 1, 2, and 3 for groups IV, III-V, II-VI, and I-VII, respectively. At the same time, the proportion coefficient (⁠$C B V fit$⁠) as a function of 1000/V₀ was obtained as 21.3, 20.4, 10.1, and 5.24 for groups IV, III-V, II-VI, and I-VII, respectively. This observation is very different from the inference of the bond charge model (model No. 2 in Table I), which was simultaneously associated with bandgap and bond charge in addition to the volume of the unit cell as variables. In a well-accepted phenomenological model (model No. 3 in Table I), following the assumption of cohesion energy, the bulk modulus is derived as B₀ ∝ 1/V₀ on thermodynamic equilibrium, of which the proportion coefficient is determined by the constants of m, n, A, and B. Although B₀ is often interpreted as quantitatively dependent on volumetric energy density connected with compression,⁵ this analytic relation implies those trends observed in Fig. 3 when the constants of m, n, A, and B are kept for every class of crystals marked by lattice symmetry and constituent elements as showed in Table II. In each crystallographic lattice, atoms occupy the same space in the same way, so the Coulomb interaction between atoms and relevant cohesion energy is thus determined by constituent elements. From both crystallographic analyses and data-mining, V₀ is found to be linearly proportional to d³ for zinc blende, rocksalt, cesium chloride, corundum (A₂O₃), and rutile (AO₂)-structured crystals.³⁵ Therefore, B₀ ∝ 1/V₀ ∝ d⁻³ occurs universally for covalent and ionic solids with various crystallographic structures.

On hardness (H_v), the proportion coefficient (⁠$C H d fit$⁠) as a function of d⁻³ was obtained as 382, 282, 31.7, and 9.65, while the proportion coefficient (⁠$C H V fit$⁠) as a function of 1000/V₀ was 4.43, 3.48, 0.398, and 0.080 for groups IV, III-V, II-VI, and I-VII, respectively. In comparison with empirical models No. 4–6 listed in Table I, a linear relation of H_v ∝ 1/d³ ∝ 1/V₀ is much more concise to depict the structure–property relation for covalent and ionic crystals. On one hand, those empirical analytic relations shown in Table I contain some parameters on both atomic (e.g., Z_i, R_i) and unit cell (e.g., V₀, d_ij, the coordination number of each atom, and the number of atoms in a unit cell) hierarchical levels.^6–8 In quantum mechanics, the parameters of a unit cell describe quantitatively the state of a solid and are determined by the features of the constituent atoms. On the other hand, each parameter at the unit cell level is a function of those at the atomic level in mathematics. For instance, crystallographic analyses have established well-established analytic relations between lattice constants and volume of unit cell (bond length) for each Bravais lattice, e.g., V₀ = a³ = (2d)³ and $V 0 = a 3 = ( 4 d / 3 ) 3$ (The format of equation should be modified for it looks like subscript.), obtained theoretically for rocksalt and diamond structures, respectively.³¹ Through data-mining extant experimental data, as shown in Fig. 3, V₀ = 8.4d³ and V₀ = 13.2d³ were obtained for rocksalt and diamond structures, respectively. So far, it is seen from Fig. 3 and Table II that there is a simple and concise analytic relation between each couple among structural parameters and mechanical properties, which is not like the complex empirical relations shown in Table I.

Although hardness is different from bulk modulus, they are considered correlative to each other for covalent and ionic crystals because both compressibility and hardness are mechanical responses of the atom close packing system to pressure. Till now, the nature of both of them has been widely understood microscopically on the basis of valence electron behaviors (bond length, charge density, ionicity, strength, and others) and chemical bonding.^2,6,7,36,37 As illustrated by Fig. 3, a linear correlation of H_v ∝ B₀ with a proportion coefficient (⁠$C H B fit$⁠) of 0.208, 0.172, 0.019, and 0.015 was obtained for groups IV, III-V, II-VI, and I-VII, respectively. In a similar way, a linear relation of H_k ∝ B₀ with a proportion coefficient of 0.228, 0.133, and 0.028 was also obtained for diamond-structured solids of groups IV, III-V, and II-VI, respectively. Hardness is a complex property related to the extent to which solids resist both elastic and plastic deformations, of which the value is defined experimentally by pressing an indenter into the surface of the material and measuring the size of the impression.⁷ Knoop’s scale uses a sharper diamond wedge than Vicker’s one. As shown by $C H B fit$ in Table II, the ratio of H_k to H_v is varied, not like the expectation of H_k being lower than H_v, which may need to take the microstructure and impurity effects into account further. So far, it has been seen that a linear correlation occurs between hardness and bulk modulus, which is not complex like those empirical models shown in Table I. In Griffith’s crack theory, both crack strength and crack toughness are determined predominantly by the material’s parameters of elastic modulus and surface energy,¹⁴ which also supports the fact that a low compressibility with a large bulk modulus is promising for hard materials.

At the moment, it should be noted that more correlative analytic relations between structure and property, between different state parameters, and between different properties are still waiting for data-mining, as illustrated in Fig. 3. However, data-mining analytic relations between composition and structure/property should be carried out priorly because they provide causal emergence and regulation principles for material science and engineering.

In Fig. 4, the room temperature volume of the unit cell (V₀), bond length (d), bulk modulus (B₀), hardness (H_v), and entropy (S⁰) were plotted with respect to $μ× r 1 r 2 / a B$ or $Z 1 * Z 2 * / Z H *$ for diamond-, zinc-blende-, and rocksalt-structured solids. Through data-mining with least squares regression, it was found an exponential relation of y = A exp(−x/B) + C (A, B, and C are constants) that depicts quantitatively well the trend of structural parameters (1/d³ and 1/V₀) and mechanical response parameters (B₀ and H_v) with respect to $μ× r 1 r 2 / a B$⁠, of which the constants (A, B, and C) were classified into three classes: diamond/zinc-blende-type solids of group IV, III-V, and II-VI compounds, rocksalt-type solids of I-VI excepting Li halides, and rocksalt-type solids of Li halides. In comparison with Fig. 3, it is noted that, for diamondlike solids of groups IV, III-V, and II-VI, the causal relations of composition–structural/mechanical parameters belong to one class, while the slope of the linear relation of structure–mechanical properties varies with changing the group of constituent elements, of which the latter incorporates the microstructure effect of solids. On rocksalt-type halides, such two classifications can be reasonably attributed to the feature of atom size for the covalent radius ratio of Li/Na (0.69) being much smaller than that of Na/K (0.84), K/Rb (0.90), and Rb/Cs (0.88). On thermodynamic parameter S⁰, rocksalt-type halides keep the same classification as structural and mechanical parameters, while diamond-type solids (group IV) and zinc-blende-type solids (groups III-V and II-VI) are separated into two classes, of which the latter is attributed to the factor of two kinds of constituent elements increasing the number of possible microscopic arrangements (configurational entropy). With respect to $Z 1 * Z 2 * / Z H *$⁠, both 1/V₀ and H_v do not exhibit a monotonic trend with a maximum response at the position of the diamond, on both sides of which an exponential trend for diamond/zinc-blende-type solids was observed. So far, the digital ensemble descriptors constructed here are seen as chemically intuitive and physically meaningful to extract major chemical compositions and structural trends of crystalline solids, which affords the mechanism of those learned relations between state/response parameters and ensemble descriptors.

In chemistry, the Periodic Table brings ordering into the chemical and physical behaviors of elements, of which the period and group are classified in accordance with an electronic configuration that is dictated by quantum mechanics and the rotation symmetry of the Coulomb force.^14,29,40 In quantum mechanics for solids, it is popularly assumed that those atomic core electrons are unchanged when a solid is formed, which reduces the problem to a description of a model solid in which a sea of valence electrons interacts with a periodic array of positive cores.⁴¹ Thereafter, the problem of predicting the existence of crystalline solids is reduced, in principle, to checking the total energy of a system of valence electrons interacting with different structural arrays of cores and with each other. Within the paradigm of computation simulation, the total energy of the atom close packing system is dictated by the lattice symmetry of the point group, i.e., by the geometric arrangement of the atom. Following the above-mentioned ideas, it is observed from Figs. 3 and 4 that the features of period and group of constituent elements and geometric arrangements of atoms combinatorially provide natural criteria for the classification of crystalline solids, which means those crystalline materials formed by constituent elements of the same group with the same geometric arrangement configuration exhibiting a similar behavior except in quantity for their state/response parameters.

Chemical bonding is an alternative language for quantum chemistry.^42,43 Within the energetic scenario of bond formation, the lowering of electron kinetic energy upon the initial encounter of atoms has long been cited as the primary origin of covalent chemical bonds. Slater identified the role of potential energy and the accumulation of electron density in the internuclear region of bonding, and Feynman developed an image of bonding in terms of electrostatic forces. Within the state scenario of a chemical bond, quantum mechanical wavefunction interference drives the lowering of the kinetic energy of two approaching atoms. Within the real space scenario, chemical bonds occur between interacting atoms as objects in the from of electron-pair bonding after final stable or metastable atom arrangements, which are divided into various kinds of covalent, ionic, metallic, and molecular bonds. Although chemical bonds have been widely argued, their descriptions remain on a qualitative level. In contrast to the electronegativity of valence electrons, the effective nuclear charge (⁠$Z i *$⁠) of an atom is a causal factor in determining bond strength and degree of ionicity and, therefore, $Z 1 * Z 2 * / Z H *$ is thought to be a good quantitative descriptor as a whole, affording the bonding strength of solids on valence electrons, regardless of the kind of chemical bond and the difference in ionicity. Simultaneously, $μ× r 1 r 2 / a B$ is closely related to the volume effect of the atom close packing system, taking the kinetic energy of relative motion into account, which is a causal factor of chemical pressure. From the above discussion, it is seen that both ensemble descriptors constructed here imply information about chemical interaction on the atomic level.

In Fig. 5, the relation between $μ× r 1 r 2 / a B$ and $Z 1 * Z 2 * / Z H *$ was checked furthermore, and an exponential analytic relation in the same form of y = A exp(−x/B) + C was also obtained. Through adjusting the scale range for alignment as illustrated in Fig. 5, the inversed volume of the unit cell (1/V₀), hardness (H_v), and fracture toughness (K_IC) on a logarithmic scale follow the same trend as the $Z 1 * Z 2 * / Z H *$ ensemble descriptor on a linear scale. It means that the physical prediction model studied here has upgraded from data regression to data generation using causal analytic relations. Of importance, it should be noted that the state/response parameter experimentally recorded for crystalline solids is a comprehensive result of chemical composition and microstructure, including grain size, pores, impurity phase, and so on, of which the latter factors change from sample to sample due to the unavoidable distribution of processing parameters in practice. In Figs. 4 and 5, it is just taking the predominant factor of chemical composition into account so that the experimental value scatters around the line of y = A exp(−x/B) + C. Keeping the microstructure effect on the state/response parameters in mind and adopting the idea of principal component analysis (PCA), the trend of mechanical properties such as hardness and fracture toughness for diamondlike solids is able to be quantitatively predicted using correlative $μ× r 1 r 2 / a B$ and $Z 1 * Z 2 * / Z H *$ descriptors, and their intrinsic value determined by chemical composition can be obtained through physics-guided small-data mining out of extant experimental data. So far, it is seen that an exponential model describes the causality of structure/thermodynamic state and mechanical properties and the correlation between chemical bonding and chemical pressure, of which the coefficient A stands for a pre-exponential factor of every dependent parameter, B for characteristic quantity with respect to each argument for crystalline solids, and C for the baseline of each class of atom close packing system, respectively, owing to both ensemble descriptors being dimensionless.

As one of their practical applications, superhard materials are useful in a wide range of industrial applications, from scratch-resistant coatings to polishing and cutting tools. Previously, theoretical and experimental efforts focused on the possibility of finding new low compressibility materials with a hardness comparable to diamonds through the correlation between valence electrons and hardness, because the strength and compressibility of a bond were thought to determine a solid’s ability to resist deformation.^2,33 Theoretical arguments have also yielded a simple scenario of covalent bonding in diamond- and zinc-blende-type crystals, promising ultrahigh bulk modulus and hardness. From Figs. 3–5 and Table II, it is observed that bulk modulus and hardness are strongly dependent on the feature of elements marked by position (period and group) in the Periodic Table of elements, and that diamond is the hardest solid with the largest bulk modulus. From the viewpoint of ensemble descriptors, lower $μ× r 1 r 2 / a B$ coupled with suitable $Z 1 * Z 2 * / Z H *$ implies higher hardness and lower compressibility for diamond/zinc blende structured solids, both of which are relevant to smaller covalent radii of elements. On one hand, the elements B, C, and N have smaller covalent radii of 0.88, 0.77, and 0.70 Å, respectively, and superhard diamond and cubic BN solids are yielded in practice.^2,41 On the other hand, for those compounds in the center of the Periodic Table of elements, cubic BC₂N was proposed as a candidate for superhard material due to its shorter bond length.⁴¹ Experimentally, cubic BC₂N was synthesized with the high-pressure synthesis method, while the bulk modulus was 282 GPa and the Vicker hardness was 76 GPa.³⁸ Through calculating its $μ× r 1 r 2 / a B$ and $Z 1 * Z 2 * / Z H *$⁠, as demonstrated in Figs. 4 and 5, it is found that the predicted H_v for cubic BC₂N agrees well with the experimental one. The hardness of cubic BC₂N is higher than that of cubic BN but less than that of diamond, which coincides with the trend predicted by chemical pressure and chemical bonding combinatorially. Moreover, cubic BC₅ was also proposed and synthesized using the high-pressure synthesis method.³⁹ As demonstrated by Figs. 4 and 5, the trend for the hardness of diamond, cubic BC₅, BC₂N, and BN solids coincides with that of $Z 1 * Z 2 * / Z H *$ with respect to $μ× r 1 r 2 / a B$⁠. So far, it has been seen that the physical prediction model established here on the basis of causal relations with respect to chemical pressure and chemical bonding ensemble descriptors provides a feasible way to predict mechanical properties from chemical formulas and provides essential guidance for searching for new materials like superhard for cutting tools and bearings.

In summary, an intelligent prediction architecture synergized by small data and physicochemical theories of materials was proposed after adopting the data science paradigm to survey the connotation and denotation of material data. This intelligent prediction model contains three key components: First, material data are classified into various classes and hierarchies in accordance with system theory, while digital ensemble descriptors, including chemical formulas and chemical bonds, are constructed using the atomic features contained in the Periodic Table of elements within the framework of quantum mechanics, the Schrödinger equation, and the symmetry principle. Second, causal analytic relations between state/response parameters and digital ensemble descriptors hidden behind extant experimental data on materials, for instance, an exponential function for diamond, zinc blende, and rocksalt-structured solids, are mined out through physics-guided small-data intelligent algorithms. Therefore, causal principles for the emergence and regulation of every target property are unveiled to understand the underlying physicochemical mechanism and to improve the performance of crystalline materials. Third, an end-to-end physical prediction model based on causal analytic relations is thus established to connect multi-objectively comprehensive performance with those features of chemical elements through digital ensemble descriptors, which imply chemical composition and atom geometric arrangement in the unit cell of crystalline materials. Our essay opens a new era other than trial-and-error and screening for both forward prediction of state/response parameters and inverse designing of chemical composition for new materials in a predictive way, e.g., superhard materials.

See the supplementary material for more details of the data and calculations of ensemble descriptors used in this study, which includes an Excel file.

The authors acknowledge the National Natural Science Foundation of China (Grant Nos. 12274066 and 61771122), the Research Project of Zhejiang Lab (Grant No. 2021PE0AC02), and the Department of Science and Technology of Zhejiang Province (Grant No. 2023C01182) for supporting this research.

The authors have no conflicts to disclose.

Zhijie Hu: Data curation (equal); Formal analysis (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Jian Yu: Conceptualization (lead); Funding acquisition (lead); Data curation (lead); Formal analysis (lead); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article and as supplementary material in an Excel file.