Transitions between Metastable States in a Simplified Model for the Thermohaline Circulation under Random Fluctuations

TL;DR

This paper investigates how non-Gaussian Levy noise influences transitions between stable states in a simplified thermohaline circulation model, revealing that small noise can weaken circulation and larger noise accelerates this process.

Contribution

It introduces a novel analysis of non-Gaussian Levy noise effects on metastable state transitions in a thermohaline circulation model, using maximal likely pathways.

Findings

Small noise can weaken THC within a certain stability index range.

Larger noise intensity shortens the bifurcation time for weakened THC.

Non-Gaussian Levy noise significantly impacts transition dynamics in the model.

Abstract

In this work, we study the impact of non-Gaussian alpha-stable Levy motion on transitions between metastable equilibrium states (or attractors) in a stochastic Stommel two-box model for thermohaline circulation (THC). By maximizing the probability density of the solution process associated with a nonlocal Fokker-Planck equation, we compute maximal likely pathways and identify corresponding maximal likely stable equilibrium states. Our numerical results indicate weakened THC may be induced by perturbation with very small noise intensity in a certain range of stability index. Moreover, larger noise intensity and larger stability index induce weakened THC within shorter bifurcation time.

Equations

$$d{Y}_t=f(Y_t)dt+\sigma d{L}_t^{\alpha },Y(0)=Y_0,$$

$$\dot{T_e}$$

$$\dot{T_p}$$

$$\dot{S_e}$$

$$\dot{S_p}$$

$$\frac{\rho }{{\rho }_0}=1+{\beta }_S(S-S_0)-{\beta }_T(T-T_0),$$

$$Q(\Delta \rho )=\frac{1}{t_d}+\frac{q}{V}(\frac{\Delta \rho }{{\rho }_0}{)}^2,$$

$$\frac{d\Delta T}{dt}$$

$$\frac{d\Delta S}{dt}$$

$$dx$$

$$dy$$

$$dy=(F-y[1+{\mu }^2(1-y{)}^2])dt.$$

$$d{Y}_t=(\bar{F}-Y_t[1+{\mu }^2(1-Y_t{)}^2])dt+\sigma d{L}_t.$$

$${\nu }_{\alpha }(dz)=C_{\alpha }|z{|}^{-(1+\alpha )}dz,$$

$$d{Y}_t=f(Y_t)dt+\sigma d{L}_t^{\alpha },Y(0)=y_0,$$

$$A\phi (y,t)=f(y){\phi }^{\prime }(y,t)+{\int }_{{\mathbb{R}}^1\backslash \{0\}}[\phi (y+\sigma z,t)-\phi (y,t)]{\nu }_{\alpha }(dz)$$

$$=f(y){\phi }^{\prime }(y,t)+{\sigma }^{\alpha }{\int }_{{\mathbb{R}}^1\backslash \{0\}}[\phi (y+z,t)-\phi (y,t)]{\nu }_{\alpha }(dz),$$

$$p_t(y,t)=A^{\ast }p(y,t),p(y,0)=\delta (y-y_0),$$

$${\int }_{{\mathbb{R}}^1\backslash \{0\}}A\phi (y)\upsilon (y)dy={\int }_{{\mathbb{R}}^1\backslash \{0\}}\phi (y)A^{\ast }\upsilon (y)dy,$$

$$A^{\ast }\upsilon (y)={\sigma }^{\alpha }{\int }_{{\mathbb{R}}^1\backslash \{0\}}[\phi (y+z)-\phi (y)]{\nu }_{\alpha }(dz).$$

$$p_t=-(f(y)p(y,t){)}_y+{\sigma }^{\alpha }{\int }_{{\mathbb{R}}^1\backslash \{0\}}[\phi (y+z)-\phi (y)]{\nu }_{\alpha }(dz).$$

$$d{Y}_t=f(Y_t)dt+\sigma d{B}_t,Y(0)=y_0$$

$$p_t=-(f(y)p(y,t){)}_y+\frac{1}{2}{\sigma }^{2}p(y,t{)}_{yy},p(y,0)=\delta (y-y_0).$$

Full text

††journal: Applied Mathematics and Computation

1 Introduction

