Nonlinear wave-mechanical effects in Korteweg fluid magma transport

TL;DR

This paper explores wave-mechanical effects in magma transport using a nonlinear wave equation, revealing phenomena like symmetry breaking, phase transitions, and solitary wave solutions in volcanic conduits.

Contribution

It introduces a novel nonlinear wave equation with logarithmic nonlinearity to model magma flow, capturing phase transitions and inhomogeneities in volcanic conduits.

Findings

Identification of wave-mechanical effects in magma transport.

Existence of solitary wave solutions with Gaussian profiles.

Demonstration of topological kink solitons in foam phases.

Abstract

Statistical mechanics arguments and Madelung hydrodynamical presentation are applied to the transport of magma in volcanic conduits. An effective wave equation with logarithmic nonlinearity becomes apparent in systems of this kind, which describes the flow of a two-phase barotropic Korteweg fluid with capillarity, and allows multiple eigensolutions thus leading to wave-mechanical effects. We study spontaneous symmetry breaking in the erupting lava which flows up the conduit, so that fluid fragmentation and nucleation of density inhomogeneities occur; therefore, changing temperature can trigger a transition between the "magma-dissolved gas" fluid and magmatic foam phases. This phase structure is studied by both analytical and numerical methods. In the fluid phase, cell-like inhomogeneities occur which are described by solitary wave solutions with a Gaussian density profile; we derive the many-body interaction potential for these inhomogeneities. For the foam phase, we demonstrate existence of topological kink solitons which describe bubble-type inhomogeneities; their stability is ensured by the conservation of a topological charge.

Equations

$$P\propto \exp (-\mathcal{E}/T)\approx \exp \left(-\mathrm{U}/T\right),$$

$$\Psi =\sqrt{\rho }\exp (iS),$$

$${\int }_V|\Psi {|}^{2}dV={\int }_V\rho dV=M,$$

$$i{\partial }_t\Psi =\left[-\frac{\mathcal{D}}{2}{\nabla }^2-b\ln \left(|\Psi {|}^2/{\rho }_0\right)\right]\Psi ,$$

$${\partial }_t\rho +\nabla \cdot (\rho \text{u})=0,$$

$${\partial }_t\text{u}+\text{u}\cdot \nabla \text{u}-\frac{1}{\rho }\nabla \cdot \mathbb{T}=0,$$

$$\mathbb{T}=-\frac{{\mathcal{D}}^2}{4\rho }\nabla \rho \otimes \nabla \rho -\tilde{p}\mathbb{I},$$

$$\nabla \mu =-\frac{\eta }{\mathcal{D}\rho }\nabla \cdot \mathbb{T}=-\eta \nabla \left[\frac{\mathcal{D}}{2}\frac{{\nabla }^2\sqrt{\rho }}{\sqrt{\rho }}+b\ln \left(\frac{\rho }{{\rho }_0}\right)\right],$$

$$\mathcal{U}=⟨\hat{\mathrm{H}}⟩=\mathcal{F}+T\mathcal{S},$$

$$T\sim b\sim {\left(\frac{\partial \mathcal{U}}{\partial \mathcal{S}}\right)}_V,$$

$$\mathcal{L}=\frac{i}{2}(\Psi {\partial }_t{\Psi }^{\ast }-{\Psi }^{\ast }{\partial }_t\Psi )+\frac{\mathcal{D}}{2}|\nabla \Psi {|}^2+\mathcal{V}(|\Psi {|}^2),$$

$$\mathcal{V}(\rho )=-b\rho \left[\ln (\rho /{\rho }_0)-1\right]+{\mathcal{V}}_0,$$

$${\Psi }_{(g)}(\text{x},t)=\pm \sqrt{{\rho }_g(\text{x})}\exp \left(-i\omega t+i\text{k}\cdot \text{x}\right),$$

$${\rho }_g(\text{x})$$

$${\omega }_g$$

$${ϵ}_g(R)\propto \frac{1}{\ell }(R-{\ell }_0)\exp \left[-(R/\ell {)}^2\right]\left[1+\mathcal{O}\left(R-\ell \right)\right],$$

$$\mathbb{U}\left(r\right)\propto {ϵ}_g\left(r\right),$$

$$\mathbb{U}(r)=\frac{U_0}{\ell }(r-{\ell }_0)\exp \left[-(r/\ell {)}^2\right]+\mathcal{O}\left(r-\ell \right),$$

$$\Psi (\text{x},t)=\prod _{j=1}^{\bar{d}}{\psi }_j(x_j,t),$$

$${\int }_{X_j}|{\psi }_j{|}^{2}d{x}_j={\int }_{X_j}{\varrho }_{j}d{x}_j=M^{1/\bar{d}},$$

