Effect of recording layer thickness on reducing switching current in double MgO/CoFeB interfaces pMTJ

One of the most essential properties of electrons, “spin,” is extensively explored in the industrial and scientific research communities and further emerges a burgeoning subject of “spintronics,” possessing a great quantity of valuable phenomena, e.g., tunneling magneto-resistance (TMR),¹ spin transfer torque (STT),^2,3 spin–orbit torque (SOT),⁴ voltage control magnetic anisotropy (VCMA),^5,6 etc. As the utmost importance for spintronic devices, magnetic tunnel junction (MTJ) is the first commercial application of read sensor or non-volatile random access memory,^7,8 in which information storage, write approaches, and read manner are commonly derived from bi-state resistances, Oersted field/STT/SOT, and TMR effect, respectively. Other potential realms of neuron network,⁹ true random number generator,^10,11 and processing-in-memory¹² are being progressively developed by MTJ's stochastic switching behavior.

Knowing that MTJ with a key ferromagnet/oxide/ferromagnet component is mainly divided into two categories according to magnetic anisotropy, i.e., perpendicular MTJ (pMTJ) and in-plane MTJ. The former produces more advantages for overcoming the physical scaling limit (≤10 nm node), maintaining high thermal stability (10 years retention), and relieving power consumption than its counterpart. In particular, the high performance of pMTJ with the double MgO/CoFeB interfaces has been demonstrated, for example, no perpendicular magnetic anisotropy (PMA) degradation even after annealing at 400 °C, the increased thermal stability (approximately × 2), and no significant difference in the intrinsic critical current compared to a single-interface frame.^13–16 However, compared with the traditional CMOS memory cell (∼1 fJ),^17,18 the switching energy (⁠ $= I C 2 R τ P$⁠, I_C: switching current; R: device resistance; and τ_P: pulse width) of commercial MTJ is still around two or three orders of magnitudes higher because of its high switching current/voltage and, thus, affects energy-saving efficiency and device's endurance.¹⁹ Additionally, it also leads to the need for large-sized transistors to provide this driven current, in turn, being detrimental to improve chip integration. To mitigate these issues, vigorous efforts have been dedicated to explore new writing techniques. SOT-driven pMTJ has shown a reliable switching down to 210 ps and ultralow power of 300 pJ, but a giant current density (∼2 × 10⁸ A/cm²) in the heavy metal layer may bring a non-negligible Joule heating and electromigration.²⁰ Voltage-controlled pMTJs can be a simple structure of CoFeB/MgO/CoFeB without the synthetic antiferromagnetic (SAF) structure and achieve a low switching current density of ∼10⁴ A cm² combined STT with VCMA by unipolar voltage pulses, whereas their diameter sizes (400 nm–50 μm) seem to be too large for practical applications.²¹ Although a pMTJ driven by a synergy between SOT, STT, and VCMA abbreviates the cumulative time (<200 ps) of the magnetization switching, three-terminal devices increase the difficulty of the integrated circuit design and manufacture procedure.²² Notably, these systems mainly face the fabrication complexity and scalability challenge since the additional wires are needed to implement an in-plane auxiliary field. In addition, the ferromagnetic thickness and fabrication conditions have an immediate impact on the dynamic and quiescent magnetic properties of MTJ, that is, Gilbert damping, effective magnetic anisotropy, spontaneous magnetization, and antiferromagnetic coupling strength.^23–26 Therefore, another efficient way of tuning MTJ's write current/voltage without sacrificing the tape-out period and integration density is to make a customized and simple design for the MTJ's interlayer structure,^27,28 which has not been yet investigated systematically, especially for nm-scale double-interfaces pMTJ device, not just its blanket film.

In our work, the advanced double-interfaces pMTJs with the bottom CoFeB/inserted Ta/upper CoFeB/top-MgO recording frame and [Co/Pt]_n-based SAF structures are applied to investigate the tuning effect of the recording layer thickness on reducing the pMTJ's switching current. A distinct reduction in the switching currents of cylinder-shaped-like pMTJ devices is observed when the bottom CoFeB/inserted Ta thickness increases or the upper CoFeB thickness decreases. Consequently, the considerable decline of the write energy reaches a maximum value of 41.9% and 43.4% for the bi-directional switching. Notably, the key parameters, that is, the ferromagnetic volume of the recording layer, its effective anisotropy energy density, and resistance-area product (RA) for different pMTJ stacks, play the decisive roles in the write current, whereas the TMR, Gilbert damping, interface magnetic anisotropy, gyromagnetic ratio, and spontaneous magnetization have a weak response on the fine-tuned recording frames. These analyses form a cornerstone of ultralow-power nonvolatile memory and computing chips.

