Skip to main content


Front. Appl. Math. Stat., 23 June 2023
Sec. Dynamical Systems
Volume 9 - 2023 |

Bifurcation of solutions through a contact manifold in bidisperse models

  • 1Department of Information Systems, Universidad del Bío-Bío, Concepción, Chile
  • 2Academic Department of Mathematics, ITAM, Ciudad de México, Mexico

This research focuses on a hyperbolic system that describes bidisperse suspensions, consisting of two types of small particles dispersed in a viscous fluid. The dependence of solutions on the relative position of contact manifolds in the phase space is examined. The wave curve method serves as the basis for the first and second analyses. The former involves the classification of elementary waves that emerge from the origin of the phase space. Analytical solutions to prototypical Riemann problems connecting the origin with any point in the state space are provided. The latter focuses on semi-analytical solutions for Riemann problems connecting any state in the phase space with the maximum packing concentration line, as observed in standard batch sedimentation tests. When the initial condition crosses the first contact manifold, a bifurcation occurs. As the initial condition approaches the second manifold, another structure appears to undergo bifurcation, although it does not represent an actual bifurcation according to the triple shock rule. The study reveals important insights into the behavior of solutions in relation to these contact manifolds. This research sheds light on the existence of emerging quasi-umbilic points within the system, which can potentially lead to new types of bifurcations as crucial elements of the elliptic/hyperbolic boundary in the system of partial differential equations. The implications of these findings and their significance are discussed.

1. Introduction

Polydisperse suspensions can be effectively modeled using balance equations that consider N continuous phases, where particles of species i with a volume fraction ϕi possess distinct properties such as size, density, and viscosity [1, 2]. A similar solution structure is applicable for modeling traffic and pedestrian flows as well [3, 4].

In the case of the bidisperse suspension model under consideration, the solution structure of the initial value problem for standard batch settling tests has been studied for situations where strict hyperbolicity is ensured [5], as well as when elliptic regions exist within the phase space [6]. This contribution focuses on the influence of a contact manifold within the phase space, which emerges under specific parameter settings. This contact manifold possesses the physical characteristic of coinciding particles settling velocities, which is of practical relevance since one common objective in the process control of solid-liquid separation processes is to minimize segregation effects [7]. At these optimal states, the phenomenon of solution bifurcation emerges.

The general form of kinematic sedimentation models for polydisperse suspensions consists of a system of N first-order hyperbolic equations. In this work, we consider a bidisperse suspension (N = 2) which exhibits a contact manifold capable of producing solution bifurcations. Specifically, the 2 × 2 system is given by:

tϕi+xfi(Φ)=0,fi(Φ)=ϕivi(Φ),i=1,2,    (1)

where t represents time, x denotes depth, and the velocity components vi(Φ) depend on the concentration vector Φ=(ϕ1,ϕ2)T. The vector Φ represents the volume fractions of the solid phases and is constrained within the phase space of physically relevant concentrations

DΦ:={Φ=(ϕ1,ϕ2)T2|ϕ1,ϕ20,ϕ:=ϕ1+ϕ2ϕ},    (2)

where the total concentration ϕ: = ϕ1 + ϕ2 is upper-bounded by the maximum packing concentration ϕ. The maximum packing manifold, denoted by

:={Φ=(ϕ1,ϕ2)T|ϕ:=ϕ1+ϕ2=ϕ},    (3)

represents the set of all maximal states.

The resulting system of conservation laws represents a mass balance for different solid species, where the non-linear flux function


can be derived from the corresponding momentum balances [1, 2, 8]. These components describe the flow process of dispersed solid phases within a liquid, considering the dispersed phases as a continuum.

The solid-fluid relative (“slip”) velocities

ui(Φ):=vi(Φ)-vf,    (4)

measure the absolute velocities vi = vi(Φ) relative to the fluid velocity vf. The total volume displacement through horizontal sections is null:


By factoring out the slip velocities ui, we can rewrite this mass balance equation as an equivalence, stating that if the solid volume is moving downward, then the fluid must move upward:

vf(Φ)+ϕ1(v1(Φ)-vf)+ϕ2(v2(Φ)-vf)=0vf(Φ)=-ΦTu(Φ)    (5)

The absolute velocities vi = vi(Φ) of representative solid particles depend on a linear combination of the relative velocities. This leads to the flux function with components

fi(Φ)=ϕivi(Φ),  vi(Φ)=ui(Φ)-ΦTu,i=1,2,    (6)

with u=(u1(Φ),u2(Φ))T. To close this flux function model, the slip velocities are further specified as

ui(Φ)=viVi(Φ),    (7)

where vi represents the Stokes velocity, quantifying the settling velocity of an individual particle in the fluid. Additionally, Vi(Φ) denotes the hindered-settling velocity, which is a non-increasing function of the components of Φ as discussed in [9]. Following Richardson and Zaki [10], the hindered-settling velocity is defined as

