Abstract

A generalization of a Burridge-Knopoff spring-block model is investigated to illustrate the dynamics of transform faults. The model can undergo Hopf bifurcation and fold bifurcation of limit cycles. Considering the cyclical nature of the spring stiffness, the model with periodic perturbation is further explored via a continuation technique and numerical bifurcation analysis. It is shown that the periodic perturbation induces abundant dynamics, the existence, the switch, and the coexistence of multiple attractors including periodic solutions with various periods, quasiperiodic solutions, chaotic solutions through torus destruction, or cascade of period doublings. Throughout the results obtained, one can see that the system manifests complex dynamical behaviors such as chaos, self-organized criticality, and the transition of dynamical behaviors when it comes to periodic perturbations. Even very small variation of a parameter can result in radical changes of the dynamics, which provides a new insight into the fault dynamics.

1. Introduction

The system contains a series of spring-blocks manifests complicated nonlinear dynamic behavior which can be related to the phenomenon in earthquake. This kind of dynamics has attracted the attention of many scientists during the last several decades. Some characteristics of earthquakes are detected, such as the periodicity, chaos, and self-organized criticality [14]. In general, earthquakes are mainly considered due to the slip of tectonic plates [5]. Burridge and Knopoff proposed a physical model whose dynamics obey Gutenberg-Richter and Omori-Utsu law for certain parameters [6]. Afterwards, in 1989, based on the Burridge-Knopoff model, Carlson and Langer further presented another version of BK model for fault dynamics consisting of a uniform chain of blocks and springs pulled slowly across a rough surface [7]. This model can give rise to events of various sizes, and the numerical evaluation of the distribution of earthquake magnitudes results in a power-law spectrum similar to the earthquake in nature. Since then, abundant spring-block models were studied by excellent numerical simulations and statistical methods [810]. Time delay is introduced in one-block system transforming the system into infinite-dimensional, which makes it possible for the occurrence of deterministic chaos [11]. When it comes to periodic parameter perturbations, Kostić et al. showed that minor perturbation to the parameters such as the oscillation amplitudes and the angular frequencies are sufficient to change the original behavior, leading to the onset of chaos [12]. The effect of a transient periodic external force was introduced to describe the dynamic triggering impact of nearby or distant earthquakes [13].

There is still an important open question in seismology. Is the temporal or spatial complexity of earthquakes generated by the nonlinearities of the equations or due to geometrical heterogeneities? In [14], Kelemen et al. noticed a periodic shear-heating for intermediate-depth earthquakes, which produces periodic oscillations in the spring stiffness during the driving faults, suggesting a form of heterogeneity in the system. In this paper, our aim of the research is to study both the unforced and periodically forced spring-block system theoretically and give bifurcation analysis in the whole parameter space. For the Burridge-Knopoff (BK) model, the nonlinear dynamical behavior is mainly induced by the friction between the blocks and the rough surface of the lower plate. In fact, there will be disturbance in the stiffness of the springs due to the environmental influence, such as the temperature changing during the plate movement at different seasons. Naturally, a new question arises: how periodic parameter disturbance works on the motion of the spring-block system? To answer this question, we focus on the spring-block model with the interaction between the blocks, which corresponds to the motion induced by the slipping of the faults [15].

The periodic forced system discussed in the present work is a nonautonomous system, which makes the problem much more difficult. Considering the classical methods cannot be applied directly, a continuation technique and numerical bifurcation analysis are applied to investigate the complex dynamical behavior. We will show that this kind of perturbation would generate much more complex dynamics. It gives rise to the onset of deterministic chaos and illustrates the switch of dynamics, the coexistence of multiple attractors. There are periodic solutions with various periods, quasiperiodic solutions. It can interpret the periodicity in deep-focus earthquake [16], the unpredictability of earthquake, and the system is self-organized to a critical state with an avalanche followed [17]. Therefore, we explore a model with periodic perturbation which can give rise to much more complex dynamical behavior in the spring-block system.

2. The Model and Bifurcations

Firstly, we start from a multiple slider-block model including a square array of slider blocks [3], which is illustrated in Figure 1. The blocks are attached to the driver plate with a driver spring, with spring constant . Each block is attached to its adjacent block with connector springs, and the spring constant is . When the block begins to slip, its motion is given bywhere is the position of block in the array, is the mass of each block, presents the dynamic frictional force between the block and the surface when the block is moving, and is time.

