Improving plasma uniformity in the inductively coupled plasma by external magnetic field Open Access

In recent years, radio frequency (rf) plasma sources have become extensively utilized in material etching, thin film deposition, and surface treatment within the microelectronics industry.^1–5 Enhancing etching efficiency and uniformity has become a critical factor in addressing this challenge to accommodate the gradual reduction in device feature sizes and the demand for high aspect ratios in etched trenches during process production.^6,7 Inductively coupled plasma (ICP) sources often generate plasma as a dry etching device.^8,9 However, traditional ICP sources often exhibit a center-high distribution of plasma density, which leads to a large radial density gradient in the chamber. This non-uniform plasma distribution often fails to meet process production requirements.¹⁰ Therefore, to improve plasma uniformity, modifications to traditional inductively coupled plasma source devices have become a focus of interest.^11–13 

The introduction of external magnetic fields into inductively coupled plasmas has been proven to be an effective method for improving plasma uniformity.^14,15 Lee et al.¹⁶ found that introducing an external magnetic field can restrict the radial transport of electrons under low-pressure conditions, resulting in a restricted radial density distribution with enhanced density levels. Another study by Lee and Kim also showed that magnetized ICP has a more uniform plasma density distribution than unmagnetized ICP.¹⁷ Briefi et al.¹⁸ utilized optical emission spectroscopy to investigate the plasma characteristics of weakly magnetized hydrogen discharge. Their results showed that the electron density on the generator axis decreases with an increase in the magnetic field, causing the electron density to shift from radially uniform to edge-high distribution. In addition, Lee et al.¹⁹ proposed various magnet combinations to enhance plasma density. Their study further demonstrated that applying electromagnetic fields increased the etching rate but decreased the uniformity. Kim et al.²⁰ discovered that using an axial dc magnetic field increases the plasma density without affecting the axial density distribution. Furthermore, Tsakadze et al.²¹ also studied plasma discharges with cross-internal oscillation currents. Results indicated that the external magnetic field enhanced the uniformity and depth of rf power deposition, thereby leading to a more uniform plasma distribution.

Previous research has experimentally demonstrated that plasma density and uniformity were improved in ICP discharges through external magnetic fields using diagnostic methods such as magnetic induction probes, spectroscopy, and Langmuir probes. However, the experimental results obtained by these methods are often relatively simplistic, and their research is not comprehensive and specific enough. Moreover, studies on the plasma dynamic behavior in fluid models have yet to be reported. Specifically, the modulation effect of magnetic coil current on plasma uniformity at different process parameters requires further exploration. These explorations are crucial for improving plasma etching and optimizing procedures. Furthermore, plasma dynamic behavior can be affected by these external conditions, which leads to different etching results. It is known that different process parameters, including gas pressure and oxygen ratio, play an essential role in inductively coupled argon–oxygen plasma discharge.^22,23 Therefore, in this work, a two-dimensional fluid model was employed to investigate the modulation effect of an external magnetic field on the plasma density, ionization rate, ion flux, and electron temperature. In addition, the impact of various process parameters, such as gas pressure and oxygen ratio, on magnetized plasma discharge was investigated. This work aims to introduce an external magnetic field into the traditional ICP source to improve plasma uniformity, which will help to deeply understand the plasma discharge behavior and improve plasma etching technology.

This paper is organized as follows: Sec. II introduces the two-dimensional fluid model and initial discharge conditions of the argon–oxygen ICP reactor. The effect of the magnetic coil current on the distribution of plasma parameters is discussed in Sec. III, especially the modulation effect on the radial uniformity of the plasma. In addition, the functional relationship between the magnetic coil current and the gas pressure and oxygen ratio is given. Finally, the conclusions of this study are presented in Sec. IV.

Given the cylindrical symmetry of the ICP reactor, a two-dimensional axisymmetric fluid model was established based on the COMSOL multi-physics simulation software. Figure 1 shows the schematic of the ICP reactor configuration. Specifically, the plasma chemical reaction occurs within a discharge chamber with a radius of 20 cm and a height of 12 cm. A 10-turn magnetic field coil located 3 cm on the right side of the chamber is used to generate a dc magnetic field. The driving power of the inductor coil above the dielectric window is fixed at 200 W, and the rf frequency is fixed at 13.56 MHz. The working gas is an argon–oxygen mixed gas. Moreover, variables within the simulation include the magnetic coil current, gas pressure, and the O₂ ratio, with their respective ranges specified as 0–15 A, 10–20 mTorr, and 0.1–0.9. A more detailed description of the ICP reactor configuration can be found in previous work.²⁴ 

