Fokker–Planck dynamics of the El Niño-Southern Oscillation

The asymmetric nature of the El Niño-Southern Oscillation (ENSO) is explored by using a probabilistic model (PROM) for ENSO. Based on a Fokker–Planck Equation (FPE), PROM describes the dynamics of a nonlinear stochastic ENSO recharge oscillator model for eastern equatorial Pacific temperature anomalies and equatorial Pacific basin-averaged thermocline depth changes. Eigen analyses of PROM provide new insights into the stationary and oscillatory solutions of the stochastic dynamical system. The first probabilistic eigenmode represents a stationary mode, which exhibits the asymmetric features of ENSO, in case deterministic nonlinearities or multiplicative noises are included. The second mode is linked to the oscillatory nature of ENSO and represents a cyclic asymmetric probability distribution, which emerges from the key dynamical processes. Other eigenmodes are associated with the temporal evolution of higher order statistical moments of the ENSO system. The model solutions demonstrate that the deterministic nonlinearity plays a stronger role in establishing the observed asymmetry of ENSO as compared to the multiplicative stochastic part.

. Observed SST-thermocline depth phase diagram for equatorial Pacific. Joint probability density distributions (shading) of monthly-mean Niño-3 SST anomalies (5° S-5° N, 150°-90° W) and equatorial Pacific thermocline depth anomalies, and phase diagram (gray line) depicted by time sequential tracks of these two variables obtained for the period 1958-2016. Thermocline depth anomalies indicated the 17 °C isotherm depth anomalies averaged over the equatorial band between 5° S-5° N, 120° E-80° W. Both variables are smoothed using a 1-2-1 filter. Contour units for probability density are in %.