The behavior of the system is determined by the two parts: the frictional force and the stiffness of the system . For very soft systems , the blocks exhibit stick-slip behavior independently. For very stiff systems, , the array of blocks behaves as a single block [3]. Linear array of multiple slider-block models and two-dimensional array of slider-block model are simulated by using the cellular-automata approach [18, 19]. Huang et al. carried out simulations on a square array of blocks using static dynamic friction and a cellular-automata approach [20]. The statistics for the slip events manifest power-law distribution. For larger stiff systems, the entire grid of slider blocks is strongly correlated and large slip events including all the blocks occur regularly. For soft systems, the stiffness of the system relatively small, there is no large event during the motion. So far most of the literatures reported the simulation results based on the slider-block model. One extension was to consider a linear chain of slider blocks pulled at one end which is a model for a propagating rupture on a fault [21]. In the following, we present theoretical analysis and simulation based on the model containing a chain of blocks.

Secondly, we explore the model containing a chain of blocks attached through a harmonic spring to a driving plate (Figure 2), which causes the stick-slip motion of the block along the rough surface of the lower plate. The block is dragged (using a spring of stiffness ) on a frictional floor by a driver moving at a rate in the direction. A chain of blocks are coupled to each other by harmonic springs with stiffness of . Our research begins with the following model:where , , and are the mass, the position coordinate, and the frictional coefficient of the ith block, respectively. The normal force is applied to the block. Because of the continuity of the earth’s medium, we will discuss the problem in the framework of a continuous system. Applying dynamical analysis, we find that there is Hopf bifurcation and fold bifurcation of limit cycles for the undisturbed system.

We discuss the discrete model (2) in the framework of a continuous system. The expression was adopted as a dimensionless variable representing the position of the ith block, where is the length between each block. Using (2), it can be deduced that

The frictional stress at the base of the block obeys a rate- and state-dependent friction law [22]where is a state variable; is a characteristic slip distance. and are constants that represent the rate and time dependence of the friction, respectively. and are the steady-state values at a reference velocity , which is chosen to be in the present study. is the slipping velocity of the block, which is equal to . From (4), we can get

In the current study, the earthquake is considered a wave propagating in the plates. The interactions between the blocks resemble the transmission of waves in a homogeneous isotropic elastic body [23]. Thus, to solve the partial differential (3), the strain evolution is traced during the deformation by analyzing the traveling wave transformation. According to the traveling wave transformation, , where is the speed of the traveling wave, , , . Equation (3) can be transformed into the following ordinary differential equation, which enables obtaining solutions in the form of propagating kinks:Correspondingly (4) and (5) are rewritten as

Set , , and , and let denote the derivative with respect to . We can deduce the system of the motion as follows:where are all positive constants. We nondimensionalize the quantities by takingSetting , , , , and , letting denote the derivative with respect to in the following, and dropping all the bars, then system (10) becomesClearly, system (12) has only one equilibrium . The Jacobian matrix at takes the formThe corresponding characteristic equation is

Then for system (12), we have the following results.

Theorem 1. (1)The equilibrium is asymptotically stable if(2) is unstable ifor .(3)System (12) undergoes a Hopf bifurcation if

Proof. DenoteThe first two items can be easily obtained from the Routh-Hurwitz criterion by regarding or as bifurcation parameters, respectively. We focus on the third case. Consider as a bifurcation parameter. Substituting () into (14) yieldsThe three roots areOn the other hand, from (14) we can obtainThus, for Similarly, if we take as a bifurcation parameter,andfor . This completes the proof.

The first Lyapunov coefficient that determines the direction of the Hopf bifurcation can be calculated directly by the formula in [24] after reducing system (12). And we choose not to present here due to the complexity for we cannot give the concise explicit condition to determine the sign of it. Therefore we show that can have alternate sign or vanish by simulation as follows. In our simulation below, we will present Hopf bifurcation curve on the plane. So we must ensure that the transversality condition holds whenever or is taken as the bifurcation parameter. That is why we choose different parameters to unfold the Hopf bifurcation in Theorem 1.

