Mathematical modelling of proton migration in Earth mantle

In the study, we address the mathematical problem of proton migration in the Earth's mantle and suggest a prototype for exploring the Earth's interior to map the effects of superionic proton conduction. The problem can be mathematically solved by deriving the self-consistent electromagnetic field potential U(x,t) and then reconstructing the distribution function f(x, v, t). Reducing the Vlasov-Maxwell system of equations to non-linear sh-Gordon hyperbolic and transport equations, the propagation of a non-linear wavefront within the domain, and transport of the boundary conditions in the form of a non-linear wave are examined. By computing a 3D model and through Fourier-analysis, the spatial and electrical characteristics of potential U(x, t) are investigated. The numerical results are compared to the Fourier transformed quantities of the potential (V) obtained through field observations of the electric potential (Kuznetsov method). The non-stationary solutions for the forced oscillation of a two-component system, and therefore, the oscillatory strengths of two types of charged particles can be usefully addressed by the proposed mathematical model. Moreover, the model, along with data analysis of the electric potential observations and probabilistic seismic hazard maps, can be used to develop an advanced seismic risk metric.


Introduction
Recent observations of prominent electrical conductivity due to superionic proton conduction in the host hydrous mineral (iron oxide-hydroxide ) [14,16] in the deep Earth mantle have brought to concrete realization a paradigm and prediction by Vernadski 100 years ago [19,33]. It is widely accepted that water in the mantle is stored in hydrous minerals. Recently, it has been reported that hydrous minerals can bring water into depths greater than 1250 km (lower mantle) [24]. In the deep Earth mantle, hydrogen is stored as hydroxyl (OH − ) in hydrous or other non-hydrous minerals. However, the amount of hydrogen in the Earth's deep interior and its form (molecular or atomic) are unclear. Hou et al. [16] demonstrated that the protons in the FeOOH mineral structure moved freely through a lattice framework of oxygen and other cations. The free movement of hydrogen within water ice crystal has been recognized as superionic proton conduction [22,30]. Furthermore, it has been experimentally confirmed that superionicity may occur in hydrous minerals in the Earth's deep interior.
Here, the H atom constitutes the unit of mass and energy transport, and the migration of electrons induced by hydrogen is more profound. High-electrical conductivity offers a new perspective for using electromagnetic data to explore the Earth's interior [16].
For the last 100 years, dipole electric lines ("antenna crosses") were utilized for measuring the electromagnetic field at the surface, induced by the flow of telluric currents. Most academic research investigates the induction of electro-telluric anomalies [31] where a unit of the domain in the subsurface is modelled as an antenna under the influence of the Earth-ionosphere (or magnetosphere) resonator. The passive telluric antenna has a range of application, from geological prospecting [23] to investigating the non-seismic forecasting of earthquakes [25,32]. It is general consensus that detection of the repetitive electric signals originating from the Earth's interior using electric line "crosses" is a highly complex task, and a global monitoring network is also not available.. This study examines the topic from two research perspectives: addressing the mathematical problem of proton migration in the Earth's mantle and suggesting a prototype for exploring the Earth's interior to map the effects of superionic proton conduction.
The development of a mathematical model for proton migration in the superionic phases [10,14,16] is highly challenging. Through literature survey including experimental results [16], we concluded that mathematical modelling in superionicity research must endeavour to capture the collective motion of charged particles (protons). The primary variable of interest is the self-consistent electromagnetic field potential created by the charges themselves.
The collective motion of charged particles is a manifestation of phase boundary motion in the bulk of the host material. Phase transition cannot occur simultaneously throughout the entire domain. Therefore, the self-consistent electromagnetic field potential changes with time in the same direction in which the phase boundary moves through the material. This significantly affects electron transport as a function of the concentration [13,15,16]. Ensemble of charged particles in time-dependent domains exhibits oscillatory or wave behaviour [9]. Key research and results focus on the properties of wave equations that enhance the understanding of practical motivating problems, such as the propagation of electromagnetic waves with a periodically moving boundary. On the theoretical side, hyperbolic partial differential equations (PDE) are mainly explored for solving the mixed initial-boundary value problem concerned with time evolution in cylindrical regions [20]. The construction of analytical solutions remains challenging for the classical, periodic wave because of the possible instability of the solutions of wave equations. In this case, the main theoretical focus is on spectral methods for studying the driven pattern formation during resonance. In the study, we extend the classical mixed initial-boundary value problem for second-order non-linear hyperbolic equations. The theoretical scope includes the development of a mathematical model of the propagation of a self-consistent electromagnetic field potential U (x, t) in a timedependent domain with moving walls, where self-consistent electromagnetic field potential U (x, t) is created by the charges of two types (H + and OH − ).
We refer to the evidence-based description of physical phenomena and Vlasov's theory [34] for developing the mathematical model of the self-consistent electromagnetic field potential U (x, t) induced in a truncated cone domain in the Earth's mantle through superionic proton conduction. The source of the collective motion of charged particles is localized at the core-mantle boundary (see figure 1). The mathematical description addresses the following question: What happens when the potential is treated as a particle distribution function of two types of charges in a phase space where both the position and velocity approach unstable stationary states? This amounts to the study of a solution for Vlasov-Maxwell (VM) equations concerned with the de-stabilization of a homogeneous stationary state and the non-linear development of the wave governed by res-onances, in general. The rationale is that these distribution functions should be i)relevant to two types of charged particles in a many-particle system, ii) governed by collision-less physics where the dominant physical interactions are caused by energy-induced Coulomb force for electrostatic interaction in an ensemble of charged particles, and iii) could provide basic building blocks to describe the quantitative properties of potential U (x, t) and its complex dynamics.
In the next section, we formulate the problem statement based on experimental results and a survey of the data in literature. Based on the problem statement, we treat the collective motion of charged particles on a mesoscopic level with distribution functions in the form of the VM system of equations. Time-dependence of the spatial domain implies that the problem is distributed in space and is therefore, infinite-dimensional. We supplement the VM system of equations with initial data and boundary conditions.
For a special ansatz of the distribution functions involving the self-consistent electromagnetic field potential, we reduce the VM system of equations to nonlinear sh-Gordon hyperbolic and transport equations. We solve this overdetermined system using the method of characteristics. From these results, we obtain the solution in the form of a travelling wave. Further, we solve the initial-boundary value problem for a non-linear hyperbolic equation using the lower-upper solutions method under the condition that the self-consistent electromagnetic field potential U (x, t) is in the form of a travelling wave.
Using the sh-Gordon equation derived from the 3D wave equation, we numerically simulate a two-component system. We verify through computer simulations, using a complex 3D model, that the self-consistent electromagnetic field potential propagates within the spatial domain as a non-linear traveling wave, and that the boundary conditions move in time as a non-linear wave. Depending on the initial data selection, the non-linear wave takes the form of either nested hierarchies on a wavefront or waves mixing with resonant frequencies. The obtained numerical results are then compared with field observations (Kuznetsov method) [18] through Fourier transform. The Kuznetsov method, which involves field observation of the electric potential, has been investigated by the global station network (Cosmetecor) in Kamchatka, Italy, Altai, Fiji, and Sakhalin [4][5][6] over last 30 years. Non-stationary changes of electric potential and their relation to the major seismic events that accompany subduction of the Pacific plate under Kamchatka are taken as an illustrative example. Incorporating the proposed PDE model and statistical models, future exploration of realistic phenomena in Earth science becomes feasible. With the proposed mathematical (PDE) model, experimental results on the phase transition (charge-transfer state) in dense hydrogen obtained 30 years ago are discussed [13,15].
In Section 3, the VM system is reduced to non-linear hyperbolic and transport equations. In Section 4, the boundary value problem is investigated. Section 5 reports and discusses the non-linear numerical simulations, and the concluding remarks are finally presented. 4 Figure 2: Illustration of the problem statement where domain Ω is limited to the space between L 1 and L 4 (upper-mantle). The source of collective motion of the charged particles is localized at z=0 (core-mantle boundary). The double-dotted arrows denote vector d whose direction coincides with the Z-axis of symmetry of the ansatz distribution function of the velocity and is positive vertically upwards. The directions of the time-dependent electric field are also denoted by double-dotted arrows in the plot. The solid-arrows represent the time-dependent self-consistent magnetic field B(x, t). The single dotted-line and the waves represent the physical processes occurring in the Earth's asthenosphere and crust. The rectangle denotes the location of the multi-electrode station in the Earth's surface boundary.