Results
Observed probability density distribution of ENSO. We estimate the probability density distribution of the observed ENSO system by calculating the 2-dimensional normalized histogram of the Niño-3 (5° S-5° N, 150°-90° W) SST anomaly index ( T E ) and the equatorial Pacific thermocline depth anomalies ( [h] ) (Fig. 1). Here the 17 °C isotherm depth is used as a proxy of thermocline depth instead of 20 °C isotherm in order to avoid an outcropping of the 20 °C isotherm to the surface. The phase-space trajectories of ( T E (t), [h](t) ) (see also Fig. 2 of Timmermann et al. 5 ) revealed the apparent amplitude-and phase-asymmetries between El Niño and La Niña, respectively. The probability density distribution is not symmetric but skewed toward to positive T E and negative [h] . Similarly, the trajectories of ( T E (t), [h](t) ) tend to exhibit a far-reaching shape toward positive T E and negative [h] , indicating a stronger warm phase and a greater oceanic heat content discharge (i.e., shallow thermocline depth) compared to its counterpart 10 . This asymmetric feature provides an easier transition from warm to cold events, compared to the opposite case 27 , as observed. This is because the larger heat discharge associated with the large positive T E produces a strong negative SST tendency through vertical thermal advection; thus, this initiates a cold event more effectively 4 . This occurrence is therefore also directly related to the transition asymmetry 28,29 . Furthermore, the length between two adjacent dots during the warm phase ( T E > 0), that indicates the rate of change ( T E , [h] ), was overall longer than that during the cold phase ( T E < 0), which represents the duration asymmetry 12 .
It has been proposed that the positively skewed T E is caused by a combination of oceanic nonlinear thermal advection due to NDH (i.e., oceanic eddy thermal flux) 15,30,31 , an asymmetric atmospheric response to the ENSOinduced SSTA 25,28,32,33 , and a state-dependent noise-induced instability 19,20,24 . The asymmetry in the duration/ transition of El Niño and La Niña has been linked to the asymmetric impact of the Indian Ocean heat capacitor effect, where warming of the Indian Ocean accompanied by El Niño reinforces a quick termination (although this rarely occurs during a La Niña event) 12,34,35 . Other hypotheses put forward to explain the duration-asymmetry include (i) nonlinear air-sea coupling 32 , (ii) asymmetry in the delayed negative feedback with respect to the reflective oceanic Kelvin wave 29,36 , or (iii) a strong meridional seasonal marching of the western Pacific surface winds (which is a quick terminator of El Niño but not of La Niña) 37 . nonlinear recharge oscillator model. Following previous studies on the role of NDH 15,30,31 in setting amplitude asymmetry 15 , we implemented NDH in the ROM, especially the dominant zonal and vertical NDHs 15,31 . Using the formulation of Jin et al. 6 , the zonal and vertical components of NDH can be represented by combinations of T E and [h] , as follows: where µ a is a sensitivity factor of the zonal surface wind stress against SSTA; β us , β uh ,β w , and β h are the sensitivity factors of zonal current against zonal wind stress and thermocline depth, anomalous upwelling against zonal wind stress, and the equatorial zonal slope of thermocline depth anomalies against zonal wind stress anomaly, respectively; α h is the sensitivity factor of the anomalous subsurface temperature anomaly against thermocline depth change; L and H s are the ocean basin length and ocean surface layer depth, respectively.; and δ is a scaling factor for the zonal gradient of mean SST within the eastern box. Derivations of Eq. (2) are presented in the "Methods" section and also found in Kim et al. 38 . State-dependent (multiplicative) noise can serve as a stochastic nonlinearity, as can be easily demonstrated by introducing logarithmic variables. Here we implement multiplicative noise, in which the stochastic noise amplitude is enhanced only when T E > 0 , which is qualitatively consistent with observations 20 . Finally, the nonlinear ROM can be represented as where ξ 1 , and ξ 2 are independent Gaussian white noise with zero mean and unit variance, and σ 1 and σ 2 are the amplitudes of noise. The noise has a random part and a state-dependent part, and H(x) is a Heaviside step function, which has a value of one if the argument positive and zero otherwise. From the observed time series of T E and [h] , we estimated ROM coefficients ( I BJ , F , ε , α, β 1 , β 2 , B ) via a least-square fitting of T E and [h] tendencies against T E and [h] 39 . Although NDH provides a physical meaning on the quadrature terms in (3a), we do not advocate it as the only nonlinear process in ENSO system. Moreover, as ROM coefficients were statistically computed from the observations, other deterministic nonlinear processes were empirically fitted on quadrature terms of either β 1 or β 2 . The amplitudes of noise ( σ 1 , σ 2 ) and B were estimated as the standard deviation of the residual from ROM coefficient fitting. Thus, all coefficients (see "Methods" for their values) were empirically computed and were not retrieved from a separated physical quantity.
probabilistic nonlinear model (pRoM). All possible solutions of ROM (Eq. 3) can be expressed as probability density functions (PDFs). To describe the evolution of the probability distribution the 2-dimensional Langevin equation can be translated into a Fokker-Planck recharge oscillator model (PROM). Unlike the ensemble method in previous studies 40,41 , the PROM directly solves for the PDF of the stochastic dynamical system, and it has thus been used to explore noise-perturbed ROM system in recent years 19,21,23,25 (see "Methods" for the details). First, the two-dimensional Langevin equation in (T, h) space was derived from the nonlinear ROM in Eq. (3) as follows, is a drift term and the noise amplitude D i (i = 1, 2) takes on the form of a diffusion term in the probability. To simplify the notation we replaced parameters used in Eq. (3), ( I BJ , F , − ε , − α ) by ( a 11 , a 12 , a 21 , a 22 ) and T E and [h] by T and h , respectively. We adopt the Stratonovich formulation for the stochastic system. Assuming that stochastic processes can be represented by Gaussian white noise, the FPE is a considerably effective tool to describe the temporal evolution of PDF of the original Langevin system. Therefore, Eq. (4) is converted to a 2-dimensional FPE with state-dependent diffusion (i.e., PROM) for the probability density P , as follows, ; and k SN = BH(T)T . The subscript DN and SN denote the deterministic and stochastic nonlinearity (multiplicative noise contributions), respectively, and δ(T) is the Dirac delta function. In the second line of Eq. (5), the first term describes the damping (growth) rate of P, and the second and third terms represent the advection toward the T and h directions with velocities (or nonlinear force) f 1 (T, h) and f 2 (T, h) , respectively. The fourth and fifth terms capture the diffusion of the PDF induced by noise. As seen in Eq. (5), the damping rate and advection velocity become asymmetric in (T, h) space due to the nonlinear parameters, β 1 and β 2 , and state-dependent noise B. (5)] and the nonlinear version were computed using two different initial conditions: an El Niño onset case (T = 0K, h = 10m) and a La Niña onset case (T = 0K, h = −10m) . The corresponding initial PDFs were set as a 2-dimensional Dirac delta function centered on the El Niño and La Niña onset cases, respectively. The time integration of Eq. (5) was then performed (see "Methods" for integration method).