Bifurcation curves of system (12) are given in Figure 3 for two different sets of parameters, where we show that the first Lyapunov coefficient can have alternate sign or vanish at some parameter values. And near the generalized Hopf bifurcation fold bifurcation of limit cycles is detected. and represent the supercritical and subcritical Hopf bifurcation, respectively. is the fold bifurcation curve of limit cycles. The triangle is the generalized Hopf bifurcation labeled with . In Figure 3(a), the equilibrium is stable in region 1. When is crossed from region 1 to region 2, becomes unstable and a stable limit cycle arises. If crossing from region 2 to region 3 (see the enlargement of Figure 3(a) in Figure 3(b)), becomes stable again and an unstable limit cycle appears. That is, in region 3 there are two limit cycles: one unstable; one stable. On curve the two limit cycles in region 3 collide leading to one semistable limit cycle which disappears in region 1.

Here we present another case in Figure 3(c) corresponding to the bifurcation diagram for the periodically forced system in Figure 6. In Figure 3(c), is unstable in region 1; it becomes stable and an unstable limit cycle arises through from region 1 to region 2. Crossing from region 2 to region 3, becomes unstable again and a stable limit cycle appears. The two limit cycles in region 3 collide on curve yielding one semistable limit cycle which disappears in region 1.

3. The Model with Periodic Perturbation

After traveling wave transformation and nondimensionalize, we introduce the perturbation with time dependent, which will present that even a small oscillation in the spring stiffness can sufficiently change the original behavior, leading to the burst of chaos or other complex dynamical behaviors.

The parameter is generated by a linear transformation from the parameter . Therefore, the parameter reflects the spring stiffness between the blocks and the driving plate. During the sliding of the blocks, considering the influence of external parameter (such as the temperature, viscous properties, or the sliding velocity), the parameter will not always keep a constant; it fluctuates around a certain constant in the realistic motion. Set it fluctuate between a minimum and a maximum value. There exists periodic shear-heating for intermediate-depth earthquakes [14], which can produce periodic oscillations in the spring stiffness during the driving faults. And the slow rupture along the faults or some nonnatural source of vibrations, the earth tides, and reservoir effects also may induce some periodic perturbation [25, 26]. So we choose a sinusoidal function to present this perturbation (see Figure 4) and then analyze how such periodic disturbance in the parameter influences the dynamical behavior of the system. As a matter of fact, although the realistic seismic waves do not present as idealistic sinusoidal perturbations, this kind of periodic perturbation is also meaningful note that a real wave pattern can be generated by the superposition of sine waves.

The periodic perturbation in which reflects the stiffness of the spring is written aswhere the parameter is the value of without perturbation or, equivalently, the average value of with perturbation; is the amplitude of the perturbation. In this setting, note that is positive; we have , while for a given , the parameter determines the oscillation amplitude. For a given , the parameter is the ordinate cut-off. We investigate the dynamical behavior in parameter space.

To simplify the analysis, we fix . System (12) can be transformed intoIn order to study the dynamics of the periodically forced system (27), we transform the system into a higher dimensional autonomous system by augmenting the forced system with a decoupled pair of ODEs which have an oscillatory solution of the form required for the forcing. That is,System (28d)-(28e) has an asymptotically stable periodic solutionIn this way the equilibrium of the two-dimensional unforced system corresponds to the periodic solution of the five-dimensional system with . We can study whether the periodic solution can survive and whether the bifurcations of the equilibrium of the unforced system can be extended to the forced system as bifurcations of periodic solutions when .

The standard procedure for identifying and locating bifurcations of a periodic system such as (27) uses the Poincaré map that transforms the continuous system into a discrete one by sampling the solution once in each forcing period. For our system the forcing period is one year, so the Poincaré map isThe phase portrait of the Poincaré map is composed of fixed points, regular and irregular invariant sets, and all other orbits. Particularly, we will consider (i)the fixed points of the th iteration of the map, which correspond to the subharmonical periodic solutions of period ;(ii)the closed and regular invariant sets of the map, which correspond to the quasiperiodic solutions (invariant tori);(iii)the irregular invariant sets, which correspond to chaotic solutions (strange attractors).

Fixed points of the map bifurcate at some parameter values where the fixed points change stability. For a periodically forced system of ODEs, the corresponding Poincaré return map undergoes bifurcations including NS (Neimark-Sacker), flip (period doubling), and tangent (fold) bifurcations. Following the convention in [27], we use , , and to denote the above bifurcation curves for the fixed points of the th iteration of the Poincaré map, respectively.