In the tremendous variety of nonlinear complex dynamical systems, among which the thermohaline circulation (THC) is one to mention, Gaussian as well as non-Gaussian noise play a pivotal role on setting up the dynamical behaviour of the system. This large-scale circulation in the ocean is decisive in the global climate change. The THC is driven by fluxes of heat and freshwater across the sea surface and posterior interior mixing of heat and salt [1]. Modeling the THC is a daunting puzzle due to the lack of, for instance, universal equation of state which connects the density of water to temperature and salinity and the complicated form of the domain which is bounded by the edges of various continents. Taking the system-level approach, i.e., as a system the ocean is just a reservoir filled with salt water, where the circulation is driven by the density (heat and salinity) difference, eases the impasse. The oceanic conveyor belt (or THC) carries warm surface water, around 17 Sv in volume and 0.5-1.5$\times 10^{15}W$ in heat on annual average, from low latitude (or equatorial) to high latitude (or polar) regime where it cools down, gets denser and sinks. The cold water at deeper levels, 2-3 km in depth, back from the polar to the equator regime as deep water-currents and surface water becomes more salty at the equator as water is removed by evaporation due to the greater heat which make it heavier and less salty thereupon it upwells. To examine this circulation Stommel pioneered a simple deterministic two-box model [2]. This idea has been further extended and improved through both conceptual and numerical models [3]-[5]. Actually box models are taken as the simplest setting to investigate stochastically forced systems. Climate records reveal incidence of abrupt climatic changes that might have some link with transition between stable equilibrium states of the THC [6]. Uncertain processes interrelated with random atmospheric fluctuations, challenges corresponding to unresolved scales and strange mechanisms worth consideration since they may lead to switching between the equilibrium states. It is customary to regard noisy fluctuations as Gaussian [7]-[10]. Velez-Belchi et al. in [6] verified stochastic resonance with additive noise induces transitions between the different stable states of the THC. A pathwise analysis of slowly varying not too large a perturbation on box model for THC has shown the system spends significant amount of time in metastable equilibrium [11]. Multiple stable states and possible critical transitions is discussed in [12].

Noisy fluctuations, meanwhile, are demonstrated to be non-Gaussian in some complex systems such as spectral analysis of paleoclimatic data [13], metastability in climate systems [14], and bursty transition in gene expression [15]. The effect of Gaussian noise on THC has been well dealt with, while the influence of non-Gaussian stable Lévy noise on the dynamical behaviour of this system is scant. It is more expedient to model these random fluctuations which portray haphazard jumps, discontinuity, heavy tail distribution and bursting sample paths by a non-Gaussian ($\alpha$-stable) Lévy motion. With respect to $\alpha$-stable Lévy noise, Wang et al. in [16] determined low variability salinity difference level domain by numerical simulation of mean exit time.

In this paper, we focus on the influence of non-Gaussian $\alpha$-stable Lévy noise on stochastic box model for THC by investigating maximal likely attractors of the maximal likely pathways, i.e., evolution of the maximizer of probability density function as time increases in stochastic phase portraits. The subtlety of phase portraits for stochastic differential equations (SDE) can be supplemented by reconveying these geometrical entities in terms maximal likely pathways [17], [18]. We will consider a simple model for THC as a scalar SDE in the following form:

[TABLE]

where $f(Y_{t})$ is a deterministic vector field, $\sigma$ is the noise intensity, $\alpha$ is stability index and {$L_{t}$}, $t\geq 0$ is a scalar non-Gaussian Lévy motion.

The paper is structured as follows. Section 2 presents details of a stochastic Stommel model. Section 3 reviews maximal likely stable equilibrium states or maximal likely attractors. Section 4 discusses numerical results of computations on maximal likely pathways under different parameters, viz, stability index, noise intensity and (nondimensional) freshwater flux.

2 Model and Method

2.1 Model

An equatorial box and polar box with temperature $T_{e}$ and $T_{p}$, and salinity $S_{e}$ and $S_{p}$ respectively having the same volume $V$ are connected by advective flow and barter heat and freshwater with the atmosphere. Moreover, each container is individually forced at its boundary. The equatorial box bounces according to a time scale $t_{r}$ back to local atmospheric temperature $T_{a}$ = $T_{0}+\frac{\theta}{2}$, where $\theta$ is the equator-to-pole atmospheric temperature difference and $T_{0}$ is a reference temperature. The polar box retains local atmospheric temperature, for the sake of symmetry, $T_{a}$ = $T_{0}-\frac{\theta}{2}$. Freshwater is dispatched from equatorial into polar vessel by a laid down freshwater flux $\frac{F_{s}}{2}$ [19]. A sketch of a variant of Stommel’s model is shown in Figure 1.