time-dependent solutions of pRoM. The time evolutions of PROM for the linear version
As shown in Fig. 2a, the solution of the Fokker-Planck model of the linear ROM produces a damped oscillatory feature initially, and then it converges to a stationary Gaussian PDF, regardless of the initial conditions. The oscillatory features perturbed by the El Niño and La Niña onset conditions are symmetric. Therefore, the time evolution of the standard deviation and the kurtosis for the El Niño and La Niña onset cases are also identical. The standard deviation increases in an oscillatory regime and then stabilizes after reaching its maximum in the stationary regime, and kurtosis is strongly damped in the oscillatory regime, but it then retains its value in the stationary regime. No skewness is observed.
The nonlinear PROM produces prominent asymmetric features during the oscillatory phase (Fig. 2b). The time evolutions of the standard deviation, skewness, and kurtosis for El Niño and those for La Niña onset cases are all different, although they merge in the stationary regime without difference. Interestingly, the standard deviation in the stationary regime obtained from the linear PROM and those from the nonlinear PROM are almost identical, while skewness and kurtosis are different.
A time integration of PROM showed already some basic features. However, the solutions depended on the initial conditions 21 . Therefore, a more general framework is required to better understand the fundamental features of linear and nonlinear PROMs.

eigensolutions of pRoM.
To obtain a more generalized solution, we apply a variable separation method, where the P(T, h, t) was written as a product of a time dependent function ( e n t ), and a phase-space-dependent function ( ψ n (T, h) ), as follows: It should be noted that ψ shall not be directly considered as a probability density, but more as a structure function in phase-space that captures major characteristics of the ENSO model.
On substituting Eq. (6) into Eq. (5), we obtain the eigenvalue problem: w h e r e L FP = − a 11 + a 22 Therefore, the solutions of PROM can be represented by the eigenfunction ψ n (T, h) and the eigenvalue n as a discrete eigenvalue spectrum 42 . It should be noted that L FP is a non-Hermitian matrix, and thus the eigenfunctions are not orthogonal. If a drift potential, U , exists such that www.nature.com/scientificreports/ L FP can be transformed as a Hermitian (self-adjoint system). In this case the FPE becomes a Schrödinger-type equation 21 . However, our 2-dimensional nonlinear system has no explicit drift potential; therefore, the nonlinear PROM cannot be transformed into a self-adjoint system. The eigenmode solution was obtained by discretizing the Fokker-Planck operator ( L FP ) into the matrix form (see "Methods" for the detail). The solution of n had two cases: only real numbers or real and imaginary numbers; here the real and imaginary numbers indicate the growth rate and frequency, respectively. The growth rate and frequency are converted to an e-folding time, 1/Re( n ) and period, 2π/Im( n ) , respectively and plotted in Fig. 3a.
The growth rates obtained from both the linear and nonlinear PROMs were all negative (see Fig. 3a), with the exception of one obtained from the nonlinear PROM, which showed a positive growth rate. This positive growth mode is physically unrealistic as the corresponding fixed point (T = −3.56, h = 127.76) was beyond a realistic range. Although the nonlinear PROM may have a specific domain that shows realistic oscillatory behavior analogous to the maximum potential intensity of ENSO 17 , the positive growth rate mode was beyond of this domain.
Based on their physical characteristics, the eigensolutions were separated in three modes: "stationary", "cyclic", and strongly damped "sub" modes (see Fig. 3a). Among these, the stationary mode was the least damped mode and the eigenvalue had only real numbers ( = −4.43 × 10 −4 for the linear PROM and = −3.36 × 10 −2 for the nonlinear PROM in units of year −1 ), which indicates an absence of oscillatory features. For example, the stationary mode corresponds to the empirical probability distribution obtained from the histogram of ENSO index. The corresponding probability density distribution (PDD) from the linear PROM (Fig. 3b) shows a symmetric Gaussian distribution, of which the maximum is located at the zero point ( T = 0, h = 0 ) (first panel of Fig. 3b), where the zero point corresponds to the mean state of this system or a fixed point ∂P ∂t = 0 , and the PDD resembles a stable focus in dynamical systems' theory. However, the PDD from the nonlinear PROM shows an asymmetric feature resembling the observed PDD, where T is positively skewed and h is negatively skewed (second panel of Fig. 3b). As the highest density is found at the zeros in the linear PROM and at near zeros in the nonlinear PROM, and the PDD diffuses outward, we could argue that the damping and diffusion terms of  www.nature.com/scientificreports/ in Fig. 1. The skewness from the stochastic nonlinearity only PROM is 0.13. However, for a fair comparison, we actually re-estimated the parameters for each PROM version independently; and with these independently computed parameters, the deterministic nonlinear PROM (0.52) still produced larger skewness than the stochastic nonlinear PROM (0.34). Thus, the positively skewed SSTA can be attributed to the nonlinear process, especially by the deterministic nonlinearity. The second least damped modes of both the linear and nonlinear PROMs are a "cyclic mode", of which the eigenvalue is a pair of complex numbers: ( = r ± i i = −0.734 ± 1.726i for linear PROM and = −0.751 ± 1.720i for the nonlinear PROM; units are year −1 ), corresponding to an oscillatory mode. Also, the eigenfunction has a conjugate pair, ψ = ψ r ± iψ i , and these conjugate solution pairs can be combined as follows: where the phase is φ(T, h) = arccos ψ r √ ψ 2 r +ψ 2 i . The combined function exhibits oscillatory behavior with period = 2π i and growth rate = r . This mode could be regarded as a stochastically forced recharge oscillator mode, with a 3.6-year period and − 0.73 (− 0.75) year −1 e-folding damping for both the linear and nonlinear PROMs. It should be noted that the stationary modes are corresponding to the equilibrium state of time-dependent solutions of PROM, while the oscillatory features in the early period of time-dependent solutions are representing the cyclic mode (see Fig. 2).
The cyclic mode is obviously related to the FPE advection terms in Eq. (5). The first and second panels of Fig. 3c, referring to the linear and nonlinear PROM, respectively, show the probability density amplitude ( ψ r (T, h) 2 + ψ i (T, h) 2 ). The probability density amplitude of the linear PROM is uniformly distributed along an ellipse, which closely resembles a limit cycle. However, the probability density amplitude of nonlinear PROM has an oval shape with a small eccentricity, which represents the non-symmetric oscillation of ENSO 10 .
The phase distribution also supports the difference between the oscillatory features of the linear and nonlinear PROMs. The expected duration of the warm event ( T > 0 ) is approximately 6 months shorter than that of the cold event ( T < 0 ) in the cyclic mode of nonlinear PROM (Fig. 3c). This difference is also observed in the stationary mode of the nonlinear PROM; however, it is relatively shorter (~ 2 months) than that in the oscillatory mode (not shown). Equation (8) further confirms that the difference in the temporal behavior can be attributed to the difference in phase, as the phase is determined by the eigen function rather than the eigenvalue.
Differences in growth rate and period between the nonlinear deterministic PROM (third panel of Fig. 3c) and the nonlinear stochastic PROM (fourth panel of Fig. 3c) are marginal, but the skewness are quite distinct. The cyclic mode from the nonlinear stochastic PROM is actually similar to that from the linear PROM, while that from the nonlinear deterministic PROM is almost identical to the fully nonlinear PROM (including deterministic nonlinearities and state-dependent noise). It indicates that the asymmetric oscillatory behavior of this system is mainly driven by the deterministic nonlinearity.
Evolution of statistical moments. To link the time evolution of statistical moment of the PDF and eigenmode, we computed the projection of each eigenmode on each statistical moment. It should be noted that the original statistical moments cannot be fully reconstructed by the eigenmode projection because of the nonorthogonality between eigenmodes. The corresponding mathematical expression is, where n indicates the mode and φ n = tan −1 T k imag n �T k � real n (see "Methods" for detailed derivations). Figure 4 shows the projection of each mode, i.e., [ T k real n ] 2 + [ T k imag n ] 2 . For the linear PROM, all the eigenfunctions is either symmetric or antisymmetric (See Supplementary Fig. 1). Hence, the odd (even) statistical moment can be solely determined by the antisymmetric (symmetric) function ( T k = ∫ T k ψ(T)dT ). The amplitude of T is solely determined by the cyclic mode (n = 2), and thus the temporal evolution feature of T is determined by the cyclic mode (the fluctuation of the PDF can be seen in Fig. 2a), whereas the amplitude of T 2 is mainly determined by the stationary mode (n = 1) and sub-modes (mainly n = 3 ~ 5) and the cyclic mode does not produce it. Both the stationary mode and sub-modes (n = 3 ~ 5) have symmetric eigenfunctions for T (i.e., ψ(T) = ∫ ψ(T, h)dh ), whereas the cyclic mode has an antisymmetric eigenfunction (see Supplementary Fig. 1), which does not contribute to the amplitude of T 2 . In contrast, T 3 is determined by antisymmetric eigenfunctions: the cyclic mode (n = 2) and the sub-modes (n = 6, 7).
In the nonlinear PROM, none of the eigenfunctions are simply either symmetrical or asymmetrical function, due to the nonlinear dynamics (see Supplementary Fig. 1) which leads to the interaction of statistical moments. The broken symmetry causes entrainment of the extra eigenmodes, which is not the case in the linear PROM. Thus, each moment is determined by a complex combination of eigenfunctions. Therefore, not only the cyclic mode (n = 2), but also several sub-modes contribute to the evolutionary feature of T : the cyclic mode also contributes T 2 ; and T 3 is determined by almost all modes, even by the stationary mode.
The power spectrum of each statistical moment further demonstrates the entrainment effect of the sub-modes through nonlinear dynamics (Supplementary Fig. 2). The spectrum of T exhibits a spectral peak that is pronounced in the cyclic mode (n = 2) and smoothly decays toward a higher frequency in both linear and nonlinear models. However, there are some weak anomalies in the higher frequency band in the nonlinear model, that are www.nature.com/scientificreports/ related to sub-modes and higher order nonlinear resonances. These anomalies that feature in the nonlinear model also appear in the standard deviation and skewness spectra, particularly in the lower frequency band, mainly due to the cyclic mode (n = 2) and the sub-modes (n = 3 ~ 4). The spectral density of skewness in the nonlinear model is significantly larger than that in the linear model.