Vi(Φ):={(1ϕ)ni1if 0ϕϕ,0otherwise,i=1,2,    (8)

where ni > 1 accounts for the slowing down of the process at increasing concentrations. Combining assumptions (7) and (8) gives

ui(Φ)=vi(1-ϕ)ni-1    (9)

for ΦDΦ.

Strictly hyperbolic systems of conservation laws, characterized by real and distinct eigenvalues of the Jacobian matrix of the flux function, provide a well-established framework for solving Riemann problems [11]. The system (1) with flux function (6) has been proven to be strictly hyperbolic for N = 2 in [2] and later extended to general N in [1]. However, these results are limited to hindered-settling factors that depend solely on the total concentration ϕ. It was previously believed that this strict hyperbolicity required constant exponents n1 = n2 = ⋯ = nN, but [9] demonstrated that this restriction is unnecessary.

In [9], it was shown that strict hyperbolicity holds for general N≥2 and hindered-settling factors of the form Vi(Φ) = Vi(ϕ), as long as the inequality

ui(1-ϕ)-ui<0    (10)

holds for all i = 1, …, N, where ui represents the derivative with respect to ϕ. In this case, the only requirement on the hindered-settling function Vi(Φ) = Vi(ϕ) is that it depends on the total concentration ϕ. The inequality (10) is satisfied when the relative velocities are ordered as u1 > u2 > ⋯ > un for any ϕ. For the hindered-settling function (8), this condition holds when the parameters are ordered as

v1>v2>>vNwith n1<n2<<nN.    (11)

The strict hyperbolicity of the model can be verified using a secular equation framework introduced in [12] and adapted to size-dependent hindered-settling factors, not necessarily of the form (8), in [13]. This framework was further applied in [14] to analyze various hindered-settling functions.

Preliminary numerical simulations of Riemann problems for N = 3, with arbitrary parameter choices, indicate abrupt behavior when strict hyperbolicity fails, leading to the coincidence of eigenvalues. Across the contact manifold, there is a clear change in the solution structure with distinct bifurcations. Solutions of Riemann problems RP(0, Φ) for different Φ on opposite sides of the coincidence plane do not exhibit smooth transitions. However, when visualized in a 3D phase space, the assumed values traverse different planes and axes. The occurrence of this phenomenon of abrupt change, appearing even for N = 2, prompted our search for analytical insights, culminating in the present contribution.

This contribution focuses on the wave classification of 2 × 2 systems of conservation laws that arise as one-dimensional kinematic models for the sedimentation of bidisperse suspensions. Analyzing bidisperse suspensions provides insights into the flow properties of polydisperse suspensions, which are mixtures of small solid particles dispersed in a viscous fluid. Specifically, this study concentrates on the properties of the 2 × 2 system (1), flux function (6), and the closures (7) and (8). The model under consideration involves the following specifications in opposition to (11), that ensures the preservation of strict hyperbolicity:

(S1) v∞1 > v∞2 > 0,

(S2) n1 > n2 > 1,

(S3) ϕ ≡ 1.

For Vi(ϕ) given by (8), the flux function f(Φ) takes the form


for values ΦDΦ and f1(Φ) = f2(Φ) = 0 otherwise. A special interest lies in the classification of the solution structure of the Riemann problem

Φ(t=0,x)={Φ if x<0,Φ+ if x>0.    (12)

For convenience, we refer to the Riemann problem consisting of the system of PDEs (1) with initial condition (12) as RP(Φ, Φ+). Here, Φ and Φ+ denote the specified left and right values, respectively.

The model is applied to the batch settling process of an initially homogeneous suspension in a closed container, characterized by the initial-boundary value problem:

tϕi+xfi(Φ)=0,i=1,2,Φ(0,x)=Φ0(x),0xL,    (13)
fi(Φ)=0,x{0,L},i=1,2,    (14)

here, L is the height of the domain, and the flux-density vector components are given by (6). Since the zero-flux boundary condition (14) holds, the initial-boundary data (13) and (14) can be replaced by the Cauchy data

Φ(0,x)=Φ0(x)={Ofor x<0,Φ0for 0xL,Φfor x>L,    (15)

where O: = (0, 0)T denotes the origin, and Φ represents a state on the maximum concentration manifold (3). Consequently, the Riemann problems RP(O, Φ0) and RP(Φ0,Φ) are of particular interest.

This contribution offers analytical insights into the solution structure of the Riemann problem RP(O, Φ), which is of particular interest for describing the interactions at the upper interface between clear liquid and an initially homogeneous suspension during a batch settling process.

Previous studies have examined weakly hyperbolic systems, which are hyperbolic but not strictly hyperbolic, as models for multi-phase flow in porous media. For instance, the Corey model with convex permeability exhibits a single isolated point, known as the umbilic point, where strict hyperbolicity fails [1519]. Another example of hyperbolicity loss occurs when the phase space contains an elliptic region, as observed in the case of the Stone model [20]. Coincidence along a curve can also lead to the loss of strict hyperbolicity, as shown in [21] for a 2 × 2 system.

To solve Riemann problems for weakly hyperbolic systems, the wave-curve method has been applied, which involves connecting a sequence of elementary waves. When strict hyperbolicity fails, the method of Liu for dealing with non-convex fluxes is not sufficient, and non-local branches of the Hugoniot locus must also be considered [22, 23]. The wave-curve method has been used in various contexts, such as the injection problem, which involves injecting a gas-water mixture into a porous medium containing oil [24], foam injection [25], and gravity effects [26].

For systems with a quadratic flux function, the solution in the vicinity of the umbilic point has been classified [18, 19, 27], and four distinct types of umbilic points associated with different shapes of nearby integral curves have been identified [19]. In the Corey model with convex permeability, only two types of umbilic points occur [18]. Based on this classification, certain types of Riemann problems have been considered in [16].

For non-strictly hyperbolic systems of two conservation laws that possess an umbilic point and an identity viscosity matrix, a systematic classification of solutions to the Riemann problem has been conducted [28, 29]. Non-classical waves and transitional shocks can arise due to a non-local Hugoniot locus, as shown in [3032], and these transitional shocks are influenced by the regularization resulting from a non-identical viscosity matrix. It is worth noting that these types of shocks also occur on the same local branches.

Lastly, the concept of quasi-umbilic points has recently emerged in the context of porous media [26, 33]. it is briefly discussed in the Appendix (due to its only weak connection to RP(O, Φ), and in order to not divert the attention from the primary focus on solution construction and bifurcations).

2. Basic definitions

In this section, several definitions [26, 27, 34] are provided in order to facilitate the appropriate classification of the system under study.

Definition 1. The Hugoniot locus of a state Φ, denoted as H(Φ-), is the set of all states Φ+ that satisfy the Rankine-Hugoniot condition

f(Φ+)-f(Φ-)=σ(Φ+-Φ-),    (16)

where σ = σ, Φ+) is the propagation velocity of the discontinuity.

The shock classification according to Lax [35] is employed to refer to a subset of the Hugoniot locus that corresponds to a certain wave family. Admissible shocks within this classification are classified based on the magnitude of the shock speed with respect to the first and second eigenvalues on both sides of the discontinuity, as well as the speed of the discontinuity itself, see, e.g., [28, 35].

Definition 2. Three types of classical admissible shocks can be distinguished. The classification applies to a left state Φ and a right state Φ+ connected by the Rankine-Hugoniot condition (16) with a jump velocity σ = σ, Φ+):

1-Lax shock: λ1(Φ+)σλ1(Φ-) and σλ2(Φ+),

2-Lax shock: λ2(Φ+)<σ<λ2(Φ-) and λ1(Φ-)σ,

Over-compressive shock (OC): λ2(Φ+)<σ<λ1(Φ-).

The definition of left- and right-characteristic shocks within this classification encompasses cases where the shock speed coincides with a characteristic speed (see Figure 1). While an over-compressible shock cannot be characteristic according to this definition, the limit of the aforementioned inequalities is included for 1-Lax shocks. On the other hand, for 2-Lax shocks, strict inequalities are required to avoid any ambiguities.


Figure 1. Schematic diagram of shock classification. Each type, as specified in Definition 2, is encompassed by cases shown with light shading. The cases corresponding to classic Lax classification are depicted with dark shading. The special case of shock velocity σ = 0 is included, where the classification relies on inequalities and equalities involving λ1-=λ1(Φ-), λ2-=λ2(Φ-), λ1+=λ1(Φ+), and λ2+=λ2(Φ+). The blue lines represent slow characteristic speeds, while the red lines represent fast characteristic speeds.

It is important to note that this extension deviates from the conventional approach presented in the seminal work of [28] or summarized in Table 1 of [36]. However, we have opted for this unconventional extension in order to maintain the continuity of the 1-Lax branch in the classification of the Hugoniot locus of the origin, as depicted in Figure 4 and established in Theorem 2. This choice facilitates a smoother and more coherent discussion in the subsequent analysis.

An inflection manifold Ii is determined for all states Φ where the i-th eigenvalue achieves either a maximum or minimum value along the integral curve of the same family.

Definition 3. The i-th inflection manifold is defined as


where λi is the i-th eigenvalue of the Jacobian matrix of the flux function and ri is the non-null corresponding eigenvector.

the invariant manifolds defined in [37], we introduce the following concept:

Definition 4. An i-th contact manifold occurs when the i-th integral curve passing through a state Φo coincides with a part of the Hugoniot locus H(Φo), such that any state Φ on this intersection satisfies


for the shock speed σo, Φ).

A necessary condition for establishing an i-th contact manifold is that the integral curve is not a rarefaction in the usual sense, but rather that the characteristic speed remains fixed along the curve.

For states on a contact manifold the following transitivity rule holds: if Φ1 and Φ2 mutually belong to the Hugoniot locus of the other, connected by a shock of speed σ = σ1, Φ2), and if Φ2 is on the Hugoniot locus of a state Φ3 by a shock of the same speed σ, then Φ1 belongs to H(Φ3) and σ1, Φ3) also coincides with σ. This fundamental property is known as the triple shock rule [15, 23, 34]. Another useful version of this rule is stated as follows.

Lemma 1. Let Φ1, Φ2, Φ3 be non-collinear states such that Φ1, Φ2 belong to H(Φ3) and Φ1 belongs to H(Φ2), then σ1, Φ2) = σ2, Φ3) = σ1, Φ3).

On a contact manifold, rarefactions and shocks are indistinguishable since all speeds align with a characteristic speed [37]. Remarkably, for the Hugoniot locus of the origin, the constant characteristic speed corresponds to a straight line and coincides with the integral curves.

3. Contact manifold

If the specifications (S1) and (S2) of the considered model hold, then it is possible to identify a contact manifold within the phase space. This manifold turns out to be decisive in characterizing solutions of Riemann problems, in particular because the origin is connected to this manifold by means of a right characteristic 1-Lax shock.

Definition 5. Given a state Φ=(ϕ1,ϕ2)DΦ, such that v1(Φ)=v2(Φ), the following subset of phase space which contains Φ is defined

C(Φ):={ΦDΦ:v1(Φ)=v2(Φ)=v1(Φ)}.    (17)

The state Φ serves as representative of the set C(Φ). The definition of a set C(Φ) unifies two complementary properties:

(1) v1(Φ) = v2(Φ) for all ΦC(Φ),

(2) vi(Φ-)=vi(Φ+) for Φ-,Φ+C(Φ) and i = 1, 2.

Property (1) denotes the local coincidence of velocities between different phases within a specific state, while property (2) describes the constancy of velocities within a single family along the manifold. Actually, a set C(Φ) is a contact manifold for a broad class of flux functions that have a the considered structure. All the states belonging to this contact manifold share the same Hugoniot locus.

Lemma 2. If the flux function of the system (1) has the structure

fi(Φ)=ϕivi(Φ),    (18)

then any set C(Φ) is a contact manifold with constant shock speed

σ(Φ-,Φ+)=v1(Φ).    (19)

Moreover, if the structure of the flux function (18) is such that the absolute velocity v depends only on the total concentration ϕ, namely

vi(Φ)=vi(ϕ),    (20)

then the contact manifold C(Φ) consists of the line

C(Φ)={ΦDΦ:|ϕ1+ϕ2=ϕ},    (21)

where ϕ=ϕ1+ϕ2 is the total concentration of Φ.

Proof. According to the definition in equation (17), any states Φ-,Φ+C(Φ) satisfy that


such that one can factorize

ϕi+vi(Φ+)-ϕi-vi(Φ-)=v1(Φ)(ϕi+-ϕi-),i=1,2.    (22)

From the specific structure of the flux function (18) one recognizes that equation (22) states the Rankine-Hugoniot condition (16) with shock speed (19). Consequently, it is established that all states belonging to the set C(Φ) mutually belong to the Hugoniot locus of each other. Furthermore, any two states within a set C(Φ) can be connected by a shock of the same speed. Since shock speeds locally converge to an eigenvalue, on a manifold with constant shock speed any two states Φ and Φ+ have an eigenvalue that coincides with the shock speed,


This result confirms the definition of a contact manifold as stated in Definition 4. The line-shaped nature of the manifold C(Φ) can be deduced from property (20), which ensures that the absolute velocities remain constant along the line in equation (21), that is, v1(Φ) = v2(Φ) for all ΦC(Φ).

With this Lemma, nothing has been said yet about the existence of such a contact manifold for the considered model equation. Neither it is decided to which characteristic family the contact manifold belongs.

A further characterization of a set C(Φ) is developed in the sequel, starting from properties that can be derived from the generic structure of the model, leading to properties that depend on particular model specifications. A property that can be used in several instances is

v1(Φ)=v2(Φ)    u1(Φ)=u2(Φ),    (23)

which holds directly from model setup (6) and the combined condition (9), together with vf(Φ)=-ΦTu [Here, property (23) takes the role of a generalized condition for (20)].

So far, we have shown that any set C(Φ) is a contact manifold, which additionally takes the shape of a line. Moreover, under the considered model specifications, such manifolds effectively exist.