Following [7], we define the heat and salinity balances by the equations:

[TABLE]

[TABLE]

Table 1. Parameters of the stochastic Stommel box model [19]

Density differences are essential for the circulation to take place. The density $\rho$ of a box, which is affected by $T$ and $S$ in an opposite way, is approximated by a linear equation of state

[TABLE]

where $\beta_{T}$ and $\beta_{S}$ are (positive) constant thermal expansion and haline contraction coefficients, respectively. We may choose, for simplicity, a positive transport function which is independent of the direction of $\Delta\rho$. This is viable because in the closed advection of water between two vessels ($Q$), the quantity of equatorial water moving into the poles and vice versa keeps unaltered in both directions of circulation. Thus,

[TABLE]

where $t_{d}$ is diffusive mixing time scale between the two vessels that would occur when there is no density difference. Subtracting the equations in (2.1) and defining $\Delta S\equiv S_{e}-S_{a}$ and $\Delta T\equiv T_{e}-T_{a}$ yields a coupled pair of equations for the time evolution of temperature and salinity differences between the vessels:

[TABLE]

With the scales $\Delta T\equiv x\theta,\Delta S\equiv\frac{y\beta_{T}\theta}{\beta_{s}}$ and time scaled with $t_{d}$, we recast (2.2) as the following coupled dimensionless system of equations.

[TABLE]

Where $\beta=\frac{t_{d}}{t_{r}}$ measures temperature restoring tensility, $\mu^{2}=\frac{qt_{d}(\beta_{T}\theta)^{2}}{V}$ is strength of the buoyancy-driven convection between the two boxes relative to the diffusive mixing. The parameter $F={\frac{\beta_{S}S_{0}t_{d}}{\beta_{T}\theta H}}F_{S}$ represents dimensionless freshwater forcing. The physical meaning of $\beta\gg 1$ is that mixing between the boxes temperatures is sluggish compared to the rapid parity of each box’s temperature with its local temperature forcing. Salinity leads to non-linearity which causes the existence of multiple equilibria and thresholds in the THC.

Due to the dominance of the $\beta$ term in equations (2.3) the temperature difference $x$ changes very little, and we can thus take $x$ as $x=1+\mathcal{O}(\beta^{-1})$. This leaves us with an ODE in $y(t)$,

[TABLE]

Suppose that $F=\bar{F}$ is constant. The time evolution of $y$ can be represented by a potential function $V(y)$ as $dy=-V^{\prime}(y)dt$ where $V(y)=-{\bar{F}}y+\frac{1+\mu^{2}}{2}y^{2}-{\frac{2\mu^{2}}{3}y^{3}}+\frac{\mu^{2}}{4}y^{4}$. The strength of the transport between boxes is governed by density difference $|\Delta\rho|=\rho_{0}|\beta_{S}\Delta S-\beta_{T}\Delta T|$ which in dimensionless form is $|y-x|$. In certain parameter ranges $V(y)$ has two stable states: one with larger northward heat transport, usually referred to as the on-state, where the salinity difference $y$ is small, and other with weak (or even reversed) circulation, where the difference $y$ is large, called the off-state, depending on differences in the initial conditions only. After Cessi’s valuation of parameter values [7], [19] $\bar{F}=1.1$, $\mu^{2}=6.2$, $V(y)$ is a double-welled potential with two stable minima at $y_{-}=0.24$ and $y_{+}=1.07$ and an unstable maximum at $y_{u}=0.69$ as can be seen in Figure 1.

Numerical simulations of the THC disclose its deep transport, stability and variability delicacy , especially with respect to freshwater flux [21],[22]. A multiplicity of models has displayed that the competition between thermal and saline forcing results in multiple equilibria. Clearly, the deterministic system has two stable steady states and the random driving moves the temperature and the salinity away from initial equilibrium. So a perturbation is certainly needed to induce a transition between both states. These perturbations are provided by background noise such as the freshwater flux. The fluctuations in the freshwater flux or in the external salinity inputs lead $F$ to vary following the difference between salinity increment rate at the equator and corresponding decline at north pole. Further, $F$ may be parameterized as a sum of time independent mean component $\bar{F}$ and a noisy fluctuating process $\frac{dL_{t}}{dt}$. As a result, we come out with the SDE:

[TABLE]

2.2 $\alpha$-stable Lévy process

A stable distribution $S_{\alpha}(\mu,\beta,\gamma)$ is the distribution for a stable random variable [17], where the parameters $\alpha\in(0,2),\mu\in[0,\infty),\beta\in[-1,1]$ and $\gamma\in(-\infty,\infty)$ respectively are indices of stability, scale, skewness and shift.

An $\alpha$-stable scalar Lévy motion $L_{t}^{\alpha},\alpha\in(0,2)$ is a non-Gaussian stochastic process satisfying the conditions:[17], [23]-[25]

i) $L_{0}^{\alpha}=0$, almost surely;

ii) Independent increments: For each $n\in\mathbb{N}$ and each $0\leq t_{1}<t_{2}<\cdots<t_{n-1}<t_{n}<\infty$, the random variables $L_{t_{2}}^{\alpha}-L_{t_{1}}^{\alpha},\cdots,L_{t_{n}}^{\alpha}-L_{t_{n-1}}^{\alpha}$ are independent;

iii) Stationary increments: $L_{t}^{\alpha}-L_{s}^{\alpha}$ and $L_{t-s}^{\alpha}$ have the same distribution $S_{\alpha}((t-s)^{\frac{1}{\alpha}},0,0)$;

iv) Stochastically continuous sample paths: Sample paths are continuous in probability. For all $\delta>0$, all $s\geq 0;\mathbb{P}(|L_{t}^{\alpha}-L_{s}^{\alpha}|)>\delta)\to 0$ as $t\to s$.

The jump measure of $L_{t}^{\alpha}$ is

[TABLE]

where $C_{\alpha}=\frac{\alpha}{2^{1-\alpha}{\sqrt{\pi}}}\frac{\Gamma(\frac{1+\alpha}{2})}{\Gamma(1-\frac{\alpha}{2})}$ is stability constant.

Brownian motion, $B_{t}$, is a symmetric 2-stable process. $B_{t}$ has independent and stationary increments, continuous sample paths almost surely and has normal distribution $\mathcal{N}(0,t)$.

2.3 Method

We will consider a simple model for THC as a scalar SDE with additive noise, i.e. the intensity of the noise is independent of the state of the THC, in the form

[TABLE]

where $f(Y_{t})=\bar{F}-Y_{t}[1+\mu^{2}(1-Y_{t})^{2}]$ is a deterministic vector field, $\sigma$ is the noise intensity and {$L_{t}^{\alpha}$}, $t\geq 0$, is a symmetric scalar $\alpha$-stable non-Gaussian Lévy motion.

The generator A of the solution process $Y_{t}$ for the SDE (2.6) with triplet $(0,0,\nu_{\alpha})$ is

[TABLE]

where $\varphi(y)\in C(\mathbb{R}^{1})$ and the Fokker-Planck equation of (2.6) in terms of the probability density function $p(y,t)$ for the solution process $Y_{t}$ given initial condition $Y_{0}=y_{0}$ is [17]:

[TABLE]

here, $\delta$ represents the dirac function and $A^{*}$, the adjoint operator of $A$ in the Hilbert space $L^{2}(\mathbb{R}^{1})$, can be computed by solving

[TABLE]

for $\varphi,\upsilon$ in the domains of definition for the operators $A$ and $A^{*}$, respectively,

[TABLE]

Consequently, we obtain the following nonlocal Fokker-Planck equation

[TABLE]

The SDE

[TABLE]

is brought in if Brownian motion, which has a Fokker-Planck equation having the form of local partial differential equation (2.11), substitutes the Lévy motion.

[TABLE]

Numerical finite difference method as in [26] and standard difference method are used in simulating equations (2.9) and (2.11) respectively.

3 Maximal Likely Trajectories

Phase portraits provide geometric pictures for lower-dimensional deterministic dynamical systems. At a given time instant $t$, the maximizer of $y_{m}(t)$ of the probability density function $p(y,t)$ yields the most probable location of this trajectory at time $t$. The deterministic entity $y_{m}(t)$ follows the uppermost crest of the surface in the $(y,t,p)$-space leaving behind trace out of the so called the maximal likely trajectory starting at $y_{0}$, as time goes on.