All plasma chemical reactions and rate coefficients involved in the fluid model are shown in Table I. For Ar–O₂ plasma, this work only considers common and important chemical reactions (e.g., ionization collision, elastic, excitation, dissociation, dissociative attachment) and ignores some reactions that can be ignored.^25,26 Specifically, the main species considered in this model are as follows: molecules (O₂, Ar), charged species (e, O₂⁺, O⁻, O⁺, Ar⁺), excited species (O(¹D), O₂(a¹Δ_g), Ar^*). Since O⁻ is the main negative species in low-pressure discharge, other negative species (O₂⁻ and O₃⁻) are ignored in this model.²⁷ In addition, three-body reactions and the existence of vibrating molecular oxygen are not considered.²⁸ It is assumed that charged and excited species [Ar⁺, O₂⁺, O⁺, O⁻, O₂(a¹Δ_g)] are neutralized or de-excited on the wall, and the adhesion coefficient is set to 1. Neutral particles [O and O(¹D)] are reflected by the wall, and the adhesion coefficient is set to 0.2. This work assumes that the electron energy distribution conforms to the Maxwell distribution,²⁹ and on this basis, the rate coefficients of various reactions are calculated.

In addition, the setting of meshing and boundary conditions is crucial in simulation calculations. To improve the convergence of the model and the calculation accuracy, the entire calculation domain of the fluid model is divided into 18 386 mesh units. Most of them are free triangle meshes. The average mesh unit quality is 0.88. The induction and magnetic coils use mapped meshes with 25 units and a unit size ratio 20. In addition, the boundary layer mesh is used at the chamber wall to solve the gradient problem accurately. Moreover, the initial boundary conditions are set as follows: the initial electron density is 1 × 10¹⁵ m⁻³, the initial electron average energy is 3 eV, the gas temperature is 300 K, the axial electron flux and electron energy flux are set to 0, and the ion flux and ion energy are continuous at the chamber wall (i.e., $∇ni=0$⁠, $∇⋅Γi=0$⁠). The reaction chamber wall is grounded.

In this subsection, in order to study the magnetizing effect of the magnetic coil current on the plasma, the basic discharge conditions are fixed: the driving frequency is 13.56 MHz, the gas pressure is 20 mTorr, the oxygen ratio is 0.5, and the magnetic coil current changes range from 0 to 15 A.

Figure 2 shows the distribution of various plasma parameters obtained through fluid simulation at different magnetic coil currents. Where, Figs. 2(a1)–2(a5) depict the spatial distribution of electron density, while Figs. 2(b1)–2(b5) illustrate the spatial distribution of the ionization rate. Moreover, the spatial distribution of the electric field is shown in Fig. 2(c). Specifically, Fig. 2(a1) illustrates the center-high distribution of electron density at the unmagnetized stage (i.e., the magnetic coil current is 0 A), which results in a large radial density gradient. In the weak magnetization stage (i.e., the magnetic coil current ranges from 1 to 8 A), the magnetic field generated by the magnetic coil enters the main plasma region from the right side of the chamber wall. Currently, the magnetic confinement effect from the external magnetic field is weak. During this process, electrons are heated, and the peak electron density region gradually increases from the chamber center to the chamber edge, resulting in a more uniform distribution, as depicted in Figs. 2(a2) and 2(a3). However, the magnetic coil enters the strong magnetization stage (i.e., the magnetic coil current ranges from 9 to 15 A) when its current rises to 9 A. Since the magnetic field at the chamber edge is much higher than that at the chamber center, electrons are mainly heated locally at the chamber edge, causing the peak electron density to shift from center-high to edge-high distribution in Figs. 2(a4) and 2(a5), which further worsens plasma uniformity. Moreover, the optimal plasma density distribution obtained at the magnetic coil current is 10 A [see Fig. 2(a4)]. These phenomena indicate that introducing an external magnetic field can significantly control the plasma density distribution.

