Symmetry, Integrability and Geometry: Methods and Applications Multi-Well Potentials in Quantum Mechanics and Stochastic Processes ⋆

Using the formalism of extended N=4 supersymmetric quantum mechanics we consider the procedure of the construction of multi-well potentials. We demonstrate the form-invariance of Hamiltonians entering the supermultiplet, using the presented relation for integrals, which contain fundamental solutions. The possibility of partial N=4 supersymmetry breaking is determined. We also obtain exact forms of multi-well potentials, both symmetric and asymmetric, using the Hamiltonian of harmonic oscillator as initial. The modification of the shape of potentials due to variation of parameters is also discussed, as well as application of the obtained results to the study of tunneling processes. We consider the case of exact, as well as partially broken N=4 supersymmetry. The distinctive feature of obtained probability densities and potentials is a parametric freedom, which allows to substantially modify their shape. We obtain the expressions for probability densities under the generalization of the Ornstein-Uhlenbeck process.


Introduction
Hamiltonians for systems with multi-well potentials are in the focus of classical and quantum dynamics. One of the most fascinating applications of multi-well potentials in classical dynamics is the description of particles transitions from one local minimum to another under the influence of different types of noise. This problem is general for many physical systems from Universe to microphysics. Calculations of the transition rates from one local minimum to another were initiated by Kramers [1] and after 70 years of development the problem is far from complete [2]. From the modern perspective another effect of amplification of transitions over the barrier under the action of weak time dependent periodic signal, which is called stochastic resonance [3,4] attracts the special attention.
The quantum mechanical dynamics in multi-well potentials is interesting for its own reasons. We just recall the tunneling effect, which is in focus again in respect to the wide-ranging research of trapping of atoms of alkali metals in superfluid state [5,6]. Moreover, when more than one barrier exists, the resonant amplification of tunneling rate is possible, which is known as resonant tunneling [7,8]. Experimental observation of this phenomenon in superconductive heterostructures forms a basis for construction of resonant tunneling diode. We have to point out the importance of having a mechanism of changing the parameters of the potential (locations of minima, heights of barriers etc.), to determine the conditions of resonant tunneling.
Classical and quantum systems mentioned in the above, are common in methods of their analysis. The dynamics of stochastic systems is described by the Fokker-Planck equation (FP) [9,10]. One of the methods for solving the FP equation is the eigenfunction expansion method. Due to a formal similarity of the FP and the Schrödinger equations, the eigenfunction expansion method, which is analogous to the bound state expansion of the Schrödinger equation, can be applied to the Fokker-Planck equation. For instance, in the case of bistable stochastic system, the knowledge of the full set of wave functions and eigenvalues of the corresponding quantum Hamiltonian completely determines the time evolution of the solutions to the FP equation. Moreover, it allows one to make conclusions about the dynamics of the corresponding metastable system. Analysis of processes in multi-well potentials is complicated due to the fact that the existing models deal usually with piecewise potentials (such as constructed from rectangular or parabolic wells and barriers), which are presumably far from the real potentials. That is why analytic wave functions and spectrum in such potentials are unknown, which implies the only numerical analysis of their properties.
Solving the FP equation exactly in this approach is closely related to the existence of exactly solvable quantum-mechanic problems, the number of which significantly increases during the last years [11,12]. Recall that the first meaningful example of construction of the significantly nonlinear models of diffusion in bistable system was constructed in [13,14]. The formalism of the Darboux transformation, used in these papers, is the one of the basic methods in constructing the isospectral Hamiltonians in supersymmetric quantum mechanics (SUSY QM). Further development in the construction of the exactly solvable stochastic models was achieved in [15], where the Fokker-Planck models with prescribed properties were constructed by use of the Darboux-Grum procedure. Existence of exactly solvable [16] and partially solvable [17] quantum mechanical models with multi-well potentials could fill the gap in theoretical analysis of the above mentioned processes and may be considered as a more realistic approximation in their research. Approach based on the accounting the instanton contribution in double-well-like potentials [18,19] is commonly used in the description of tunneling processes and, in particular, in studying the features of the Bose-Einstein condencates in multi-well traps. The instanton calculus is also well applied to the transition between the wells under the influence of noise, because the dynamics of these processes is mainly determined by the energy of the first excited state.
This paper is aimed at determining the relations between different types of potentials in the framework of the extended supersymmetric quantum mechanics (N = 4 SUSY QM) [20,21,22]. We also study the possibility of changing parameters of potentials in wide range and obtaining the exact expressions for spectrum and wave functions. The latter is especially important when potential is equipped with two and even more local minima, so the "resonant" tunneling may be realized. The obtained expressions are used for derivation of new exactly solvable stochastic models. Here we develop the approach of [23] to the construction of exactly solvable stochastic models for potentials with several local minima. Contrary to the case of previously established connection between the FP equation and N = 2 SUSY QM, the considered approach allows one to extend the range of exactly solvable stochastic models, since the super Hamiltonian of N = 4 SUSY QM contains large number of the isospectral Hamiltonians. A special characteristic of the approach is the existence of a parametric freedom in probability densities and new potentials entering the Langevin equation. It makes, in particular, possible to change the shape of potentials and densities without changing the time dependence of the probability density. This fact makes possible to study the characteristics of stochastic dynamics of particles in multi-well potentials with variable parameters. On the other hand, the development of new exactly solvable diffusional models leads to the construction of more adequate approximations to real processes.
In Section 2 we briefly discuss the procedure of construction of isospectral Hamiltonians with additional states below the ground state of the initial Hamiltonian in N = 4 SUSY QM. We es-pecially emphasize that multi-well potentials can arise due to general properties of the solutions to the Schrödinger type auxiliary equation. In particular, we discuss the phenomenon of partial supersymmetry breaking in N = 4 SUSY QM. In the next section we give general expressions, that allows one to analyze the obtained potentials and wave functions, independently on the concrete form of the initial Hamiltonian. It allows, for instance, to generalize the concept of form-invariant potentials [24] and to calculate the normalization constants of zero-modes wave functions in the case of exact and partially broken supersymmetry. Section 4 is devoted to the construction of stochastic models, by use of connections in the supermultiplet of isospectral Hamiltonians of N = 4 SUSY QM. We obtain the expressions for corresponding probability functions and potentials, entering the Langevin equation. Using the harmonic oscillator Hamiltonian as the initial one, we obtain the expressions of the potentials of isospectral Hamiltonians and wave functions for the case of exact and partially broken supersymmetry. The analysis of conditions for emergence of double and triple-well potentials and ways for varying of their form in wide range is provided in Section 5. The proposed scheme of construction of new stochastic models is demonstrated in Section 6 on the example of generalization of well known Ornstein-Uhlenbeck process. We give the results of the calculations of probability densities and discuss the ways of their modification when parameters of scheme vary. In conclusion we summarize the main results and their possible applications.