Our double MgO/CoFeB interfaces pMTJ stacks from the bottom to the top consisted of Ta/Si-based seed layer, [Co/Pt]_n-based pinned layer, Ru-based bridge layer, [Co/Pt]_n/CoFeB/Ta/CoFeB-based reference layer, MgO-based barrier layer, bottom CoFeB/inserted Ta/upper CoFeB/top MgO-based recording layer, and Ru/Ta-based capping layer, as shown in Fig. 1(a). The [Co/Pt]_n-based SAF structure existed in the pinned layer/bridge layer/reference layer by Ruderman–Kittel–Kasuya–Yosida (RKKY) effect of Ru. These films were deposited on the thermally oxidized Si substrates by magnetron sputtering and then in situ annealed at the same temperature. The four series of comparative experiments for six different recording frames are conducted, as illustrated in Table I. Only the thicknesses of the recording layers were changed with the wafer, while the constituents and thicknesses of other pMTJ's stacks remained unchanged. The thicknesses (t) of the bottom CoFeB, inserted Ta, and upper CoFeB in the recording layer were varied from 1.30 to 0.95 nm, from 0.40 to 0.30 nm and from 0.75 to 0.55 nm, respectively. The cylinder-shaped-like pMTJ devices with an 80 nm diameter pattern were prepared by electron beam lithography, reactive ion etching/Ar ion-milling from the corresponding blanket pMTJ, oxide refill, and chemical-mechanical planarization.

All electrical properties were measured by utilizing a room temperature test station with a 40 GHz GSG probe and the sourcemeter instrument. All data of the electrical properties shown in this paper were the average value of ten different devices of a batch and measured at ambient temperature. The transmission electron microscopy (TEM), energy disperse spectroscopy (EDS), vibrating sample magnetometer (VSM), and vector network analyzer ferromagnetic resonance (VNA-FMR) spectrometer were performed to characterize the microstructure of pMTJ, chemical compositions, and static and dynamic magnetic properties of the recording frames at ambient temperature, respectively.

Figures 1(a) and 1(b) show the representative EDS line profile and TEM image of a blanket pMTJ sample with a recording frame of top MgO(0.30 nm)/CoFeB(0.55 nm)/Ta(0.30 nm)/CoFeB(0.95 nm), respectively. The MgO barrier is composed of the well fcc-crystallized parts and a part of amorphous zones and has a significantly better crystalline quality than the top MgO in the recording layer. Despite the recording layer appears to be an amorphous layer and its ultrathin inserted layer (∼3 Å) sandwiched with CoFeB is difficult to be observed clearly, an existence of Ta has been proved by a visible peak in the EDS line profile at a middle region of the bottom CoFeB/inserted Ta/upper CoFeB/top-MgO recording layer. Nevertheless, the actual thicknesses of the barrier-MgO/bottom CoFeB/inserted Ta/upper CoFeB in the blanket film are approximately comparable to their nominal thicknesses although a slight diffusion of top MgO at the top-MgO/capping layer interface exists. Thus, the nominal thicknesses are used to mark the pMTJs' recording stacks in subsequent discussion. Based on the etching process, the diameter (D) of cylinder-shaped-like pMTJ devices is fixed at ∼92.5 nm, as depicted in Figs. 1(c) and 1(d). Note that the TEM results from the blanket pMTJs also hold for the pMTJ devices, because the etching process mainly affects the devices' edge within a 2–3 nm damage region.²⁹ 