By means of continuation, the above bifurcation curves can be numerically computed. In the bifurcation diagrams presented below, we use the following codimension 2 bifurcation points of the Poincaré map as the organizing centers:(i), or : 1:1 resonance;(ii), or : 1:2 resonance;(iii), or : 1:3 resonance;(iv), or : 1:4 resonance.

In the following, we will explore and present bifurcation diagrams of the Poincaré map corresponding to two different cases of the unforced system. The correspondence of solutions for the original continuous system to those for the Poincaré map has been given above.

Case 1. The unforced system (12) undergoes a supercritical Hopf bifurcation when is varied with other parameters fixed.

Taking , system (12) undergoes a supercritical Hopf bifurcation when is varied (Figure 3(a)). For these values, the unforced system oscillates on a stable limit cycle. The asymptotically period of the cycle (evaluated numerically) is when approaches . The bifurcation diagram of system (28a), (28b), (28c), (28d), and (28e) on the plane is illustrated in Figure 5.

In Figure 5(a), the point () on the -axis corresponds to the supercritical Hopf bifurcation in the unforced system. At point the period of the stable limit cycle of the unforced system is 2. In region 1, system (28a), (28b), (28c), (28d), and (28e) has a stable period one solution. Actually, it is . The periodic solution loses its stability in region 2 (below ) and a stable quasiperiodic solution appears which may be destroyed in some subregion of 2 by homoclinic structure near strong resonance leading to chaos. When or is crossed to the right (in region 3), two periodic solutions with period two arise: one is stable; the other one is saddle type. In regions 4 and 5, only one period two solution remains which is stable in 4 and unstable in 5. That is, the saddle type period two solution disappears through . On the other hand, the period one solution becomes saddle in regions 4 and 5 (crossing ). In region 6, the period two solution in regions 4 and 5 becomes saddle type and a period four solution arises.

Case 2. The unforced system (12) undergoes a subcritical Hopf bifurcation when is varied with other parameters fixed.

Taking , system (12) undergoes a subcritical Hopf bifurcation when is varied (Figure 3(a)). For these values, the unforced system has an unstable limit cycle. The asymptotically period of the cycle (evaluated numerically) is when approaches . The bifurcation diagram of system (28a), (28b), (28c), (28d), and (28e) on the plane is illustrated in Figure 6. Note that Figure 6(e) is almost the same as Figure 6(c). That is, the part of Figure 6(a) above is almost the same as the part between and . We must mention that, above , there are similar structures as well. We choose not to present more here. And in the following we only analyze the part shown in Figure 6(c).

In Figure 6(a), the point () on the -axis corresponds to the subcritical Hopf bifurcation in the unforced system. We will understand Figure 6(a) from its enlargements in Figures 6(b), 6(c), 6(d), and 6(e). In Figure 6(b), the system has an unstable period one solution () in region 1; it becomes stable in region 2 and an unstable quasiperiodic solution arises, which may be destroyed by some homoclinic structures. In region 3, the period one solution becomes saddle-type and a stable period two solution appears which becomes unstable in regions 4 and 5. Moreover in region 4, there is a periodic solution of period four.

In Figure 6(c), the system has a stable period one solution in region 1; it becomes unstable in regions 2 and 3 and 4. In region 3, there is a saddle-type periodic solution of period 2. In region 4, the system has a stable period two solution which turns unstable in regions 5, 6, and 7. In regions 5 and 7 there is a periodic solution of period four. Furthermore, there is a very close to . Between and (region 1 in Figure 6(d)), there are two period two solutions: one stable; one saddle type.

4. Complex Dynamic Behaviors

We have explored the dynamics of the spring-block model considering the interaction between a series of blocks. Comparing bifurcations in above sections, we can see that periodic perturbation can generate abundant dynamics.