In addition, the ionization rate moves from the center to the edge of the chamber with rises in the magnetic coil current, as depicted in Figs. 2(b1)–2(b5). At this time, due to the increase in magnetic field, the charge distribution within the plasma region becomes uneven, which triggers the electric field reversal phenomenon when the magnetic coil current is 8 A [see Fig. 2(c2)]. Interestingly, different from the electric field reversal in CCP discharge,^32–37 the electric field reversal here is manifested as the bipolar field lines change from right-convex distribution to left-convex distribution. The emergence of the reversing electric field causes the ionization rate to shift entirely away from the chamber center and toward the chamber edge. Additionally, under the influence of a strong magnetic field, electrons spiral along the magnetic field lines and are gradually drawn toward the magnetic field source. The probability of collisions between electrons and gas molecules within a unit volume also increases, further enhancing the ionization rate and the trend of shifting toward the chamber edge.

At the same time, the variation of various plasma parameters with the magnetic coil current is given in Fig. 3. Specifically, when the magnetic coil current is 10 A, the peak electron density is 5.42 × 10¹⁶ m⁻³, representing a 2.44 × 10¹⁶ m⁻³ improvement over the initial value [see Fig. 3(a)]. The uniformity degree is 94%, which means a 39% improvement over the initial value [see Fig. 3(b)]. In addition, it can be seen from Fig. 3(b) that the plasma uniformity degree improves from 55% to 83%, going from 0 to 8 A in the magnetic coil current. However, the uniformity decreases to 82% with further rises to 15 A in the magnetic coil current. This indicates that the weak magnetization stage positively impacts the uniformity distribution of plasma density. Conversely, during the strong magnetization stage, the confinement effect of the strong magnetic field on electrons leads to the tendency of electron density to distribute toward the chamber edge, thereby breaking the radial uniformity distribution of the plasma. Furthermore, when the magnetic coil current increases from 0 to 15 A, the ionization rate increases from 3.84 × 10²¹ to 2.22 × 10²² m⁻³ s⁻¹, which also means that the magnetic coil current can significantly improve the ionization rate [see Fig. 3(a)]. On the other hand, a higher ionization rate implies that more electrons and ions are produced in the plasma region, leading to an increase in the ion flux reaching the substrate. As a result, the O⁻ and O₂⁺ ion flux increases with the magnetic coil current, as shown in Figs. 3(c) and 3(d). However, due to the magnetic confinement effect, the electron flux reaching the substrate is limited. Therefore, we can consider the ions to be unmagnetized.³⁸ That is, the external magnetic field does not directly affect the distribution of the ion flux. Thus, the ion flux does not show a distribution similar to the electron density (i.e., edge-high distribution).

It is worth noting that in a strong magnetic field, the movement of electrons is constrained, leading to a significant increase in the ionization rate and electron density near the chamber edge. Furthermore, electron density escalates with increases in the magnetic coil current, that is, in the magnetization stage, including weak and strong magnetization, with the peak electron density continuously rising [see Fig. 3(a)]. This phenomenon is attributed to the magnetic confinement effect exerted by an external magnetic field on electrons. As the external magnetic field increases, the electron cyclotron radius becomes smaller. Consequently, higher magnetic field intensity restricts the diffusive movement of electrons, thereby reducing electron losses on the chamber walls and, in turn, enhancing the electron density. Moreover, with the increase in the magnetic field, there is a tendency for electron density and ionization rate to transition from the center to the edge of the chamber, a trend that is accentuated under the influence of a strong magnetic field [see Figs. 2(a5) and 2(b5)].

Process parameters are also crucial for modulating plasma density distribution. Therefore, to deeply explore the influence mechanism of the external magnetic field on plasma uniformity, we further studied the modulation effect of the external magnetic field under different gas pressures and O₂ ratios. Specifically, Fig. 4 shows how different magnetic coil currents change the distribution of plasma parameters, including electron density [see Figs. 4(a1)–4(b5)], peak electron density [see Fig. 4(c)], and uniformity degree [see Fig. 4(d)], at gas pressures of 10 and 15 mTorr. It is observed that with the increase in magnetic coil current, the electron density radial distribution becomes increasingly uniform in Figs. 4(a1)–4(a5). This indicates the plasma discharge is mainly in the weak magnetization stage within 0–15 A. This phenomenon is attributed to the effects under low-pressure conditions (10 mTorr). As the magnetic field increases, the electrons are more strongly magnetized. Due to the lower collision frequency between particles and the lower electron absorption power, which results in the radial diffusion movement of electrons in the plasma region slows down, and the occurrence of the reversal electric field is also suppressed. Moreover, as the magnetic coil current increases, the plasma density distribution becomes more uniform. When the magnetic coil current reaches 15 A, the peak electron density increases from 1.54 × 10¹⁶ to 3.81 × 10¹⁶ m⁻³ [see Fig. 4(c)]. At this point, the plasma uniformity reaches the best, with a uniformity of 85%, marking a 29% increase from the initial value, as shown in Fig. 4(d).