Hereafter, the magnetization switching characteristics of the double-interfaces pMTJ devices are investigated. Figure 2(a) presents the resistances (R) of ten pMTJ devices at wafer 5 vs applied voltage (V_DC) without resorting to the external magnetic field. The R–V_DC curves of other pMTJs at different wafers exhibit the same trend as Fig. 2(a) (not shown). This V_DC with a millisecond-scale pulse width allows us to acquire devices' resistance simultaneously when a fixed V_DC is applied to the upper and bottom electrodes of the device, which can be regarded as a quasi-static test process. The positive V_DC generates a charge current from the recording layer to the reference layer, and, thus, the magnetization directions of the two layers favor the parallel (P) configuration, vice versa, the negative V_DC is inclined to cause their antiparallel (AP) configuration. When a sufficient amount of the spin-polarized current passes through a pMTJ device, the magnetization orientation of the recording layer is able to switch; if not, the torque is insufficient to make it switch. It can be seen from Fig. 2(b) that the average TMR and its deviation for our ten pMTJ devices at wafer 1–wafer 6 are 91% (standard deviation σ = 1%), 91% (σ = 1%), 80% (σ = 1%), 90% (σ = 1%), 85% (σ = 1%), and 89% (σ = 2%), respectively. The TMR around 88% responds weakly to six different recording structures (σ = 5%), whereas the effective resistance-area product (RA) calculated by the product of P-state resistance (R_P) at a 20 mV bias voltage and circular-shaped sectional area (=πD²/4) appears to fluctuate obviously in range of 9.2–18.0 Ω μm² (mean = 13.0 Ω μm² and σ = 23%). The resistance and write voltage of the pMTJs at magnetization reversal for the P-to-AP (AP-to-P) direction denoted as R_P@write (R_AP@write) and V_P-to-AP (V_AP-to-P) are shown in Figs. 2(c) and 2(d), respectively. The R_P@write (R_AP@write) for our ten pMTJ devices at wafer 1–wafer 6 are 1981 Ω (2907 Ω), 1313 Ω (1974 Ω), 1487 Ω (2187 Ω), 1841 Ω (2641 Ω), 2503 Ω (3741 Ω), and 1809 Ω (2803 Ω), and the maximum of their deviations is 4%. The V_P-to-AP (V_AP-to-P) for our ten pMTJ devices at wafer 1–wafer 6 are −0.377 V (0.472 V), −0.313 V (0.448 V), −0.354 V (0.459 V), −0.357 V (0.506 V), −0.350 V (0.452 V), and −0.282 V (0.375 V), respectively, as well as their deviations are no more than 10%. The V_AP-to-P is larger than $| V P - to - AP|$ at a specific stack structure resulted from the STT asymmetry,³⁰ while the trend of V_{P-to-AP(AP-to-P)} with different interlayered structures (wafer 1–6) is resemble to that of R_P(AP)@write or effective RA and its change is no more than 0.131 V.

We mainly focus on the current-driven switching current of the double-interfaces pMTJs in the thermal activation regime so as to seek its tuning mechanism. The corresponding currents (I_{P-to-AP(AP-to-P)}) for the bi-directional switching and their average values (I_C) are obtained from the formulas of V_{P-to-AP(AP-to-P)}/R_P(AP)@write and $( I AP - to - P+ | I P - to - AP | ) / 2$⁠, respectively. The average switching current for the P-to-AP (AP-to-P) switching I_P-to-AP (I_AP-to-P) for our ten pMTJ devices at wafer 1–wafer 6 is −190.7 μA (162.8 μA), −238.7 μA (227.4 μA), −238.2 μA (210.2 μA), −193.9 μA (191.7 μA), −140.0 μA (120.9 μA), and −155.9 μA (134.2 μA), respectively. From the foregoing, the $| I P - to - AP|$ is slightly larger than I_AP-to-P and their asymmetry is less than 28 μA. Therefore, the I_C of the pMTJ devices as a function of the six recording frames can be obtained, as depicted in Fig. 3. It is found that (i) with decreasing the bottom CoFeB thickness from 1.30 to 0.95 nm, a significant increase of the I_C is observed at first (31.8%), but soon it persists essentially unchanged (3.8%); (ii) with decreasing the inserted Ta thickness from 0.40 to 0.30 nm, the I_C increases by a factor of 20.9%; (iii) with decreasing the upper CoFeB thickness from 0.75 to 0.55 nm, the I_C both drops rapidly in two cases, i.e., 41.8% for t_Ta = 0.30 nm and 24.8% for t_Ta = 0.40 nm. The average switching current densities [J_C = (J_AP-to-P+|J_P-to-AP|)/2, here, J_AP-to-P and J_P-to-AP are the switching current for the AP-to-P and P-to-AP switching, respectively] for our ten pMTJ devices at wafer 1–wafer 6 are 2.63, 3.47, 3.34, 2.87, 1.94, and 2.16 MA/cm², respectively. Consequently, the decline of their write energies can achieve a maximum value of 41.9% (43.4%) for the P-to-AP (AP-to-P) switching by tuning the recording layer thicknesses, and the magnitude of the STT-pMTJ's write energy is expected to drop down to picojoule or even femtojoule if the write speed is as high as several megahertz³¹ due to the fact that the switching current/voltage depends weakly on pulse width for long pulse cases in a range of 50 ns–2 μs, and it rapidly rises as the pulse widths continue to decrease.^32,33 Take one pMTJ device with a recording frame of top MgO(0.30 nm)/CoFeB(0.55 nm)/Ta(0.30 nm)/CoFeB(0.95 nm) as an example, this storage cell of MRAM can achieve quickly writes (500 MHz), and the estimated write energy for the P-to-AP switching is about 3 pJ with a write pulse width of 2 ns, as shown in Fig. 4.