Isospectral Hamiltonians of N = 4 SUSY QM with additional states
Extended N = 4 SUSY QM [20,21,22] is equivalent to the second-order polynomial SUSY QM (reducible case) [25,26,27,28] and assumes the existence of complex operators of sypersymmetries Q 1 (Q 1 ) and Q 2 (Q 2 ), through which the Hamiltonians H σ 2 σ 1 can be expressed. Hamiltonian of N = 4 SUSY QM has a form ( = m = 1): where W (x) is a superpotential and σ (i) 3 are matrices, which commute with each other and have eigenvalues ±1, σ 3 = 1 ⊗ σ 3 . Supercharges Q i of extended supersymmetric quantum mechanics form the algebra: Hamiltonian and supercharges act on four-dimensional internal space and Hamiltonian is diagonal on vectors ψ σ 2 3 . Supercharges Q i (Q i ) act as lowering (raising) operators for indices σ 1 , σ 2 . It is convenient to represent the Hamiltonian structure and the connection between wave functions in diagram form (see below).
Obviously, due to commutativity of operators Q i andQ i with Hamiltonian, all ψ σ 2 σ 1 (x, E) are eigenfunctions of Hamiltonian with the same eigenvalue E. The only exception is the case, when the wave functions turn to 0 under the action of generators of supersymmetry.
Construction of isospectral Hamiltonians within N = 4 SUSY QM is based on the fact that four Hamiltonians are combined into the supermultiplet H σ 2 σ 1 . Nevertheless, it has to be ≡ H σ 2 −σ 1 only three of them are nontrivial. Let's note, that such a relation in the case of higher-derivative SUSY [26] establishes the correspondence between a quasi-Hamiltonian and operators of the Schrödinger type, and is identical for any superpotentials. The procedure of construction of isospectral Hamiltonians, when ground state is removed from the initial Hamiltonian, is considered in [22] in detail. We will consider the construction of isospectral Hamiltonians by adding the states above the ground state of initial Hamiltonian. Similar procedure was already performed (e.g. in [23]), but the distinctive feature of the present research is to obtain of general results without specification of concrete form of the initial Hamiltonian. Let's consider the auxiliary equation: Let's take one of the Hamiltonians as the initial one where ε is the so-called factorization energy. Hereafter the energy is measured from ε. Strictly speaking, the supersymmetry relations with supercharges (2.1) expressed through superpotential W (x) are satisfied for shifted by ε Hamiltonian H − ε. However, due to commutativity of supercharges with constant ε, the relations between wave functions of H and H − ε remain the same. When fixing operator H σ 2 σ 1 , the form of W (x) depends on the choice of factorization energy ε, hence the Hamiltonians H σ 2 −σ 1 , H −σ 2 σ 1 , H −σ 2 −σ 1 also have nontrivial dependence on ε. When ε < E 0 (where E 0 is the ground state energy of initial Hamiltonian), the auxiliary equation (2.2) has two linear independent solutions ϕ i (x, ε), i = 1, 2, which are nonnegative and have the following asymptotics [15]: the general solution has the form ϕ(x, ε, c) = N (ϕ 1 (x, ε) + cϕ 2 (x, ε)) with appropriately chosen constants (N is the normalization constant) and has no zeros. Thus, the functionφ(x, ε, c) = N −1 ϕ(x,ε,c) is finite and can be normalized at every concrete choice of ε and c. Let's note that, for some values of ε and c, ϕ(x, ε, c) can have local extrema. In this case the natural choice of the initial Hamiltonian is H + − or H − + (which are identical due to the symmetry of H σ 2 σ 1 under σ 1 ↔ σ 2 ). Then the superpotential has the form: with two new arbitrary parameters λ, x i , but one of them is inessential, because it gives an additional contribution to W (x). All the Hamiltonians forming the supermultiplet have nontrivial dependence on these parameters. To consider the connection between Hamiltonians from the supermultiplet, let's take H − + as the initial one. Denoting the solution to as ψ − + (x, E) and using the first representation of the Hamiltonian H σ 2 σ 1 (the l.h.s. of (2.3)), we obtain the following relation between H − − , ψ − − (x, E) and the initial expressions: The new state with E = 0 and by definition normalized wave function (remind that energies are measured from ε) appears in the Hamiltonian. For the discrete spectrum the normalization of the excited states wave functions is conserved. Using the second representation (the r.h.s. of (2.3)) and initial Hamiltonian and wave functions can be obtained: It is worth mentioning that the normalization of the wave function ψ + + (x, E = 0), as in the case of one-well potentials, can always be performed at anyφ(x, ε, c) and λ by use of the following expression in the normalization condition: .
From this relation it is easy to derive the relation between the normalization constants: . The usage of superpotential (2.4) withφ(x, ε, c) corresponds to exact supersymmetry, which leads to existence of zero-modes in H − − and H + + . Existence of two zero-modes in super Hamiltonian of N = 4 SUSY QM is caused by the fact that the Witten index theorem has to be modified when intertwining conditions are nonlinear, as discussed in detail in [25].
Let's consider the case when expression (2.4) contains one of the particular solutions, e.g. ϕ 1 (x, ε), instead of ϕ(x, ε, c). If one of the particular solutions of second order differential equation is known, the second solution can be obtained from the relation .
Thus, the superpotential (2.4) becomes It is easy to show, that the state with the energy E = 0 in the spectrum of H − − is absent (i.e., the wave function ψ − − (x, E = 0) is nonnormalizable) hence the spontaneously broken N = 2 supersymmetry exists. On the one hand, the zero energy state has the wave function normalizable for certain values of λ > 0, appears in the spectrum of H + + . Here the exact N = 2 supersymmetry takes place. Furthermore, it is known [29,30], that the partial supersymmetry breaking is impossible in N = 4 SUSY QM without central charges. This contradiction resolves with taking into consideration that the employment of a factorization energy ε in the construction of isospectral Hamiltonians is the simplest way of incorporation of central charges in N = 4 SUSY QM. The similar situation occurs in consideration of form-invariant potentials [31]. More complete and consistent consideration of N = 4 SUSY QM with central charges is given in [32], but this point is over the scope of the current paper.