A state that attracts all nearby states is called a maximal likely stable equilibrium state. Maximal likely stable equilibrium states are dependent on non-Gaussianity index $\alpha$, noise intensity $\sigma$ and freshwater flux $\bar{F}$. The most probable phase portrait [17], [18], an extension of phase portraits to envision stochastic dynamics, is the state space with representative maximal likely pathway including stable equilibrium states. Most probable phase portraits like the phase portraits are geometric entities. They depict sample pathways and stable equilibrium states that are maximal likely.

We examine the number and location of maximal likely stable equilibrium salinity difference states as a parameter varies. We may define bifurcation time as the time between the change in number of maximal likely equilibrium states. It is a time scale for the birth of a new most probable stable equilibrium state.

4 Result and Discussion

In this section, we will discuss the number and value of maximal likely stable equilibrium states for the stochastic THC model (2.6). For simplicity, we consider four maximum likely trajectories with initial points $y_{0}=0.2,0.4,0.8,1.5$, to both sides of the equilibria in the deterministic dynamical system. Due to the lengthy computation process time is capped to $T=50$.

Next we examine maximal likely trajectories when system parameters change.

As Lévy noise parameters vary:

The deterministic dynamical system is bistable in some range of freshwater flux parameter that contains $\bar{F}=1.1$. A perturbation with very small noise intensity on the deterministic system invokes an interesting phenomenon as the Lévy noise index varies.

For $\sigma=0.001$, in line with the existence of a transient state the interval of alpha stability can be partitioned in two.

i) $0.1\lessapprox\alpha<0.7$: In this interval of stability index, we have low salinity difference transient maximal likely stable equilibrium state. The transient state terminates since maximal likely salinity difference trajectories in the stochastic THC model merge by jumping to maximal likely metastable equilibrium state as time counts on. Bifurcation time is a time interval for the emergence of a new maximal likely stable equilibrium state. This time scale predicts transition from low salinity difference maximal likely state (or strong THC) to high salinity difference maximal likely state (or weak THC). Jump from low salinity

difference to high salinity difference maximal likely equilibrium state comes out at different bifurcation time as index of stability in the interval varies. As can be seen from figure 3, bifurcation time attains its minimum at $t\thickapprox 15.4683$ when $\alpha\thickapprox 0.15$. This duration should be multiplied by the diffusion time, taken here to be $t_{d}=180$ years as in Table (1). Corresponding dimensional value of this minimum time is therefore about 2,784 years. As shown in figure 4, bifurcation time decreases, but increases after it attains minimum as $\alpha$ increases in this interval. In figure 3 bifurcation time is at $t\thickapprox 22.5279$ (about 4,055 years) when $\alpha=0.3$. The transient maximum likely equilibrium state lasts for relatively short time interval before it ends jumping up to the maximal likely metastable equilibrium state. Next, we turn our focus on discussing why this jump takes place. An $\alpha$-stable process moves predominantly by big jumps as $\alpha$ nears to 0, and small jumps override as $\alpha$ is approaches 2. The jumps that occur in the maximal likely trajectories due to the Lévy noise perturbation in this interval are seldom but big in size, as a result the salinity difference trajectories originating from a small neighbourhood of the deep well have a substantial probability of bypassing the barrier wall to land in the shallow well. It can also be clearly seen that the bifurcation time increases or equivalently transient maximal likely stable equilibrium state lasts longer, with the increase in $\alpha$ bound to the interval under consideration. Since $y$ is dimensionless salinity difference, it is proportional to the salinity difference between the equatorial and polar regions, higher $y$ values imply reduced density difference.

Thus, higher values of $y$ may consequently lead to diminished THC. This implies noise with such a small noise intensity and some range of values of stability indices may invoke weakened THC.

For the purpose of comparison, the influence of Gaussian noise with the same noise intensity is shown in figure 3. The stochastic system (2.10) is bistable. The reason for transition in the Lévy noise case may be attributed to the jump in the evolution of salinity difference trajectories.