To understand the tuning mechanism on the I_C reduction, the Oersted-field controlled switching properties of the double MgO/CoFeB interface pMTJs and their dynamic/static magnetic properties are also explored. Figure 5(a) reveals the resistance vs magnetic field (R–H) loops of ten randomly selected pMTJ devices at wafer 5. The well-defined reversal behavior of two-level states indicates that the ferromagnetic coupling strength between the bottom CoFeB and upper CoFeB is strong enough to make their magnetization keep coincident switching, and the TMR originated from the R–H loops is consistent with that from the R–V_DC data. Furthermore, the offset field H_O [=(H_C+ + H_C−)/2] and coercivity H_C [=(H_C+ − H_C−)/2] are shown in Figs. 5(b) and 5(c), respectively. The H_O, coming from the stray field of the SAF structure or orange-peel coupling between the recording layer and reference layer, is around 64 Oe that favors the P state. The small H_O can ensure that an external magnetic field is not required nor enforced during the write operation of MTJ. It is believed that an increase/decrease of the H_C is most probably connected with the change of the magnetic anisotropy field H_K (Refs. 21 and 34) for the coherent rotation of the magnetization reversal process (H_C = 2 K/M_S ∝ H_K; K: anisotropy constant; and M_S: spontaneous magnetization³⁵), suggesting that the H_K of the recording layers monotonically increases with shrinking the bottom CoFeB thickness and persists essentially unchanged with fine tuning the inserted Ta thickness, nonetheless, it markedly decreases with a thinner thickness of the upper CoFeB in the two cases.

The magnetic moments of the upper CoFeB and bottom CoFeB are ferromagnetic coupling based on the fact that the strong antiferromagnetic coupling at the first oscillation peak occurs at a Ta thickness of 0.7 nm,³⁶ which is significantly thicker than the inserted Ta thickness employed in our work (t_Ta = 0.3 or 0.4 nm). Therefore, the two magnetic moments of the upper CoFeB and bottom CoFeB are considered as a single entity for analyzing the dynamic/static magnetic properties of the recording layer. Figure 6(a) shows minor the magnetic hysteresis (M–H) loops of the blanket pMTJs with different recording structures in the presence of an external field perpendicular to the sample surface. The magnetization of the recording layers possesses an apparent out-of-plane easy axis except for the thickest case of the bottom CoFeB, confirming a low H_keff of the recording layer with a top MgO(0.30 nm)/CoFeB(0.55 nm)/Ta(0.30 nm)/CoFeB(1.30 nm) frame. Moreover, the total dead layer thickness t_dead and M_S of ferromagnet including the upper CoFeB and bottom CoFeB in the recording layer can be estimated by the x axis intercept and slope of a linear fitting based on the formula of m = M_S(t_CoFeB − t_dead) (m: magnetic moment of recording layer per unit area and t_CoFeB: total CoFeB thickness in the recording layer without consideration to the dead layer, i.e., nominal CoFeB thickness), as shown in Fig. 6(b). It can be found that the non-negligible t_dead obviously changes from 0.53 to 0.68 nm in a thicker t_Ta case, and the effect of the inserted Ta thickness on its M_S turns out to be weak, not exceeding 4% [Fig. 6(c)]. Their full magnetic hysteresis loops are also provided in Fig. 6(d), indicating a quick flip for the magnetization switching of the reference layer and pinned layer resulted from their good PMA [M_S(reference):∼975.6 emu/cm³, σ = 3.9%; M_S(pinned):∼1098.6 emu/cm³, σ = 4.1%], and a large exchange coupling field H_AFC [=(H_AFC+ + |H_AFC−|)/2] of [Co/Pt]_n-based SAF structures with a magnitude of ∼4.9 kOe (σ = 1.0%), as exhibited in Fig. 6(e). The similar static magnetic properties of the pMTJ multilayers with the parts of the same architecture design (e.g., reference layer, bridge layer, pinned layer, etc.) in different wafers also demonstrate the good repeatability and conformability of their manufacture process.