General properties of multi-well potentials
We get started the consideration of properties of isospectral Hamiltanians, derived in previous section without appealing to the form of the initial Hamiltonian, from the following usefull expressions. Let y 1 and y 2 be two linear independent solutions of a homogeneous second order differential equation. Then the following expression holds [33]: Here W {y 1 , y 2 } = y 1 y ′ 2 − y ′ 1 y 2 is the Wronskian, which for the second order differential equation, reduced to canonical form, as the Schrödinger equation, is independent on x and thus could be taken out of the integral. This expression is very useful for calculations of integrals in expressions (2.4)-(2.6). First of all, it is natural to set tends to 0 in the limit and the function with asymptotic Let's fix c = 1 in (2.6). This choice allows one to simplify the consideration, but nevertheless, reveals the fundamental features of Hamiltonians and wave functions of N = 4 SUSY QM, using only general properties of the solutions to the auxiliary equation under ε < E 0 . Using relations (3.1), it is easy to obtain the following modulo constant term expression for the superpotential: Let's note that ∆(±∞, ε, λ), entering Λ(ε, λ), are determined by the asymptotics of solutions to the auxiliary equation. Due to this fact, and since for H − − the potential is determined by the symmetric combination ϕ 1 (x, ε) + ϕ 2 (x, ε), while in the case of H + + -by the asymmetric combination ϕ 1 (x, ε) + Λ(ε, λ)ϕ 2 (x, ε), we get: In some sense, the potentialsŪ − − (x, ε) andŪ + + (x, ε, λ) are form-invariant [15], i.e. potentials and wave functions transform to each other by changing of parameters and their spectra are identical and this holds independently on the choice of the initial Hamiltonian. As it will be shown in the next section, if the potential in H − − is a multi-well symmetrical potential, then in H + + it should be asymmetrical. Moreover, varying ε, as well as λ, the form ofŪ + + (x, ε, λ) can change. Relations (3.1) and (3.2) are also useful to derive the exact form of the wave functions ψ − − (x, E) and ψ + + (x, E). Thus, the expression for the wave functions ψ + + (x, E) comes from the similar expression for ψ − − (x, E) by the substitution ϕ 1 (x, ε) + ϕ 2 (x, ε) → (ϕ 1 (x, ε) + Λ(ε, λ)ϕ 2 (x, ε)). In particular, the normalization constant for ψ + + (x, E = 0) can be obtained from the corresponding expression for ψ − − (x, E = 0). To demonstrate the possibility of partial supersymmetry breaking in N = 4 SUSY QM we have to show that the wave function ψ + + (x, E = 0) ∼ 1 ϕ(x,ε,λ) is normalizable. The value of the normalization constant can be derived from (3.1) under c = λ > 0. We will return in what follows to the calculation of this constant for the concrete form of the initial Hamiltonian. In its turn, the potential in H + + has the same form as the corresponding potentials in the case of exact supersymmetry under replacement Λ(ε, λ) → λ.
As it follows from (2.6) and (3.3), the parametric dependence on λ is only present in H + + and ψ + + (x, E). It should be noted that the value of λ is not to be fixed by the normalization condition for the wave function. Existence of a parametric freedom is common in building the isospectral Hamiltonians, which is based on different versions of the inverse scattering problem [34,35,36]. It is caused by the ambiguity in reconstructing quantum-mechanical potentials from the spectral data. We will show latter that similar situation exists also in stochastic models. This is unexpectedly enough, since the potentials of stochastic models have direct physical meaning in contrast to quantum mechanical potentials.