However, in Figs. 4(b1)–4(b5), the weak magnetization stage ranges from 0 to 10 A at a pressure of 15 mTorr, which is delayed by 2 A compared with Fig. 2(a2). In addition, interestingly, with a magnetic coil current of 10 A, the E × B drift resulting from the magnetic field leads to uneven charge distribution within the plasma, thereby triggering the electric field reversal phenomenon. The reversal electric field caused by the magnetic field significantly enhances the ionization rate and electron absorption power, leading to an increase in peak electron density from 3.65 × 10¹⁶ to 8.94 × 10¹⁶ m⁻³ [see Fig. 4(c)]. In short, as the gas pressure rises from low to high, due to the increase in the frequency of collisions between particles, the magnetization behavior of the plasma caused by the external magnetic field will be delayed. Furthermore, as illustrated in Figs. 2(a5) and 4(a5), the plasma uniformity degree decreases from 85% to 82% [see Figs. 2(d) and 4(d)]. This suggests that higher gas pressures undermine the modulating effect of the magnetic field on plasma density distribution, thereby leading to worse plasma uniformity. On the other hand, when the magnetic coil current is 0 A in Fig. 4(c), the gas pressure rises from 10 to 15 mTorr, and the peak electron density increases from 1.54 × 10¹⁶ to 2.36 × 10¹⁶ m⁻³. This suggests that electron density escalates with increasing gas pressure, though plasma uniformity remains relatively unchanged.

In addition, the sharp increase in the curve of the peak electron density at a gas pressure of 15 mTorr in Fig. 4(c) is the result of the reverse electric field. As shown in Fig. 5(a5), when the gas pressure is 10 mTorr, the electric field reversal appears when the magnetic coil current is 15 A. Currently, the electron density growth trend is relatively slow [see Fig. 4(c)]. However, when the gas pressure is 15 mTorr, the electric field reversal appears when the magnetic coil current is 10 A, as shown in Fig. 5(b4). After that, as the magnetic coil current further increases, the electron density growth curve gradually steepens [see Fig. 4(c)]. This means that the electric field reversal can significantly increase the electron density. In addition, a lower gas pressure (10 mTorr) inhibits the appearance of the electric field reversal. Furthermore, the variation trend of the peak ion flux with the magnetic coil current is given in Figs. 5(c) and 5(d). It can be seen that the O⁻ ion flux increases with the increase in gas pressure and magnetic coil current [see Fig. 5(c)]. The O₂⁺ ion flux decreases with the rise of gas pressure in the magnetic coil current range from 0 to 12 A. Interestingly, it increases with the increase in gas pressure as the magnetic coil current further increases [see Fig. 5(d)]. This results from the electric field reversal caused by the magnetic field at higher gas pressure.

In order to explore the evolution of electron density under different magnetic coil currents, the spatial distribution of electron temperature is given, as shown in Fig. 6. It can be clearly seen that in the low-pressure (10 mTorr) unmagnetized stage [see Fig. 6(a1)], the peak of the electron temperature is mainly located below the dielectric window. At this time, the electrons are mainly heated below the coil, and due to the obvious power deposition results in a chamber center-high distribution of electron density. As the magnetic coil current increases to 10 A [see Fig. 6(a4)], the peak electron temperature below the coil decreases from 4.34 to 4.14 eV, as shown in Fig. 6(c). At this time, the electron temperature appears two peaks below the coil and near the chamber wall, and the peak electron temperature exhibits a gradual transition toward the chamber wall. This directly causes the center-high distribution to gradually diffuse toward the chamber wall of electron density, thus improving the radial uniformity of the plasma. When the magnetic coil current further increases to 15 A, the peak electron density increases to 4.59 eV. Meanwhile, the high distribution under the coil is completely transformed into a high distribution at the chamber wall of peak electron density. It is obvious that the control effect of the magnetic coil on the plasma density is more obvious, and due to the reduction of the collision frequency in the plasma under low-pressure conditions, the electron energy loss is reduced, which further leads to the peak electron density being enhanced. At the same time, the best plasma radial distribution is obtained when the magnetic coil current is 15 A [see Fig. 4(a5)].