In the following, we take the case in Figure 5 as an example to understand the complex dynamics. Firstly, periodic perturbation supports the existence of multiple attractors. Indeed, the unforced system () has always a unique attractor; it is an equilibrium above and a limit cycle below . For the periodically forced system (), there exist multiple attractors. For example, in some subregion of 3 (below ) a stable cycle of period 2 first coexists with a quasiperiodic solution and then with a chaotic solution; see Figures 7, 8, and 9. Figures 7(a) and 7(b) display solutions for , where we see a periodic solution of period two and a quasiperiodic solution (or a very high period solution), respectively. Figures 7(c) and 7(d) are solutions for . The solution in Figure 7(c) is of period two. The solution in Figure 7(d) is considered as a chaotic solution, which will be explained later. What makes differences between Figures 7(a), 7(b), 7(c), and 7(d) is the different initial value. To make types of solutions more readable, we give their phase portraits and Poincaré sections in Figures 8 and 9. Specially, the quasiperiodic solution in Figure 7(b) and the chaotic solution in Figure 7(d) are too similar to distinguish whether they are quasiperiodic or chaotic. Even their Poincaré sections look similar. However, from the partial enlargement of Figure 9(d) on the top right corner we can conclude that the Poincaré section of Figure 7(d) is not a regular closed curve like the one for Figure 7(b). Therefore we consider the solution in Figure 7(d) as chaotic and will calculate the maximal Lyapunov exponent to verify it later.

The deterministic chaos in Figure 7(d) is further corroborated by calculating the positive value of the maximal Lyapunov exponent. We use Wolf method to calculate the largest Lyapunov exponents [28]. First of all, we reconstruct phase space based on the time series in Figure 7(d), where is the number of data points. Before reconstructing the phase space, the time delay is determined by using mutual information method [29], and the embedding dimension is calculated by the Cao-method [30]. Define , where is the i-th evolution time; the set of constitutes the reconstructed attractor. Start from an initial point, , and find its nearest neighbor point, . Set the distance between the two points to be ; we track the evolution of the two points, after a time of at which the distance of the two points evolves to be , where is a given small constant which is slightly larger than the minimum distance of each two points. Then find the nearest neighbor point of , and write it as , for which the angular separation between the evolved and replacement elements is as small as possible. The distance between the two points is . Track the evolution to get and repeat the above steps until the end of the time series. Set the number of the iterations during the evolution , then the largest Lyapunov exponent is . In this case we calculate that the largest Lyapunov exponent is , where the time delay , and the embedding dimension . The positive largest Lyapunov exponent suggests that the solution in Figure 7(d) is chaotic.

Secondly, in some region of Figure 5 even very small change of a parameter can lead to radical changes of the behavior of the system. For example, in region 1 the system has only one attractor, namely, a stable cycle of period 1. If is slowly increased with fixed, the stable cycle varies smoothly and loses stability through . After is crossed, the system oscillates on another attractor, i.e., a stable cycle of period two. On the other hand, if is now slowly decreased, after crossing , the stable period two cycle disappears; then the system oscillates on a stable period one cycle.

Moreover, in Figure 5, we can find cascade of period doublings , which can lead to chaos in some subregion of region 6 like the one shown in Figure 10. Quasiperiodic solutions can be destroyed by some homoclinic structures yielding chaos like the one shown in Figures 9(c) and 9(d).

In addition, a stick-slip solution of system (27) is given in Figure 11, suggesting that the stick-slip vibrations exist in the system. There are shorter slip phases between longer stick phases, and the stick-slip motion crosses the zero velocity point with a stop. The stick phases correspond to the accumulation of energy and the slip phases are corresponding to the release of energy. The time of energy accumulation (stick phase) is longer than the relaxation time (slip phase). It means that after longer time for the energy accumulation, the system reaches a critical state and then the energy is fully released, which is a characteristic of self-organized critical state.

To further investigate the influence of a periodic perturbation on the spring-block system, we present the numerical simulation of the discrete system (2). Note that is a linear transformation of , for the periodic perturbation , the discrete system (2) with a rate- and state-dependent friction law (4), (5) is solved numerically for different parameters under the periodic boundary conditions of , ,. Considering, at , the sliding speed of each blocks is unknown, so set as a random number which is smaller than . There are different sizes of sliding events, and some of them are shown in Figures 12(a)12(c) in the form of the sliding velocity as a function of positions and time . The parameters used in the numerical simulations are shown in Figure 12. For a given , the sliding speed of the blocks is lower for a smaller (Figures 12(a) and 12(b)). For a given , it demonstrates that the blocks slide with a lower sliding speed for a smaller (Figures 12(a) and 12(c)). The most distinctive feature in Figures 12(a)12(c) is that there is oscillation in sliding speed which is corresponding to the stick-slip vibrations of the system. From the numerical simulation, there are large sliding events when the parameters or are large, which means the sliding speed can be decreased by changing the parameter perturbation.