Discussion
In this study, the original linear ROM was expanded to a nonlinear PROM by including deterministic and stochastic nonlinearities. This nonlinear PROM is different from other previously modified ROMs or delayed oscillator models presented in previous studies 19,25,29,32,43 , where empirical parameterizations rather than dynamical parameterizations have been adopted. Furthermore, to provide a more complete phase-space perspective, rather than a local or initial-condition dependent picture, the probability density was directly computed using an eigen analysis; thus it was possible to compare all possible probability distributions of the ENSO system derived from the nonlinear ROM to the observed probability distribution. The nonlinear PROM well mimics the observed nonlinear features of ENSO and shows asymmetries in amplitude, duration, and oscillatory behavior. Moreover, the eigenanalysis of the nonlinear PROM elegantly reveals features about the evolution of the probability distribution and its moments that are not as evident by just analyzing the Langevin equation of the nonlinear stochastic ROM. With the PROM, the relative role of deterministic and stochastic nonlinearities in driving asymmetric features of ENSO can be assessed both qualitatively and quantitatively, in terms of not only basic statistics (variance, skewness) but also temporal behaviors in a phase space. Our results indicate that the deterministic nonlinearity is more important to induce the observed ENSO asymmetry than the stochastic nonlinearity related to the state-dependent noise. In summary, the decomposition of intrinsic probability density distribution features of ENSO physics, including stationary, cyclic, and other sub-modes, provide elegant and new insights into the nature of ENSO's temporal complexity, and the method used in this study could thus become a useful diagnostic tool to further study the complexity of ENSO in complex climate model simulations.
Because of two spatial degrees of freedom, PROM could not address the spatial complexity of ENSO such as the spatial difference in the maximum loading of SST anomaly between El Niño and La Niña 44 , and two types of El Niño (i.e., Central and Eastern Pacific types 45,46 ). Some of the other processes driving the temporal ENSO complexity are also missing, such as the annual cycle-ENSO interaction 5 .    6 , where the first and second terms represent the wind-driven current and the geostrophic current, respectively, and β us and β uh are the corresponding sensitivity parameters. As τ x = µ a T E , which indicates the wind stress response against SSTA, and µ a is a sensitivity parameter, the zonal component of NDH becomes −u ∂T ∂x → µ a β us L/2 δ�T� 2 E + β uh L/2 δ[h]�T� E . The vertical component of NDH −w dT dz over the equatorial eastern Pacific can be written as −�w� E