However, at higher gas pressure (20 mTorr), this transition behavior of electron temperature is advanced, and an obvious double-peak electron temperature distribution appears in the main plasma region when the magnetic coil current is 8 A, as shown in Fig. 6(b2). This also means that the plasma is more evenly distributed. In addition, the average peak electron temperature (3.65 eV) is also lower than the low-pressure case (4.59 eV) due to the increased electron energy loss due to collisions in high-pressure conditions [see Fig. 6(c)].

In this subsection, Fig. 7 shows how different magnetic coil currents change the distribution of plasma parameters, including electron density [see Figs. 7(a1)–7(b5)], peak electron density [see Fig. 7(c)], and uniformity degree [see Fig. 7(d)], at O₂ ratios of 0.3 and 0.7. The peak electron density declines from 5.18 × 10¹⁶ to 1.86 × 10¹⁶ m⁻³ in Figs. 7(a1) and 7(b1), indicating a decrease in electron density with rising O₂ ratios. Moreover, the electron density transitions from center-high to uniform to edge-high distribution with increases in the magnetic coil current, as depicted in Figs. 7(a1)–7(a5). This evolution trend is similar to Figs. 2(a1)–2(a5). Furthermore, in the weak magnetization stage (0–8 A), electrons in the chamber center are distant from the external magnetic field source and experience minimal magnetic confinement effects, leading to a high and uniform electron density. However, at the strong magnetization stage (9–15 A), the electron density exhibits an edge-high distribution due to the strong magnetic field at the outer wall of the chamber. At this time, the plasma uniformity degree is increasingly lower. It can also be seen from the peak electron density growth trend in Fig. 7(c) that the growth trend of the electron density is slow in the weak magnetization stage, and the growth trend of the electron density is significantly intensified in the strong magnetization stage. The electric field reversal phenomenon is observed when the magnetic coil current reaches 8 A. This reversal electric field accelerates electron movement and imparts additional energy, thereby enhancing electron absorption power and the ionization rate, which in turn increases electron density. Interestingly, Figs. 7(a1)–7(a5) and 7(b1)–7(b5) exhibit identical change patterns, which indicates that the O₂ ratio does not influence the external magnetic field effect. Furthermore, different O₂ ratios include 0.3, 0.5, and 0.7, all of which achieve the optimal plasma density distribution when the magnetic coil current is 10 A, and the uniformity degree is 90%, 95%, and 94%, respectively, as shown in Figs. 2(d) and 7(d).

In addition, the effect of the magnetic coil current on other plasma parameters, such as the electric field and ion flux, is shown in Fig. 8. As mentioned above, the electric field reversal caused by the external magnetic field is the main factor causing the sharp change in electron density. Specifically, when the O₂ ratio is 0.3, an obvious electric field reversal can be observed at a magnetic coil current of 8 A [see Fig. 8(a2)]. Similarly, the electric field change when the O₂ ratio is 0.7 is similar to the former [see Fig. 8(b2)]. This shows that the oxygen content will not affect the electric field reversal. In addition, the changing trend of the ion flux is shown in Figs. 8(c) and 8(d). The ion flux increases with the increase in the magnetic coil current. In addition, the O⁻ ion flux decreases with the increase in the O₂ ratio [see Fig. 8(c)], and the O₂⁺ ion flux increases with the increase in the O₂ ratio [see Fig. 8(d)]. At the same time, a significant growth trend of the ion flux is observed in the strong magnetization stage. This change trend is the result of the reversed electric field. In summary, in addition to their inherent impact on plasma discharge, varying oxygen contents do not further affect the control effect of the external magnetic field on the plasma density distribution. On the contrary, different magnetic induction strengths will positively or negatively impact the argon–oxygen plasma discharge process.