Probability densities functions
The FP equation is equivalent to the Langevin equation, however its application in physics is more wide, since it is formulated in more appropriate language of the probability densities P ± (x, t; x 0 , t 0 ). According to [9,12] the FP equation has the form: where U ± (x) is the potential entering the Langevin equation. Let's set t 0 = 0. The Fokker-Planck equation describes the stochastic dynamics of particles in potentials U + (x) and U − (x) = −U + (x). Substituting the FP equation transforms into the imaginary time Schrödinger equation: in which the diffusion constant D can be treated as the "Planck constant", while H ± has the form of N = 2 SUSY QM Hamiltonian. Equations (4.1) are basic to apply of the eigenfunctions expansion method for construct exactly solvable stochastic models by use of the results of corresponding quantum mechanical problems. Details of this method in the framework of N = 2 SUSY QM can be found in [9,12]. Utilization of the extended N = 4 SUSY QM formalism gives additional possibilities to construct new exactly-solvable models of stochastic processes, since the Hamiltonian H σ 2 σ 1 of N = 4 SUSY QM includes four isospectral Hamiltonians. It allows to obtain new K σ 2 σ 1 (x, t) and U σ 2 σ 1 (x), hence P σ 2 σ 1 (x, t; x 0 ). To obtain the expressions for probability densities P σ 2 σ 1 (x, t; x 0 ) we should to know not only the wave functions and the spectrum of H σ 2 σ 1 , but those of the corresponding U σ 2 σ 1 (x). It is easy to see that Further consideration is based on account of the symmetry of N = 4 SUSY QM H σ 2 σ 1 under σ 1 ↔ σ 2 , that leads to the relation: First equality indicates the existence of the expression of H (+) in terms of Q 1 (Q 1 ), while the second equality implies that H (−) is expressed in terms of Q 2 (Q 2 ). At the same time, the supercharges entering H (+) and H (−) are substantially different: According to (4.3), quantum H + − and H − + have the same spectrum and the wave functions, thus the corresponding K(x, t) are identical. But the corresponding stochastic models are described by essentially different potentials, therefore their P (x, t; x 0 ) are different. Moreover, the potential U + − (x, ε, Λ) = −U + + (x, ε, Λ) = −D ln |ϕ 1 (x, ε) + Λ(ε, λ) ϕ 2 (x, ε)| has non-trivial parametric dependence on λ, which implies the existence of the family of stochastic models with the same time dependence of the probability density, but substantially different coordinates dependence. Surprisingly, the parametric freedom allows for a modification of the shape of potential.
Recall that physical quantities, e.g. time of passing the potential maximum, considerably depend on the local modification of the potential form [37,38,39].
Denoting ϕ(x, ε, c = 1) = ϕ 1 (x, ε) + ϕ 2 (x, ε), the expression for calculation of distribution function P − − (x, t; x 0 ) takes the form: For simplicity let's assume that the spectrum of the initial Hamiltonian has only the discrete states. When the continuous spectrum exists the summation should be replaced with integration over the corresponding density of states. Substituting the expressions for wave functions ψ − − (x, E n ) (2.5) into (4.5) we obtain: The expression for P + + (x, t; x 0 ) can be obtained from equation (4.6) by substituting ϕ(x, ε, 1) → ϕ(x, ε, Λ(λ, ε)), N −2 → N −2 Λ . Both of P − − (x, t; x 0 ) and P + + (x, t; x 0 ) have the equilibrium distributions at t → ∞, which coincide with square of the zero modes wave functions in Hamiltonians H − − and H + + . Moreover, it is easy to see that terms from the excited states of H − − and H + + do not contribute to the normalization condition of these density functions, i.e. the normalization of these functions is guaranteed by their equilibrium values. Dependence on parameter λ in P + + (x, t; x 0 ) is not eliminated by the normalization condition. The distribution P − + (x, t; x 0 ) (P + − (x, t; x 0 )) describes the stochastic dynamics of a particle in the metastable state of the potential, inverted to U − + (x, ε) (−U + + (x, ε, Λ(ε, λ)). Thus, it has not the equilibrium limit and is determined by the expression: Similarly to the previous case P + − (x, t; x 0 ) can be obtained from P − + (x, t; x 0 ) by substitution ϕ(x, ε, 1) → ϕ(x, ε, Λ(λ, ε)). The parametric dependence on λ can be eliminated by the normalization condition for P + − (x, t; x 0 ), but as it will be demonstrated below, choosing the initial Hamiltonian of harmonic oscillator, the parametric arbitrariness in P + − (x, t; x 0 ) remains. To summarize, we should note that in the framework of N = 4 SUSY QM four models of stochastic dynamics emerge. Two of them correspond to the motion of particle in double-well potentials (symmetric and asymmetric) under the action of the Gaussian noise, the other ones describe the dynamics of the particle in the metastable state, when N = 4 supersymmetry is exact. For the case of the partial supersymmetry breaking P − − (x, t; x 0 ), P − + (x, t; x 0 ) and P + − (x, t; x 0 ) describe stochastic processes in the metastable state, while P + + (x, t; x 0 ) -in bistable. Some of them possess the parametric freedom, that allows to change the parameters of external field.