Problem statement
Let Ω be an anisotropic and multi-layered domain in the Earth's lower-and upper-mantle (see figure 2), where Ω = L 1 ∪ L 2 ∪ L 3 ∪ L 4 . Under certain temperature and pressure, it comprises a fixed ratio of atoms in a defined spatial homogeneous arrangement. During a short period, the arrangement can be approximated to thread-like and plate-like crystalline structures. The arrangement migrates in a stepped manner from L 1 to L 4 sequentially. The quantitative determination of the full distribution of the "atomic-level anisotropy" in domain Ω is beyond the scope of this study.
Based on experimental results, we consider the positive ions as proton (H + , atomic) condensates, analogous to electron condensates, and the negative ions as OH − (molecular) condensates. We introduce a distribution function f i that describes the behaviour of the proton and molecular condensates interacting through the Coulomb force, ignoring collisions. In the Vlasov theory [34], the distribution functions evolve in time as an entity and are not constructed from particle statistics. The main advantage is that the distribution function becomes noiseless. Therefore, distribution function f i (x, v, t) characterizes domain Ω.
Consider the self-consistent electromagnetic field potential U (x, t) induced in domain Ω through superionic proton conduction. The source of collective motion of the charged particles is localized at the core-mantle boundary (z=0, see figure 2). The superionic phase is not constant in time and is impulsive by nature. Therefore, self-consistent electromagnetic field potential is changed in time propagating in time-dependent domain with moving walls, where selfconsistent electromagnetic field potential is created by the charges of two types (H + and OH − ).
The self-consistent electromagnetic field, E(x, t) and B(x, t), is time dependent on the spatial boundary; the self-consistent electric field is aligned with the unit outward normal vector, the self-consistent magnetic field is orthogonal to the unit outward normal vector, and the distribution functions satisfy the reflection conditions with respect to the particle velocity.
We will derive the self-consistent electromagnetic field potential U (x, t) in domain Ω with the boundary condition on the boundary of domain ∂Ω and then we will reconstruct the distribution functions f i (x, v, t).
We treat distribution function f i (x, v, t) only for one selected truncated cone in Figure 2 because all the truncated cones are symmetrical about the centre, which is the Earth's core (see Figure 1). We consider a sub-vertical cylindrical truncated cone to extend from z=0 (core-mantle boundary) to z=L, i.e., approximately 70 km. Figure 2 shows the physical processes manifesting in the uppermost part of the cone (in the Earth's asthenosphere and the crust) as a single dotted line and waves, addressable by observing the electric potential. The multi-electrode Cosmetecor station (Kuznetsov method) records the electric potential in the sub-surface of the Earth's surface boundary.