ii) $0.7\lessapprox\alpha<2$: Our numerical results display the stochastic Stommel two-box model for THC in the SDE (2.6) has two maximal likely stable equilibrium states. The stable equilibrium states in the absence of noise match with these maximal likely stable states after perturbation. The stochastic system in this interval of stability index is, therefore, not affected by the noise. In figure 3, for instance, the two most likely stable salinity difference are at 0.24 and 1.07. An implication of this phenomenon in the perturbed THC model is the possibility that there is neither a shift in potential wells from the original location nor a state transition. Intensity of the noise is not large enough to stir a slight change in the maximal likely equilibrium states and the jumps occurring in the maximal likely trajectories are lower than the height of the barrier wall between the deep potential well and the shallow potential well (see figure 2).

Numerical experiments reveal the range of stability index where maximum likely transient equilibrium state prevails gets wider as noise intensity increases from $0.001$ to $0.1$ as shown in 4. In figure 4, the system is monostable when stability index is below 0.7 whereas the system is bistable when stability index exceeds 1.9 in figure 4. Mean while, the bifurcation time is affected adversely by this noise intensity increment. The minimum bifurcation time decreases from 15.4683 (about 2,784 years) in figure 4 to 7.8275 (about 1,409 years) in figure 4. Notice that the stability index as well increases from 0.15 to 0.5. This result hints that the transition from weak to strong THC may be enhanced as the Lévy noise parameters increase.

As value of $\bar{F}$ varies:

$\bar{F}$ is of particular interest as it monitors the strength of the freshwater forcing and determines the stable states and transitions between stable states of THC. In the absence of noise, the model is bistable in the interval $0.96\lessapprox\bar{F}\lessapprox 1.3$, otherwise there is only one stable state: on-state for $\bar{F}\lessapprox 0.96$ or off-state when $\bar{F}\gtrapprox 1.3$. It is desirable to examine how this stochastic THC model is influenced as freshwater flux parameter varies. We consider very small values of the parameters $\alpha=0.05$ and $\sigma=0.01$ and the same initial conditions. When the value of the freshwater flux parameter varies while the system is under the influence of non-Gaussian Lévy noise, our numerical results show that, for $\bar{F}=1.1$ there are two maximal likely stable equilibrium states at 0.24 and 1.07.

It is worth noting that maximal likely stable equilibrium states of the stochastic system under the given freshwater flux parameter, stability constant and noise intensity remain the same as that of the deterministic system. When $\bar{F}=0.9$, as in figure 5, the stochastic system has only one maximal likely stable equilibrium state at 0.17. Similarly, at 0.185 for $\bar{F}=0.95$. Actually in these, the stochastic system is monostable and these maximal likely low salinity difference states correspond to strong THC. Existence of transient maximal likely equilibrium state is an interesting event that takes place when $\bar{F}=1$, and $\bar{F}=1.15$. In figure 5, the transient maximum likely high salinity difference equilibrium state 1 merges to the maximal likely low salinity difference stable equilibrium state 0.205 at the bifurcation time $t=34.9$ (about 6,282 years) when $\bar{F}=1$. For $\bar{F}=1.15$, the transient maximal likely equilibrium state is at 0.27 and the maximal likely metastable equilibrium state is at 1.09. The jump is from high to low salinity difference state and it takes place at bifurcation time $t=21.8606$ (about 3,934 years) as in figure 5. When there is a transient state, numerical simulation results show bifurcation time increases and the transient maximal likely stable equilibrium state is low salinity difference state or strong THC which finally merges to high salinity difference state or weak THC when $\bar{F}<1.1$, while if $\bar{F}>1.1$, bifurcation time decreases and the transition is reversed. The system has maximum likely high salinity difference equilibrium state (or weak THC) 1.17 if $\bar{F}=1.4$ as shown in 5.

Abrupt climate changes may have some degree of relationship with sudden fluctuations of the THC due to vast freshwater flux caused by break up and melting of iceberg and glaciers. This conveyor belt transfers huge amount of heat from low to high latitudes regulating the global climate by redistributing heat across the planet. Weak THC may have far reaching consequences in water sources, precipitation for farming and well-being of low latitude ecosystems. On the other hand, it may lead to cooler ocean temperature in certain pockets of high latitudes and warmer temperature in other pockets with potential to extra sea level rise along coasts and ascends sea temperature adversely affecting the oceanic ecology.

Acknowledgements

The authors would like to thank Xiujun Chen, Hui Wang, Xiaoli Chen and Ziying He for helpful discussions. This work was partially supported by NSFC grant 11801192.