The distribution of electron temperature is shown in Fig. 9. In the weak magnetization stage, the electron temperature exhibits a double-peak distribution, as shown in Figs. 9(a2) and 9(b2). This phenomenon is more pronounced in the case of a high O₂ ratio (0.7). In the strong magnetization stage, i.e., when the magnetic coil current is 9 A [see Figs. 9(a3) and 9(b3)], since the electrons are primarily heated by the magnetic coil, the peak electron temperature exhibits a high distribution at the chamber wall. This behavior is more pronounced at lower ratios (0.3). Consequently, when the O₂ ratio was 0.7, the best plasma uniformity (95%) was achieved at a magnetic coil current of 10 A, as depicted in Fig. 7(b4). Furthermore, as the O₂ ratio increases from 0.3 to 0.7, the peak electron temperature rises from 3.27 to 4.11 eV [see Figs. 9(a5) and 9(b5)]. This indicates that the electron temperature increases with the oxygen ratio.

On the other hand, as the magnetic coil current increases from 0 to 15 A, the electron temperature shows an evolution that first decreases and then increases, as shown in Fig. 9(c). For example, in the case of an oxygen ratio of 0.7, the peak electron temperature decreases from 3.64 to 3.45 eV and then increases to 4.11 eV. This is because the energy coupled into the plasma by the magnetic coil power in the weak magnetization stage is small. The electron energy is lost faster due to the increase in collision frequency, which causes the electron temperature to decrease. When the magnetic coil current further increases to the strong magnetization stage, the enhanced magnetic confinement effect reduces the collision frequency between electrons and other particles, thereby slowing down electron energy loss. At the same time, due to the appearance of electric field reversal, the electron absorption power is enhanced, resulting in electrons gaining more energy from the external magnetic field, thereby increasing the electron temperature in the plasma.

In this work, a two-dimensional fluid model is used to study the modulation effect of an external magnetic field on plasma uniformity in argon–oxygen plasma discharges. By discussing the distribution of plasma parameters, including electron density, ionization rate, electric field, ion flux, and electron temperature, the effects of different magnetic coil currents, gas pressures, and oxygen ratios on the plasma were studied.

The plasma density increases as the magnetic coil current increases at fixed gas pressure and oxygen ratio. In the weak magnetization stage, the external magnetic field can significantly improve the plasma density and uniformity. The best plasma uniformity degree of 94% was obtained when the magnetic coil current was 10 A. Instead, in the strong magnetization stage, high magnetic induction intensity will cause the plasma density to change from center-high to edge-high distribution. At this stage, the electric field reversal caused by the magnetic field is the main factor causing the sharp increase in electron density. As the magnetic coil current further rises, the plasma uniformity deteriorates again. Moreover, the ionization rate and the ion flux increase with the increase in the magnetic coil current, in which the ionization rate changes from the center distribution to the chamber edge distribution. In addition, under low-pressure conditions, plasma uniformity is positively correlated with magnetic coil current. The best plasma uniformity degree of 91% was obtained at 15 mTorr and 15 A. In contrast, high-pressure conditions will destroy the modulation effect of strong magnetic fields on plasma density. On the other hand, as the magnetic coil current increases, the electron temperature in the plasma shows a double-peak phenomenon, and its peak value has an evolution trend of first decreasing and then increasing. Furthermore, the magnetization effect of the external magnetic field on the plasma is not affected by the oxygen content.

In summary, the plasma density and uniformity can be effectively improved by introducing an external magnetic field. In addition, process parameters such as gas pressure and oxygen ratio can also have a positive or negative impact on the behavior of the magnetized plasma by an external magnetic field. These findings indicate that introducing an external magnetic field into a traditional ICP source can help improve etching uniformity and optimize the production process.

This work was financially supported by the Opening Project of Science and Technology on Reliability Physics and Application Technology of Electronic Component Laboratory (No. ZHD201701), the Natural Science Basic Research Plan in Shaanxi Province of China (No. 2024JC-YBMS-342), and the Youth Innovation Team of Shaanxi Universities.

The authors have no conflicts to disclose.

Yang Zhao: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Xiaohua Zhou: Conceptualization (equal); Investigation (equal); Supervision (equal). Jianxiang Zhang: Data curation (equal); Investigation (equal); Writing – review & editing (equal). Shasha Song: Data curation (equal); Investigation (equal); Writing – review & editing (equal). Yuzhen Zhao: Project administration (equal); Supervision (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.