Isospectral Hamiltonians with almost equidistant spectrum
To construct the explicit expressions for potentials and wave functions we choose the initial Hamiltonian with the harmonic oscillator (HO) potential. Let's consider the solution to the auxiliary equation for ε < E 0 = ω 2 , ( = m = 1): Introducing dimensionless variables ξ = √ 2ωx we obtain the equation for ϕ(ξ,ε), whereε = ε ω : Now we consider two different cases. a) Exact N = 4 SUSY. This equation has two linear independent solutions: parabolic ). According to terminology in the above, we denote [40], with gamma-function Γ(−ν). Following the procedure from the previous section, the general solution to the auxiliary equation is chosen to be: As one can see from (5.1), ϕ(x,ε, 1) is an even function of ξ. To obtain the exact form of the superpotential, the integral, entering the definition of the W (x, ε, λ) as well as the normalization constant N −2 have to be calculated from (3.1), (3.2). Due to the symmetry of ϕ(ξ,ε), expression of the integral simplifies and takes the form: Hence, the superpotential is: .
Using the asymptotic value of the parabolic cylinder function, we have ∆(+∞,ε, 1) = −1 and the expression for the superpotential simplifies: .
Analysis of (5.4) reveals that there exist several local minima only for 0 <ε < 1 2 . With values ofε to be close to the right boundary, the third local minimum appears (Fig. 1) and depth of the outside minima increases. It should be noted that in terms of the dimensionless variable ξ the only way to vary the form of the potential is by varyingε and λ. In the case of natural units, additionally, the form of the potential (in particular, positions of the local minima) can be changed by variation of ω.
b) Partial N = 4 SUSY breaking. The utilization of the particular solution of auxiliary equation for derivation of superpotential (5.2) realizes the situation of partial supersymmetry breaking. As it was mentioned above, the spectrum of H − − does not contain states with E = 0, because their wave function is nonnormalizable. The shape ofŪ − − (ξ,ε) is presented in Fig. 3. At the same time, zero state appears in the spectrum of H + + with the wave function: ε)) . The normalization constant N −2 λ is calculated by use of (5.4) and has the form N −2 λ = 2λ √ π Γ(−ν) . Thus, it is restricted to λ > 0. The Hamiltonian becomes H + ,ε, λ). The potential shape (ε = 0.47) is presented in Fig. 2 with λ > 0.
, with Φ(ξ) to be the error function), potentials U − − (x,ε) correspond to that of [41] (the only difference is in additional constant term) This potential is single-well, as well as correspondingŪ + + (ξ,ε = − 1 2 , λ), because existence of several local minima is possible only for 0 <ε < 1 2 . Using the form-invariance property ofŪ + + : The spectrum of the Hamiltonian with (5.6) contains, in contrast to (5.5), the state with E = 0 and the wave function of the form: ) .