$$i{\partial }_t{\psi }_j=\left[-\frac{\mathcal{D}}{2}{\partial }_{x_j{x}_j}^2+|b|\ln \left(|{\psi }_j{|}^2/{\varrho }_0\right)\right]{\psi }_j,$$

$$\frac{\mathcal{D}}{2}\frac{d^2\psi }{d{x}^2}=\frac{d\bar{\mathcal{V}}(|\psi {|}^2)}{d\psi },$$

$$\frac{d\psi }{dx}=\pm \sqrt{2\tilde{U}(\psi )},$$

$$Q={\varrho }_0^{-1/2}[\psi (+\infty )-\psi (-\infty )],$$

Full text

††preprint: Europhys. Lett. (EPL) 122, 39001 (2018)  [DOI: 10.1209/0295-5075/122/39001]

An active volcano can be modeled as a hydrodynamic system consisting of a magma reservoir or chamber which is connected to a long conduit. Before eruption, magma is accumulated inside pockets under the Earth’s crust. Changes in pressure trigger various chemical and physical processes, such as cavitation and gas diffusion, which fragment magma, form a foam, eventually leading to breaking of the crust, followed by a volcanic eruption. The ascent of magma is usually described as a steady one-dimensional flow of an isothermal compressible fluid in a porous medium driven by pore compaction and dilation fo84 ; mc84 ; ss84 ; ssw85 ; toc90 .

The use of barotropic Korteweg capillary fluid as a model of magma transport in a conduit was recently proposed in works dmf03 ; gl08 . In this class of theories, commonly referred as diffuse interface models, the capillary interface is viewed as a diffusion transition domain of rapid smooth variation of density, while surface tension is intrinsically incorporated ds85 ; dk93 ; amw98 . This allows us to describe flows with spontaneous nucleation, coalescence and breakdown of density inhomogeneities in two-phase fluids an96 .

Wave equation. – We consider a magmatic fluid flowing in a cylindrical channel, which is in thermal contact with a reservoir of infinitely large heat capacity so as to maintain constant temperature. We assume also that this fluid is formed by partial melting of solids with a mostly low thermal conductivity, such as silica. Therefore, its microscopic structure can be regarded as a many-body system of particles, atoms and molecules, whose average potential energy is larger than kinetic. Then the probability $P$ of a microstate in this system is given by a Boltzmann-type rule, in which kinetic energy can be neglected compared to potential:

[TABLE]

where $T$, ${\cal E}$ and ${\rm U}$ are, respectively, the temperature, energy and potential energy of a many-body system; here we work in units where the Boltzmann constant $k_{B}=1$.

At a microscopic level, this many-body system is described by a rather large set of particles’ evolution equations, but employing the fluid approximation used in continuum mechanics, we can make a transition to a single equation for the fluid wavefunction in a Madelung form:

[TABLE]

where $\rho=\rho(\textbf{x},t)$ is a fluid density, and $S=S(\textbf{x},t)$ is a phase which is related to fluid velocity via the gradient: $\textbf{u}\sim\boldsymbol{\nabla}S$ (we assume irrotational flow). This wavefunction should not be confused with particle wavefunctions in quantum mechanics, but rather it is a complex-valued function which stores macroscopic information about the underlying many-body system, such as the density and velocity of a fluid element or parcel ry99 (another wave-mechanical analogy can be found in, e.g., classical electromagnetic wave theory sybook09 ; z16 ). Nevertheless, some mathematical similarities between these categories of functions do exist and will be used here. Specifically, the function $\Psi$ obeys a normalization condition

[TABLE]

where $M$ and $V$ are the total mass and volume of the fluid. If we consider that this function is also a solution of a wave equation, this poses restrictions, which are somewhat similar to those of a quantum-mechanical case: a set of all normalizable fluid wave functions must constitute a Hilbert space, such as $L^{2}(\mathbb{R}^{{\bar{d}}})$ where ${\bar{d}}$ is the number of spatial dimensions of the fluid.

To derive an equation for $\Psi$, we correlate the probability sample space of an underlying microscopic many-body system with the configuration space of fluid’s degrees of freedom, i.e., one can assume a correspondence between probability (1) and fluid density: $P\sim\rho$ hence $|\Psi|^{2}\sim\exp{\left(-{\rm U}/T\right)}$. From the latter formula, a general expression follows for an operator of the potential ${\rm U}$: $\hat{\rm U}=-(T-T_{0})\ln{(\rho/\rho_{0})},$ where $\rho_{0}$ and $T_{0}$ are reference values of, respectively, fluid density and temperature.