Moreover, based on the numerical results, the statistics of the sliding speed of the -th block shows a power-law distribution with the fitting exponent (Figure 12(d)). Using enables deducing that the stress drops, , also obey a power-law distribution, which means the self-organized critical behavior in the spring-block system. The system evolves to a self-organized critical state and then a collapse bursts with large sliding speed, which can be obtained in the real observed earthquakes based on the data analysis [3].

5. Conclusion

In the present work, we investigated a generalization of a Burridge-Knopoff spring-block model. It undergoes Hopf bifurcation and fold bifurcation of limit cycles for the system without parameter perturbation. Considering a small periodic fluctuation in the spring stiffness, the model with periodic perturbation is further investigated via a continuation technique and numerical bifurcation analysis. It suggests that periodic perturbation induces transition of the dynamics, the existence, the switch, and the coexistence of multiple attractors including periodic solutions with various periods, quasiperiodic solutions, and chaotic solutions through torus destruction or cascade of period doublings, which are the mechanism of the oscillations in the spring-block model. In addition, based on the numerical simulation, the statistics of the sliding speed of the i-th blocks shows a power-law distribution suggesting the self-organized critical state, which is a feature in the real observed earthquakes. It is significant to illustrate the dynamic behaviors such as chaos and self-organized state. When the system is chaos, the evolution of the system is unpredicted due to the solution’s sensitivity with respect to the initial condition, and rupture is unpredicted theoretically. For the self-organized critical state, the system is in a stable state and relatively controllable when the accumulated energy is less than the threshold. The system will collapse once the energy accumulated in the system exceeds a certain threshold

The present work gives a new approach to chaos, quasiperiodicity, and periodicity for the spring-block model. In [31], David R. Shelly examined an 8.5-year sequence of more than 900 recurring low-frequency earthquakes. These events exhibit tightly clustered recurrence, containing periodic, chaotic, and doubled earthquakes. Large nonrandom deviations from periodic behavior may interrupt sequences of nearly periodic events. In particular, even very small variation of the parameter can produce radical changes of the dynamics. The observed phenomenon is corresponding to the results in the present work. In addition, the numerical simulation shows that the sliding speed of the blocks can be deceased by reducing the parameter perturbation, which will facilitate the reduction of large seismic activity.

It showed that the burst of deterministic chaos could be observed, even for a very small perturbation, . This result is in contrast to the previous research in [32], which showed the occurrence of deterministic chaos in Burridge-Knopoff single-block model for a control parameter . In addition, Erickson et al. showed that chaos bursts as the number of blocks increases from 20 to 21 for [1]. Therefore, the present work also implies that the burst of chaos does not have to depend on the number of blocks, which is consistent with the previous research on the dynamics of spring-block model with time delay [11].

The results obtained in this study could be used in the development of landslide model when the system is subject to small external perturbation [33, 34], considering the sliding process could be described by a single sliding block moving along a rough surface. Kostić, S. et al. analyzed the dynamics of a single-block model and found that small perturbation of initial shear stress induces deterministic chaos [35]. When it comes to the periodic disturbances produced by the heterogeneity of the internal structure, an interesting problem arises, what is the difference between the chaos induced by the external condition or internal structures? The results in present work may provide a theoretical reference for this question.

This spring-block model could be related to transform faults, and it is also meaningful note that the research on such kind of spring-block Burridge-Knopoff model can provide useful clues for observed phenomena in seismology [31], which is the key to solve the scientific problem of whether earthquakes can be predicted or transformed. Meanwhile, the present work also provides a reference for the problem of plastic deformation in disordered materials, which is related to a spring-block model [3639]. Some experts studied the damage process and discussed a discrete model based on the potential energy of the block system [4042]. The models of point-to-face contact and edge-to-edge contact are presented which form a part of the contact theory [43, 44]. In the present model, the tangential spring is considered. If normal spring is also introduced, the system becomes much more complex. It remains an open and challenging question which needs further investigation.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research is supported by the National Natural Science Foundation of China Grants (nos. 11771407 and 11701521) and China Postdoctoral Science Foundation (no. 2018M632790).