Vlasov-Maxwell system: Reduction to non-linear hyperbolic and transport equations
We consider the classical solutions of the VM equations in Ω × R 3 × (0, T ) with boundary conditions on the electromagnetic field on ∂Ω × (0, T ), where Ω = G × R is an infinite cylinder and G ⊂ R 2 is a bounded domain with boundary ∂G ∈ C 1 . Here, are the position and velocity, respectively, of a particle; E = E(x, t) and B = B(x, t) are self-consistent electric and magnetic fields, respectively; q 1 > 0 and m 1 are the charge and mass, respectively, of a positive proton; q 2 < 0 and m 2 are the charge and mass, respectively, of a negative ion.
The charge and current densities are defined as follows: We search for the distribution functions of protons and negative ions in the form and the corresponding fields [0, ∞[×Ω → R is an electromagnetic field potential that has to be determined. Constant vector d i defines the symmetry of the problem in the domain. We reduce the VM system for the ansatz of distribution function (8) to a hyperbolic equation and transport equation Based on (9) and (10), we formulate the following theorem: Theorem 1.Let f i (s) be arbitrary differential functions; moreover, Then, every solution U (x, t) of hyperbolic equation (9) with condition (10) corresponds to a solution of the system (1)- (5) in the form where ϕ 0 (x) is an arbitrary function satisfying ϕ 0 = 0, ∇ϕ 0 ⊥ d.
The proof of the theorem is given in [27].
Corollary 1 In a stationary case, (9) is transformed to the form Remark 1 If In this case, (9) can be expressed as follows: For case N = 2 (2D case), this equation is transformed to When l = −1, (15) becomes an sh-Gordon equation: Remark 2 In accordance with the conditions of Theorem 1, the scalar Φ and vector A potentials are defined as follows: where and p(x) is an arbitrary harmonic function. As function U (x, t) satisfies (9), potentials Φ and A are connected by Lorentz calibration For analysing (16), we direct a constant vector d ∈ R 3 along the Z axis; i.e., we Solution (19) describes the wave spreading velocity in the positive direction along the Z axis at constant velocity d/2α (direction of symmetry of domain (16) to where Moreover, introducing a new variable η = (4α 2 c 2 /(4α 2 c 2 − d 2 )) 1/2 ξ, we transform (20)

Boundary value problem
Here, we consider the classical solutions (f 1 , . . . , f n , E, B) of the VM system in a special form ((11), (12)), which we express as follows: where functionsf i : R → [0, ∞[ and vector d ∈ R 3 \{0} are given, and function U : [0, ∞[×Ω → R has to be defined. Assuming ∂Ω ∈ C 1 , we introduce the boundary conditions for the electromagnetic field and the specular reflection condition for the distribution function on the boundary where n Ω is a unit vector normal to ∂Ω.
satisfies the boundary value problem (29), and E, B, f i are expressed in terms of U (x, t) in (22)- (24). Adams and Fournier [1] introduced definitions and theorems for Sobolev spaces and other related spaces of the function, which are reflected in our study.
To prove the existence of classical solutions for (1)- (5) and (25)-(26), we apply the lower-upper solutions method developed for non-linear elliptic systems. In contrast to the stationary problem, the non-stationary one is more complicated because we need to add an equation of the first order (10) to non-linear wave equation (9). Hence, the problem is not "strongly" elliptic, and further development of the method of lower-upper solutions is needed.
is a monotonically increasing function for every x ∈ Ω. Then, boundary value problem Proof. Due to the monotonicity of h, it is easy to check whether there exist for all x ∈Ω. Let u 01 = min(u 0 , 0) and u 02 = max(u 0 , 0). Let u k ∈ C 2,α (Ω) be a solution of the linear boundary value problem for k ∈ (1, 2) According to the maximum principle, u 1 ≤ 0 ≤ u 2 inΩ. From the last one, it follows that u 1 is a lower solution and u 2 is an upper solution for (27). Then, from the theorem of existence (see [26], Theorem 7.1), it follows that (27) has a unique solution u ∈ C 2,α (Ω).
Remark 3 Lemma 1 is a well-known statement and does not require additional comments. We only remark that the monotonicity condition for function h(x, ·) for the VP system is first applied by [28].
We introduce the following conditions for functionf : Then, the following claims hold: 1. Assume conditions (f2) and (f3). Then, h f : R → R is continuous and non-negative 3. Assume conditions (f4) and |d| < 1. Then, the following conditions (f2), (f3), h f are continuously differentiable and h f (u) ≤ Ce −u for all u ∈ R are satisfied.
Proof. According to Lemma 2, h f is a continuous and continuously differentiable function. Function U satisfies (9). Therefore, it follows from Theorem 1 that f 1 , . . . , f n is a solution of the Vlasov equation, and E, B is a solution of the Maxwell system. As U vanishes on ∂Ω, from the definition of U and the translation invariance Ω in d, we obtain that U and ∂ t U vanish on [0, ∞[×∂Ω.

Numerical illustration
In this section, we demonstrate the propagation of a non-linear travelling wavefront of the self-consistent electromagnetic field potential U (x, t) within domain Ω. Mathematically, the transport equation describes the non-stationary transport of charged particles in the domain Ω, and the boundary conditions move in time similar to a non-linear wave. We verify this through computer simulation within a complex 3D model, and compute the numerical solutions of the 3D hyperbolic sh-Gordon equation for the self-consistent electromagnetic field potential U (x, t) derived in the previous sections.
For this numerical illustration, we solve ∆u = sinh(u) in Ω with the Dirichlet boundary conditions given by u(x) = x 1 and also u(x) = sin(x 1 + x 2 + x 3 ). Domain Ω has been defined in the previous sections. We run the simulation through the truncated cone z=0 to Z=L. For the numerical computations, we utilize the finite element library FEniCS. For the domain geometry and mesh, we use FreeCAD and Gmsh software. See [2,11,12,21]. Visualization is realized with ParaView, see [29] The Newton method convergence of the FEniCS automatic differentiation is indicated in Table 1, and the solution is visualization in Figures 3 and 4.

Second-order hyperbolic sh-Gordon equation
In this section, we present a numerical example for the solution of where Ω is the previously introduced truncated cone. We divide the boundary of the cone in ∂Ω = ∂ T ∪ ∂ H ∪ ∂ B for the top, hull (body of the cone), and bottom of the cone. For the boundary data, we impose Additionally, we impose homogeneous Neumann data on ∂ T . We also impose the initial potential and the initial velocity We only impose special boundary data u ∂ that satisfies the transport equation along the boundary: where u B is defined at ∂ B and for x ∈ ∂ H ; z(x, t) ∈ ∂ B represents the intersection of the bottom-circle of the cone with a line in the cone that is also contained in the tangent plane of the cone hull at x. We follow [3, step-23] for the solution of the wave equation. We introduce v = du dt and represent the first order (in time) system as with the same boundary conditions described above. We then discretize time interval [0, T ] into T /δt time intervals t 0 = 0, t 1 = δt, . . . . Potential u is denoted as u n := u(t n , ·), and velocity v is denoted likewise. We then express the 0.5-Crank-Nicolson scheme in time for the linear hyperbolic part and set the nonlinear term in explicit form. We obtain, .
We then perform Galerkin projection of these equations into finite element space considering the previously explained boundary conditions.

Numerical illustration (Case I)
Throughout the numerical simulation, we consider potential u ∂ = 1 at the contact Ω 0 of domain Ω. This value is transported along the boundary with a velocity that is positive vertically upwards (along the y−axis), as shown in Figure 5. The solution of the wave equation is depicted in Figure 6. As a result, the travelling wave in domain Ω becomes decoupled from the potential at Ω 0 , allowing the dynamics in the domain Ω to propagate. This approach corresponds to the commonly used experimental conditions. When the travelling wave moves in time and approaches the top (∂ T ) of the truncated cone as shown in Figures6a, 6b, and 6c, the charged particle density along the vertical axis increases, generating a precursory filamentary structure of the enhanced self-consistent electromagnetic field potential U (x, t). Its typical form is shown in Figure 6d domain Ω, the value of the potential strongly increases in the wavefront. Further increase of the potential results in complete reconfiguration of the spatial modulated filaments as observed in Figure 6e. Non-stationary transport of the charged particles continues to occur synchronously at the top (∂ T ) of domain Ω.
The charged particle density increases, and a stable stationary structure with distinct boundaries is formed. At the top (∂ T ) of domain Ω, nested coaxial-like multiple-walled tubes emerge strictly along the vertical symmetry axis. Each ring-like tube is characterized by the value of the self-consistent electromagnetic field potential. The potential decreases in the radial direction, with the least value farthest from the symmetry axis; potential U (x, t) decays in the direction of the electric field. In this regime, no filamentary structures exist in the lower-part of domain Ω. Its typical form is shown in Figure 6e. In a longer time-interval, the periodic nested hierarchies transform into an apparent chaotic spatial scenario. Its typical form is shown in Figure 6f. This regime is characterized by an apparent maximum of the self-consistent electromagnetic field potential U (x, t) along the vertical axis of symmetry in the middle of domain Ω. In this regime, no periodicity exists in domain Ω.

Numerical illustration (Case II)
Throughout the numerical simulation, domain Ω is controlled by the potential at contact space Ω 0 : This value is transported along the boundary with a velocity that is positive vertically upwards (along the y−axis), as shown in Figure 7. The solution of the wave equation is depicted in Figure 8. An illustrative video is presented in the Supplementary material.
It can be shown that the charged particle density and consequently, the value of potential U (x, t) are maximum at the wavefront, when the travelling wave propagates through domain Ω, as depicted in Figures 8a and 8b top (∂ T ), the enhanced self-consistent electromagnetic field potential U (x, t) becomes non-linear. As shown in Figures 8c, d, and e, the non-linear travelling wavefront transforms into a spot on the vertical axis. The circular rings evolve concurrently with radial periodicity to the central spot. Each ring is characterized by the value U (x, t) of the self-consistent electromagnetic field. During this transformation, the charged particle density along the vertical axis retains the precursory filamentary structure of the enhanced self-consistent electromagnetic field potential U (x, t) in the rest of domain Ω. It is to be noted that no regular modulation of the longitudinal and radial spatial structures exists in the small space between the top and upper-middle part of domain Ω. Its typical form is shown in Figures 8c, d, and e. As depicted in Figure 8f, the propagation of the non-linear travelling wave (case II) results in nested hierarchies, and is bent in the direction of the Lorentz force. The final spatial structure becomes fixed in the domain, and forms two completely distinctive configurations. The coupling between the two regions is most intriguing. It can be shown that the formation of the two configurations and the mostlikely formation of nested hierarchies satisfy the phase-matching resonance, i.e. , coupling occurs through the phase and amplitude of the oscillating waves. Through discrete Fourier transform (DFT) of the numerical data of the amplitude of potential U (x, t) , the values of the frequencies can be determined.
To assess the efficiency of the wave mixing process, a set of finite signals are extracted from selected points along the vertical axis of the cone's symmetry as shown in Figure 9a. As the time series of potential U (x, t) is in real space, the DFT [8] results have conjugate symmetry; thus, the imaginary part of the data is the complex conjugate of the real frequency data. Therefore, the DFT between 0 and N/2 represents the positive frequency data for N time-steps. Figures 9b and 9c display the numerical data of the potential of the frequency component and the frequency of the velocity component, respectively. Analytically, the signal domain is first transformed from time-based to frequency-based. The numerical solver generates the time series of potential U(x, t); the data set includes discrete samples in both time and space. DFT analysis of the discrete signal is performed over k set of points: Thereby, the peaks of the spectrum evaluated with a confidence interval of ±2 · σ(u k (f )) are located(shaded region in Figures 9b and c). To increase the accuracy, Fourier transform is performed in the vicinity of the peaks (continuous lines) with the application of a high-frequency filter. A threshold of 50 (U) is applied for the resulting set of potentials of the frequency components.  The amplitudes of the components at frequencies f are shown in Figure 9d.
The frequencies are f 1 = 2 Hz, f 2 = 13 Hz, f 3 = 17 Hz, and f 4 = 23 Hz. Certain parameter values at which the frequencies satisfy the resonance condition need to be investigated further. In this study, we compare the field measurements of the electric potential (Kuznetsov method) [6,18] and the numerical results.
Note that although there are the four frequency components at frequencies f 1 − f 4 , the two peaks at frequencies f 1 and f 2 represent the spectral line series. These are the main and second spectral lines, responsible for exchange between the non-linear wavefront propagating within the domain and the transport of the boundary conditions in the form of a non-linear wave that varies with time. The two peaks at frequencies f 3 and f 4 describe the wave-mixing effects, and are in resonance with any two components at frequencies f 1 and f 2 . Figure 10 compares the numerical results with the Fourier transform of the potential (V) obtained through field observations of the electric potential (Kuznetsov method). This method, which involves field observation of the electric potential, has been investigated by the extended station network (Cosmetecor) in Kamchatka, Italy, Altai, Fiji, and Sakhalin [4,5]. This can be considered as a prototype for exploring the Earth's interior to map the effects of superionic proton conduction [16].
We consider the large mantle earthquake M8.3, which occurred on 24-th May 2013 as an illustrative example because its occurrence depth is 598 km and Cosmetecor station No4-IMFSET is located in Esso (55.927N 158.703E), which is 362 km from the epi-centre of the earthquake. Technical details on the measurement stations are available in [6]. Quantity U (f ) was first computed through fast-Fourier transform and quantity U (t) was computed in a two-day sliding window. These computations were repeated for six channels in total for each 10-day interval between 10-th April and 20-th June 2013. In the resultant spectrum, three peaks could be identified. The three frequency components appeared before earthquake M8.3 and decayed from the moment of its occurrence.
As seen in Figure 11, the computation for a group of four stations located 0-5-10-25 km from reference point 53.05°158.65°can reveal the spatial localization of the electric potential at the Earth's surface boundary. M7.2, the large earthquake on 30-th January 2016, was the single large event that occurred 100km from reference point 53.05°158.65 in a 10-year interval from 01-st July 2006 to 30-th January 2016. The accuracy of the electric potential measurements (Kuznetsov method) can be improved by considering the non-linear effect. As can be seen, the field measurements of the electric potential and the numerical results are in good agreement.

Concluding remarks
We reduced the VM system of equations to non-linear sh-Gordon hyperbolic and transport equations, which is a solution for the system of differential equations describing the dynamics over time of the distribution function consisting of two types of charged particles. In the computed 3D model, we cast a nonlinear wavefront propagating within a domain and transport of the boundary conditions in the form of a non-linear wave varying with time.
Depending on the initial potential u ∂ at the bottom of domain Ω, the abovementioned non-linear wave takes the form of either nested hierarchies on a wavefront describing the process of two-wave interaction in the case where damping occurs rapidly or a wave mixing process with resonant frequencies if the phasematching resonance is satisfied, i.e. , coupling occurs through the phase and amplitude of the oscillating waves. Certain parameter values at which the frequencies satisfy the resonance condition in case II need to be investigated further. In case II, spectrum u(f ) of amplitude u(t) describes the spectral line series, where frequency components f 3 and f 4 are in resonance with any two components at frequencies f 1 and f 2 .
A mathematical model of proton migration in the mantle was developed, and it was concluded that the electric potential measurements (Kuznetsov method) were in good agreement with the numerical simulation results of the 3D-model. Using the sh-Gordon equation derived from the 3D wave equation, numerical simulation of a two-particle system correctly predicted the temporal features (anomalies) from the time series data (electric potential) through FFT. Large mantle (depth=598 km) earthquake M8.3 on 24-th May, 2013 was considered as an example. The numerical results showed that number of frequency components closely matched those of the field measurements.
In case II, the final spatial structure was fixed in the domain, with two completely distinctive configurations; the nested coaxial-like multiple-walled tubes at the top of the domain showed bending in the direction of the Lorentz force. Such non-linear behaviour of the self-consistent electromagnetic field potential was in the agreement with the spatial localization of the aggregated sum of the electric potential obtained at a reference point on the Earth's surface boundary. Large earthquake M7.2 on 30-th January 2016 was the single large event that occurred 100 km from reference point 53.05°158.65 during a 10-year interval from 01 July 2006 to 30 January 2016. The Kuznetsov method of electric potential observation and the existing station network (Cosmetecor) in Kamchatka, Italy, Altai, Fiji, and Sakhalin [4,5] can be considered as a prototype for the exploration of free hydrogen movement in the Earth's interior.
It is worthwhile to mention that free hydrogen (atom) migration and its collective motion can be viewed as a source that can drastically enhance the oscillator strength (f ), which can be regarded as the effective number of electrons active in the transition [17]. The classical theory of plasma oscillations (referred to as Langmuir waves or electrostatic oscillations) dictates that harmonic oscillators (the electrons in hydrogen atoms) are permitted not only to emit electromagnetic waves but also to absorb them. If the Coulomb force that quantifies the electrostatic interaction between the traverse wave and oscillator is Fourier transformed in a harmonic series, the resulting equation is a differential equation of motion of a one-dimensional ideal atomic oscillator: where ω 0 -is the cyclic frequency of free oscillation, F ω is the frequencydependent Coulomb force for electrostatic interaction between the traverse wave and oscillator, ω is the frequency of the traverse wave, and δ ω is the initial phase.
Absorption of the traverse wave will cause the ideal atomic oscillator to oscillate only if the circular frequency of the wave matches the circular frequency ω 0 of the oscillator. According to classical mechanics, the frequency of oscillation of the Coulomb force [F] must coincide with the frequency of oscillation of the oscillator strength, where the latter is the force required to maintain the equilibrium position of the harmonic oscillator. When the circular frequencies of these forces are slightly diverse, broadening of the spectral absorption lines occurs.
In quantum mechanics, the oscillator strength is a dimensionless quantity that links the theories of emission and absorption to the experiment. Formula 42 of the oscillator strength for absorption is applied in the form: where n and m are quantum numbers, h is Planck's constant, g m is the statistical weight of the absorption levels of the atomic oscillators, and αv are the sum of the quantum indices of the occupied levels of the atomic oscillators. The energy difference between the two states results in broadening of the spectral line.
To summarize, as plasma carries significant free energy to excite collective oscillation, such as the well-known electrostatic Langmuir oscillation, the plasma frequency is w 0 ∼ √ n. The classical Langmuir theory offers a general interpretation of a single spectral line for the collective oscillation of charged particles, of which there is only one type in the system under consideration. In classical terms, a single harmonic oscillator model (one-component system) can be developed.
However, the experimental results for phase transition (charge-transfer state) in dense hydrogen [13,15], and the recent experimental results for superionic proton conduction [14,16] combined with the numerical results presented in this study suggest further investigation of a model of a collective of coupled harmonic oscillators (two-component system).
The existence of non-stationary solutions for VM system of equations has been proved [27]. Formula (43) in [27] describes the behaviour of two types of charged particles in accordance with Vlasov's theory [34].
where v is the velocity vector and r is the radius vector of the charges, F is the force acting on a point charge, and α i >0 and d i ∈ R3 are free parameters.
We express the first VM equation [27] in the form (44) of an equation representing forces. The form is close to the definition of force in Newton's second law of mechanics.
where m is the mass of a unit charge, and F C and F L are the Coulomb and Lorentz forces applied to the unit charge. The first VM equation in the form (44) enables the determination of nonstationary solutions for describing the dynamics of the q 1 (H + ) condensate and q 2 (OH − ) condensate. More specifically, the first VM equation in the form 44 enables the determination of non-stationary solutions for the forced oscillations of two-component system.
The non-stationary solutions for the forced oscillation of two-component system, and therefore, the oscillatory strengths of two types of charged particles can be usefully addressed by the proposed PDE model. We intend to explore the experimental results [13][14][15][16] in future studies.
Appropriate treatment of the oscillatory strengths of two-component system is fundamental to correctly predict the critical pressure for superionic phase transition or the frequencies of hydrogen-dominated phonon modes in multicomponent environment. The increased value of the potential can be considered as a control parameter at which resonance can occur in two-component system. This part is a research question for further investigation.
The second important contribution of the mathematical model is that it enables better understanding of the seismic risk probabilities in a seismic-prone region in combination with the two-layer Bayesian stochastic and Neural Network model for the analysis of electric potential measurements at the Earth's surface boundary (Bobrovskiy, paper under review) and seismic probability risk maps. The time series includes 20 years of continuous observations. For the Kamchatka simulation case, this results in the recognition of two cases with high-probability (the occurrence/non-occurrence of a large earthquake), and the probability distribution over various time frames (from months to weeks). A shortcoming of the existing approach is that most of the data has been collected in locations limited to certain cities in Kamchatka, Altai, Italy, Fiji, and Sakhalin. Some locations in the Pacific Ring such as Indonesia are more seismically active than these locations. Therefore, coordinated effort by the scientific community is required for the continuous collection of data related to free hydrogen transport in the Earth's interior by installing stations (Kuznetsov method) in such locations. The second shortcoming is related to the existing earthquake early warning approach. In this approach, earthquake alerts are activated according to a performance-based engineering framework, and there is an engineering-based risk model to determine whether the alerts are triggered. Therefore, the current efforts of our team focus on real-time seismic analysis to include the developed mathematical model, and the two-layer Bayesian Stochastic and neural network model for various earthquake scenarios. The uncertainties that are an integral part of an earthquake warning system as well as the complex interdependence of urban lifelines and daily life need to be considered. Such interdependence necessitates appropriate seismic risk information to make decisions, against a risk metric, on the mobilization of a wide range of capacities to mitigate earthquake consequences.

Data availability
The datasets generated and analysed during the current study are available in the FigShare repository. Bobrovskiy, Vadim (2022): Data Time Dependent Wave. figshare. Dataset.