Generalization of the Ornstein-Uhlenbeck process
Let us give an example of the construction of new stochastic models. We choose the harmonic oscillator (HO) Hamiltonian as the initial one. Then, we consider the auxiliary equation with ε < E 0 = Dω 2 (recall, that D plays the role of the "Planck constant" in the FP equation and ε = ε Dω ). Here we consider two cases. a) Exact N = 4 SUSY. For many problems of stochastic dynamics it is reasonable to use symmetric double-well potentials, so we consider first the case c = 1. As it can be seen from (5.1), ϕ(ξ,ε, 1) is even function, and the normalization constant N −2 can be calculated from (5.2). One of the main characteristics of the stochastic process is the potential entering the substantially influence on many characteristics, for example, on the rates of the potential barrier crossing [37,38,39]. U + − (x, Λ) = −U + + (x, Λ) = −D ln |ϕ 1 (x,ε) + (λ + 1)ϕ 2 (x,ε)|, obtained within the framework of proposed model, is presented in Fig. 4. Changes of the U + + (ξ, λ) shape due to varying λ become more and more sharp under λ → −1. The analysis of the expressions of U + + (ξ, λ) shows, that several local minima are possible only when 0 <ε < 1 2 . Variation ofε inside this range changes the height of the barrier, which significantly increases withε → 1 2 . Moreover, the shape of the potential can be changed (especially the minima locations) by varying ω, when switching to natural variables. To summarize, the existence of (ω,ε, λ) significantly change the shape of the U + + (x, λ). Substituting the HO wave functions to (4.7) and using the Mehler formula [40], we get P − + (ξ, z; ξ 0 ): where z = e −ωt . As it was noted in the above P − + (ξ, z; ξ 0 ) does not have the equilibrium value and could be normalized (the proof of this claim is given in Appendix), i.e. +∞ −∞ dξ P − + (ξ, z; ξ 0 ) = 1.
The probability density P − − (ξ, z; ξ 0 ) is characterized by the existence of the equilibrium value, which is determined by the zero energy wave function of H − − , i.e. .
After some calculations (see Appendix) we obtain the expression for P − − (ξ, z; ξ 0 ): The probability density P + − (ξ, z; ξ 0 ) can be obtained from P − + (ξ, z; ξ 0 ) (6.1) by substituting ϕ(x,ε, c = 1) → ϕ(x,ε, λ + 1). Note that the normalization condition does not eliminate the λ-freedom (see Appendix). As it has been noted in previous section, the expression for P + + (ξ, z; ξ 0 ) is obtained from P − − (ξ, z; ξ 0 ) with replacing ϕ(x,ε, c = 1) → ϕ(x,ε, λ + 1), . The dependencies of P + + (ξ, z; ξ 0 ) on spatial and time variables, presented in Fig. 5, demonstrate sharp modifications (with fixedε, ω, D) at different values of λ, especially with λ → −1, when the considerable assymmetry of P + + (ξ, z; ξ 0 ) as well as decreasing the time of passing to bimodality are observed. Moreover, modifications of the shape of potential U + + (ξ, λ), as well as P + + (ξ, z; ξ 0 ), can take place with variation of ω, when we pass from dimensionless variables ξ to physical x. This leads to the shift of the local minima and change the height of the barriers. Whenε → 1 2 , the barrier between local minima significantly increases, that allows to study the transition time from one local minimum to another as a function of ∆U D . b) Partial N = 4 SUSY breaking. Let's briefly discuss partial N = 4 SUSY breaking and its consequences to stochastic dynamics. In this case three models for the description of the stochastic dynamics of particles in metastable states (U − + (ξ,ε), U − − (ξ,ε), U + − (ξ,ε, λ)) and one in the bistable state (U + + (ξ,ε, λ) ) arise in the considered approach. The corresponding potentials have the form: The range of λ variation in U + + (ξ,ε, λ) and U + − (ξ,ε, λ) is restricted to be λ > 0. Fig. 4b presents the dependencies of the potentials on spatial variables at fixedε and λ. Stochastic dynamics in mentioned metastable states satisfies to different boundary conditions: reflective boundary at left (right) and absorbing at right (left) for U − − (ξ,ε) (U − + (ξ,ε)) and absorbing boundaries at both sides for U + − (ξ,ε, λ). The corresponding probability densities can be calculated with taking into account previously obtained expressions with appropriate modifications.

