Nonlinear Helmholtz equations with sign-changing diffusion coefficient

In this paper we study nonlinear Helmholtz equations with sign-changing diffusion coefficients on bounded domains. The existence of an orthonormal basis of eigenfunctions is established making use of weak T-coercivity theory. All eigenvalues are proved to be bifurcation points and the bifurcating branches are investigated both theoretically and numerically. In a one-dimensional model example we obtain the existence of infinitely many bifurcating branches that are mutually disjoint, unbounded, and consist of solutions with a fixed nodal pattern.


Introduction
In this paper, we are interested in nonlinear Helmholtz equations of the form where Ω ⊂ R N is a bounded domain and the diffusion coefficient σ is sign-changing. As we will explain in Section 1.1, such problems occur in the study of time-harmonic wave propagation through metamaterials with negative permeability and nonlinear Kerrtype permittivity. Up to now, the linear theory dealing with the well-posedness of such problems for right-hand sides f (x) instead of κ(x) u 3 has been studied to some extent both analytically and numerically [6,5,9,3,8,4]. Here, the main difficulty is that the differential operator u → − div(σ(x) ∇u) is not elliptic on the whole domain Ω. Accordingly, the standard theory for elliptic boundary value problems based on the Lax-Milgram Lemma does not apply. In the papers [6,5] the (weak) T-coercivity approach was introduced to develop a solution theory for such linear problems. The guiding idea of this method is to require that the strongly indefinite bilinear form satisfies the assumptions of the Lax-Milgram Lemma up to some compact perturbation once v is replaced by Tv for some isomorphism T : H 1 0 (Ω) → H 1 0 (Ω). Our intention is to combine this approach with methods from nonlinear analysis to study the nonlinear Helmholtz equation Eq. (1).
Under reasonable assumptions on Ω, σ, and c our main contributions are the following: (i) There is an orthonormal basis (φ j ) j∈Z of L 2 (Ω) that consists of eigenfunctions of the linear differential operator u → −c(x) −1 div(σ(x) ∇u). The corresponding eigenvalue sequence (λ j ) j∈Z is unbounded from above and from below. Moreover, the eigenfunctions are dense in H 1 0 (Ω). (ii) If κ ∈ L ∞ (Ω) then each of the eigenvalues λ j is a bifurcation point for Eq. (1) with respect to the trivial solution family {(0, λ) : λ ∈ R}. By definition, this means that for all j ∈ Z there is a sequence u n j , λ n j n∈N ⊂ H 1 0 (Ω) \ {0} × R of nontrivial solutions such that (u n j , λ n j ) → (0, λ j ) in H 1 0 (Ω) × R as n → ∞. Our numerical illustrations in Section 3 indicate that these solutions may be located on a smooth and unbounded curve in H 1 0 (Ω) × R going through the point (0, λ j ). We can also prove this for some one-dimensional model problem.
(iii) In a one-dimensional case, for any given λ ∈ R there are infinitely many nontrivial solutions for Eq. (1) provided that κ ∈ L ∞ (Ω) is uniformly positive or uniformly negative. This result is obtained using variational methods instead of bifurcation theory.
A few comments regarding Items (i) to (iii) are in order. As to (i), the existence of an orthonormal basis of eigenfunctions in L 2 (Ω) may be considered as a known fact in view of [8,Section 1]. However, we could not find a reference in the literature that covers our setting, so we briefly review this in Section 4. On the other hand, the Weyl law asymptotics seem to be new. Our proofs rely on the weak T-coercivity approach developed in [6,5]. Under this assumption, the linear theory turns out to be analogous to the linear (Fredholm) Theory for elliptic boundary value problems. The construction of isomorphisms T : H 1 0 (Ω) → H 1 0 (Ω), however, is a research topic on its own and depends on the precise setting, notably on the nature of the interface where the sign of σ jumps, see, e.g. [5,4].
Our main bifurcation theoretical results from (ii) also rely on the functional analytical framework given by the weak T-coercivity approach. The task is to detect nontrivial solutions of Eq. (1) that bifurcate from the trivial solution family. Being given the linear theory and the Implicit Function Theorem, one knows that such bifurcations can only occur at λ = λ j for some j ∈ Z. To prove the occurrence of bifurcations we resort to variational bifurcation theory (see below for references). Here the main difficulty comes from the fact that the associated energy functional is strongly indefinite due to the sign change of σ. Strong indefiniteness means that the quadratic part of Ψ λ is positive definite on an infinite-dimensional subspace of H 1 0 (Ω), and it is negative definite on another infinite-dimensional subspace. As a consequence, standard results in this area going back to Böhme [2, Satz II.1], Marino [18], and Rabinowitz [23,Theorem 11.4] do not apply. Instead, we demonstrate how to apply the more recent variational bifurcation theory developed by Fitzpatrick, Pejsachowicz, Recht, Waterstraat [14,21]. Better results are obtained for eigenvalues with odd geometric multiplicity, which is based on Rabinowitz' Global Bifurcation Theorem [22]. In our onedimensional model case we significantly improve and numerically illustrate our bifurcation results with the aid of bifurcation diagrams (Section 3). The latter provide a qualitative picture of the bifurcation scenario given that each point in such a diagram corresponds to (λ, u ) where (u, λ) solves Eq. (1).
As to (iii), the strong indefiniteness of Ψ λ also makes it harder to prove the existence of critical points. Note that critical points u ∈ H 1 0 (Ω) satisfy Ψ λ (u) = 0, which is equivalent to Eq. (1). In the case of positive diffusion coefficients, the Symmetric Mountain Pass Theorem [1,Theorem 10.18] applies and yields infinitely many nontrivial solutions. In the context of strongly indefinite functionals, an analogous result was established only recently by Szulkin and Weth [24,Section 5]. We will show how to apply their abstract results under reasonable extra assumptions in order to obtain infinitely many nontrivial solutions of Eq. (1) for any given λ ∈ R assuming κ ≥ α > 0 or κ ≤ −α < 0.
Before commenting on the physical background of Eq. (1) we wish to emphasize that our main goal is to bring (weak) T-coercivity theory and nonlinear analysis together. Accordingly, we do not aim for the most general assumptions for our results to hold true. For instance, we avoid technicalities related to the regularity of the interface where the sign of σ jumps. Similarly, we content ourselves with the special nonlinearity κ(x) u 3 . Only little effort is needed to generalize our bifurcation results as well as our variational existence results to more general nonlinearities in all space dimensions N ≥ 1.
1.1. Physical motivation. We comment on the physical background of Eq. (1). The propagation of electromagnetic waves with a fixed temporal frequency parameter ω ∈ R is governed by the time-harmonic Maxwell's equations Here, charges and currents are assumed to be absent. The symbols E, D : R 3 → C 3 denote the electric field and the electric induction and H, B : R 3 → C 3 represent the magnetic field and the magnetic induction, respectively. In nonlinear Kerr media the constitutive relations between these fields are given by where ε, µ, χ are real-valued, see [7,Chapter 4]. In physics, these quantities are called permittivity, permeability and third-order susceptibility of the given medium, respectively. Plugging in this ansatz into Eq. (3) one finds Now we assume that the propagation of electromagnetic waves is considered in a closed waveguide Ω × R having a bounded cross-section Ω ⊂ R 2 and all material parameters only depend on the cross-section variable. With an abuse of notation, this means µ(x) = µ(y), ε(x) = ε(y), and χ(x) = χ(y) where x = (y, z) ∈ Ω × R. This is a natural assumption for the modelling of layered cylindrical waveguides. If then the electric field is of the special form E(x) = (0, 0, u(y)) with u real-valued, we infer − div(µ(y) −1 ∇u) − ω 2 ε(y) u = ω 2 χ(y) u 3 in Ω, u ∈ H 1 0 (Ω), which corresponds to Eq. (1) in the two-dimensional case. For complex-valued u the nonlinearity is given by |u| 2 u. Sign-changing diffusion coefficients σ(y) := µ(y) −1 may occur if one of the layers of the waveguide is filled with a negative-index metamaterial (NIM) where the permeability µ may be negative, see, e.g., [20]. Note that any such solution u determines E, D, B, H via Eqs. (4) and (5). We mention that one may equally solve for the magnetic field, which leads to Helmholtz-type problems with nonlinear (i.e. solution-dependent) and possibly sign-changing diffusion coefficients.
1.2. Notation. In the following, we equip the Hilbert spaces H 1 0 (Ω) and L 2 (Ω) with the inner products respectively. Our assumptions on Ω, c, and σ will imply that the associated norms · |σ| and · c are equivalent to the standard norms on these spaces. Moreover, we introduce the bilinear form 1.3. Outline. The paper is organized as follows. Section 2 contains a mathematically rigorous statement of our main results dealing with bifurcations for Eq. (1) from the trivial solution family. These results are illustrated numerically in Section 3 with the aid of bifurcation diagrams. Those illustrate the evolution of solutions along the branches as well the global behavior of the latter. In Section 4, we set up the linear theory that we need to prove our main bifurcation theoretical results in Section 5. Finally, Section 6 contains further existence results for nontrivial solutions of Eq. (1) obtained by a variational approach. The proof is based on the Critical Point Theory from [24,Chapter 4].