In addition, some other key parameters (i.e., resonance field μ₀H_res and full width at half maximum μ₀ΔH) of the recording layer are obtained from fitting the VNA-FMR data with the Lorentzian curves [Figs. 7(a) and 7(b)]. Therefore, the effective magnetic anisotropy field μ₀H_keff of the recording layer in the blanket pMTJs can be calculated by a linear fitting based on Kittel's formula f = γμ₀(H_res + H_keff)/2π = γ(μ₀H_res − δNM_S + 2μ₀K_i/M_St_CoFeB′)/2π (f: microwave frequency; γ: gyromagnetic ratio; μ₀: permeability; δN: the difference of demagnetization factor projected onto z and x directions; K_i: interface magnetic anisotropy; and t_CoFeB′: total CoFeB thickness in the recording layer with consideration to the dead layer),^37,38 as illustrated in Fig. 6(c). Considering the discrepancy of δN between the blanket film and pMTJ device (i.e., δN = 1 for the blanket film and δN ≈ 0.96 for pMTJ devices), the μ₀H_keff of the recording layer in a pMTJ device grows nearly double times from 1276 to 2700 Oe although the bottom CoFeB thickness just has a small change (∼0.15 nm), almost in agreement with the variation of their H_C (∼1.6 times) and the minor out-of-plane M–H loops [Fig. 6(a)]. In contrast, the K_i and γ are nearly independent of the bottom CoFeB thickness: the K_i (γ) is as high as 3.1727 mJ/m² (186.75 GHz/T) for the pMTJ devices at wafer 1 and 3.0469 mJ/m² (187.29 GHz/T) for the pMTJ devices at wafer 2, as depicted in Table II. Note that this calculation overlooks the change of the M_S and K_i before and after pattern etching. On the other hand, the Gilbert damping α of the recording layer can be calculated by a linear fitting to the μ₀ΔH vs f relationships through exploiting the formula of μ₀ΔH = μ₀ΔH₀ + 4παf/γ (ΔH₀: inhomogeneous linewidth broadening),³⁷ indicating that the α as low as 5 × 10⁻³ makes a weak response to a fine tuning of the CoFeB thickness. The very faint VNA-FMR signals of other samples are hard to our measurement due to the fact that their effective CoFeB thickness gets thinner.