Assuming that our flow is Hamiltonian, we consider an operator form for its total energy, which can be written as an identity: $\hat{\rm H}=\hat{\rm T}+\hat{\rm U}$, where $\hat{\rm H}=i\eta\,\partial_{t}$ is a Liouville-type generator of time translations, $\eta$ being a scaling dimensionful parameter (we can assume it positive for definiteness), $\hat{\rm T}\propto\frac{1}{2}\hat{\textbf{p}}\cdot\hat{\textbf{p}}$ is a kinetic energy operator, $\hat{\textbf{p}}\propto-i\boldsymbol{\nabla}$ is a generator of spatial translations, and $\hat{\rm U}$ is a potential energy operator derived above. These operators act in the above-mentioned Hilbert space of states described by rays $|\Psi\rangle$, therefore, this identity can be recast in the form $\hat{\rm H}|\Psi\rangle=(\hat{\rm T}+\hat{\rm U})|\Psi\rangle$, which brings us to the logarithmic wave equation of a Schrödinger  type (where the role of a Planck constant is played by $\eta$):

[TABLE]

where $b=(T-T_{0})/\eta$ and ${\cal D}$ are real constants. The latter is a material parameter of the theory related to surface tension; in wave-mechanical context it acts as a measure of the inertia of a capillary flow. It becomes apparent that $T_{0}$ is the critical temperature at which the coupling $b$ changes its sign, so it is related to the phase transition discussed below.

The wave equation must be supplemented with the normalization condition (3). Besides, if the fluid flows inside some kind of vessel or external potential $V_{\text{ext}}(\textbf{x},t)$, such as Earth’s gravity, then a corresponding term $\eta^{-1}V_{\text{ext}}(\textbf{x},t)\,\Psi$ must be added to a right-hand side of eq. (4). Since we are assuming a leading-order approximation here, we can neglect geometrical constraint’s effects and consider a trapless flow, i.e., $V_{\text{ext}}\equiv 0$ for now.

One can easily verify that by substituting eq. (2) into (4) one recovers hydrodynamic laws for mass and momentum conservation for a two-phase compressible inviscid fluid with internal capillarity whose flow is irrotational and isothermal dmf03 ; gl08 :

[TABLE]

with $\textbf{u}={\cal D}\boldsymbol{\nabla}S$ and the stress tensor $\mathbb{T}$ of the Korteweg form with capillarity ds85 , which is used to model fluid mixtures with phase changes and diffuse interfaces amw98 ; an96 :

[TABLE]

where $\mathbb{I}$ is the identity matrix, $\tilde{p}=p(\rho)-\frac{1}{4}{\cal D}^{2}\boldsymbol{\nabla}^{2}\rho$ is capillary pressure, and $p(\rho)=-{\cal D}b\rho$ is a barotropic equation of state for fluid pressure $p$. The stress tensor is also related to the chemical potential $\mu$ through the formula:

[TABLE]

which indicates that the chemical potential can be directly determined from eq. (4) using a stationary state ansatz $\Psi(\textbf{x},t)=\exp{\left(-i\mu t/\eta\right)}\Psi(\textbf{x})$.

By using a standard averaging procedure with respect to the inner product in the Hilbert space of the fluid wavefunctions $\Psi$, one can show that an averaged form of eq. (4) can be written as a formula for a wave-mechanical internal energy of the fluid:

[TABLE]

where ${\cal F}=\langle\hat{{\rm T}}\rangle=-\frac{{\cal D}}{2M}\int_{V}\Psi^{*}\boldsymbol{\nabla}^{2}\Psi dV$ is a wave-mechanical free energy, the temperature is counted with respect to a reference value $T_{0}$, and ${\cal S}=-\frac{1}{M}\int_{V}\rho\ln{\!(\rho/\rho_{0})}\,dV$ is an entropy per mass; according to aforesaid ${\cal S}\sim-\int_{V}P\ln P\,dV$. Thus, $T$ and ${\cal S}$ must be thermodynamically conjugated:

[TABLE]

if one regards ${\cal U}$ as a thermodynamical potential.

In summary, eq. (4) is a concise form of writing hydrodynamic equations for Korteweg materials, which makes it a useful model in studies of the Korteweg-type magmas, especially considering the substantial amount of information accumulated so far about the properties of wave equations with logarithmic nonlinearity in the different branches of physics, see works ros68 ; bb76 ; z10gc ; dz11 ; dmz15 ; szm16 ; zrz17 and references therein, including the theory of logarithmic Bose-Einstein condensates gg15 ; az11 ; z12eb ; bo15 ; z17zna .