Main results
We now come to the precise formulation of our main results for Eq.
Here, Ω and the coefficient functions σ, c, κ will be chosen as follows: Assumption (B). There is a bounded linear invertible operator T : Later, in Corollary 4.5, we show that Assumptions (A) and (B) ensure the existence of an orthonormal basis (φ j ) j∈Z of (L 2 (Ω), (·, ·) c ) consisting of eigenfunctions associated to the linear differential operator u → −c(x) −1 div(σ(x) ∇u) appearing in Eq. (1). Due to the sign-change of σ the corresponding sequence of eigenvalues (λ j ) j∈Z can be indexed in such a way that λ j → ±∞ holds as j → ±∞. In Theorem 2.1 below we show that nontrivial solutions of Eq. (1) bifurcate from the trivial branch {(0, λ) : λ ∈ R} at any of these eigenvalues. If the eigenvalue comes with an odd-dimensional eigenspace, we even find that the bifurcating nontrivial solutions lie on connected sets C j ⊂ H 1 0 (Ω) × R that are unbounded or return to the trivial solution branch {(0, λ) : λ ∈ R} at some other bifurcation point. As in [22], for any given j ∈ Z, the set C j is defined as the connected component of (0, λ j ) in S, which in turn is defined as the closure of all nontrivial solutions of Eq. (1) in the space H 1 0 (Ω) × R. Our first main result reads as follows: Theorem 2.1. Assume (A) and (B). Let (λ j ) j∈Z denote the unbounded sequence of eigenvalues from Corollary 4.5. Then each (0, λ j ) is a bifurcation point for Eq. (1). If λ j has odd geometric multiplicity, then the connected component C j in S containing (0, λ j ) satisfies Rabinowitz' alternative: (I) C j is unbounded in H 1 0 (Ω) × R or (II) C j contains another trivial solution (0, λ k ) with k = j.
(a) In Assumption (A) we may as well assume c(x) ≤ −α < 0; it suffices to replace (c, λ) by (−c, −λ). On the other hand we cannot assume c to be sign-changing since we will need that (·, ·) c is an inner product on L 2 (Ω). In Remark 4.6 (c), we show that one cannot expect our results to hold for general sign-changing c ∈ L ∞ (Ω).
We strengthen our result in some 1D model example where we can show the following: • Assumptions (A) and (B) are satisfied.
• The eigenvalues (λ j ) are simple and in particular have odd geometric multiplicity.
• The case (II) in Rabinowitz' Alternative is ruled out, hence all C j are unbounded.
The seemingly complicated ordering of the eigenvalues is exclusively motivated by the nodal patterns given by Items (i) to (iii). Here, |u | > 0 on Ω ± means that the continuous extension of |u | : Ω ± → R to Ω ± is positive. We stress that nontrivial solutions u are smooth away from the interface x = 0 and continuous at x = 0, but they are not continuously differentiable at this point. In fact, σu is continuous on Ω so that u (0) does not exist in the classical sense. In the following Section 3 our results are illustrated with the aid of bifurcation diagrams.

