Rogue waves, self-similar statistics, and self-similar intermediate asymptotics

TL;DR

This paper develops a statistical theory explaining how rogue waves and extreme events emerge as self-similar giant fluctuations in nonlinear wave systems, supported by analytical and numerical evidence.

Contribution

It introduces a novel self-similar statistical framework for rogue wave emergence in nonlinear systems with gain, advancing understanding of extreme event formation.

Findings

Rogue waves appear as parabolic giant fluctuations in self-similar regimes.

The non-Gaussian statistics of rogue waves are analytically characterized.

Numerical simulations validate the self-similar structure of extreme events.

Abstract

We advance a statistical theory of extreme event emergence in random nonlinear wave systems with self-similar intermediate asymptotics. We show, within the framework of a generic (1 + 1)D nonlinear Schrodinger equation with linear gain, that extreme events and even rogue waves in weakly nonlinear, statistical open systems emerge as parabolic-shape giant fluctuations in the self-similar asymptotic propagation regime. We analytically demonstrate the self-similar structure of the non-Gaussian statistics of emergent rogue waves and validate our results with numerical simulations. Our results shed new light on generic statistical features of extreme events in nonlinear open systems with self-similar intermediate asymptotics.

Equations

$$i{\partial }_{Z}U-\frac{1}{2}{\beta }_2{\partial }_{TT}^{2}U+\gamma |U{|}^{2}U-\frac{i}{2}GU=0,$$

$$\Psi =U/\sqrt{⟨P_0⟩};t=T/T_p,z=Z/L.$$

$$i\sigma {\partial }_z\Psi -\frac{1}{2}{\sigma }^2{\partial }_{tt}^2\Psi +|\Psi {|}^2\Psi -\frac{i}{2}\sigma g\Psi =0,$$

$$\Psi =A{e}^{i\varphi /\sigma },$$

$$A(t,z)=a(z)f[t/t_p(z)],$$

$$\varphi (t,z)={\varphi }_0(z)+c(z)t^2,$$

$$a(z)=a_0{e}^{gz/3},$$

$$t_p(z)=a_0^{-2}{e}^{gz/3},$$

$$a_0=\frac{1}{2}(\sqrt{2}g{w}_0{)}^{1/3},w_0=W_0/⟨W_0⟩,$$

$$⟨W_0⟩={\int }_{-\infty }^{\infty }dT⟨|U(T,0){|}^2⟩.$$