Symmetry breaking and phase transitions. – Using a variational approach, eq. (4) can be derived as an Euler-Lagrange equation for Galilean-invariant Lagrangian density:

[TABLE]

with potential density given by

[TABLE]

where ${\cal V}_{0}={\cal V}(0)$ is a constant whose value can be chosen to ensure that the potential equals zero at the local minima; in our case we can set ${\cal V}_{0}=\rho_{0}(|b|-b)/2$. For complex values of $\Psi$, this potential has nontrivial local extrema at $|\Psi_{e}|=\sqrt{\rho_{0}}$, which can be minima or maxima depending on a sign of $b$, cf. Fig. 1. This indicates a possibility of symmetry breaking (or restoration), and the existence of at least two phases in the fluid which correspond to the positive and negative values of $b$.

It should be noted that for the same values of $b$ and $\rho_{0}$, and for the same set of boundary conditions, eq. (4) allows multiple normalized eigensolutions. This can be observed by studying stationary solutions, because propagating solutions can be always generated by means of the Galilean transformation $\textbf{x}\to\textbf{x}-\textbf{u}\,t$. The equation for a stationary fluid wavefunction, $\Psi(\textbf{x},t)=\exp{\left(-i\omega t\right)}\Psi(\textbf{x})$, indicates the existence of multiple states in the Hilbert space of the system, which correspond to different eigenvalues of the wave frequency $\omega=\mu/\eta$. As discussed after eq. (2), this means that the Korteweg-type magmatic fluid spontaneously selects one of the possible eigensolutions, which makes it somewhat similar to quantum liquids. However, the ground state is still the preferred one, as it corresponds to a minimum of $\omega$ which is also a minimum of chemical potential.

These aspects will be further discussed in the next two sections, where we consider the phase structure of the model (4) in details.

Cellular phase. – If temperature of Korteweg magma fluid satisfies the condition $T>T_{0}$ then the nonlinear coupling $b$ is positive, hence the potential density ${\cal V}$ has an upside-down Mexican-hat shape, with local degenerate maxima at $|\Psi_{e}|=\sqrt{\rho_{0}}$, cf. solid and dashed curves in Fig. 1.

In this case, one solitary wave solution of eq. (4) can be found analytically. It corresponds to the ground state (i.e., the one with a lowest eigenvalue of the frequency $\omega$) and can be written in the form of a Gaussian parcel modulated by a plane wave:

[TABLE]

where k is a constant vector, and the density’s normalized eigenfunction and frequency eigenvalue are, respectively:

[TABLE]

where $\tilde{\rho}=M/\tilde{V}$, $\tilde{V}=\pi^{{\bar{d}}/2}\ell^{\bar{d}}$ and $\ell=\sqrt{|{\cal D}/(2b)|}$ are the density peak value, effective volume and Gaussian width, respectively; $\textbf{k}=0$ if $b\not=0$.

Therefore, in this phase, Korteweg-type magma tends to fragment into clusters of density inhomogeneities with a Gaussian profile, referred here as cells. One can effectively treat these cells as classical point particles, while encoding their nonzero size in the two-body interaction potential $\mathbb{U}\left(|\textbf{x}-\textbf{x}^{\prime}|\right)$. The latter can be approximately derived in the following way.

A spherical cell of a size $R\sim\ell$ stores an amount of internal bulk mass-energy $\epsilon_{g}(R)\propto\int_{0}^{R}\rho_{g}(r^{\prime})r^{\prime 2}dr^{\prime}$. By direct computation, we obtain

[TABLE]

where ${\ell_{0}}=\ell\left[1/2+1/(\text{e}\sqrt{\pi}\,\text{erf}(1))\right]\approx 0.75\,\ell$. Because each cell tries to maintain its size and mass when interacting with its environment, to alter these values an amount of energy must be supplied which is proportional to $\epsilon_{g}$. This energy can transport only through the interactions with other cells, therefore

[TABLE]

where $r=|\textbf{x}-\textbf{x}^{\prime}|$ is a distance between centers of mass of cells. We introduce the proportionality factor $U_{0}$, assume it to be constant in the leading-order approximation, and obtain the cells’ two-body interaction potential function:

[TABLE]

where ${\cal O}$-terms can be discarded unless interactions deform cells so strongly that their interior structure should be taken into consideration. The potential (18) can be further used in many-body simulations of the Korteweg-type materials including magmas in the “liquid-dissolved gas” phase, as discussed in the conclusion.

Foam phase. – If temperature $T<T_{0}$ then $b$ turns negative, and the potential density ${\cal V}$ acquires a conventional Mexican-hat shape, with local degenerate minima at $|\Psi_{e}|=\sqrt{\rho_{0}}$, cf. dotted and dash-dotted curves in Fig. 1. Therefore, topologically nontrivial solitons must exist which interpolate between the local minima and maximum at $\Psi=0$.