Lemma 3. If conditions (S1), (S2), and (S3) are satisfied, then two distinct contact manifolds exist in the domain DΦ, namely

C(Φx):={Φ=(ϕ1,ϕ2)T|ϕ1+ϕ2=ϕx},    (24)
C(Φ):=={Φ=(ϕ1,ϕ2)T|ϕ1+ϕ2=ϕ}.    (25)

The manifold C(Φ) is represented by any state Φ ∈ ∂, that for handiness can be specified as Φ = (0, 1)T. A representative state Φx that defines the contact manifold C(Φx) that is distinct to the manifold C(Φ) can be identified as

Φx=(ϕx,0)T,ϕx=1-(v2/v1)1/(n1-n2).    (26)

In a geometrical interpretation, the plane C(Φx) corresponds to a crossing of speeds, and C(Φ) to a final encounter.

Proof. First, it is shown that the boundary ∂ is a contact manifold. For any state Φ ∈ ∂, if (S3) then one has u1(Φ)=u2(Φ)=0. By property (23) and since (5) holds, one obtains v1(Φ)=v2(Φ) for any state Φ ∈ ∂, satisfying the Definition (17) of the set C(Φ), which, according to Lemma 2, is a contact manifold.

To show that there is an additional contact manifold C(Φx), it is necessary to find a set of states ΦxC(Φ) such that v1(Φx)=v2(Φx). Because of property (23), this is equivalent to finding a Φx such that u1(Φx)=u2(Φx). Solving


with respect to ϕx, which is the total concentration of any representative Φx, gives (26). The conditions (S1) and (S2) on the parameters guarantee the existence of ϕx ∈ (0, 1) (Actually, from (S2) is only used n1n2).

Throughout this work, the notation C(Φ) is used to refer to a generic contact manifold, whereas C(Φx) and C(Φ) refer to specific contact manifolds with assigned representative values Φx and Φ, respectively. It turns out that there exists a peculiar Hugoniot locus in the model for the origin O as reference point. There are two local branches intersecting at the primary bifurcation given in O, these principal branches are the axes. However, the locus has four secondary bifurcation points with transverse lines (as non-local branches) joining them two by two in conjunction with these lines given by contact manifolds C(Φx) and C(Φ).

Lemma 4 (Hugoniot locus of origin). The Hugoniot locus H(O) of the origin O = (0, 0)T consists of four branches: the two coordinate axes as local branches,

1:={Φ=(ϕ,0)T,ϕ[0,1]},2:={Φ=(0,ϕ)T,ϕ[0,1]},    (27)

with variable shock speed

σ(O,Φ)=vi(Φ)for Φi,i=1,2,    (28)

and the two contact manifolds C(Φx) and C(Φ), identified by (25) and (24), as transversal branches with constant shock speed

σ(O,Φ)=v1(Φ)=v2(Φ)for ΦC(Φ).    (29)

Proof. With the structure of the flux function (6), the Rankine-Hugoniot condition (16) connecting the origin O with any state Φ takes the form


This system of two equations has two possible kinds of solution: On the local branches on the axes, ∂1 and ∂2, one of the two equations becomes obsolete, since both sides vanish on the considered axis. Therefore, the velocity of the remaining equation determines the shock speed as (28). On the transversal branches, C(Φx) and C(Φ), no such cancellation occurs. Therefore, for a solution to exist, the velocities must be equal, resulting in the speed (29).

The triple shock rule as stated in Lemma 1 applies to the connection of the origin to any state ΦC(Φ) having speed (29) with any middle state ΦMC) having speed (19). Indeed, O, Φ, ΦM are not collinear, thus we have σ(O, ΦM) = σM, Φ) = σ(O, Φ). This means that any (shock) solution O → Φ can be constructed by a first shock O → ΦM followed by a second shock ΦM → Φ of same speed; both solutions determine the same wave pattern, thus the same solution (This structural property plays a central role in the continuity of the bifurcation of solutions).

Another useful property is the convertibility of the relative velocity ui(Φ) to the absolute velocity vi(Φ) on the edges. A state Φ ∈ ∂i, i = 1, 2, on an axis has the representation Φ = ϕδi1e1 + ϕδi2e2, where ek, k = 1, 2, are unit basis vectors and δik is the Kronecker symbol. With this notation the definition of the relative velocity vi(Φ) reduces to

vi(Φ)=ui(Φ)-ΦTu=ui(Φ)-ϕui(Φ)=(1-ϕ)ui(Φ).    (30)

4. Characteristic speeds

System (1) can be written in quasi-linear form as


where RP(Φ) is the Jacobian matrix of the vector valued flux function f(Φ):=(f1(Φ),f2(Φ))T. The structure of the Jacobian matrix is examined for general N in [9, 14]; for our case it becomes

J(Φ)=(J11(Φ)J12(Φ)J21(Φ)J22(Φ))=(v1+u11ϕ1u12ϕ1u21ϕ2v2+u22ϕ2),    (31)

or componentwise as


where δij is the Kronecker symbol, vi is the absolute velocity (4), and uij is specified as

uij=uij(Φ)=ui(Φ)-ΦTu(Φ)-uj(Φ),i,j=1,2,    (32)

where u(Φ)=u1(Φ)u2(Φ)T. The short notation ui(Φ)=viVi(Φ), i=1,2 is used for the derivative of (7) [see, e.g., [9]], which, in a standard model like (9), acts as a derivative with respect to ϕ.

By definition a system is strictly hyperbolic if the Jacobian matrix of the flux function has distinct real eigenvalues. Thus, the system of conservation laws (1) is strictly hyperbolic if the discriminant


of the Jacobian matrix (31) is positive. For a hindered settling factor given by (8), strict hyperbolicity holds for identical exponents n1 = n2, see [2], which is proofed by showing algebraically that ΔΦ > 0 for Φ in the interior of the phase space DΦ under the specification (S3) with ϕ = 1. Moreover, strict hyperbolicity also holds for the case with different exponents n1n2. This can be shown by a straightforward calculation [9], which yields the discriminant composed by a sum of a square and a positive term

ΔΦ=[(n1ϕ1-1)u1(Φ)-(n2ϕ2-1)u2(Φ)]2      +4n1n2ϕ1ϕ2u1(Φ)u2(Φ).    (33)

This term is positive because of conditions (S1) and (S2) together with the bounds ϕ1, ϕ2 > 0, ϕ = ϕ1 + ϕ2 < 1 (given by definition of the invariance domain DΦ) and the definition of ui(Φ). The positivity of the discriminant ΔΦ > 0 indicates strict hyperbolicity of the system (1).

The structure of the eigensystem of the Jacobian matrix RP of the flux function (6) is derived as outlined below. The eigenvalue calculation is tedious but straightforward, leading to the following result.

Lemma 5. The eigenvalues of the Jacobian matrix (31) are calculated as


where ΔΦ is the discriminant that takes the form (33).

With the help of this eigenvalue specification, the characterization of the contact manifold can be completed, which is formulated for any set C(Φ), making it particularly applicable to the sets C(Φx) and C(Φ).

Theorem 1. The eigenvalues for any state ΦC(Φ) can be specified as

λ1(Φ)=v1(Φ)-ΔΦand λ2(Φ)=v1(Φ).    (34)

Therefore, any set C(Φ) is a contact manifold with respect to the second characteristic family.

Proof. Since any state Φ in the set C(Φ) is characterized by definition by the property v1(Φ)=v2(Φ)=v1(Φ), the eigenvalues become

λ1,2(Φ)=v1(Φ)-12[n1ϕ1u1(Φ)+n2ϕ2u2(Φ)]±12ΔΦ                =v1(Φ)-12ΔΦ±12ΔΦ,

giving the pair of eigenvalues (34). Here, the last step is justified by property (23), i.e., the equivalence of the equalities v1(Φ) = v2(Φ) and u1(Φ) = u2(Φ). Specifically, one can relate the term within the brackets to the discriminant (33) as

[n1ϕ1u1+n2ϕ2u2]2=[n1ϕ1u1-n2ϕ2u2]2+4n1n2ϕ1ϕ2u1u2                                              =ΔΦ,

which is valid in the special case when u1(Φ) = u2(Φ).

The association of the contact manifold with the second family can be seen from the fact that the eigenvalue λ2, according to the established values in (34), is constant on C(Φ), λ2(Φ)v1(Φ) for all ΦC(Φ). With the shock speed established in (19) in Lemma 2, for the connection of any two states on C(Φ), one gets


for all Φ-,Φ+C(Φ) establishing that C(Φ) is a 2-th contact manifold.

Theorem 1 applies to both contact manifold C(Φx) and C(Φ). For the line C(Φ), the contact manifold is additionally characterized with respect to the first characteristic family.

For ΦC(Φx), the smaller eigenvalue λ1 depends on the discriminant, while only the larger eigenvalue λ2(Φ)v1(Φx) is constant, Thus, a contact manifold C(Φ) is an integral curve of the second family. On the maximum packing manifold for Φ=C(Φ), where ϕ1 + ϕ2 = 1 holds, both eigenvalues vanish; λ1(Φ) = λ2(Φ) = 0.

Theorem 1 states explicit expressions for the eigenvalues on the contact manifold C(Φx) and C(Φ), which are part of the Hugoniot locus of the origin as identified in Lemma 4. The remaining eigenvalues along the Hugoniot locus of the origin can be directly evaluated from the general eigenvalues according to Lemma 5. A more comprehensive presentation of these findings can be found in [38].