$$f(\eta)=\left\{\begin{array}[]{cc}\sqrt{1-\textstyle\frac{g^{2}}{18a_{0}^{6}}\eta^{2}},&|\eta|\leq 3\sqrt{2}a_{0}^{3}/g,\\ 0,&|\eta|\geq 3\sqrt{2}a_{0}^{3}/g.\end{array}\right.$$

$${\Gamma }_0(t_1,t_2)=⟨{\Psi }^{\ast }(t_1,0)\Psi (t_2,0)⟩,$$

$${\Gamma }_0(t_1,t_2)=\frac{1}{\sqrt{\pi }}\exp \left(-\frac{t_1^2+t_2^2}{2}\right)\exp \left[-\frac{(t_1-t_2{)}^2}{2{\sigma }_c^2}\right].$$

$$\Psi (t,0)=\sum _{n=0}^N{c}_n{\psi }_n(t),$$

$${\psi }_n(t)={\left(\frac{2\xi }{\pi }\right)}^{1/4}{\left(\frac{1}{2^{n}n!}\right)}^{1/2}{H}_n(\sqrt{2\xi }t)e^{-\xi t^2}.$$

$$\alpha =1/2,\beta =1/(2{\sigma }_c^2),$$

$$\xi =\sqrt{{\alpha }^2+2\alpha \beta }.$$

$${\int }_{-\infty }^{\infty }dt{\psi }_m^{\ast }(t){\psi }_n(t)={\delta }_{mn}.$$

$$⟨c_n^{\ast }{c}_m⟩={\lambda }_n{\delta }_{mn},$$

$${\lambda }_n=A{\left(\frac{\beta }{\alpha +\beta +\xi }\right)}^n.$$

$$1={\int }_{-\infty }^{\infty }dt⟨|\Psi (t,0){|}^2⟩=\sum _{n=0}^N⟨|c_n{|}^2⟩=\sum _{n=0}^N{\lambda }_n,$$

$${\lambda }_n={\left(\frac{\beta }{\alpha +\beta +\xi }\right)}^n\frac{1-\frac{\beta }{\alpha +\beta +\xi }}{1-{\left(\frac{\beta }{\alpha +\beta +\xi }\right)}^{N+1}}.$$

$$c_n=\sqrt{E_n}{e}^{i{\varphi }_n},$$

$$\mathcal{P}(E_n,{\varphi }_n)=\frac{\theta (E_n)}{2\pi {\lambda }_n}{e}^{-E_n/{\lambda }_n}.$$

$$w_0={\int }_{-\infty }^{\infty }dt|\Psi (t,0){|}^2=\sum _{n=0}^N|c_n{|}^2=\sum _{n=0}^N{E}_n.$$

$${\chi }_{{\rho }_n}(s)=⟨e^{is{E}_n}⟩=\frac{1}{1-is{\lambda }_n}.$$

$${\chi }_{w_0}(s)=\prod _{n=0}^N\frac{1}{1-is{\lambda }_n}.$$

$$\mathcal{P}(w_0)={\int }_{-\infty }^{\infty }\frac{ds}{2\pi }{e}^{-is{w}_0}\prod _{n=0}^N\frac{1}{1-is{\lambda }_n}.$$

$$\mathcal{P}(w_0)={\int }_{-\infty }^{\infty }\frac{ds}{2\pi }{e}^{-is{w}_0}\prod _{n=0}^N\frac{i/{\lambda }_n}{s+i/{\lambda }_n}.$$

$$\mathcal{P}(w_0)=\sum _{m=0}^N\frac{1}{{\lambda }_m}{e}^{-w_0/{\lambda }_m}\prod _{\underset{n\ne m}{n=0}}^N\frac{{\lambda }_n^{-1}}{{\lambda }_n^{-1}-{\lambda }_m^{-1}}.$$

Full text

I Introduction

Rogue waves (RW), extremely rare, giant-amplitude waves obeying non-Gaussian statistics, were originally discussed in the oceanographic context Ono ; Zakh ; RW-rev . The concept has been quickly recognized as germane to generic wave supporting physics settings and RWs have been discovered, among others, in supercontinuum generating optical fibers Solli ; SCdud1 , optical cavities Res1 , Bose-Einstein condensates Blud , Raman fiber amplifiers Ham1 ; Ham2 , fiber lasers RWlas1 ; RWlas2 , laser filamentation Kasp , plasmas Mos , stimulated Raman scattering Yash1 ; Yash2 , discrete nonlinear lattices latt , and even in the multimode optical fibers and microwave transport in the linear propagation regime RWf ; Hell .

Although there apparently exists no universal mechanism describing RW generation in any physical system Dud-rev , many weakly nonlinear and dispersive statistical wave systems are governed by a generic $(1+1)$D nonlinear Schrödinger equation (NLSE) Agra . The RW excitation in the NLSE model with random input wave fields has been studied in the anomalous dispersion regime both numerically Pic ; Dud ; Aga ; Akh1 ; Akh2 ; Sur1 ; Sur2 and experimentally Sur1 ; Sur2 . These studies revealed modulation instability driven RW excitation scenarios in weakly nonlinear conservative systems and elucidated the respective roles of spontaneous Peregrine-like breather excitation from a noisy environment and of random soliton collisions in triggering the emergence of heavy-tailed probability density distributions (PDF) of wave intensities. Such heavy-tailed PDFs herald the RW generation in the system Dud ; Sur1 ; Sur2 ; Akh1 ; Akh2 .

At the same time, open physical systems often cannot support either solitons or breathers, at least as their long-term asymptotic states, because the energy supply from—or loss to—the environment precludes the establishment of precise balance between the dispersion and nonlinearity necessary for soliton formation. Although dissipative soliton formation is possible in some open systems if, on the one hand, the nonlinearity is balanced by dispersion/diffraction and, on the other hand, gain is balanced by loss Dissol , a multitude of open systems across physics disciplines exhibit self-similar dynamics instead. Examples range from blast waves in gas dynamics and turbulent bursts in fluids Bar to nonlinear waves in optical fiber Krug1 ; Krug2 ; Ser and graded-index waveguide PSA2 ; PSA3 amplifiers, saturable two-level absorbers PSA4 , and growing Bose-Einstein condensates Drum ; Vog . Moreover, there exists a wide class of open wave systems displaying self-similar evolution in the intermediate range of parameters such that particular initial conditions at the source no longer play any role, though the system has not yet reached its steady-state Bar ; Krug1 ; Krug2 ; PSA4 ; Drum . This observation prompts a natural question: Is there an universal scenario of extreme event generation in the statistical wave systems with self-similar intermediate asymptotics? A related fundamental issue has to do with the influence of self-similar dynamics on the wave ensemble statistics in the self-similar evolution regime.

In this work, we take the first step toward addressing these fundamental topics by advancing a statistical theory of extreme events in weakly nonlinear random wave systems with unsaturated gain, described within the framework of a generic NLSE modified by a linear gain term. The modified NLSE possesses self-similar intermediate asymptotics with a parabolic intensity profile in the normal dispersion regime Krug2 . We develop a statistical theory of RW generation in the system by studying random input wave propagation there. We analytically derive and numerically verify the PDF of a wave peak power ensemble, establish its non-Gaussian statistics, and demonstrate the self-similar evolution of the ensemble statistics on random pulse propagation in the intermediate regime. We stress that, to our knowledge, this is the first demonstration of the RW statistics self-similarity in any nonlinear random wave system. Our analytical results are independent of a particular source ensemble model. In addition, they are in excellent agreement with numerical simulations of the modified NLSE.

As the NLSE with linear gain captures salient features of optical wave propagation in any weakly nonlinear amplifying media Krug1 ; Krug2 , our findings are expected to be generic. We also note that since the NLSE with a linear gain term and harmonic trapping potential attains a self-similar asymptotics as well Drum , our results apply, at least qualitatively, to optical waves in graded-index waveguide amplifiers and matter waves in growing Bose-Einstein condensates.

This work is organized as follows. In Section I, we introduce a generic dimensionless NLSE with linear gain as our mathematical model of self-similar dynamics in a wide variety of nonlinear wave systems and briefly review its parabolic self-similar solutions. In Section II, we formulate our statistical source ensemble model and present numerical results for a typical ensemble member evolution. In Section III, we present an analytical derivation of the wave peak power PDF in the self-similar regime and verify our findings with comprehensive numerical simulations of the modified NLSE. We present our conclusions in Section IV.

II Modified dimensionless nonlinear Schrödinger equation and its parabolic self-similar solutions

We consider statistical pulse propagation in a nonlinear fiber amplifier in the normal dispersion regime governed by the modified NLSE in the form

[TABLE]

where $G$ is linear gain; $\beta_{2}$ is a group-velocity dispersion, and $\gamma$ is a Kerr nonlinearity coefficients. To explore generic features of random waves with self-similar intermediate asymptotics, which are independent of the source and medium particulars, it will prove convenient to work with dimensionless variables, defined as

[TABLE]

Here $T_{p}$ is an average pulse width and ${\langle}P_{0}{\rangle}$ is an average peak power of a statistical input pulse ensemble; $L=\sqrt{L_{NL}L_{D}}$, where $L_{NL}=(\gamma{\langle}P_{0}{\rangle})^{-1}$ and $L_{D}=T_{p}^{2}/\beta_{2}$ are the usual nonlinear and dispersion lengths. In the dimensionless variables, the modified NLSE reads

[TABLE]

where we introduced a dimensionless gain $g=GL$ and soliton $\sigma=\sqrt{L_{\mathrm{NL}}/L_{\mathrm{D}}}$ parameters which entirely determine the system dynamics given the source coherence state. We note that the smaller the soliton parameter, the faster the self-similarity is attained Krug1 .

We now present a parabolic self-similar solution to Eq. (3) by transforming the original results of Krug1 to a dimensionless form. To this end, we express the parabolic self-similar solution (similariton) in the polar form

[TABLE]

where the amplitude and phase are given by the expressions

[TABLE]

and

[TABLE]

respectively. Henceforth we are mainly focusing on the similariton amplitude as we will derive the similariton peak power PDF; thus we ignore the phase. In Eq. (5), $a(z)$ and $t_{p}(z)$ are dimensionless amplitude and width of the similariton, given by

[TABLE]

and

[TABLE]

respectively. Here

[TABLE]

where the dimensionless input pulse energy $w_{0}$ is normalized to the average energy of the input pulse ensemble,

[TABLE]

Here the angle brackets denote ensemble averaging. Further, introducing the dimensionless similarity variable, $\eta=t/t_{p}(z)$, we can express the similariton profile as

[TABLE]

It follows at once from Eqs. (5) through (11) that the similariton power profile has a parabolic shape which, together with its amplitude and width, is completely determined by the (scaled) input pulse energy $w_{0}$ as well as the medium gain $g$. To study the individual pulse evolution numerically, we have to specify a random pulse ensemble at the source.

III Statistical ensemble formulation of input pulses

We now describe the input pulse ensemble in terms of a generic Gaussian Schell model  MW previously employed in extreme event studies Yash1 ; Yash2 . The(GSM) ensemble has a Gaussian average intensity and Gaussian degree of the second-order temporal coherence MW ; Laleh . The GSM ensemble mutual intensity, defined as

[TABLE]

reads then

[TABLE]

Here $\sigma_{c}=T_{c}/T_{p}$ is a source coherence parameter. Given $\sigma_{c}$, we can approximate any realistic GSM source with a finite number $N=N(\sigma_{c})$ of (uncorrelated) excited coherent modes via the Karhunen-Loève expansion,

[TABLE]

where $\{c_{n}\}$ are complex random amplitudes and $\{\psi_{n}(t)\}$ are coherent mode functions, known for a GSM source to be MW

[TABLE]

Here $H_{n}(x)$ is a Hermite polynomial of the order $n$, and we introduced the notations

[TABLE]

and

[TABLE]

We note that the mode functions are orthonormal such that

[TABLE]

The second-order statistics of the (uncorrelated) random amplitudes $\{c_{n}\}$ are determined by

[TABLE]

with the modal weight distribution,

[TABLE]

The set of $\{\lambda_{n}\}$ determines average energies of the excited coherent modes representing the source. The constant $A$ is determined by a relevant physical normalization; in our case, the (scaled) average energy ${\langle}w_{0}{\rangle}$ of the pulse must be normalized to unity:

[TABLE]

implying that

[TABLE]

Note, the case $N=0$ implies there are no excited modes, yielding $\lambda_{n}=\lambda_{0}=1$ as expected; that is, only a ground mode of the source is excited—this is an example of a fully (temporarily) coherent source at the second order.

To specify the source ensemble beyond the second order, we express the set of $\{c_{n}\}$ in the polar form

[TABLE]

and assume the following joint PDF of the phases and energies

[TABLE]

Here $\theta(x)$ is a Heaviside step function; in other words, the phases are uniformly distributed in the interval, $-\pi\leq\phi_{n}\leq\pi$, and the mode powers obey the Rayleigh distribution. Eq. (24) guarantees Gaussian statistics of the input pulse ensemble $\{\Psi(0,t)\}$. We stress that coherent mode amplitude fluctuations are crucial as they determine the total source energy fluctuations, which, in turn, shape the RW statistics in the self-similar intermediate propagation regime.

We performed numerical simulations of a GSM ensemble of $10^{4}$ random pulses. In our numerical simulations we take $T_{p}=0.2$ ps, ${\langle}P_{0}{\rangle}=2.8$ W, $\gamma=5.8$ W*-1* m*-1*, $\beta_{2}=0.025$ ps2m*-1*, and $G=1.9$m*-1*, implying that the system is in the nonlinearity dominated regime, $\sigma\ll 1$ corresponding to the experiment Krug1 , except the average input pulse ensemble power and fiber nonlinearity are scaled by three orders of magnitude down and up, respectively. Thus, the nonlinear length remains the same and the parabolic self-similar regime is reached at sufficiently short amplifier lengths. The reason behind the scaling is practical difficulty to generate high-power (kW) ultrashort pulse sources with thermal-like power fluctuations, whereas engineering such a statistical source with the power around 1 W has been recently reported Sur1 . On the other hand, fiber amplifiers with such high nonlinearities can, in principle, be realized, for instance, by dye doping high nonlinear refractive index liquids, filling the cores of specially designed photonic crystal fibers of the kind reported in c.f. Giess ; Yash3 . Alternatively, one can attain the self-similar regime with statistical beams propagating in highly nonlinear chalcogenide planar waveguide amplifiers with defocusing nonlinearities of the order of $\gamma\simeq 10$ W*-1m-1* Eggl .

A particular ensemble realization evolution is illustrated in Fig.1. We can infer from the figure that the power profile of the realization attains a parabolic shape at the distance $Z=5.5$m; this conclusion holds for any realization, albeit the propagation distance over which the parabolic profile is reached varies from realization to realization. Most importantly, as soon as the parabolic profile is reached, it remains unchanged, up to scaling, indicating the self-similar regime ensues.

IV Self-similar statistics of peak power PDF of the pulse ensemble

The self-similar dynamics of each realization results in a remarkably simple statistical evolution of the ensemble as a whole which can be revealed by analytically deriving its peak power PDF. To this end, we can find the total energy distribution of the source pulse ensemble. As the system dynamics is entirely determined by the source energy $w_{0}$, we derive a generic source energy PDF without committing to a specific source model under the only assumptions that (i) the mode powers have thermal-like distributions of Eq. (24) and (ii) we know the set of average mode energies $\{\lambda_{n}\}$. The source energy is then given by the expression

[TABLE]

As all modes are uncorrelated, we can use characteristic functions to determine $\mathcal{P}(w_{0})$. To this end, we first determine a characteristic function of $E_{n}$,

[TABLE]

It follows at once that

[TABLE]

Taking an inverse Fourier transform of Eq. (27), we obtain the source energy PDF as

[TABLE]

Eq. (28) can be cast into the form

[TABLE]

Assuming the source modes have no degeneracy, which is certainly true of the GSM source modes, all $\lambda^{\prime}s$ are distinct. It follows that the integral on the r.h.s of Eq. (29) is straightforward to do in the complex plane—the result is a sum of residues at simple poles, $s_{m}=-i/\lambda_{m}$, $0\leq m\leq N$. The result then reads

[TABLE]

Note the source energy PDF is a weighted superposition of exponential distributions of the modes carrying the average energies $\{\lambda_{n}\}$.

Next, let us write down the pulse peak power as

[TABLE]

where we made use of Eqs. (4) through (9). The peak power PDF of the pulse ensemble is defined as

[TABLE]

where the angle brackets denote averaging over an ensemble of random incident pulses. The averaging in Eq. (32) can be carried out using the delta function property,

[TABLE]

Here $w_{0\pm}$ are the two roots of the equation

[TABLE]

which are written explicitly as

[TABLE]

However, only the positive root is physical because $w_{0}\geq 0$. We can then drop the other root, perform a trivial integration with the delta function and, collecting all terms, we obtain a properly normalized PDF as

[TABLE]

where the set of normalization constants $\{b_{m}\}$ reads

[TABLE]

Here $N$ coherent modes can correspond to physical modes of a multimode fluctuating source. The normalization constants are found from the condition $\int_{0}^{\infty}dP_{\ast}\mathcal{P}(P_{\ast},z)=1$. Eq. (38) happens to represent a weighted superposition of Weibull distributions with the index $3/2$ Wei . Finally, we observe that the PDF can be expressed in a manifestly self-similar form as

[TABLE]

in terms of the similarity variable $P_{\ast}e^{-2gz/3}$. We display the PDF in Fig 2 (left panel) at three propagation distances: $Z=4.5$, $Z=5$, $Z=5.5$ meters in black solid, red dashed and blue dash-dotted curves, respectively. To visualize the self-similarity, we plot the PDF in the scaled variables, explicitly demonstrating in Fig. 2 (right panel) that all three curves coalesce into one. We stress that the derived self-similar structure of the PDF is independent of either a specific source model or the source coherence level, although the source coherence parameter $\sigma_{c}$ affects the overall PDF shape as we will show below. Thus, our fundings apply to any statistical pulse ensemble with thermal-like power fluctuations in the self-similar evolution regime.

To verify our analytical results and ascertain the RW existence among the extreme events in the self-similar intermediate propagation regime, we performed numerical simulations of the modified NLSE, Eq. (3). The results are shown in Fig. 3 for the same propagation distances as in Fig.2. As is evidenced by the figure, numerical PDFs, shown by green histograms, are in excellent agreement with the theoretical curves, confirming that the system is indeed in the self-similar regime at these propagation distances. The RW emergence is marked by the abnormality index (AI) greater than two Dud-rev . The AI is defined as the ratio of an RW intensity $I_{RW}$ to that of a “significant intensity” $I_{1/3}$,

[TABLE]

where the “significant intensity” is defined as the mean intensity of one third of the highest peak intensity events. We can see that at all selected distances $AI\geq 2$, indicating RW emergence in the system. We stress that we limit ourselves to statistical signatures of RW emergence; the examination of the other characteristic of RWs as the waves appearing from nowhere and disappearing without a trace Taki ; Akh3 lies outside the scope of this work. In the inset panels to Fig. 3 we exhibit the corresponding RWs in red curves which clearly acquire a parabolic shape. Moreover, we display the average intensity distribution of the pulse ensemble at the corresponding propagation distance in the same inset panel in a black dashed curve to show that it has a Gaussian-like rather than the parabolic shape. This is because each ensemble realization attains a parabolic shape with a different pulse width and peak pulse power, given by Eqs. (8) and (31), respectively, such that the ensemble average over these parabolic pulse profiles yields a Gaussian-like pulse. Thus, we conclude that extreme events in general—and RWs, in particular—emerge as giant parabolic shape fluctuations away from the average in the self-similar intermediate asymptotic regime. Further, we note that the AI should remain the same for a given pulse ensemble in this regime because $I_{RW}$ and $I_{1/3}$ scale the same way with the propagation distance there. This observation is borne out by our numerical simulations displayed in Fig. 3: We can infer from the figure that $AI\simeq 2.4$ up to numerical round-off errors.

Finally, we examine the ensemble peak power PDF dependence on the source coherence state. To this end, we performed simulations for a very coherent, $\sigma_{c}=10$, and rather less coherent $\sigma_{c}=0.5$ pulse ensembles and compared their PDFs in the self-similar asymptotic propagation regime in Fig. 4 using Eqs. (36) and (37) as well as our numerical data. We can infer from Fig. 4 by comparing either analytical PDF curves or the histograms on the left and right panels that as the source coherence increases—and so does $\sigma_{c}$—the PDF tail stretches as well. Indeed, $AI\simeq 3$ for the very coherent source, while $AI\simeq 2$ for the less coherent one in Fig. 4. This conclusion appears to be at odds with our previous results on RW generation in stimulated Raman scattering with a noisy Stokes pulse input ensemble Yash2 . To reconcile the two observations, we recall that as $\sigma_{c}$ increases, so does the effective number of excited coherent modes of the source MW . Further, the Raman nonlinearity has very long memory implying that the greater the number of coherent modes, the greater the chances for a (giant power) “champion” mode to emerge within a Stokes pulse ensemble. This is because Raman medium memory reinforces unequal energy redistribution from a pump pulse among the Stokes pulse coherent modes throughout multiple Raman scattering cycles. In the case of weakly nonlinear amplifying media, however, the instantaneous Kerr nonlinearity has no coherent memory. Consequently, the lack of cumulative reinforcement of unequal power gain among the pulse ensemble modes favours the chances of giant power mode emergence for sources with a few excited coherent modes.

V Conclusions

In conclusion, we have elucidated the emergence scenario and salient statistical properties of extreme events in weakly nonlinear, statistical open systems exhibiting self-similar intermediate asymptotic evolution. We have demonstrated that rogue waves manifest themselves as giant self-similar fluctuations away from the average and they acquire self-similar, non-Gaussian statistics in such systems. We stress that our generic results hold in the intermediate evolution regime where gain saturation and ensuing amplified spontaneous emission noise are negligible. As the pulse intensity will have sufficiently grown up, we can no longer neglect photon emission by the excited levels of the medium atoms, stimulated by the pulse field; the stimulated emission, in turn, causes gain saturation. At the same time, spontaneously emitted photons stimulate the emission of more random photons, adding amplified spontaneous emission noise to the system. The latter, which can be the dominant noise contribution at high pulse amplification levels, is expected to ultimately lead to the break-down of the system self-similarity, cause the eventual destruction of statistically self-similar RWs and hence invalidate the proposed generic RW excitation mechanism beyond the intermediate evolution regime. The RW nature and statistical properties near gain saturation in presence of pronounced amplified spontaneous emission is a challenging open problem which we plan to address in the future.

S.A.P. acknowledges financial support from Natural Science and Engineering Research Council of Canada, (RGPIN-2018-05497); F.W. acknowledges financial support from National Natural Science Foundation of China, (11874046); Y.C. acknowledges financial support from National Natural Science Foundation of China, (91750201, 11525418).