The Slonczewski model that reflects the thermally activated characteristics of the magnetization switching can be adapted to qualitatively explain our observations of the I_C reduction, that is, I_C = I_C0{1 − (k_BT/E)⋅ln(τ_P/τ₀)}; here, I_C0 = 2αγeE/μ_Bg, E = H_kM_SV_CoFeB′/2, g = [−4 + (1 + P)³⋅(3 + cosθ)/4P^1.5,^2,39 [I_C0: critical current; E: energy barrier; k_B: Boltzmann constant; T: device's temperature; τ_P: pulse width; τ₀: inverse of attempt frequency (=1 ns); e: electron charge, μ_B: Bohr magnetron, g: spin polarization-related parameter (∼0.195 at parallel state; ∼1.795 at antiparallel state); θ: included angle of the magnetic moment between the recording layer and reference layer; $P (= TMR / ( TMR + 2)$⁠: spin polarization determined by Julliere's formula⁴⁰]. Furthermore, the dead layer in the ferromagnet consisting of the upper CoFeB and bottom CoFeB makes no contributions to the CoFeB volume in the recording layer; hence, the effective CoFeB volume (V_CoFeB′) is taken into consideration in the following discussions.

If the RA-dependent effect as well as the slight changes of the α and g as a function of the layered structures are ignored, the devices' temperatures are approximately the same as ambient temperature, and, thus, the magnitude of E in the pMTJ devices plays an essential role in I_C. Under these circumstances, two essential parameters, just as the ferromagnetic volume in the recording layer and PMA's energy density (=H_kM_S/2), have a combined effect on the variation of the I_C with different recording frames as Table III shows. (i) For tuning the bottom CoFeB thickness case, a decline of the effective CoFeB volume (V_CoFeB′) is helpful for a low I_C, whereas an increase pseudo energy density of PMA (K_eff′ = H_CM_S/2) that reflects the change of PMA's energy density with the recording frames is not favorable to reduce the energy dissipation on recording binary information. This competition makes the lowest I_C of the pMTJ devices at wafer 1, in which the product of the K_eff′ and V_CoFeB′ is much lower than that at the other two cases, even reaching a reduction of ∼38%. (ii) For tuning the inserted Ta thickness case, the I_C for the specific t_Ta (=0.30 nm) samples at wafer 4 accelerates to 20.9% at the t_Ta = 0.4 nm circumstance, seemingly in contradiction with the nominal CoFeB volume and K_eff′ remaining roughly static at different t_Ta cases. However, the V_CoFeB′ of the pMTJ devices at wafer 4 (8181 nm³) is smaller than that at wafer 5 (9227 nm³), likely causing the lower E and I_C. (iii) For tuning the upper CoFeB thickness cases with the two series of the comparative experiments, the positive effects of the V_CoFeB′ and K_eff′ are superimposed on reducing the I_C as the upper t_CoFeB decreases, which provides an effective approach to conserve the write energy and improve memory density.

In addition, the effective spin maxing conductance, Joule heating, or VCMA are used frequently to explain the RA-dependent switching current density in prior studies.^41,42 Therefore, the discrepancy of the RA tightly tied to the crystallization of MgO and CoFeB, element diffusion, and edge damage/re-deposition effect from the etch process^43–46 may also attach great importance to tuning the I_C. It can be observed that a larger RA is likely to result in a reduction of the I_C for a specific cross-sectional area. For example, as the upper CoFeB thickness decreases, the I_C in the t_Ta = 0.3 nm case falls faster than that in the t_Ta = 0.4 nm case, as explained below. On one hand, the drop of the pseudo E (=K_eff′V_CoFeB′) in the former case is just slightly smaller (44.7%) compared to its counterparts (54.0%). On the other hand, the effective RA of pMTJ devices at wafer 3 (10.6 Ω μm²) is far below than that at wafer 5 (18.0 Ω μm²), whereas the effective RAs for wafer 4 and wafer 6 are comparative (∼12.8 Ω μm²), implying that the former (wafer 3 and wafer 5) has significant difference in the RA-dependent additional effect compared to the latter (wafer 4 and wafer 6). The foregoing discussion is aimed at the qualitative analysis of the combined effect of the E and RA on the current-induced magnetization switching. The future work will be concentrated on a quantitative analysis by simulation and the optimization of the effective RA of our pMTJ devices by tuning the preparation conditions (e.g., etching time, annealing time, annealing temperature, etc.).

In our work, an advanced double MgO/CoFeB interfaces pMTJs with the [Co/Pt]_n-based SAF structures are employed to investigate their quasi-static switching characteristics by fine tuning the thickness of each recording layer, but otherwise, the thickness and composition of the reference/pinned layer and sizes of the pMTJ devices are almost the same. We observe a distinct reduction in the switching current of the pMTJ devices with different recording frames when the bottom CoFeB/inserted Ta thickness increases or the upper CoFeB thickness decreases, achieving a maximum drop of 41.8%. The competition/complementation between the V_CoFeB′ and K_eff′ of the recording layer and the change of RA are pointed out as the main tuning mechanism of the write current. Note that the strength of PMA in the MgO/CoFeB/Ta/CoFeB/MgO structure seems to fade out in the thicker bottom CoFeB cases (>1.30 nm), and an increase of the inserted Ta thickness is detrimental to acquire an ultrathin dead layer in the recording layer.

This work was supported by the Shanghai Sailing Program (No. 22YF1456400).

The authors have no conflicts to disclose.

Lili Lang: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Yujie Jiang: Data curation (equal); Supervision (equal); Validation (equal). Cailu Wang: Project administration (supporting); Supervision (equal); Validation (equal). Yemin Dong: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Validation (equal).

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