Visualization of bifurcation results via PDE2path
In this section, we illustrate our theoretical results of Theorem 2.1 and Corollary 2.3 with numerical bifurcation diagrams. These diagrams show the value of λ on the x-axis and the L 2 -norm of solutions u for that λ on the y-axis. Thereby, (numerical) bifurcation diagrams allow to get an overview of the "structure" of solutions and, in particular, to visualize the connected components C j of Theorem 2.1. The results were obtained with the package pde2path [25,11], version 2.9b and using Matlab 2018b. The code to reproduce the numerical results is available on Zenodo with DOI 10.5281/zenodo.5707422.
3.1. One-dimensional example. We consider Ω = (−5, 5) with Ω − = (−5, 0), Ω + = (0, 5) and c ≡ 1. The diffusion coefficient σ is chosen piecewise constant, set σ + = 1 and compare two different values for σ − , namely σ − ∈ {−2, −1.005}. We consider Eq. (6) in this special case, i.e., We choose a tailored finite element mesh which is refined close to Γ = {0} in the following way. We start with an equidistant mesh with h = 2 −9 , i.e., Ω is divided into 5120 equal subintervals. Then, we refine all intervals which are closer than 0.1 to Γ five times by halving them. This finally means that intervals close to Γ are only 2 −14 long. We point out that this finely resolved mesh is required to faithfully represent the interface behavior at Γ = {0}, especially for σ − = −1.005. An insufficient mesh resolution does not only influence the numerical quality of the eigenfunctions or solutions along the branches, but also the (qualitative picture) of the bifurcation diagram. We validated our results by assuring that a further refinement of the mesh (halving all intervals) leads to the same results and conclusions.
Qualitatively, they are quite similar with clearly separated, apparently unbounded branches without secondary bifurcations. Note that the bending direction of the branches to the left is determined by the sign of the nonlinear term and can be predicted by the bifurcation formulae (I.6.11) in [15]. The first striking phenomenon due to the sign-changing coefficient is the occurrence of eigenvalues and, hence, bifurcation points, with negative value. In fact, for sign-changing σ, there are two families of eigenvalues diverging to ±∞, see Theorem 2.1. We use the following labeling of branches (cf. Fig. 1): The branch starting closest to zero is labeled as C 0 and the branches for negative and positive bifurcation points are labeled as C −i and C i with i ∈ N, respectively. The absolute value of i increases as |λ| → ∞. In our setting this labeling of the branches is consistent with the notation introduced in Section 4.
Besides the eigenvalues, we also study the eigenfunctions by considering the solutions at the first point of each branch in Fig. 2. We display the branch name according to Fig. 1 as well as the value of λ at the bifurcation point.
As (partly) expected from [8], we make the following observations. Firstly, the solutions are concentrated (w.r.t. the L 2 -norm) on the "oscillatory part", which is Ω − for negative eigenvalues (left column of Fig. 2) and Ω + for positive eigenvalues (right column of Fig. 2). The eigenvalue closest to zero (from which C 0 emanates) plays a special role (middle column of Fig. 2). Secondly, with increasing |λ|, the number of maxima and minima increases as one observes also for the eigenfunctions of the Laplacian. Thirdly, the transmission condition at Γ requires the (normal) derivative of u to change sign, such that the solutions have a "tip" at the interface. Taking a closer look at the bifurcation values and the corresponding solutions in Fig. 2, we note that C 0 starts much closer to zero for σ − = −1.005 than for σ − = −2. This illustrates the theoretical expectation that due to the symmetry of the domain Ω, we have an eigenvalue approaching zero for the contrast going to −1. Moreover, we observe a certain shrinking of the negative bifurcation values towards zero when the contrast approaches −1.