For Φ on an edge ∂i, i = 1, 2, i.e., ϕi = 0 for some i = 1, 2, as specified in (27), the discriminant reduces to ΔΦ=[u3-i(Φ)(n3-iϕ3-i-1)]2, where the subindex 3−i becomes 2 on the axis ∂1 and 1 on the axis ∂2. A simple calculation shows that for Φ = (ϕ, 0)T ∈ ∂1, the eigenvalues are

λ1,2(Φ)=12{[(11)(2+(11)n1)ϕ]u1(Φ)+(1±1)u2(Φ)}                 = {u2(Φ)ϕu1(Φ),(1(1+n1)ϕ)u1(Φ),

and for Φ = (0, ϕ)T ∈ ∂2 the eigenvalues are

λ1,2(Φ) =12{(11)u1(Φ)+[(1±1)(2+(1±1)n2)ϕ]u2(Φ)}                ={(1(1+n2)ϕ)u2(Φ),u1(Φ)ϕu2(Φ).

Within these eigenvalue characterizations we can distinguish the two types

{λa(Φ) : =(1(1+ni)ϕ)ui(Φ), for  Φi, i=1,2,λb(Φ) : =u3i(Φ)ϕui(Φ), for  Φi, i=1,2.    (35)

The notation λa(Φ), λb(Φ) with new subindices instead of λ1(Φ), λ2(Φ) is used since the order

λa(Φ)=λ1(Φ)λ2(Φ)=λb(Φ)    (36)

cannot be always guaranteed. However, the notation λ1, λ2 is used whenever the order (36) can be assured, e.g., at the origin, where


Generally, order changes can occur on both axes, as can be seen from Figure 3B. The order (36) holds only in the absence of such order changes [cf., [38]]. However, the occurrence of coincidence points does not affect the subsequent classification. In the sequel, it is shown how the eigenvalues λa(Φ) and λb(Φ) assume the role of the eigenvalue of the first or the second family in dependence of the value ϕ.

4.1. Inflection curves

The highly non-linear structure of the flux function makes it difficult –if not impossible– to obtain an explicit formula for inflection manifolds [27]. However, a characterization for states in the phase space can be obtained, since the inflection points there correspond to those of scalar equations. Namely, the eigenvalue λa(Φ) has the same form as the first derivative of the scalar flux function, whereas there is no correspondence for the additional eigenvalue λb(Φ), which emerges for the 2 × 2 system. The derivative

λa(ϕ)=ni((1+ni)ϕ-2)ui(Φ)/(1-ϕ),i=1,2    (37)

is well-defined for ϕ ∈ [0, 1). From (37) we have that λa(ϕ) may vanish at ϕ = 1 and at ϕ=2/(1+ni). The first case is of less significance because it represents a state at a vertex of the phase space. We have that ϕ ∈ (0, 1) since ni is larger than one due to condition (S2). Thus, the state ϕ represents an inflection point and the states Φi such that ϕi=ϕ and ϕj = 0 for ij ∈ {1, 2} belong to an inflection manifold (On edge ∂1, there is a unique inflection point, see Lemma A1). The pertinence of the inflection points to the first characteristic family can be shown for general parameter settings.

Lemma 6. There is at least one inflection point located on each axis of the phase space DΦ, namely Φ1:=(2/(1+n1),0)T and Φ2:=(0,2/(1+n2))T. Moreover, the corresponding eigenvalues belong to the first characteristic family.

Proof. The location of the inflection points on the axes is determined by finding the zeros of the corresponding eigenvalue derivatives (37).

The attribution of the inflection points to a certain characteristic family is done by comparing the magnitudes of the eigenvalues λa(Φi) and λb(Φi), both given by (35): Substituting ϕ=2/(1+ni) one observes that


Since ui, uj are non-negative and (1−ni)/(1+ni) is always negative by the parameter setting (S2), we have that


Therefore, in the inflection points on the axes, the eigenvalue λa corresponds to the first (smaller) eigenvalue λ1 such that the eigenvalue order (36) holds.

The first inflection manifold, which has Φ1 and Φ2 as its endpoints, was numerically calculated and plotted in Figure 2A. Its location will be important for constructing the solution of RP(Φ, Φ).


Figure 2. Phase space of Example 1 with parameters (44). (A) Integral curves in the ϕ1ϕ2−coordinate plane plotted as continuous curves. The arrows point into the direction of increasing eigenvalues. The first family (blue) crosses all lines with constant ϕ and the second family (red) connects the axes. Additionally, the contact manifold (solid dark curve) and the inflection manifold of the first family (dotted-dashed curve) are plotted. (B) Solution of Riemann problems RP[O,(0.2, 0.25)], RP[O,(0.2, 0.35)] using a finite difference method (low resolution).

4.2. Eigenvectors

The Jacobian matrix (31) with components (32) can be written as a rank two modification

J=D+k=12akbkT=D+BAT    (38)

of the diagonal matrix D=diag(v1,v2)2×2, where one introduces column vectors ak, bk ∈ ℝ2, k ∈ {1, 2}, in matrices

A=(1u1+ΦTu(Φ)1u2+ΦTu(Φ))and B=(ϕ1u1(Φ)-ϕ1ϕ2u2(Φ)-ϕ2),    (39)

which both have rank two. The formula for the components of a right eigenvector can be deduced from secular equations and is given as a function of eigenvalues [14],

rij(λ)=1vjλ[bj1k=12ak1bk2vkλbj2(1+k=1N ak1bk1vkλ)],    (40)

where the parameters

aj1=1,aj2=uj(Φ)+ΦTu,bj1=ϕjuj,bj2=-ϕj,    (41)

correspond to the entries of the matrices A and B in (39). Whereas formula (40) holds for general rank two modifications of form (38), which are valid for systems of arbitrary size. By the abbreviation


one can rewrite the eigenvector (40) with parameters (41) compactly as


with the components


Multiplying by (v1 − λ)(v2 − λ), then the right eigenvectors get the form

r(Φ,λi(Φ))=(ϕ1(v2(Φ)-λi(Φ))+ϕ1ϕ2(u2(Φ)-u1(Φ))ϕ2(v1(Φ)-λi(Φ))-ϕ1ϕ2(u2(Φ)-u1(Φ))),    (42)

for the corresponding eigenvalues λi(Φ), i = 1, 2.

Special eigenvectors, in particular those for the values of the Hugoniot locus of the origin, can be obtained from either exploiting the structure of the Jacobian matrix (31) or from using the analytical form of the eigenvectors (42). For instance, the eigenvectors on a set C(Φ) that correspond to the second eigenvalue λ2(Φ)=v1(Φ) have the form

r2(Φ)=(11).    (43)

In the origin one has r1(O)=(0,1)T and r2(O)=(1,0)T.

4.3. Illustration of a benchmark example

For illustration, a benchmark example (Example 1) is considered with parameter setting

v1=1,v2=1/2,  and  n1=4,n2=3,    (44)

that satisfies specifications (S1)–(S3). The integral curves in the ϕ1ϕ2−coordinate plane are shown in Figure 2A. The arrows point into the direction of increasing eigenvalues. The first characteristic family crosses all lines where ϕ is constant (except ϕ = 1) and the second family is connecting the axes. It can be recognized that the second family is genuinely non-linear out of Φx∪Φ, whereas genuine non-linearity of the first family is lost at an inflection manifold, where the corresponding eigenvalues take their minimum.

The direction of increasing eigenvalues of the second characteristic family switches at the contact manifold. This direction switch impacts on the solution structure of the Riemann problem RP(O, Φ+), depending on the position of the right state Φ+ relative to the contact manifold, see Figure 2B for results of simulations by a finite difference method, namely a second-order central-upwind scheme where the numerical flux is a generalization of the Lax-Friedrichs flux, where the constant weight (Δt)/(2Δx) at the smoothing term is replaced by a locally adaptive weight, that take an upper estimate of the local velocities. Here we can observe that only a slight change of the Riemann data provokes a fundamental change of the solution path in the phase space.

5. Shock classification of Hugoniot locus of origin

In this Section the types of shocks that are connected to the origin are classified. The shock classification of Definition 2 distinguishes three different types of admissible shocks, namely 1-Lax, 2-Lax and overcompressive shocks. The locations of the shocks on the Hugoniot H(O) of the origin O = (0, 0)T are identified in Lemma 4, which includes two contact manifolds. The shock classification builds on the eigenvalue analysis of Section 4 by comparing the shock speed and the eigenvalues in each state on the Hugoniot locus [see [38] for further graphical details].

A key feature in the determination of the shock type involves the switching of the order of the relative velocities at the threshold concentration ϕx,

{u1(Φ)>u2(Φ) for ϕ[0, ϕx),u1(Φ)=u2(Φ) for ϕ{ϕx, ϕ},u1(Φ)<u2(Φ) for ϕ(ϕx, ϕ).    (45)

This order switch is a direct consequence of the definition of the relative velocity (7), (8), together with the specifications (S1) and (S2).

A visual guide for the shock classification of shocks occurring between the origin and states located on the edges ∂1 and ∂2 is given in Figure 3, where the shock speeds are compared to the characteristic speeds. The shock classification is essentially obtained by speed comparisons. In the following result, the shock characterization is established for general parameter choices.

In Figure 3B, two important values appear:

(1) The threshold value ϕσ is defined as

ϕσ=1-n1v2/v1.    (46)

The value ϕσ characterizes a state Φσ = (ϕσ, 0)T where σ(O,Φσ)=λ1(O);

(2) The value ϕu is defined as the unique zero root of

(1-n1ϕ)(1-ϕ)n1-n2-v2/v1=0.    (47)

The value ϕu characterizes a state Φu = (ϕu, 0)T where λ1(Φu)=λ2(Φu); a coincidence point over ∂1-edge, a fact that is only proven in Lemma A1, which is of vital importance for the description of quasi-umbilic points.


Figure 3. Eigenvalues and shock speeds of shocks from the origin to (A) the states Φ = (0, ϕ) ∈ ∂2, ϕ ∈ [0, 1], (B) the states Φ = (ϕ, 0) ∈ ∂1, ϕ ∈ [0, 1]. Bold solid curves represent the shock speed σ(O, Φ). The eigenvalues (characteristic speeds) given in (35) are represented as light solid curves for λa(Φ) and dashed curves for λb(Φ). The dotted horizontal lines represent heights corresponding to ϕx, ϕσ, and ϕu, ordered from top to bottom, the dotted vertical lines represent λ1(O) and λ2(O); hence the position of ϕσ.

Theorem 2. The Riemann problem


connecting the origin to a state of its Hugoniot locus is solved by a single shock of speed σ = σ(O, Φ), which can be classified according to Definition 2 as stated below.

Classification for Φ=(ϕ1,0)T1_:

2-Lax:   λ1,2(Φ)σ<λ2(O) and λ1(O)<σfor ϕ1(0,ϕσ),               λ1,2(Φ)σ<λ2(O) and λ1(O)=σfor ϕ1=ϕσ.    OC:                   λ1,2(Φ)<σ<λ1,2(O)for ϕ1(ϕσ,ϕx).1-Lax:   λ1(Φ)<σ<λ1,2(O) and σ=λ2(Φ)for ϕ1{ϕx,ϕ},               λ1(Φ)<σ<λ1,2(O) and σ<λ2(Φ)for ϕ1(ϕx,ϕ).    (48)

The threshold value ϕσ given by (46).

Classification for Φ=(0,ϕ2)T2_:

1-Lax:   λ1(Φ)<σ<λ1,2(O) and σ<λ2(Φ)for ϕ2(0,ϕx),               λ1(Φ)<σ<λ1,2(O) and σ=λ2(Φ)for ϕ2{ϕx,ϕ}.     OC:                 λ1,2(Φ)<σ<λ1,2(O)for ϕ2(ϕx,ϕ).    (49)

Classification for ΦC(Φx)_:

1-Lax:λ1(Φ)<σ<λ1,2(O) and σ=λ2(Φ).    (50)

Classification for ΦC(Φ)_:

1-Lax:λ1(Φ)=σ<λ1,2(O) and σ=λ2(Φ).    (51)

Strict inequalities in the shock type classification (48)–(51) mean that the corresponding shocks are not characteristic shocks. For instance, for Φ ∈ ∂1 there are no 1-Lax shocks that are characteristic at the left datum O.

Proof. The proof is done in two main steps:

1.- The first step is to relate the characteristic speeds λ1, 2(O) and λa, b(Φ) to the shock speed σ(O, Φ) depending on the position of Φ.

2.- The second step consists in determining which shock classification applies according to Definition 2.

Table 1 provides an overview on the speed magnitudes for right states on the axes ∂1 and ∂2.


Table 1. Shock classification on Hugoniot locus of origin for states on the axes ∂1 and ∂2.

λ2(O) compared to σ(O, Φ) on ∂1 and ∂2

On the axes Φ ∈ ∂i\O the shock speed is limited as

σ(O,Φ)=vi(Φ)<vi={v1=λ2(O), if i=1,v2=λ1(O)<λ2(O), if i=2,    (52)

since the shock speed


is monotonically decreasing on both axes ∂1 and ∂2 [Note that (52) excludes 2-Lax shocks on edge ∂2].

λa(Φ) compared to σ(O, Φ) on ∂1 and ∂2

The eigenvalue λa(Φ) as specified in (35) leads to common properties for both axes. For all points Φ ∈ ∂i\O, i ∈ {1, 2}, along the axes one has

λa(Φ)=[1-(1+ni)ϕi]ui(Φ)            <(1-ϕi)ui(Φ)=vi(Φ)=σ(O,Φ),i=1,2,    (53)

due to properties that apply on the axes, namely, the speed convertibility (30) and the shock speed (28).

λb(Φ) compared to σ(O, Φ) on ∂1 and ∂2

The composition of the eigenvalue λb(Φ) depends on the inequality (45) between the relative velocities u1(Φ) and u2(Φ), and induces quite different behaviors on the edges. On the axis ∂1 the eigenvalue λb(Φ) behaves as

{λb(Φ)<σ(O,Φ) for ϕ1(0,ϕx),λb(Φ)=σ(O,Φ) for ϕ1{ϕx,ϕ},λb(Φ)>σ(O,Φ) for ϕ1(ϕx,ϕ),    (54)

because on ∂1, the inequalities (45) between u1(Φ) and u2(Φ) implies that for the interval ϕ1[0,ϕx) the eigenvalue (35) satisfies


and for the other interval, where ϕ1(ϕx,ϕ], the inequality switches sign.

On the axis ∂2 one has

{λb(Φ)>σ(O,Φ) for ϕ2[0,ϕx),λb(Φ)=σ(O,Φ) for ϕ2{ϕx,ϕ},λb(Φ)<σ(O,Φ) for ϕ2(ϕx,ϕ),    (55)

since, from the inequalities (45) for ϕ2[0,ϕx), one obtains


and for ϕ ∈ (ϕx, ϕ] the inequality switches sign.

On the edge ∂2 the value of λb(Φ) determines whether the shock is 1-Lax or over-compressive. The shock type is decided by the relative magnitudes of u1(Φ) contra u2(Φ): The equality u1(Φ) = u2(Φ) only holds for ϕ2=ϕx (and ϕ2=ϕ). On edge ∂1 there is a coincidence of eigenvalues, see also Lemma A1.

λ1(O) compared to σ(O, Φ) on ∂1 and ∂2

Note that, along both axes σ(O, Φ) is a monotonically decreasing function with σ(O, Φ) = 0. By (52) it is assured that σ(O, Φ) < λ1(O) holds for all states on the axis ∂2\O. However, this does not hold for all states on the axis ∂1\O.

Since limε0σ(O,(ε,0))=λ2(O)>λ1(O)>0 on the axis ∂1, there exists a state Φσ = (ϕσ, 0)T such that σ(O,Φσ)=λ1(O). In this state Φσ the equality


holds, from which the characterization of ϕσ by (46) can be deduced. The shock speed σ(O, Φ) relates to λ1(O) in dependence of ϕσ as

{σ(O,Φ)>λ1(O) for ϕ1[0,ϕσ),σ(O,Φ)=λ1(O) for ϕ1=ϕσ,σ(O,Φ)<λ1(O) for ϕ1(ϕσ,ϕ].    (56)

Combining the previously established (52), (53), and (54), with the discrimination (56) implies that, on the axis ∂1, for ϕ < ϕσ there is a 2-Lax shock, whereas for ϕσ < ϕ < ϕx the shock is over-compressive.

According to (53) and (54), there is a clear separation between the shock speeds for ϕ ∈ (ϕx, ϕ) on the axis ∂1:

λ1(Φ)=λa(Φ)<σ(O,Φ)<λb(Φ)=λ2(Φ).    (57)

This separation allows an association of the eigenvalues according to (36). Similarly, according to (53) and (55), there is a clear separation between the shock speeds (57) for ϕ ∈ (0, ϕx) on ∂2, which establishes the order of the eigenvalues as (36). For ϕ ∈ (ϕx, ϕ) there is no clear eigenvalue separation, which however does not affect the fact that the shock is over-compressive.

Now we are able to conclude the shock classification on the axes: Putting the inequalities (52), (53), (55) together gives the shock classification (49) for right states on the axis ∂2. Putting the inequalities (52), (53), (54), (56) together gives the shock classification (48) for right states on the axis ∂1. See the resume in Table 1.

Contact manifold C(Φx)_

The shock classification for states on the contact manifold ΦC(Φx) assures that all states are connected to the origin by a 1-Lax shock: Since the functions v1(Φ) and v2(Φ) take their maximum in the origin O and the origin is not part of the contact manifold, one has the strict inequality

σ(O,Φ)=v1(Φx)<v1=λ1(O)<v2=λ2(O).    (58)

By the eigenvalue characterization (34) for states on C(Φx), and noting that ΔΦ > 0 for all ΦC(Φx), one gets

λ1(Φ)<σ(O,Φ)=v1(Φ)=λ2(Φ)    (59)

for all ΦC(Φx). The properties (58) and (59) are summarized in the shock classification (50).

Contact manifold =C(Φ)_

For states on the maximum packing manifold Φ ∈ ∂ one has vanishing eigenvalues,


which leads to classification (51).

Figure 4 displays the different shock types for states connected with the origin, as elaborated in Theorem 2. The 1-Lax states of H(O) are shown as blue dashed lines. States on the contact manifold C(Φx) satisfy the 1-Lax conditions being characteristic in the second family at the right state. Similarly, the states on the maximum packing C(Φ)= satisfy the conditions to be 1-Lax shocks, being characteristic at the right state with respect to both families.


Figure 4. Different shock types for states connected to the origin. Continuous curves are 2-Lax shocks, dashed curves are 1-Lax shocks and dash-dotted are over-compressive shocks. The Hugoniot locus of the origin comprises the two edges and the contact manifolds C(Φx), C(Φ). Additionally, some 2-Lax shock curves from arbitrarily chosen generic points Φ1M, Φ2M, Φ-M, and Φ+M are indicated. Reference values on the horizontal axis ∂1 are ϕσ, ϕx, and ϕ, on the vertical axis ∂2 are ϕx and ϕ. Recently, we became aware that the 2-Lax branch of H(Φ2M) crossing the OC branch of H(O) corresponds to the fishbone structure that was studied in [36].

On the edge ∂2, at Φx = (0,ϕx)T and Φ = (0,ϕ)T the shocks are 1-Lax and characteristic at the right states, say σ(O,Φx)=λ2(Φx) and σ(O,Φ)=λ2(Φ). On the edge ∂1, at Φσ = (ϕσ, 0)T, and Φx = (ϕx, 0)T and Φ = (ϕ, 0)T the shocks are 2-Lax (red and continuous) and 1-Lax, respectively, being characteristic at left or right states, say σ(O,Φσ)=λ1(O), σ(O,Φx)=λ2(Φx) and σ(O,Φ)=λ2(Φ). The over-compressive shocks (green and dot-dashed) are restricted in between Φσ and Φx on ∂1 and in between Φx and Φ on ∂2. Several states located in Figure 3B have characteristic shocks: If ΦC(Φx), then the shock is right characteristic for the second family. If ϕ = ϕ, then it is characteristic for both families. If Φ = (0,ϕσ)T, then the shock is 2-Lax but left characteristic in the first family.

6. Riemann problems

In this section, we construct the solution of the Riemann problems RP(O, Φ) and RP(Φ, Φ), where O = (0, 0)T, Φ ∈ ∂, and ΦDΦ is a generic right or left state in the phase space. These Riemann problems are derived from the standard initial condition (15).

6.1. The Riemann problem RP (O, Φ)

The behavior of solutions to the Riemann problem RP(O, Φ) depends on the position of Φ in the interior of DΦ with respect to the contact manifold C(Φx), which divides the domain into two regions, see Figure 4:


where C(Φx)=cl(Dx)cl(Dx+) marks the intersection of the closures of those two regions. The main difference between both domains is the opposite direction of the characteristic speeds on the integral curves. In Dx- the second eigenvalue increases from edge ∂1 to edge ∂2, and, reversely, in Dx+ it increases from edge ∂2 to edge ∂1. See the orientation of the second rarefactions in Figure 2A, where the arrows point into the directions of increasing eigenvalues.

A Riemann solution from O to any state Φ in DΦ generally consists of a 1-Lax shock followed by a 2-Lax shock. For Φ belonging to Dx- the middle state that connects the 1-wave with the 2-wave is denoted by Φ-M, and for Φ belonging to Dx+ the middle state is denoted by Φ+M. In both cases the middle state ΦM belongs to the 1-Lax locus of the origin O. Since both waves are shocks any middle state Φ+M or Φ-M is located at H(Φ)H(O). It turns out that Φ+M belongs to the edge ∂1 and Φ-M belongs to the edge ∂2.

For a right state ΦDx- the solution of RP(O, Φ) comprises the 1-Lax shock from O to a state ΦM = (0,ϕM)T such that ϕM ∈ (0, ϕx) and the 2-Lax shock from ΦM to Φ. Similarly, for a right state ΦDx+ the solution of RP(O, Φ) comprises the 1-Lax shock from O to a state ΦM = (ϕM, 0)T such that ϕM ∈ (ϕx, ϕ) and the 2-Lax shock from ΦM to Φ. Such cases are depicted in Figure 4; the former as Φ-Dx- and the latter as Φ+Dx+.

We noticed a change in the nature of solution behavior when Φ crosses the contact manifold C(Φx). Such a bifurcation clearly occurs in the phase space, see also Figure 2B. However, in the solution profile the changes are smooth through this manifold. To do this, we need to understand the solution of the RP(O, Φ) with Φ in C(Φx). In [38], there is a detailed explanation of the L1-continuity of solutions as Φ changes along DΦ.

The solution of the Riemann problem RP(O, Φ) with Φ on the Hugoniot locus of O comprises a single shock from O to Φ with a classification that depends on the position of Φ and is elaborated in Theorem 2:

1. For Φ = (ϕ, 0)T with ϕ ∈ (ϕσ, ϕx) the shock is over-compressive,

2. For Φ = (0, ϕ)T with ϕ ∈ (ϕx, ϕ) the shock is over-compressive,

3. For Φ = (ϕ, 0)T with ϕ ∈ (0, ϕσ) the shock is 2-Lax,

4. For any other Φ in H(O) the shock is 1-Lax.

For a state ΦC(Φx), the solution of the Riemann problem RP(O, Φ) consists of a single 1-Lax shock from O to Φ. The two consecutive shocks from O to (0,ϕx)T and from (0,ϕx)T to Φ have both the same speed v1(ϕx)=v2(ϕx), which in turn is the same speed of the direct shock from O to Φ. Since all these shocks have the same speed, the middle state (0,ϕx)T is “invisible” in the solution profile in the physical space; this is because of Lemma 1. Therefore, the solution structure for ΦDΦ is represented as


For ΦC(Φx) the middle state ΦM may assume any other state on the same contact manifold and even collapse with Φ. This is the basis of the L1-continuity in the solution profiles.

Please refer to Figure 4 and notice that as any Φ-Dx- tends to a ΦC(Φx) the middle state Φ-M tends to (0,ϕx)T. Similarly, as Φ+Dx+ tends to a ΦC(Φx), the state Φ+M tends to (ϕx, 0)T. Thus, we notice a continuous dependence of solutions to the Riemann problem RP(O, Φ) on the right datum, even when it has a vivid bifurcation in the phase space.

6.2. The Riemann problem RP (Φ, Φ)

In this Section, the dependence of the solution structure on the exponents n1 and n2 is illustrated by two examples: Example 1 and Example 2, with the corresponding parameter settings (44) and (60), respectively.

In Example 1, the solution consists of a simple 1-wave comprising a shock followed by a rarefaction or a single rarefaction. The solution structure of Example 2 depends on the existence of one or two detached inflection curves for the first characteristic family, which seem to emerge within the quasi-umbilic points Q2 and Q3 described in Appendix 1 (Supplementary material). As C(Φ) is a contact manifold, the goal is to consider a generic state ΦDΦ from which waves reaching any state in Φ∈∂ are constructed. Remarkably, the solution structure for RP(O, Φ) remains the same for both examples.

Construction for Example 1. The construction for Example 1 can be orientated by the integral curves and the inflection manifold in Figure 2. In Figure 2A, the inflection manifold I1 of the first family is shown as a dash-dotted curve, which is almost a straight line, connecting the states Φ1=(0.4,0)T and Φ2=(0,0.5)T. This can be compared to Lemma 6. Therefore, a wave curve starting at a state on the upper-right hand side of I1 follows a centered rarefaction of the first family until the maximum packing manifold ∂ is reached. The final rarefaction point is (0, 1)T, unless the starting state is on the ∂1 axis, in such a case, the final rarefaction point is (1, 0)T. Again, according to the triple shock rule over a contact manifold C(Φ), the latter may also be seen as finishing at (1, 0)T. The clear bifurcation in the phase space does not exist in the solution profile.

The characteristic velocities near the contact manifold C(Φ) are close to zero, and the characteristic directions are close to r2=r1=(1,-1)T, see (43). If the left state Φ is on the lower-left side of the inflection I1 then the wave is obtained by a backward 1-wave construction. Namely, all shocks from a state Φ at the lower-left hand side of I1 are connected to a state ΦMH(Φ) on the right of I1 satisfying σ(Φ,ΦM)=λ1(ΦM). Therefore, the solution for such a state Φ comprises the 1-Lax shock from Φ to ΦM which is characteristic at ΦM, and the rarefaction curve from ΦM to (0, 1)T, or (1, 0)T if Φ belongs to ∂1.

Construction of Example 2. In this example the parameters are specified as

v1=1,v2=1/2,  n1=4.6,n2=1.5.    (60)

The construction of Example 2 is visualized in Figure 5. For the study of the solution profiles, the growth direction of the slow eigenvalues along the integral curves is key in their structure for both examples RP(Φ, Φ). Figure 5A shows the I1 inflection manifold in two detached parts called T and B; we use the latter for the bottom part and the former for the top part. In T, the inflection consists only of minimum eigenvalues. Point B possesses certain qualities, it is the limit of the integral curves that do not cross B, from the edge ∂1 to B, the eigenvalues are minimum, and from B to the maximum concentration, they are maximum.


Figure 5. Curves of Example 2. In (A), the first family rarefaction curves are plotted as continuous curves, the arrows point into the direction of increasing eigenvalues. The dash-dotted curves are the inflection manifold I1 splitted into the two branches: T and B. In (B), the solid curve C1,1 is the double-contact of the first family (notice the two components, one in the light region and other in the white region); A is the state where a limit shock curve from T does not cross B; it is tangent to it at B. In (C), the presentation of the wave curve solutions includes continuous curves corresponding to rarefaction fans and dashed curves corresponding to shock waves.

In Figure 5B the light (blue) and dark (brown) shaded regions represent right states Φ for which the solutions are similar to that of Example 1:

1. When the states belong to the upper corner above T within the light shaded region, then the wave curve comprises a single 1-rarefaction connecting the state Φ moving along increasing eigenvalues λ1 toward Φ = (1, 0)T;

2. For a left state Φ in the dark shaded region, a characteristic 1-Lax shock to a middle state in the light shaded region crosses T, from this middle state a first family rarefaction follows to the maximum package concentration.

The construction of RP(Φ, Φ) solutions for Φ in the white region of Figure 5B is outlined below. All solutions start with a right characteristic 1-Lax shock to a state ΦMH(Φ) such that σ(Φ,ΦM)=λ1(ΦM). We distinguish three cases according to the number of flights over B or T:

(1) The shock curve crosses B once;

(2) The shock curve crosses T once;

(3) The shock curve crosses B twice, and T once.

In Figure 5C, the left states of case (1) are represented by Φ1 and Φ2, case (2) are represented by left states Φ in the dark region of panel (b) (see point 2 above) and, case (3) are represented by Φ3. The wave curves starting from Φ1, Φ2 and Φ3 have the following structure:

Φ11-Lax|Φ1M1-rarΦ1B|1-Lax|Φ1C1-rar Φ,    (61)
Φ21-Lax|Φ2M1-rarΦ2B|1-Lax|Φ2C1-rarΦ,    (62)
Φ31-Lax|Φ3M1-rarΦ,    (63)

where the symbol | indicates where a shock is characteristic. For case (3) the construction of the shock curves proceeds as before, namely a 1-Lax shock followed by a 1-rarefaction connecting a middle state Φ3M in the light shade region in Figure 5B to Φ. For cases (1) and (2), the characteristic shock (Φi,ΦiM) precedes a 1-rarefaction from ΦiM toward B at a state ΦiB on C1, 1, from there another 1-Lax left and right characteristic shock connects to a state ΦiC on the other side of C1, 1. From ΦiC the wave curve terminates with a 1-rarefaction to Φ.

If the state Φ2 is approximated to a state Φ3 then a continuous change in the solution profile is observed such that the wave group (62) collapses to the wave group (63). Indeed, if Φ2 comes closer to Φ3, then Φ2M comes closer to Φ2B and Φ2C to Φ3M in such a way that the shock speeds σ(Φ2,Φ2M) and σ(Φ2B,Φ2C) approximate σ(Φ3,Φ3M), until, in the limit, the 1-rarefaction wave from Φ2M to Φ2B vanishes.

In Figure 6, the wave group for sedimentation with initial concentration Φ0=(0.1,0.2)T is shown as an example of solutions of the first type, which turns out to be most complete by an additional sequence of shocks upon crossing the double contact C1, 1. The concentration profiles are depicted along the vertical direction, denoted by the dimensionless variable ξ = x/t, where the species experience gravitational effects to reach the base at x = L. In Figure 6A, a single shock for the setting of Example 1 can be observed, similarly as it occurs for Φ3 of Example 2. Figure 6C shows the concentrations for the setting of Example 2. In this case, we can appreciate two 1-Lax shocks. The first shock, with the initial datum on the left and Φ1M on the right, has velocity σ1M. The second shock, denoted with a larger velocity σBC, crosses the double contact C1, 1 from Φ1B to Φ1C.


Figure 6. (A) Typical solution profile for Example 1, with parameters v∞1 = 1, v∞2 = 1/2, n1 = 4, n2 = 3; solid curves represent rarefaction fans, dashed curves correspond to shock waves; blue curves correspond to phase 1, red to phase 2; (B) path of the wave group in phase space plots: the blue curve represents the wave curve associated with Example 1, the red curve represents the path for the solution in Example 2; (C) typical solution profile for Example 2, with parameters v1∞ = 1, v2∞ = 1/2, n1 = 4.6, n2 = 1.5.

Why does the sediment rise from the very bottom with dominance of particles of phase 2 marked by red lines? A physically motivated explanation is that, although phase 2 by design has a smaller terminal velocity (v∞2 < v∞1), the smaller hindrance exponent (n2 < n1) compensates this by an overall larger settling velocity at large concentrations. If particle shape and density can be assumed to be the same, then the smaller terminal velocity, i.e., the Stokes velocity at vanishing concentrations, comes from a smaller particles size. However, at large concentrations approaching the maximum packing ∂, the hindrance factor (8) becomes dominant in determining the settling velocity. Consequently, the smaller particles get attached to the bottom, whereas the bigger particles marked in blue color tend to settle upon them. Toward the sediment, lighter particles float above denser particles. The state Φ = (0, 1)T at the very bottom means that the hindrance factor is the dominant and definitive criterion of the horizontal order of the phases.

Notably, the behavior of the sedimentation process is similar from Φ0M to Φ in Example 1 and from Φ1M to Φ1B in Example 2. The rarefaction from Φ1C to Φ has no counterpart in Example 1. However, this only occurs in the solution profiles, as noted in the central panel where both curves reach Φ with similar slopes. The change in convexity in the behavior of the sediment near the base is due to the different exponents in the hindered-settling velocities, as observed in the settings for each example given in Eqs. (44) and (60), respectively.

7. Discussion

We have successfully analyzed Riemann problem solutions related to particle sedimentation testing. The semi-analytical analysis for the solution of the RP(Φ, Φ) shows that extending it to systems with more components is complicated, and its development requires a handcrafted approach as all the manifolds through which the solutions can bifurcate must be studied. Example 2 is an instance of this situation, say with the introduction of a double contact manifold C1, 1. In several phase space dimensions more complicated situations can occur.

On the other hand, the study of the solution for RP(O, Φ) is much simpler and can be done analytically. This situation can be scaled to systems with several different particles as long as the conditions (S1) and (S2) are consistently scaled, namely v∞1 > v∞2 > ⋯ > vN > 0 and n1 > n2 > ⋯ > nN > 1, respectively. In these situations, the contact manifold has codimension one and is identified by a specific ϕx as in the Defs. 4 and 5 with the respective extension of Lemma 2. In particular, the family associated with the manifold depends on the type of codimension of the particular state in it. The inner points of this surface belong to a N-th contact manifold, its boundary of codimension 2 is a (N−1)-th contact manifold, this reduction of family continues until reaching the edges that represent a 2-th contact manifold. Still, C(Φx) is part of the Hugoniot locus of origin.

The solution is constructed by concatenating Lax shocks, whose position depends solely on the relative location of Φ in relation to C(Φx). When ϕ < ϕx, the 1-Lax shock usually occurs on the corresponding axis ∂N, and for ϕ > ϕx, it typically occurs on the axis ∂1. It is important to note that this may vary if Φ lies on the axes or, more generally, if ϕk = 0 for some k ∈ 1, …, N. Preliminary examples show that for N = 3, the construction can be achieved using a simple algorithm. Starting from Φ, we can calculate its Hugoniot locus and obtain the ΦM′ state: If ϕ < ϕx is satisfied, then ϕ3M is set to 0, and if ϕ > ϕx, then ϕ1M is set to 0; this yields a 3-Lax shock. From ΦM′, we can find ΦM at the corresponding boundary, as discussed above. Thus, the sought-after solution is given by the sequence


There is currently limited knowledge regarding the analysis of umbilic points for systems with more than two equations. The investigation of quasi-umbilic and DRS points (see [39, 40]) presents great potential for advancing our understanding in this research field. The initial identification of DRS points on the elliptic/hyperbolic boundary is provided in Čanić [41]. A comprehensive discussion on the additional connections between these points can be found in the Supplementary material.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.


This work was supported by Conacyt (National Council for Science and Technology of Mexico) through the Basic Science program with Grant A1-S-26012. We would also like to thank Conicyt (Chile) for their support through Fondecyt project #1120587. SB acknowledges the support of Universidad del Bío-Bío through project INES I+D 22-14 (Co-construyendo un ecosistema para la transformación digital y el desarrollo territorial de las Regiones de Biobío y Ñuble) and ANID (Chile) for support through Fondecyt project #1230560.


We would like to extend our sincere appreciation to the editor for their professionalism and valuable input in shaping this work. We are also grateful to the reviewers whose feedback greatly improved the presentation and content of the manuscript. Additionally, we thank ITAM for their installations and partial support from Asociación de Cultura Mexicana A.C.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at:


1. Berres S, Bürger R, Karlsen KH, Tory EM. Strongly degenerate parabolic-hyperbolic systems modeling polydisperse sedimentation with compression. SIAM J Appl Math. (2003) 64:41–80. doi: 10.1137/S0036139902408163

CrossRef Full Text | Google Scholar

2. Bürger R, Karlsen KH, Tory EM, Wendland WL. Model equations and instability regions for the sedimentation of polydisperse suspensions of spheres. Z Angew Math Mech. (2002) 82:699–722. doi: 10.1002/1521-4001(200210)82:10<699::AID-ZAMM699>3.0.CO;2-#

CrossRef Full Text | Google Scholar

3. Benzoni-Gavage S, Colombo RM. An n-populations model for traffic flow. Eur J Appl Math. (2003) 14:587–612. doi: 10.1017/S0956792503005266

CrossRef Full Text | Google Scholar

4. Berres S, Ruiz-Baier R, Schwandt H, Tory E. An adaptive finite-volume method for a model of two-phase pedestrian flow. Netw Het Media. (2011) 6:401–23. doi: 10.3934/nhm.2011.6.401

CrossRef Full Text | Google Scholar

5. Berres S, Bürger R. On Riemann problems and front tracking for a model of sedimentation of polydisperse suspensions. Z Angew Math Mech. (2007) 87:665–91. doi: 10.1002/zamm.200710343

CrossRef Full Text | Google Scholar

6. Berres S, Bürger R. On the settling of a bidisperse suspension with particles having different sizes and densities. Acta Mecanica. (2008) 201:47–62. doi: 10.1007/s00707-008-0072-0

CrossRef Full Text | Google Scholar

7. Biesheuvel PM. Particle segregation during pressure filtration for cast formation. Chem Eng Sci. (2000) 55:2595–606. doi: 10.1016/S0009-2509(99)00536-9

CrossRef Full Text | Google Scholar

8. Masliyah JH. Hindered settling in a multiple-species particle system. Chem Eng Sci. (1979) 34:1166–8.

Google Scholar

9. Basson DK, Berres S, Bürger R. On models of polydisperse sedimentation with particle-size-specific hindered-settling factors. Appl Math Modell. (2009) 33:1815–35. doi: 10.1016/j.apm.2008.03.021

CrossRef Full Text | Google Scholar

10. Richardson JF, Zaki WN. Sedimentation and fluidization: part I. Trans Inst Chem Engrs. (1954) 32:35–53.

Google Scholar

11. Dafermos CM. Hyperbolic Conservation Laws in Continuum Physics. Berlin: Springer Verlag (2000).

Google Scholar

12. Donat R, Mulet P. A secular equation for the Jacobian matrix of certain multi-species kinematic flow models. Num Methods Partial Diff Equat. (2010) 26:159–75. doi: 10.1002/num.20423

CrossRef Full Text | Google Scholar

13. Berres S, Voitovich T. On the spectrum of a rank two modification of a diagonal matrix for linearized fluxes modelling polydisperse sedimentation. In:Tadmor E, , editor. Hyperbolic Problems. Theory, Numerics and Applications. College Park, MD: American Mathematical Society (AMS) (2009). p. 409–18.

Google Scholar

14. Bürger R, Donat R, Mulet P, Vega CA. Hyperbolicity analysis of polydisperse sedimentation models via a secular equation for the flux Jacobian. SIAM J Appl Math. (2010) 70:2186–213. doi: 10.1137/09076163X

CrossRef Full Text | Google Scholar

15. Castañeda P, Furtado F, Marchesin D. On singular points for convex permeability models. In:Ancona F, Bressan A, Marcati P, Marson A, , editors. Hyperbolic Problems: Theory, Numerics, Applications. Springfield, MO: American Mathematical Society (AMS) (2014). p. 415–422.

Google Scholar

16. Isaacson E, Marchesin D, Plohr B, Temple B. The Riemann problem near a hyperbolic singularity: the classification of solutions of quadratic Riemann problems. SIAM J Appl Math. (1988) 48:1009–32.

Google Scholar

17. Isaacson E, Marchesin D, Plohr B, Temple B. Multiphase flow models with singular Riemann problems. Mat Apl Comput. (1992) 11:147–66.

Google Scholar

18. Matos V, Castañeda P, Marchesin D. Classification of the umbilic point in immiscible three-phase flow in porous media. In:Ancona F, Bressan A, Marcati P, Marson A, , editors. Hyperbolic Problems: Theory, Numerics, Applications. Springfield, MO: American Mathematical Society (AMS) (2014). p. 791–9.

Google Scholar

19. Schaeffer DG, Shearer M. The classification of 2 x 2 systems of non-strictly hyperbolic conservation laws with application to oil recovery. Comm Pure and Appl Math. (1987) 40:141–78.

Google Scholar

20. Fayers FJ. Extension of Stone's method 1 and conditions for real characteristic three-phase flow. SPE Reservoir Eng. (1989) 4:437–45.

Google Scholar

21. Isaacson E, Temple B. Analysis of a singular hyperbolic system of conservation laws. J Differ Equ. (1986) 65:250–68.

Google Scholar

22. Abreu E, Matos V, Perez J, Rodriguez-Bermadez P. A class of Lagrangian-Eulerian shock-capturing schemes for first-order hyperbolic problems with forcing terms. J Sci Comput. (2021) 86:14. doi: 10.1007/s10915-020-01392-w

CrossRef Full Text | Google Scholar

23. Castañeda P, Abreu E, Furtado F, Marchesin D. On a universal structure for immiscible three-phase flow in virgin reservoirs. Comput Geosci. (2016) 20:171–85. doi: 10.1007/s10596-016-9556-5

CrossRef Full Text | Google Scholar

24. Azevedo AV, de Souza AJ, Furtado F, Marchesin D, Plohr B. The solution by the wave curve method of three-phase flow in virgin reservoirs. Transp Porous Media. (2010) 83:99–125. doi: 10.1007/s11242-009-9508-9

CrossRef Full Text | Google Scholar

25. Tang J, Castañeda P, Rossen WR, Marchesin D. Three-phase fractional-flow theory of foam-oil displacement in porous media with multiple steady states. Water Resour Res. (2019) 55:10319–39. doi: 10.1029/2019WR025264

CrossRef Full Text | Google Scholar

26. Rodríguez-Bermúdez P, Marchesin D. Riemann solutions for vertical flow of three phases in porous media: simple cases. J Hyperbolic Differ Equ. (2013) 10:335–70. doi: 10.1142/S0219891613500124

CrossRef Full Text | Google Scholar

27. de Souza AJ. Stability of singular fundamental solutions under perturbations for flow in porous media. Mat Aplic Comp. (1992) 11:73–115.

Google Scholar

28. Schecter S, Marchesin D, Plohr B. Structurally stable Riemann solutions. J Differ Equ. (1996) 126:303–54.

Google Scholar

29. Schecter S, Plohr BJ, Marchesin D. Classification of codimension-one Riemann solution. J Dyn Diff Equ. (2001) 13:523–88. doi: 10.1023/A:1016634307145

CrossRef Full Text | Google Scholar

30. Azevedo AV, Marchesin D, Plohr B, Zumbrun K. Capillary instability in models for three-phase flow. Zeitschrift für Angewandte Mathematik und Physik. (2002) 53:713–46. doi: 10.1007/s00033-002-8180-5

CrossRef Full Text | Google Scholar

31. Lozano Guerrero LF, Marchesin D. Diffusive Riemann solutions for 3-phase flow in porous media. In: CNMAC 2019 - XXXIX Congresso Nacional de Matematica Aplicada e Computacional. Uberlândia (2020). p. 7.

Google Scholar

32. Majda A, Pego RL. Stable viscosity matrices for systems of conservation laws. J Diff Equ. (1985) 56:229–62.

Google Scholar

33. Rodríguez-Bermúdez P, Marchesin D. Loss of strict hyperbolicity for vertical three-phase flow in porous media. In:Ancona F, Bressan A, Marcati P, Marson A, , editors. Hyperbolic Problems: Theory, Numerics, Applications. Springfield, MO: American Mathematical Society (AMS) (2014). p. 881–8.

Google Scholar

34. Furtado F. Structural stability of nonlinear waves for conservation laws. PhD thesis. New York University, New York, NY, United States (1989).

Google Scholar

35. Lax P. Hyperbolic systems of conservation laws II. Commun Pure Appl Math. (1957) 10:537–66.

Google Scholar

36. Matos V, Azevedo AV, Mota JCD, Marchesin D. Bifurcation under parameter change of Riemann solutions for nonstrictly hyperbolic systems. Z Angew Math Phys. (2015) 66:1413–52. doi: 10.1007/s00033-014-0469-7

CrossRef Full Text | Google Scholar

37. Temple B. Systems of conservation laws with invariant submanifolds. Trans Am Math Soc. (1983) 280:781–95.

Google Scholar

38. Berres S, Castañeda P. Identification of shock profile solutions for bidisperse suspensions. Bull Braz Math Soc. (2016) 47:105–15. doi: 10.1007/s00574-016-0125-2

CrossRef Full Text | Google Scholar

39. Dumortier F, Roussarie R, Sotomayor J. Generic 3-Parameter Families of Vector Fields: Unfoldings of Saddle, Focus, and Elliptic Singularities with Nilpotent Linear Parts. Berlin; Heidelberg; New York, NY: Springer-Verlag (1991).

Google Scholar

40. Azevedo AV, Marchesin D, Plohr B, Zumbrun K. Bifurcation of nonclassical viscous shock profiles from the constant state. Commun Math Phys. (1999) 202:267–90.

Google Scholar

41. Čanić S. On the influence of viscosity on Riemann solutions. J Dyn Differ Equ. (1998) 10:109–49.

Google Scholar

Keywords: system of non-linear conservation laws, quasi-umbilic point, contact manifold, Hugoniot locus, Riemann problem

Citation: Berres S and Castañeda P (2023) Bifurcation of solutions through a contact manifold in bidisperse models. Front. Appl. Math. Stat. 9:1199011. doi: 10.3389/fams.2023.1199011

Received: 02 April 2023; Accepted: 05 June 2023;
Published: 23 June 2023.

Edited by:

Pablo Aguirre, Federico Santa María Technical University, Chile

Reviewed by:

Vítor Matos, University of Porto, Portugal
Hongjun Cheng, Yunnan University, China

Copyright © 2023 Berres and Castañeda. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Pablo Castañeda,