Conclusions
In the paper we consider the construction and the study, within the framework of N = 4 SUSY QM, of the general properties of Hamiltonians with multi-well potentials. These results have been used to obtain exactly solvable models of stochastic dynamics. The special attention was made to studies of features of these Hamiltonians without concretization of the form of initial Hamiltonian. Relations for certain type of integrals, containing the fundamental solutions to the Schrödinger type equations, allows one to show the form-invariance of the isospectral Hamiltonians with multi-well potentials, obtained within of N = 4 SUSY QM. In other words, having the identical spectra of H − − and H + + , the corresponding potentials can be obtained from each other through replacing parameters of the problem. Moreover, this relation allows to analytically calculate the normalization constants of zero modes, using only asymptotic values of fundamental solutions. The construction of isospectral Hamiltonians within N = 4 SUSY QM also implies the partial supersymmetry breaking. Taking the model of harmonic oscillator as an example we derive the exact form of the isospectral Hamiltonians with multi-well potentials and the corresponding wave functions. The existence of a parameteric freedom gives a possibility to vary the shape of the potential in the wide range. This becomes important in the research of different phenomena, such as tunneling processes, which are sensitive to the structure of the multi-well potentials.
The obtained results have been used for the construction of exactly-solvable stochastic models. Obtained potentials, entering the Langevin equations, and probability functions are characterized by the parametric dependence, which allows sufficiently modify their shape. A parametric freedom is typical for quantum mechanical isospectral Hamiltonians, where different versions of the inverse scattering problem are used [34,35,36]. Such a freedom arise in view of an ambiguity reconstructing the potential from the spectral data. On the other hand, the parametric freedom in stochastic models is quite surprising, since their potentials, in contrast to the quantum mechanics, have direct physical meaning. Many characteristics of stochastic dynamics, such as potential peak passage times and the metastable state lifetime [37,38,39], substantially depend on the shape of the potential. Therefore, it is interesting to investigate the characteristics of stochastic processes in multi-well potentials, such as the Kramers problem, the stochastic resonance etc., and their dependence on a modification of the shape of the potential. We note that our results admit a generalization to non-Markovian processes, by use of the approach of [42]. Furthermore, the obtained results can be generalized to the stochastic models of polymer dynamics [43]. We have also considered non-trivial consequences of the partial supersymmetry breaking in N = 4 SUSY QM to the description of stochastic dynamics. In this case three probability density functions with the same time dependence (i.e. with the identical spectrum of the Fokker-Planck operator) describe stochastic dynamics of particles in metastable states, and one function corresponds to the dynamics in the bistable state.