3.1.2.
Patterns of solutions along branches. We now take a closer look at how solutions evolve along branches -depending on whether the corresponding bifurcation value is negative, close to zero or positive. According to the previous discussion, we focus on σ − = −1.005 in the following because it shows the phenomena in a particularly pronounced form and is close to the interesting "critical" contrast of −1. In general, we observe that a certain limit pattern or profile of the solution evolves on each branch which remains qualitatively stable (values of maxima, minima and plateaus of course change with λ). As example for a negatively indexed bifurcation branch away from zero, we consider C −2 , cf. Fig. 1. The first, 50th, and 100th solution on the branch are depicted in Fig. 3. As described above, the solution concentrates in Ω − where it oscillates, while it decays exponentially in Ω + . This profile remains stable over the branch, but we note that the maxima and minima become wider along the branch. This widening of the extrema in Ω − is also noted for the other branches emanating from a negative bifurcation point. Yet, the more oscillations occur for the branches as λ → −∞, the less pronounced the effect becomes because we have more extrema over the same interval. We emphasize that this effect of widening extrema is specific to the sign-changing case and especially to bifurcations starting at negative λ.  As an example for a positively indexed bifurcation branch away from zero, we study the branch C 5 , cf. Fig. 1. As expected, we observe in Fig. 4 that the first solution concentrates on Ω + , where it oscillates as typical for an eigenfunction of the Laplacian, and shows an exponential decay in Ω − . The oscillatory pattern in Ω + is preserved along the branch. The behavior in Ω − changes when λ gets negative: Instead of an exponential decay to zero, we now see an exponential decay to (almost) a plateau (with value ± √ −λ) and a sharp transition to the zero boundary value. Once this pattern is established, it remains stable as well. This appearance of a plateau different from zero is also a specific phenomenon of the sign-changing case. The occurrence of a plateau in Ω − is also observed in Fig. 5 for the branch C 0 closest to zero, cf. Fig. 1. While the first solution has a similar shape in Ω + and Ω − with a linear decay in each subdomain, the ensuing solutions on the branch quickly evolve a plateau in Ω − and an exponential decay in Ω + . This pattern then remains stable along the branch. All in all, we observe a certain stability of profiles along branches. The form of the profiles depends on where the bifurcation starts. Moreover, we always recognize a concentration to the oscillatory part and further the establishment of plateaus different from zero in Ω − . As already emphasized, both effects are specific to the sign-changing case. This qualitative description of solutions seems to transfer to other contrasts, but the bifurcation points closest to zero and the (quantitative) decay in Ω ± significantly depend on the contrast as already discussed above. Note that our mesh satisfies the symmetry conditions laid out in [4], which may be challenging in more complicated geometries though. We focus on the behavior of solutions in this numerical experiment and let λ vary in [− 12,15]. There are three different types of eigenfunctions either concentrated on Ω − , on Γ, or on Ω + . In contrast to the one-dimensional case, there are several different eigenfunctions concentrated on Γ. As before, the eigenfunctions concentrated on Ω − or Γ are associated with negative values of λ. In Figs. 6 to 8, we show the evolution of solutions along a branch for each of the three types described above.  Furthermore, plateaus in Ω − evolve for negative λ in Figs. 7 and 8. Due to the second space dimension in the problem, we can have two (or more) different plateaus evolving in Ω − . In Fig. 7 for a branch with concentration on Γ, we note that the plateaus and the transition between them seems to slightly change the oscillatory pattern on Γ as well. While the two maxima have almost the same height for the first and 25th point (  Finally, for Fig. 8 and a branch with concentration on Ω + , we emphasize that the solution in Ω + evolves like a solution of the standard Laplacian along a branch. In particular, the extrema become thinner, i.e., more spatially localized, which should be contrasted with solutions concentrated in Ω − in Fig. 6.

Linear Theory
In this section we want to describe the linear theory for weakly T-coercive problems. As pointed out earlier, this theory is essentially well-known [6,5,8]. Since it is short and rather self-contained, we provide the details here, which will moreover allow us to fix the required notation. Furthermore, we prove some Weyl law asymptotics that have not appeared in the literature yet. We want to deal with linear problems of the form The a priori unknown solution u is to be found in the Sobolev space H 1 0 (Ω) and the coefficient functions σ and c are assumed to satisfy the conditions (A) and (B) from Section 1. To develop a solution theory for the variational problem Eq. (7) both in H 1 0 (Ω) and L 2 (Ω) we assume F ∈ H −1 (Ω) = H 1 0 (Ω) . We may rewrite Eq. (7) as We introduce the bounded linear operator A : H 1 0 (Ω) → H 1 0 (Ω) and the compact linear operator C : This is possible by Riesz' Representation Theorem and σ, c ∈ L ∞ (Ω), see Assumption (A). We will also use that the operator C can be written as C =Cι whereC : L 2 (Ω) → H 1 0 (Ω) is a bounded linear operator and ι : H 1 0 (Ω) → L 2 (Ω) denotes the embedding operator which is compact by the Rellich-Kondrachov Theorem.
To prove the invertibility of A define the family of operators z → A z := A + zC for z ∈ C on the complex Hilbert space H 1 0 (Ω; C). The bilinear form associated with A z is given by Here, the first summand is invertible while the other two summands are compact. Therefore, So, we have ker(A z ) = {0}, which implies that A z has a bounded inverse as an injective Fredholm operator. Using the analytic Fredholm theorem on C, see [12, Theorem C.8], the set {A −1 z : z ∈ C} is a meromorphic family of operators with poles of finite rank. Therefore, the operator (A + zC) −1 exists for all z ∈ C \ Λ for a discrete set Λ ⊂ R. In particular, there exists ∈ R such that A is an invertible Fredholm operator.
Regarding Eq. (7) as an equation in H 1 0 (Ω), we thus obtain the following: Proposition 4.2. Let Assumptions (A) and (B) hold as well as F ∈ H −1 (Ω). Then Eq. (7) is equivalent to where To prove the existence of an orthonormal basis of eigenfunctions for Eq. (7) we now turn towards an alternative formulation in L 2 (Ω). From Proposition 4.1 and Proposition 4.2 we obtain that Eq. (7) is equivalent to where K c := ιA −1 C and KF := ιA −1 F.
The compact operator K c : L 2 (Ω) → L 2 (Ω) is self-adjoint with respect to (·, ·) c because of for all u, v ∈ L 2 (Ω). We have thus proved the following.
We thus conclude that there is D > 0 such that This finishes the proof of Eq. (12).
To facilitate the application of this result we add a corollary.
(a) In Corollary 4.5, the ordering of the eigenvalues (λ j ) j∈Z is fixed up to translations of the indices and permutations within eigenspaces. The former ambiguity can be removed by specifying λ 0 . A natural way to do this is to require that λ 0 has the smallest absolute value. As mentioned earlier, we do not choose such an ordering in our 1D model example from Corollary 2.3 because it is in general not consistent with the j-dependent nodal patterns.
(b) A reasonable min-max formula for the eigenvalues (λ j ) j∈Z in terms of the bilinear form a does not seem to exist despite the simple formula a(φ, φ) = j∈Z λ j c 2 j for φ = j∈Z c j φ j . In fact, the bilinear form (u, v) → a(u, v) has a totally isotropic subspace of infinite dimension, for example Section 1], the authors provide some explicit one-dimensional example showing that all statements in this section may be false when c ∈ L ∞ (Ω) is signchanging. In fact, they showed that for some tailor-made σ as in Assumption (A) and c := σ the operator u → −c(x) −1 div(σ(x) ∇u) may have the whole complex plane as spectrum. In particular, the spectral theory of (compact) self-adjoint operators does not apply in this context. (d) We mention some similarities and differences concerning the spectral properties of the differential operator u → −c(x) −1 div(σ(x)∇u) for (I) sign-changing σ and c = 1, (II) σ = 1 and sign-changing c. In the case (I) Proposition 4.4 and Corollary 4.5 show that the sequence of eigenvalues is unbounded from above and from below. This is also true for (II), see the Propositions 1.10 and 1.11 in [10]. On the other hand, there are subtle differences. As we will see in Lemma 5.5, in our one-dimensional model example for case (I) there is precisely one positive eigenfunction φ 0 with associated eigenvalue λ 0 , and that one might not have the smallest absolute value among all eigenvalues. In fact, |λ 0 | can be much larger than |λ −1 |, see Remark 5.6. In particular, there is little hope to prove the existence of positive eigenvalues via some straightforward application of the Krein-Rutman theorem. This is different for the case (II) where Manes-Micheletti [17] (see also [10, Theorem 1.13]) proved the existence of one positive and one negative principal eigenvalue, i.e., algebraically simple eigenvalues of the smallest absolute value among the positive and negative eigenvalues, respectively, coming with positive eigenfunctions. So here the two models exhibit different phenomena.