Methods
, where the subscript sub indicates the oceanic subsurface, and z (= H s ) is the ocean surface layer depth. As upwelling velocity, Coefficient values for ROM. The coefficients used in Eq. (5) are listed in Table 1 where values in (*) indicate coefficients used in the linear ROM, otherwise in the nonlinear ROM.

Langevin equation and FPE.
A Langevin equation is a stochastic differential equation that describes dynamics of a system under random forcing. Originally, it was developed to describe Brownian motion 50 , the random movement of a particle in a fluid. Generalizing the idea, it has been adopted in climate dynamics field. In this adoption the main climate variable of the Langevin equation corresponds to a slowly changing climate state (e.g., SST variability over tropical Pacific) and noise corresponding to weather, a short time scale process 51 .
Another equivalent formulation of the Langevin equation is the FPE, which describes the temporal evolution of the probability distribution of the stochastic system. It is identical to ensemble solution for Langevin equation.
The explicit formulation of the dynamics of the PDF through the FPE enables us to explore fundamentals of the stochastic system, rather than tracking individual random trajectories. = 0 , and are referred to as absorbing boundaries that resemble the observed range of (T, h) (e.g., Fig. 1). The Dirac delta function in the equation is approximated as δ(x) = 1/(2�x), for − �x < x < �x 0, for elsewhere where x is grid size. We confirmed that the total probability, p(T, h)dTdh is almost conserved as 1 over time (not shown here), indicating that the numerical method is suitable. The sensitivity test for initial condition was carried out by varying h from default initial condition ( T = 0K, h = ±10m ), and the overall features of solution ( Fig. 2) are unchanged regardless of choice of initial h (not shown here). Similarly, the eigensolution (Eq. 7) is solved by discretizing Fokker-Planck operator ( L FP ) using the central difference method for the same boundary and numerical conditions as in the time-dependent FPE case. eigenmode formulation of statistical moment. The time evolution of the k-th statistical moment ( T k (t) ) is formulated as the linear sum of eigenmodes (Eq. 9). The statistical moment of the marginal distribution for T is written as T k (t) = T k p(T, h, t)dTdh . Here, the time evolution of PDF ( p(T, h, t) ) can be expanded as p(T, h, t) = ∞ n=1 ψ n (T, h)e n t = ∞ n=1 e r n t cos i n t + isin i n t ψ r n (T, h) + iψ i n (T, h) . Upon substituting the expanded PDF with the statistical moment with respect to only the real part, we get Here, the imaginary part is intentionally neglected, as it would be eliminated by its complex conjugate. Using trigonometric identities, the time evolution of the statistical moment becomes where T k real n = T k ψ r n (T, h)dTdh, T k imag n = T k ψ r n (T, h)dTdh, φ n = tan −1  www.nature.com/scientificreports/