To find these solitons explicitly, let us look for solutions of eq. (4) in a form of the product

[TABLE]

where each of $\psi$’s being normalized separately:

[TABLE]

where $\varrho_{j}=|\psi_{j}|^{2}$ and $X_{j}$ are, respectively, linear density and extent of the fluid along the $j$th coordinate. Then, due to its separability in Cartesian coordinates, eq. (4) decomposes into ${\bar{d}}$ identical 1D equations of the form:

[TABLE]

where $\varrho_{0}=\rho_{0}^{1/{\bar{d}}}$, and $\partial_{x_{j}}$ is a spatial derivative with respect to $j$th coordinate. Therefore, in this case we deal with only one 1D differential equation, so the index $j$ can be omitted for brevity. In view of a Galilean invariance, we resort to a static case $\psi=\psi(x)$ and write this equation in the form

[TABLE]

where $\bar{\cal V}(f)=|b|f\left[\ln{(f/\varrho_{0})}-1\right]+\bar{\cal V}_{0},$ $\bar{\cal V}_{0}=|b|\varrho_{0}$, and $x$ refers to any of the coordinates $x_{1},...,x_{\bar{d}}$. According to the standard approach Rajaraman:1982is , our solitons must saturate the Bogomolny-Prasad-Sommerfield (BPS) bound, which means that they must be solutions of the first-order differential equation

[TABLE]

where $\tilde{U}(\psi)=(2/{\cal D})\bar{\cal V}(|\psi|^{2})$ is a soliton particle potential, and a sign $\pm$ refers to soliton and anti-soliton solutions. While the analytical solution for this case is unknown, numerical study reveals the existence of topological solitons of a kink type, see Fig. 2. These solutions have a nonzero topological charge,

[TABLE]

which enhances their stability against decay into a trivial state $\psi=0$. Therefore, all nonsingular finite-energy solutions of eq. (23) can be cast into four topological sectors gg15 . Two of these sectors have a nonzero topological charge which ensures stability of corresponding BPS solitons.

Since the density of each of kink solitons grows from the center of mass outwards, cf. dashed and dash-dotted curves in Fig. 2, these solitons are viewed as models of bubbles with a characteristic size $\ell$. While in a single-soliton system these solutions would extend for the whole $x$-axis, in the real fluid the kinks would match with antikinks at distances of order $\ell$. Therefore, Korteweg-type magma in this phase tends to form a foam, thus starting a process of releasing the previously dissolved gas.

Conclusions. – We have used statistical mechanics arguments to derive a flow equation for Korteweg-type fluid, which contains logarithmic nonlinearity making it similar to wave equations used in the theory of dense Bose-Einstein condensates. In the classical hydrodynamics context, this equation defines a diffuse interface model with spontaneous symmetry breaking which causes phase transitions. It allows multiple (eigen)solutions whose set forms a Hilbert space of different states, in which the Korteweg fluid can be found. Therefore, as temperature changes, fragmentation and nucleation of inhomogeneities occur in erupting magma, followed by a transition between the “magma-dissolved gas” fluid phase and the magmatic foam which flows up the conduit. We have studied this spontaneous symmetry breaking phenomenon and related phase structure, by both analytical and numerical methods.

For the cellular phase, we recover a well-known solitary wave solution with a Gaussian density profile, which describes localized cell-like density inhomogeneities in magma. For such cellular structures an effective interaction potential and many-body Hamiltonian are derived, which can be used for statistical simulations of magma and other Korteweg-type materials. The derived many-body interaction potential (18) between Gaussian-shaped density inhomogeneities in the liquid phase of Korteweg fluid suggests a direction for further statistical mechanics’ studies of cellular structures in magmas, where these inhomogeneities are effectively regarded as point particles with nonsingular interactions.

Considering the foam phase, we demonstrate existence of the topological solitons of a kink type which describe bubbles or other inhomogeneities with density increasing from the center of mass outwards, up to distances of order $\ell$. We show that these solitons are saturating BPS bound and belong to a topological sector distinct from the one containing the trivial solution. Thus, their stability is ensured by the conservation of a topological charge. The foam-like structure arises in the fluid due to the matching between kinks and antikinks at distances of order $\ell$; this structure facilitates release of any gas if it was previously dissolved in the fluid.

Furthermore, studies of the logarithmic model’s Hilbert space and corresponding eigenfunctions and eigenvalues will shed light upon transitions in Korteweg fluids in general and appropriate kinds of magma in particular.

---