Proof of Theorem 2.1 and Corollary 2.3
We now prove the theoretical bifurcation results with the aid of known bifurcation results for equations of the form F (u, λ) = 0 where F ∈ C 2 (H × R, H) for some Hilbert space H. We will consider bifurcation from the trivial solution branch, so F is supposed to satisfy F (0, λ) = 0 for all λ ∈ R. To this end we proceed as follows: First, we present two abstract bifurcation theorems that allow to detect local respectively global bifurcation from the trivial solution. Next, we show how to apply these results to prove Theorem 2.1, which is straightforward. Finally, we sharpen our results in the one-dimensional case of Eq. (6) by proving Corollary 2.3.
Proof. Our assumptions (i), (ii), (iii) imply that (Ψ λ ) λ∈R is a continuous family of C 2functionals in the sense of [21, p.537]. If λ ∈ R is as required, then Theorem 2.1(i) in [21] proves that the interval [λ − ε, λ + ε] contains a bifurcation point provided that the Hessians Ψ λ ±ε (0) are invertible and the spectral flow of this family over the interval I := [λ − ε, λ + ε] is non-zero. In fact, since L is invertible and K is compact, the linear operator Ψ λ (0) = L − λK has a nontrivial kernel only for λ belonging to a discrete subset of R. So we may choose ε > 0 so small that Ψ λ +ε (0), Ψ λ −ε (0) are invertible and λ is the only candidate for bifurcation in I by the Implicit Function Theorem. Using then the positivity of K we get from Remark (3) in [21] that the spectral flow over I is the dimension of ker(Ψ λ (0)), which is positive by assumption. So λ is a bifurcation point.
The more classical variational bifurcation theorems by Marino [18], Böhme [2, Satz II.1] and Rabinowitz [23,Theorem 11.4] apply if there is µ ∈ R such that the self-adjoint operator L + µK generates a norm. This assumption holds in the context of nonlinear elliptic boundary value problems involving divergence-form operators with diffusion coefficients σ having a fixed sign. In our setting with sign-changing σ, this is not the case. [22] states that the solutions bifurcating from an eigenvalue of odd algebraic multiplicity lie on solution continua that are unbounded or return to the trivial solution branch {(0, λ) : λ ∈ R} at some other bifurcation point. Here, a solution continuum is a closed and connected set consisting of solutions. Given that the proof of this bifurcation theorem uses Leray-Schauder degree theory, more restrictive compactness properties are required to be compared to Theorem 5.1. On the other hand, no variational structure is assumed. In order to avoid technicalities, we state a simplified variant of this result from Theorem II.3.3 in [15]. The set S ⊂ H × R denotes the closure of nontrivial solutions of F (u, λ) = 0 in H × R. In particular, the statement (0, λ ) ∈ S is equivalent to saying that (0, λ ) is a bifurcation point for F (u, λ) = 0, i.e., there are solutions (u n , λ n ) ∈ H \{0}×R such that (u n , λ n ) → (0, λ ) in H × R as n → ∞.

Global Bifurcation. Rabinowitz' Global Bifurcation Theorem
Theorem 5.2 (Rabinowitz). Suppose H is a separable real Hilbert space and that F ∈ Suppose that λ ∈ R is such that the dimension of ker(L − λ K) is odd. Then (0, λ ) ∈ S. Moreover, if C denotes the connected component of (0, λ ) in S, then (I) C is unbounded or (II) C contains a point (0, λ ) with λ = λ .
A more general version of this result holds in Banach spaces and does not involve any self-adjointness assumption. It then claims the above-mentioned properties of C assuming that λ is an eigenvalue of odd algebraic multiplicity. Under our more restrictive assumptions including self-adjointness the algebraic multiplicity of λ is equal to its geometric multiplicity and hence to the dimension of the corresponding eigenspace.
for λ j as in Corollary 4.5. So Theorem 5.1 implies that each λ j is a bifurcation point for Eq. (1). Finally, Theorem 5.2 shows that every such eigenvalue with odd-dimensional eigenspace satisfies Rabinowitz' Alternative (I) or (II) from above, which finishes the proof of Theorem 2.1.

5.3.
Proof of Corollary 2.3. We now sharpen our results from Theorem 2.1 for the one-dimensional boundary value problem (6). The assumptions on σ, c, κ and Ω = (a − , a + ) ⊂ R were specified in the Introduction. We want to verify that Assumptions (A) and (B) are satisfied in this context. While Assumption (A) is trivial, the verification Assumption (B) dealing with the weak T-coercivity of (u, v) → a(u, v) := Ω σ(x) u v dx requires some work. The following result seems to be well-known to experts, but a reference appears to be missing in the literature. Proof. Let χ ∈ C ∞ 0 (Ω) with χ(x) = 1 for x close to 0. Then define where m ∈ 0, a + |a − | will be chosen sufficiently small. Then T is a well-defined bijective operator on H 1 0 (Ω) because of T • T = id. Moreover, is a compact operator, see Eq. (8). Hence, a is weakly T-coercive.
Remark 5.4. In higher-dimensional settings the verification of the weak T-coercivity condition is much more sensitive with respect to the data. This concerns Ω + , Ω − , the geometric properties of the interface Γ = Ω + ∩ Ω − and the coefficients σ + , σ − . In particular weak Tcoercivity may break down for critical ranges of the contrast σ + (x)/σ − (x), see for instance Theorem 5.1 and Theorem 5.4 in [4].
We thus conclude that the assumptions of Theorem 2.1 are verified and hence the existence of infinitely many bifurcating branches C j is ensured.
To finish the proof of Corollary 2.3 it now remains to prove the nodal characterization of all nontrivial solutions (u, λ) ∈ C j emanating from λ j . Choosing an appropriate numbering of the eigenvalue sequence (λ j ) we need to prove that nontrivial solutions (u, λ) ∈ C j have the stated nodal pattern. To prove this we first compute and analyze the eigenpairs of the linear problem.
This sequence can be ordered in the following way: and φ j has |j| interior zeros in Ω − with φ j > 0 on Ω + .
(i) Positive eigenvalues. Solving the ODE and exploiting the continuity of eigenfunctions as well as the homogeneous Dirichlet boundary conditions, we find the following formula for eigenfunctions φ j associated with positive eigenvalues (λ j > 0) The parameter α j ∈ R \ {0} is chosen such that φ j c = 1. The equation for λ j now results from the condition that σφ j has to be continuous. This means By elementary monotonicity considerations one finds that this equation has a unique solution such that τ j k + a + ∈ jπ, (j + 1 2 )π for j ≥ 1. Moreover, it has a unique solution such that τ 0 k + a + ∈ 0, 1 2 π if and only if σ + a − σ − a + > 1. No further solutions exist. We thus obtain: • For j ≥ 1, there is a unique solution λ j in the interval j 2 π 2 and φ j has j interior zeros in Ω + with φ j > 0 on Ω − .
(ii) Negative eigenvalues. Similarly, we obtain for the negative eigenvalues (λ j < 0) and As for the positive eigenvalues one finds: and φ j has |j| interior zeros in Ω − with φ j > 0 on Ω + .
By the subsequence-of-subsequence argument, it suffices to prove that a subsequence of (ũ n ) has this property. Since (ũ n ) is bounded, there is a weakly convergent subsequence with limit φ ∈ H 1 0 (Ω). Since u n solves Eq. (6) we have by Proposition 4.2 u n = A −1 (λ n + )Cũ n + u n 2 |σ| A −1 G(ũ n ) where A −1 is bounded and C, G are compact, see the proof of Theorem 2.1. Soũ n φ and u n 2 |σ| → 0 impliesũ n → φ in H 1 0 (Ω) where φ = A −1 (λ j + )Cφ. In other words, φ is an eigenfunction of φ → −c(x) −1 div(σ(x)∇φ) associated with the eigenvalue λ j . Since the eigenspaces are one-dimensional, φ is a multiple of φ j . Moreover, from ũ n |σ| = 1 we infer φ |σ| = 1, so ũ n , φ j |σ| ≥ 0 (by choice of γ n ) implies φ = +φ j / φ j |σ| . We thus conclude thatũ n → φ j / φ j |σ| in H 1 0 (Ω) and hence uniformly on Ω as n → ∞. Integrating Eq. (6) once, one finds that the convergence even holds in C 1 (Ω + ) and C 1 (Ω − ). So if infinitely many u n had more than j interior zeros in Ω ± , then the collapse of zeros would cause at least one double zero of φ j , but this is false in view of our formulas for these eigenfunctions from above. So almost all u n have at most j interior zeros in Ω ± . Similarly, almost all u n have at least j zeros. So we conclude that the solutions close to the bifurcation point have exactly j interior zeros in Ω ± and are strictly monotone in Ω ∓ .
Step 3: Nodal characterization along the whole branch. We finally claim that this nodal property is preserved on connected subsets of S that do not contain the trivial solution. Indeed, the set of solutions on C j \ {(0, λ j )} with this property is open in S with respect to the topology of H 1 0 (Ω) × R. It is also closed in S since double zeros cannot occur (by the same arguments as above) and zeros cannot converge to the interface at x = 0 as the solutions evolve along the branch. Indeed, in the latter case the equation on the monotone part would imply that the solution has to vanish identically there, whence u ≡ 0 on Ω, which is impossible. So we conclude that all elements on C j \ {(0, λ j )} have the claimed property and the proof is finished.

Variational methods
We want to show that variational methods can be used to prove further existence and multiplicity results for Eq. (1) under stronger assumptions than before. Our aim is to prove the existence of infinitely many nontrivial solutions of for any given λ ∈ R. This means that any vertical line in a bifurcation diagram for Eq. (1), say Fig. 1, hits infinitely many nontrivial solutions of Eq. (23). To show this we apply the generalized Nehari manifold approach from [24,Chapter 4]. We start by considering the general case where Ω ⊂ R N is an arbitrary bound domain and N ∈ {1, 2, 3}. In this setting we will need more information about the orthonormal basis of eigenfunctions from Corollary 4.5, see Assumption (C) below. Since the latter can be checked in our one-dimensional model example, the general result applies and leads to infinitely many solutions for any given λ ∈ R for Eq. (6), see Corollary 6.5.
6.1. The general case. The variational approach aims at proving the existence of critical points of energy functionals. In our case such a functional is given by and λ ∈ R is fixed, see Eq. (2). Of course, Ψ λ is a well-defined smooth functional on H 1 0 (Ω), but the more natural setting for our analysis involves another Hilbert space H that will be smaller or equal to H 1 0 (Ω) under our assumptions. To define it, denote by M the subspace of H 1 0 (Ω) consisting of all finite linear combinations of the eigenfunctions (φ j ) j∈Z from Corollary 4.5. Those exist by Assumptions (A) and (B). We recall Then we define H to be the completion of M with respect to the norm · := ·, · 1/2 generated by the positive definite bilinear form We will need H ⊂ H 1 0 (Ω) in order to benefit both from Sobolev's Embedding Theorem and the Rellich-Kondrachov Theorem in our analysis. So we have to ensure that the norm · dominates the norm on H 1 0 (Ω). For this reason we add the following hypothesis. Assumption (C). There is a D > 0 such that for all sequences (c j ) j∈Z with only finitely many non-zero entries we have where (φ j ) j∈Z denotes the orthonormal basis from Corollary 4.5.
We emphasize that we do not require the opposite bound · · |σ| , which appears to be much harder to verify. Note that this would imply the opposite inclusion H ⊃ H 1 0 (Ω). We now implement the variational approach in the Hilbert space (H, ·, · ) and prove the existence of critical points of Ψ λ in this space. Given that H is dense in H 1 0 (Ω), Ψ λ (u) = 0 and u ∈ H in fact implies that u ∈ H 1 0 (Ω) is a weak solution in the H 1 0 (Ω)-sense. In other words, critical points of Ψ λ : H → R provide solutions of Eq. (23).
We first provide the functional analytical framework required by the Critical Point Theory from [24]. Being given the inner product from Eq. (25), we have an orthogonal The subspaces E + , E − are infinite-dimensional whereas E 0 is finite-dimensional, which is a consequence of |λ j | → ∞ as |j| → ∞, see Eq. (24). Here, E 0 = {0} is admissible. Let Π ± : H → E ± denote the corresponding orthogonal projectors, and we will write u ± := Π ± u in the following. The whole point about the inner product Eq. (25) is that a(u, u) = u + 2 − u − 2 , so we may rewrite the functional Ψ λ from (2) as Now, Ψ λ has the right form to apply the critical point theorem from [24,Theorem 35]. Applying this result in our setting, we get the following. Note that we can also treat κ(x) ≤ −α < 0 by considering (−σ, −λ). In that case, Eq. (23) has infinitely many solutions in H among which a maximal energy solution u * characterized by Ψ λ (u * ) = max {Ψ λ (u) : u ∈ H \ {0} solves (23)}.
Proof. We check the assumptions (i)-(v) from Theorem 6.2. First note that the functional I is continuously differentiable with Fréchet derivative at u ∈ H given by This is a consequence of H ⊂ H 1 0 (Ω) ⊂ L 4 (Ω) with continuous injections by Proposition 6.1 and Sobolev's Embedding Theorem for N ∈ {1, 2, 3}. From Eq. (26) and I (u)[u] = 4I(u) > 0 for u = 0 we infer the first part of (i). The weak lower semicontinuity of I follows from Fatou's Lemma given that u k u in H implies u k u in H 1 0 (Ω) and thus u k → u pointwise almost everywhere. Hypothesis (iii) is immediate and (iv) holds because I(su)/s 2 = s 2 I(u) where inf B I > 0 for any weakly compact subset B ⊂ H \ {0}. Assumption (v) follows from the above formula for I and the compactness of the embedding H → H 1 0 (Ω) → L 4 (Ω) due to N ∈ {1, 2, 3}. The property (ii) is more difficult to prove. We refer to [24, pp. 31-32] where this has been carried out even in the case of more general nonlinearities. 6.2. An example in 1D. We now show that the general result from above applies in the one-dimensional setting that we already discussed in our bifurcation analysis from Corollary 2.3. So we consider the problem Eq. (6), namely Lemma 6.4. Let Ω, σ, c, κ be given as in Corollary 2.3 with κ(x) ≥ α > 0 for almost all x ∈ Ω. Then Assumptions (A) to (C) hold.
Proof. The assumptions of Corollary 2.3 imply Assumption (A) and Lemma 5.3 yields Assumption (B). So it remains to verify Assumption (C), which is based on the explicit formulas for the eigenpairs from the proof of Lemma 5.5. Moreover, we use Eq. (19),Eq. (21) for τ j := |λ j |, i.e., Recall k ± = c ± /|σ ± | and that a − , σ − are negative whereas a + , σ + , c + , c − are positive. All the following computations of integrals can be checked using the Python library Sympy [19] and can be found in the Jupyter notebook symbolic_1d (with format IPYNB, HTML, and PDF) supply in the zenodo archive 10.5281/zenodo.5592105.
Using the precise expressions for α −2 j from the second step and Eq. (28) we get By the first step, there is F > 0 such that φ j 2 |σ| ≤ F (1 + |λ j |) for all j ∈ Z. 4th step: Formulas and bounds for φ i , φ j |σ| . For i = j explicit computations exploiting Eq.
To estimate the first of these terms set A(z) := z tanh(z). Using the Lipschitz continuity of A, the estimate τ i ≤ √ 1 + λ i and (28) we get for λ i , λ j > 0 for some C > 0. Using the asymptotics for (α i ) from the second step and performing analogous estimates in the other cases we find that there is G > 0 independent of i, j such that | φ i , φ j |σ| | ≤ G 1 + |λ i | 1 + |λ j | 1 + |i| + |j| whenever i = j.
Combining Lemma 6.4 and Theorem 6.3 we thus obtain: Corollary 6.5. Let Ω, σ, c, κ be given as in Corollary 2.3 and λ ∈ R. Then equation Eq. (6) has infinitely many nontrivial solutions in H, among which a least energy solution.

Open problems
We finally address some open problems that we believe to be interesting to study: (1) In Corollary 4.5 we showed that the eigenvalues satisfy the bounds 1 + |λ j | ≥ c (1 + |j|) 2/N for all j ∈ N. An upper bound for the eigenvalues is unfortunately missing, let alone precise Weyl law asymptotics for the counting function Λ → Card {λ j | |λ j | ≤ Λ} that one might compare to well-known ones for the Laplacian.
In the 1D case this may be based on Eq. (28). In the higher-dimensional case we expect new difficulties due to "plasmonic" eigenvalues coming from a concentration near the interface. (2) In view of the physical context described in Section 1.1, the case of sign-changing σ and c and even κ is relevant. As explained in Remark 4.6, there is no hope to get an analogous spectral theory for all c ∈ L ∞ (Ω), but it would be interesting to identify those functions that admit such a one. (3) Corollary 2.3 relies on the precise knowledge of eigenfunctions in the 1D case, especially regarding their nodal structure. Is there a way to find similar properties in more general one-dimensional problems or even in higher-dimensional settings? (4) In Remark 6.6 we commented on the numerical evidence that φ j / 1 + |λ j | j∈Z is a Riesz basis of H 1 0 (Ω). In our one-dimensional case, a theoretical verification seems doable by extensive analysis of explicit formulas as in Lemma 6.4, but we have not found a reasonably short and elegant proof. We believe that such a one could be of interest.
Appendix A. Riesz basis property for the 1d example We numerically illustrate that φ j 1 + |λ j | j∈Z indeed satisfies the Riesz' basis property in the one-dimensional setting studied earlier. In other words, we numerically show that both inequalities · |σ| · and · · |σ| may be expected to hold. Note that the first inequality was rigorously proven in Proposition 6.1. In the following, we use the same notation as in Section 6.
For the numerical investigation, we choose the one-dimensional setting as in Section 3: We consider the linear operator φ → − d dx σ(x) φ on Ω = (−5, 5) with σ(x) = σ + = 1 on Ω + = (0, 5) and various choices for σ(x) = σ − < 0 on Ω − = (−5, 0). We numerically compute the eigenvalues of the matrix M Λ = φ i , φ j |σ| 1 + |λ i | 1 + |λ j | |λ i |,|λ j |≤Λ where the (λ j , φ j ) denote the eigenpairs from Corollary 4.5. The code is available on Zenodo with DOI 10.5281/zenodo.5707422. Figure 9 shows the maximal and minimal eigenvalues in dependence on Λ for different choices of σ − on the left and right, respectively. For all considered choices of σ − , we clearly observe that the maximal and minimal eigenvalues asymptotically tend to a finite, non-zero value. As discussed above, this behavior was rigorously proved for the maximal eigenvalue in Proposition 6.1 so that we focus on the minimal eigenvalue. The asymptotic value is reached quickly. Note that for |σ − | → 0 the minimal eigenvalue is expected to degenerate to zero, which explains the dependence on σ − observable in Fig. 9 (right). This is in line with Remark 5.6.