Junyan Cao

Jean-Pierre Demailly

Shin

A GENERAL EXTENSION THEOREM FOR COHOMOLOGY CLASSES ON NON REDUCED ANALYTIC SUBSPACES

Keywords: April 2, 2017, version 2. 2010 Mathematics Subject Classification. Primary 32L10; Secondary 32E05 Compact Kähler manifold, singular hermitian metric, coherent sheaf cohomology, Dolbeault cohomology, plurisubharmonic function, L 2 estimates, Ohsawa-Takegoshi extension theorem, multiplier ideal sheaf

The main purpose of this paper is to generalize the celebrated L 2 extension theorem of Ohsawa-Takegoshi in several directions : the holomorphic sections to extend are taken in a possibly singular hermitian line bundle, the subvariety from which the extension is performed may be non reduced, the ambient manifold is Kähler and holomorphically convex, but not necessarily compact.

Introduction and preliminaries

The purpose of this paper is to generalize the celebrated L 2 extension theorem of Ohsawa-Takegoshi [START_REF] Ohsawa | On the extension of L 2 holomorphic functions[END_REF] under the weakest possible hypotheses, along the lines of [START_REF] Demailly | Extension of holomorphic functions defined on non reduced analytic subvarieties[END_REF] and [START_REF] Matsumura | An injectivity theorem with multiplier ideal sheaves for higher direct images under Kähler morphisms[END_REF]. Especially, the ambient complex manifold X is a Kähler manifold that is only assumed to be holomorphically convex, and is not necessarily compact; by the Remmert reduction theorem, this is the same as a Kähler manifold X that admits a proper holomorphic map π : X → S onto a Stein complex space S. This allows in particular to consider relative situations over a Stein base. We consider a holomorphic line bundle E → X equipped with a singular hermitian metric h, namely a metric which can be expressed locally as h = e -ϕ where ϕ is a quasi-psh function, i.e. a function that is locally the sum ϕ = ϕ 0 + u of a plurisubharmonic function ϕ 0 and of a smooth function u. Such a bundle admits a curvature current (1) Θ E,h := i∂∂ϕ = i∂∂ϕ 0 + i∂∂u which is locally the sum of a positive (1, 1)-current i∂∂ϕ 0 and a smooth (1, 1)-form i∂∂u. Our goal is to extend sections that are defined on a (non necessarily reduced) complex subspace Y ⊂ X, when the structure sheaf O Y := O X /I(e -ψ ) is given by the multiplier ideal sheaf of a quasi-psh function ψ with neat analytic singularities, i.e. locally on a neighborhood V of an arbitrary point x 0 ∈ X we have

(2)

ψ(z) = c log |g j (z)| 2 + v(z), g j ∈ O X (V ), v ∈ C ∞ (V ).
Let us recall that the multiplier ideal sheaf I(e -ϕ ) of a quasi-psh function ϕ is defined by

(3)

I(e -ϕ ) x 0 = f ∈ O X,x 0 ; ∃U ∋ x 0 , U |f | 2 e -ϕ dλ < +∞
with respect to the Lebesgue measure λ in some local coordinates near x 0 . As usual, we also denote by K X = Λ n T * X the canonical bundle of an n-dimensional complex manifold X. As is well known, I(e -ϕ ) ⊂ O X is a coherent ideal sheaf (see e.g. [Dem-book]). Our main result is given by the following general statement.

Theorem 1.1. Let E be a holomorphic line bundle over a holomorphically convex Kähler manifold X. Let h be a possibly singular hermitian metric on E, ψ a quasi-psh function with neat analytic singularities on X. Assume that there exists a positive continuous function δ > 0 on X such that (4) Θ E,h + (1 + αδ)i∂∂ψ ≥ 0 in the sense of currents, for all α ∈ [0, 1].

Then the morphism induced by the natural inclusion I(he -ψ ) → I(h)

(5) H q (X, K X ⊗ E ⊗ I(he -ψ )) → H q (X, K X ⊗ E ⊗ I(h))

is injective for every q ≥ 0. In other words, the morphism induced by the natural sheaf surjection I(h) → I(h)/I(he -ψ ) (6) H q (X, K X ⊗ E ⊗ I(h)) → H q (X, K X ⊗ E ⊗ I(h)/I(he -ψ ))

is surjective for every q ≥ 0.

Remark 1.2. If h is smooth, we have I(h) = O X and I(h)/I(he -ψ ) = O X /I(e -ψ ) := O Y where Y is the zero subvariety of the ideal sheaf I(e -ψ ). Then for q = 0, the surjectivity statement can be interpreted an extension theorem for holomorphic sections, with respect to the restriction morphism

(7) H 0 (X, K X ⊗ E) → H 0 (Y, (K X ⊗ E) |Y ).
In general, the quotient sheaf I(h)/I(he -ψ ) is supported in an analytic subvariety Y ⊂ X, which is the zero set of the quotient ideal

J Y := I(he -ψ ) : I(h) = f ∈ O X ; f • I(h) ⊂ I(he -ψ ) ,
and (6) can be considered as a restriction morphism to Y .

The crucial idea of the proof is to prove the results (say, in the form of the surjectivity statement), only up to approximation. This is done by solving a ∂-equation

∂u ε + w ε = v
where the right hand side v is given and w ε is an error term such that w ε = O(ε a ) as ε → 0, for some constant a > 0. A twisted Bochner-Kodaira-Nakano identity introduced by Donnelly and Fefferman [START_REF] Donnelly | L 2 -cohomology and index theorem for the Bergman metric[END_REF], and Ohsawa and Takegoshi [START_REF] Ohsawa | On the extension of L 2 holomorphic functions[END_REF] is used for that purpose, with an additional correction term. The version we need can be stated as follows.

Proposition 1.3. (see [START_REF] Demailly | Extension of holomorphic functions defined on non reduced analytic subvarieties[END_REF]Prop. 3.12]) Let X be a complete Kähler manifold equipped with a (non necessarily complete) Kähler metric ω, and let (E, h) be a Hermitian vector bundle over X. Assume that there are smooth and bounded functions η, λ > 0 on X such that the curvature operator

B = B n,q E,h,ω,η,λ = [η Θ E,h -i ∂∂η -iλ -1 d∂η ∧ ∂η, Λ ω ] ∈ C ∞ (X, Herm(Λ n,q T * X ⊗ E))
satisfies B + εI > 0 for some ε > 0 (so that B can be just semi-positive or even slightly negative; here I is the identity endomorphism). Given a section v ∈ L 2 (X, Λ n,q T * X ⊗ E) such that ∂v = 0 and

M (ε) := X (B + εI) -1 v, v dV X,ω < +∞,
there exists an approximate solution f ε ∈ L 2 (X, Λ n,q-1 T * X ⊗ E) and a correction term

w ε ∈ L 2 (X, Λ n,q T * X ⊗ E) such that ∂u ε = v -w ε and X (η + λ) -1 |u ε | 2 dV X,ω + 1 ε X |w ε | 2 dV X,ω ≤ M (ε).
Moreover, if v is smooth, then u ε and w ε can be taken smooth.

In our situation, the main part of the solution, namely u ε , may very well explode as ε → 0. In order to show that the equation ∂u = v can be solved, it is therefore needed to check that the space of coboundaries is closed in the space of cocycles in the Fréchet topology under consideration (here, the L 2 loc topology), in other words, that the related cohomology group H q (X, F) is Hausdorff. In this respect, the fact of considering ∂-cohomology of smooth forms equipped with the C ∞ topology on the one hand, or cohomology of forms u ∈ L 2 loc with ∂u ∈ L 2 loc on the other hand, yields the same topology on the resulting cohomology group H q (X, F). This comes from the fact that both complexes yield fine resolutions of the same coherent sheaf F, and the topology of H q (X, F) can also be obtained by using Čech cochains with respect to a Stein covering U of X. The required Hausdorff property then comes from the following well known fact.

Lemma 1.4. Let X be a holomorphically convex complex space and F a coherent analytic sheaf over X. Then all cohomology groups H q (X, F) are Hausdorff with respect to their natural topology (induced by the Fréchet topology of local uniform convergence of holomorphic cochains). 1

In fact, the Remmert reduction theorem implies that X admits a proper holomorphic map π : X → S onto a Stein space S, and Grauert's direct image theorem shows that all direct images R q π * F are coherent sheaves on S. Now, as S is Stein, Leray's theorem combined with Cartan's theorem B tells us that we have an isomorphism H q (X, F) ≃ H 0 (S, R q π * F). More generally, if U ⊂ S is a Stein open subset, we have

(8) H q (π -1 (U ), F) ≃ H 0 (U, R q π * F)
and when U ⋐ S is relatively compact, it is easily seen that this a topological isomorphism of Fréchet spaces since both sides are O S (U ) modules of finite type and can be seen as a Fréchet quotient of some direct sum O S (U ) ⊕N by looking at local generators and local relations of R q π * F. Therefore H q (X, F) ≃ H 0 (S, R q π * F) is a topological isomorphism and the space of sections in the right hand side is a Fréchet space. In particular, H q (X, F) is Hausdorff.

The isomorphism (8) shows that it is enough to prove Theorem 1.1 locally over X, i.e., we can replace X by X ′ = π -1 (S ′ ) ⋐ X where S ′ ⋐ S. Therefore, we can assume that δ > 0 is a constant rather than a continuous function.

Proof of the extension theorem

In this section, we give a proof of Theorem 1.1 based on a generalization of the arguments of [Dem15b, Th. 2.14]. We start by proving the special case of the extension result for holomorphic sections (q = 0). Theorem 2.1. Let (X, ω) be a holomorphically convex Kähler manifold and ψ be a quasipsh function with neat analytic singularities. Let E be a line bundle with a possibly singular metric h, and Y the support of the sheaf I(h)/I(he -ψ ), along with the structure sheaf O Y := I(he -ψ ) : I(h). Assume that there is a continuous function δ > 0 such that

iΘ E,h + (1 + αδ)i∂∂ψ ≥ 0 in the sense of currents, for all α ∈ [0, 1].
Then the restriction morphism

H 0 (X, O X (K X ⊗ E) ⊗ I(h)) → H 0 (Y, O X (K X ⊗ E) ⊗ I(h)/I(he -ψ ) |Y ) is surjective.
Proof. (a) Let us first assume for simplicity that h is smooth. We will explain the general case later. Then I(h) = O X and I(h)/I(he -ψ ) = O Y = O X /I(e -ψ ). After possibly shrinking X into a relatively compact holomorphically convex open subset X ′ = π -1 (S ′ ) ⋐ X, we can suppose

1 It was pointed out to us by Prof. Takeo Ohsawa that this result does not hold under the assumption that X is weakly pseudoconvex, i.e., if we only assume that X admits a smooth psh exhaustion. A counter-example can be derived from [START_REF] Kazama | ∂-cohomology of (H, C)-groups[END_REF]. As a consequence, it is unclear whether the results of the present paper extend to the Kähler weakly pseudoconvex case, although the main L 2 estimates are still valid in that situation.

that δ > 0 is a constant and that ψ ≤ 0, after subtracting a large constant to ψ. Also, without loss of generality, we can assume that ψ admits a discrete sequence of "jumping numbers"

(9) 0 = m 0 < m 1 < • • • < m p < • • • such that I(mψ) = I(m p ψ) for m ∈ [m p , m p+1 [.
Since ψ is assumed to have analytic singularities, this follows from using a log resolution of singularities, thanks to the Hironaka desingularization theorem (by the much deeper result of [START_REF] Qi'an Guan | A proof of Demailly's strong openness conjecture[END_REF] on the strong openness conjecture, one could even possibly eliminate the assumption that ψ has analytic singularities). We fix here p such that m p ≤ 1 < m p+1 , and in the notation of [START_REF] Demailly | Extension of holomorphic functions defined on non reduced analytic subvarieties[END_REF], we let Y = Y (mp) be defined by the non necessarily reduced structure sheaf

O Y = O X /I(e -ψ ) = O X /I(e -mpψ ).
Step 1 (Construction of a smooth extension). Take

f ∈ H 0 (Y, O X (K X ⊗ E) |Y ) = H 0 (X, O X (K X ⊗ E) ⊗ O X /I(e -mpψ )).
Let U = (U i ) be a Stein covering of X and let (ρ i ) be a partition of unity subordinate to (U i ).

Thanks to the exact sequence

(10) 0 → I(e -ψ ) → O X → O X /I(e -ψ ) → 0, we can find a f i ∈ H 0 (U i , O X (K X ⊗ E)) such that f i | Y ∩U i = f | Y ∩U i .
Then (10) implies that (11)

f i -f j ∈ H 0 (U i ∩ U j , O X (K X ⊗ E) ⊗ I(e -ψ )).
As a consequence, the smooth section f :=

i ρ i • f i is a smooth extension of f and satisfies ∂ f = i (∂ρ i ) • ( f i -f j ) on U j , hence (12)
X |∂ f | 2 ω,h e -ψ dV X,ω = X j ρ j i (∂ρ i ) • ( f i -f j ) 2 ω,h e -ψ dV X,ω < +∞.
Step 2 (L 2 -estimates). We follow here the arguments of [Dem15b, proof of th. 2.14, p. 217]. Let t ∈ Z -and let χ t be the negative convex increasing function defined in [Dem15b, (5.8 * ), p. 211]. Put

η t := 1 -δ • χ t (ψ) and λ t := 2δ (χ 2 t (ψ)) 2 χ ′′ t (ψ) . We set R t := η t (Θ E,h + i∂∂ψ) -i∂∂η t -λ -1 t i∂η t ∧ ∂η t = η t (Θ E,h + (1 + δη -1 t χ ′ t (ψ))i∂∂ψ) + δ • χ ′′ t (ψ) 2 i∂ψ ∧ ∂ψ. Note that χ ′′ t (ψ) ≥ 1 8 on W t = {t < ψ < t + 1}. The curvature assumption (4) implies Θ E,h + (1 + δη -1 t χ ′ t (ψ)) i∂∂ψ ≥ 0 on X.
As in [START_REF] Demailly | Extension of holomorphic functions defined on non reduced analytic subvarieties[END_REF], we find

(13) R t ≥ 0 on X and (14) R t ≥ δ 16 i∂ψ ∧ ∂ψ on W t = {t < ψ < t + 1}. Let θ : [-∞, +∞[ → [0, 1] be a smooth non increasing real function satisfying θ(x) = 1 for x ≤ 0, θ(x) = 0 for x ≥ 1 and |θ ′ | ≤ 2.
By applying the L 2 estimate (Proposition 1.3), for every ε > 0 we can find sections u t,ε , w t,ε satisfying (15)

∂u t,ε + w t,ε = v t := ∂ θ(ψ -t) • f and (16) X (η t +λ t ) -1 |u t,ε | 2 ω,h e -ψ dV X,ω + 1 ε X |w t,ε | 2 ω,h e -ψ dV X,ω ≤ X (R t +εI) -1 v t , v t e -ψ dV X,ω , where (17) v t = ∂ θ(ψ -t) f = θ ′ (ψ -t) ∂ψ ∧ f + θ(ψ -t) ∂ f .
Combining (13), ( 14), ( 16) and (17), we get X |u t,ε | 2 ω,h e -ψ dV X,ω < +∞ and (18)

X |w t,ε | 2 ω,h e -ψ dV X,ω ≤ 128 ε δ {t<ψ<t+1} | f | 2 ω,h e -ψ dV X,ω + 2 {ψ<t+1} |∂ f | 2 ω,h e -ψ dV X,ω .
We now estimate the right hand side of (18). Since f is smooth, we have an obvious upper bound of the first term

(19) {t<ψ<t+1} | f | 2 ω,h e -ψ dV X,ω ≤ C 1 e -t ,
where C 1 is the C 0 norm of f . For the second term, thanks to (9), ( 11) and ( 12), we have (20)

X |∂ f | 2 ω,h e -(1+α)ψ dV X,ω < +∞ for any α ∈ ]0, m p+1 -1[. As a consequence, we get (21) {ψ<t+1} |∂ f | 2 ω,h e -ψ dV X,ω ≤ C 2 e αt
for some constant C 2 depending only on α. By taking ε = e (1+α)t , ( 18), ( 19) and ( 21) imply

(22) X |w t,ε | 2 ω,h e -ψ dV X,ω ≤ C 3 e αt = O(ε α 1+α ),
for some constant C 3 , whence the error tends to 0 as t → -∞ and ε → 0.

Step 3 (Final conclusion). Putting everything together and redefining u t = u t,ε , w t = w t,ε for simplicity of notation, we get

(23) ∂(θ(ψ -t) • f -u t ) = w t , X |u t | 2 h e -ψ dV X,ω < +∞ and (24) lim t→-∞ X |w t | 2 ω,h e -ψ dV X,ω = 0.
After shrinking X, we can assume that we have a finite Stein covering U = (U i ) where the U i are biholomorphic to bounded pseudoconvex domains. The standard Hörmander L 2 estimates then provide L 2 sections s t,j on U j such that ∂s t,j = w t on U j and

(25) lim t→-∞ U j |s t,j | 2 ω,h e -ψ dV X,ω = 0. Then ∂ θ(ψ -t) • f -u t - j ρ j s t,j = - j (∂ρ j ) • s t,j on X = - j (∂ρ j ) • (s t,j -s t,i ) on U i . (26)
As ∂(s t,j -s t,i ) = 0 on U i ∩ U j , the difference is holomorphic and the right hand side of (26) is smooth. Moreover, (25) shows that these differences converge uniformly to 0, hence the right hand side of (26) converges to 0 in C ∞ topology. The left hand side implies that this is a coboundary in the C ∞ Dolbeault resolution of O X (K X ⊗ E). By applying Lemma 1.4, we conclude that there is a C ∞ section σ t of K X ⊗ E converging uniformly to 0 on compact subsets of X as t → -∞, such that ∂σ t = j (∂ρ j ) • s t,j on X. This implies that

f t := θ(ψ -t) • f -u t - j ρ j s t,j + σ t
is holomorphic on X. Hörmander's L 2 estimates also produce local smooth solutions σ t,i on U i with the additional property that lim t→-∞ U i |σ t,i | 2 ω,h e -ψ dV X,ω = 0. Therefore

f t,i := θ(ψ -t) • f -u t - j ρ j s t,j + σ t,i
is holomorphic on U i and f t -f t,i converges uniformly to 0 on compact subsets of U i . However, by construction, f t,i -f i is a holomorphic section on U i that satisfies the L 2 estimate with respect to the weight e -ψ , hence

f t,i -f i is a section of O X (K X ⊗ E) ⊗ I(e -ψ ) on U i , in other words the image of f t,i in H 0 (U i , O X (K X ⊗ E) ⊗ O X /I(e -ψ ))
coincides with f |U i . As a consequence, the image of f t in

H 0 (X, O X (K X ⊗ E) ⊗ O X /I(e -ψ )) = H 0 (Y, (K X ⊗ E) |Y )
converges to f . By the direct image argument used in the preliminary section, this density property implies the surjectivity of the restriction morphism to Y .

(b) We now prove the theorem for the general case when h = e -ϕ is not necessarily smooth. We can reduce ourselves to the case when ψ has divisorial singularities (see [START_REF] Demailly | Extension of holomorphic functions defined on non reduced analytic subvarieties[END_REF] or the next section for a more detailed argument). Let us pick a section

f ∈ H 0 (X, O X (K X ⊗ E) ⊗ I(h)/I(he -ψ )).
By using the same reasoning as in Step 1, we can find a smooth extension f

∈ C ∞ (X, K X ⊗ E) of f such that (27) X |∂ f | 2 ω,h e -ψ dV X,ω < +∞.
For every t ∈ Z -fixed, as ψ has divisorial singularities, we still have

Θ E,h + (1 + δη -1 t χ ′ t (ψ))(i∂∂ψ) ac ≥ 0 on X,
where (i∂∂ψ) ac is the absolutely continuous part of i∂∂ψ. The regularization techniques of [START_REF] Demailly | Pseudo-effective line bundles on compact Kähler manifolds[END_REF] and [Dem15a, Th. 1.7, Remark 1.11] (cf. also the next section) produce a family of singular metrics {h t,ε } +∞ k=1 which are smooth in the complement X Z t,ε of an analytic set, such that I(h t,ε ) = I(h), I(h t,ε e -ψ ) = I(he -ψ ) and

Θ E,ht,ε + (1 + δη -1 t χ ′ t (ψ)) i∂∂ψ ≥ - 1 2 εω on X.
The additional error term -1 2 εω is irrelevant when we use Proposition 1.3, as it is absorbed by taking the hermitian operator B + εI. Therefore for every t ∈ Z -, with the adjustment ε = e αt , α ∈ ]0, m k+1 -1[, we can find a singular metric h t = h t,ε which is smooth in the complement X \ Z t of an analytic set, such that I(h t ) = I(h), I(h t e -ψ ) = I(he -ψ ) and h t ↑ h as t → -∞, and approximate solutions of the ∂-equation such that

∂(θ(ψ -t) • f -u t ) = w t , X |u t | 2 ω,ht e -ψ dV X,ω < +∞ and lim t→-∞ X |w t | 2 ω,ht e -ψ dV X,ω = 0.
Proposition 1.3 can indeed be applied since X Z t is complete Kähler (at least after we shrink X a little bit as X ′ = π -1 (S ′ ), cf. [START_REF] Demailly | Estimations L 2 pour l'opérateur ∂ d'un fibré vectoriel holomorphe semi-positif au-dessus d'une variété kählérienne complète[END_REF]). The theorem is then proved by using the same argument as in Step 3; it is enough to notice that the holomorphic sections s t,j -s t,i and f t,i -f i satisfy the L 2 -estimate with respect to (h t , ψ) [instead of the expected (h, ψ)], but the multiplier ideal sheaves involved are unchanged. The Hausdorff property is applied to the cohomology group

H 1 (X, O X (K X ⊗ E) ⊗ I(h)) instead of H 1 (X, K X ⊗ E)
, and the density property to the morphism of direct image sheaves

π * O X (K X ⊗ E) ⊗ I(h) → π * O X (K X ⊗ E) ⊗ I(h)/I(he -ψ )
over the Stein space S.

Proof of the extension theorem for degree q cohomology classes. The reasoning is extremely similar, so we only explain the few additional arguments needed. In fact, Proposition 1.3 can be applied right away to arbitrary (n, q)-forms with q ≥ 1, and the twisted Bochner-Kodaira-Nakano inequality yields exactly the same estimates. Any cohomology class in

H q (Y, O X (K X ⊗ E) ⊗ I(h)/I(he -ψ ))
is represented by a holomorphic Čech q-cocycle with respect to the Stein covering U = (U i ), say

(c i 0 ...iq ), c i 0 ...iq ∈ H 0 U i 0 ∩ . . . ∩ U iq , O X (K X ⊗ E) ⊗ I(h)/I(he -ψ ) .
By the standard sheaf theoretic isomorphisms with Dolbeault cohomology (cf. e.g. [Dem-e-book]), this class is represented by a smooth (n, q)-form

f = i 0 ,...,iq c i 0 ...iq ρ i 0 ∂ρ i 1 ∧ . . . ∂ρ iq
by means of a partition of unity (ρ i ) subordinate to (U i ). This form is to be interpreted as a form on the (non reduced) analytic subvariety Y associated with the ideal sheaf J = I(he -ψ ) : I(h) and the structure sheaf O Y = O X /J . We get an extension as a smooth (no longer ∂-closed) (n, q)-form on X by taking

f = i 0 ,...,iq c i 0 ...iq ρ i 0 ∂ρ i 1 ∧ . . . ∂ρ iq where c i 0 ...iq is an extension of c i 0 ...iq from U i 0 ∩ . . . ∩ U iq ∩ Y to U i 0 ∩ . . . ∩ U iq .
Again, we can find approximate L 2 solutions of the ∂-equation such that

∂(θ(ψ -t) • f -u t ) = w t , X |u t | 2 ω,ht e -ψ dV X,ω < +∞ and lim t→-∞ X |w t | 2 ω,ht e -ψ dV X,ω = 0.
The difficulty is that L 2 sections cannot be restricted in a continuous way to a subvariety. In order to overcome this problem, we play again the game of returning to Čech cohomology by solving inductively ∂-equations for w t on U i 0 ∩ . . . ∩ U i k , until we reach an equality

(28) ∂ θ(ψ -t) • f -u t = w t := - i 0 ,...,i q-1 s t,i 0 ...iq ∂ρ i 0 ∧ ∂ρ i 1 ∧ . . . ∂ρ iq with holomorphic sections s t,I = s t,i 0 ...iq on U I = U i 0 ∩ . . . ∩ U iq , such that lim t→-∞ U I |s t,I | 2 ω,ht e -ψ dV X,ω = 0.
Then the right hand side of (28) is smooth, and more precisely has coefficients in the sheaf C ∞ ⊗ O I(he -ψ ), and w t → 0 in C ∞ topology. A priori, u t is an L 2 (n, q)-form equal to u t plus a combination ρ i s t,i of the local solutions of ∂s t,i = w t , plus ρ i s t,i,j ∧ ∂ρ j where ∂s t,i,j = s t,j -s t,i , plus etc . . . , and is such that

X | u t | 2
ω,ht e -ψ dV X,ω < +∞.

Since H q (X, O X (K X ⊗ E) ⊗ I(he -ψ )) can be computed with the L 2 loc resolution of the coherent sheaf, or alternatively with the ∂-complex of (n, •)-forms with coefficients in C ∞ ⊗ O I(he -ψ ), we may assume that u t ∈ C ∞ ⊗ O I(he -ψ ), after playing again with Čech cohomology. Lemma 1.4 yields a sequence of smooth (n, q)-forms σ t with coefficients in C ∞ ⊗ O I(h), such that ∂σ t = w t and σ t → 0 in C ∞ -topology. Then

f t = θ(ψ -t) • f -u t -σ t is a ∂-closed (n, q)-form on X with values in C ∞ ⊗ O I(h) ⊗ O X (E), whose image in H q (X, O X (K X ⊗ E) ⊗ I(h)/I(he -ψ )) converges to {f } in C ∞ Fréchet
topology. We conclude by a density argument on the Stein space S, by looking at the coherent sheaf morphism

R q π * O X (K X ⊗ E) ⊗ I(h) → R q π * O X (K X ⊗ E) ⊗ I(h)/I(he -ψ ) .

An alternative proof based on injectivity theorems

We give here an alternative proof based on injectivity theorems, in the case when X is compact Kähler. The case of a holomorphically convex manifold is entirely similar, so we will content ourselves to indicate the required additional arguments at the end.

Proof of Theorem 1.1. First of all, we reduce the proof of Theorem 1.1 to the case when ψ has divisorial singularities. Since ψ has analytic singularities, there exists a modification π : X ′ → X such that the pull-back π * ψ has divisorial singularities. For the singular hermitian line bundle (E ′ , h ′ ) := (π * E, π * h) and the quasi-psh function ψ ′ := π * ψ, we can easily check that

π * (K X ′ ⊗ E ′ ⊗ I(h ′ e -ψ ′ )) = K X ⊗ E ⊗ I(he -ψ ), π * (K X ′ ⊗ E ′ ⊗ I(h ′ )) = K X ⊗ E ⊗ I(h).
Hence we obtain the following commutative diagram :

H q (X, K X ⊗ E ⊗ I(he -ψ )) ∼ = π * f / / H q (X, K X ⊗ E ⊗ I(h)) π * H q (X ′ , K X ′ ⊗ E ′ ⊗ I(h ′ e -ψ ′ )) g / / H q (X ′ , K X ′ ⊗ E ′ ⊗ I(h ′ )),
where f , g are the morphisms induced by the natural inclusions and π * is the natural edge morphism. It follows that the left edge morphism π * is an isomorphism since the curvature of the singular hermitian metric h ′ e -ψ ′ on E ′ is semi-positive by the assumption. Indeed, even if h ′ does not have analytic singularities, we can see that R q π * (K X ′ ⊗ E ′ ⊗ I(h ′ e -ψ ′ )) = 0 for every q > 0 by [Mat16b, Corollary 1.5]. (In the case of X being a projective variety, a relatively easy proof can be found in [START_REF] Fujino | Injectivity theorem for pseudo-effective line bundles and its applications[END_REF].) If Theorem 1.1 can be proven when ψ has divisorial singularities, it follows that the morphism g in the above diagram is injective since (E ′ , h ′ ) = (π * E, π * h) and ψ ′ = π * ψ satisfy the assumptions in Theorem 1.1 and ψ ′ has divisorial singularities. Therefore the morphism f is also injective by the commutative diagram. Now we explain the idea of the proof of Theorem 1.1. If we can obtain equisingular approximations h ε of h satisfying the following properties : Θ E,hε + i∂∂ψ ≥ -εω and Θ E,hε + (1 + δ)i∂∂ψ ≥ -εω, then a proof similar to [START_REF] Fujino | Injectivity theorem for pseudo-effective line bundles and its applications[END_REF] works, where ω is a fixed Kähler form on X. In the case when ψ has divisorial singularities, we can attain either of the above curvature properties, but we do not know whether we can attain them at the same time. For this reason, we will look for an essential curvature condition arising from the assumptions on the curvatures in Theorem 1.1, in order to use the "twisted" Bochner-Kodaira-Nakano identity.

From now on, we consider a quasi-psh ψ with divisorial singularities. Then there exist an effective R-divisor D and a smooth (1, 1)-form γ on X such that

i 2π ∂∂ψ = [D] + 1 2π γ
in the sense of (1, 1)-currents, where [D] denotes the current of the integration over D. For the irreducible decomposition D = N i=1 a i D i and the defining section t i of D i , we can take a smooth hermitian metric b i on D i such that

e ψ = |s| 2 b := |t 1 | 2a 1 b 1 |t 2 | 2a 2 b 2 • • • |t N | 2a N b N and -γ = Θ b (D) := N i=1 a i Θ b i (D i ).
For a positive number 0 < c ≪ 1, we define the continuous functions σ and η on X by

σ = σ c := log(|s| 2 b + c) and η = η c := 1 c -χ(σ),
where χ(t) := t -log(-t).

Remark 3.1. (i) We may assume that |s| 2 b < 1/5 by subtracting a positive constant from ψ. Further we may assume that σ < log(1/5) and 1 < χ ′ (σ) < 7/4 by choosing a sufficiently small c > 0. (ii) Further, the function η is a continuous function on X with η > 1/c. The function η is smooth on X \ D, but it need not be smooth on X since |s| 2 b is not smooth in the case when 0 < a i < 1 for some i.

Throughout the proof, we fix a Kähler form ω on X. The following proposition gives a suitable approximation of a singular hermitian metric h on E, which enables us to use the twisted Bochner-Kodaira-Nakano identity. The proof is based on the argument in [START_REF] Ohsawa | On a curvature condition that implies a cohomology injectivity theorem of Kollár-Skoda type[END_REF], [START_REF] Fujino | A transcendental approach to Kollár's injectivity theorem II[END_REF] and the equisingular approximation theorem in [DPS01, Theorem 2.3].

Proposition 3.2. There exist singular hermitian metrics {h ε } 0<ε≪1 on E with the following properties :

(a) h ε is smooth on X \ Z ε , where Z ε is a proper subvariety on X. (b) h ε ′ ≤ h ε ′′ ≤ h holds on X for ε ′ > ε ′′ > 0. (c) I(h) = I(h ε ) and I(he -ψ ) = I(h ε e -ψ ) on X. (d) η(Θ hε (E) + γ) -i∂∂η -η -2 i∂η ∧ ∂η ≥ -εω on X \ D.
(e) For arbitrary t > 0, by taking a sufficiently small ε > 0, we have e -tφ -e -tφε < ∞, where φ (resp. φ ε ) is a local weight of h (resp. h ε ).

Proof. We fix a a sufficiently small c with 7c/4 ≤ δ. Then, by Remark 3.1, we can easily check that

χ ′ (σ)|s| 2 b η(|s| 2 b + c) ≤ 7 4η ≤ 7 4 c ≤ δ.
In particular, it follows that

Θ h + 1 + χ ′ (σ)|s| 2 b η(|s| 2 b + c) γ ≥ 0 on X
since ψ has divisorial singularities and satisfies the assumptions in Theorem 1.1. By applying the equisingular approximation theorem ([DPS01, Theorem 2.3]) to h, we can take singular hermitian metrics {h ε } 0<ε≪1 on E satisfying properties (a), (b), (e), the former conclusion of (c), and the following curvature property:

Θ hε + 1 + χ ′ (σ)|s| 2 b η(|s| 2 b + c) γ ≥ -εω on X.
Now we check property (d) from the above curvature property. The function η may not be smooth on X, but it is smooth on X \ D. Therefore the same computation as in [START_REF] Ohsawa | On a curvature condition that implies a cohomology injectivity theorem of Kollár-Skoda type[END_REF] and [START_REF] Fujino | A transcendental approach to Kollár's injectivity theorem II[END_REF] works on X \ D. In particular, from a complicated but straightforward computation, we obtain [START_REF] Fujino | A transcendental approach to Kollár's injectivity theorem II[END_REF] for the precise computation). Then, by -γ = Θ b (D) on X \ D, we can see that

-i∂∂η = - χ ′ (σ)|s| 2 b |s| 2 b + c Θ b (D) + c χ ′ (σ)|s| 2 b + χ ′′ (σ) χ ′ (σ) 2 i∂η ∧ ∂η on X \ D (see
η(Θ hε (E) + γ) -i∂∂η - 1 η 2 i∂η ∧ ∂η = c χ ′ (σ)|s| 2 b + χ ′′ (σ) χ ′ (σ) 2 - 1 η 2 i∂η ∧ ∂η + η Θ hε (E) + (1 + χ ′ (σ)|s| 2 b η(|s| 2 b + c) )γ ≥ c χ ′ (σ)|s| 2 b + χ ′′ (σ) χ ′ (σ) 2 - 1 η 2 i∂η ∧ ∂η -εηω on X \ D.
A straightforward computation yields that χ ′′ (σ)/χ ′ (σ) 2 ≥ 1/η 2 , and thus the first term is semi-positive. Since η is bounded above, we infer that property (d) holds.

Finally we check the last conclusion of property (c) by proving the following lemma, which can be obtained from the strong openness theorem (see [START_REF] Qi'an Guan | A proof of Demailly's strong openness conjecture[END_REF], [START_REF] Lempert | Modules of square integrable holomorphic germs[END_REF], [START_REF] Hiêp | The weighted log canonical threshold[END_REF]) and property (e).

Lemma 3.3. For a quasi-psh function ϕ, we have I(he -ϕ ) = I(h ε e -ϕ ). In particular, we obtain the last conclusion of property (c).

Proof. We have the inclusion I(he -ϕ ) ⊂ I(h ε e -ϕ ) by h ε ≤ h. To get the converse inclusion, we consider a local holomorphic function g such that |g| 2 e -ϕ-φε is integrable, where φ ε (resp. φ) is a local weight of h ε (resp. h). Then Hölder's inequality yields

|g| 2 e -φ-ϕ = |g| 2 e -ϕ-φε e -φ+φε ≤ |g| 2p e -p(ϕ+φε) 1/p • e -q(φ-φε) 1/q ,
where p, q are real numbers such that 1/p + 1/q = 1 and p > 1. By the strong openness theorem, the function |g| 2p e -p(ϕ+φε) is integrable when p is sufficiently close to one. On the other hand, we have e -q(φ-φε) -1 = e qφε e -qφ -e -qφε ≤ sup e qφε e -qφ -e -qφε .

The right hand side is finite for a sufficiently small ε by property (e).

This concludes the proof of Proposition 3.2.

From now on, we proceed to prove Theorem 1.1 by using Proposition 3.2. In the same way as in [FM16, Section 5], one constructs a family of complete Kähler forms {ω ε,δ } 0<δ≪1 on Y ε := X \ (Z ε ∪ D) with the following properties :

(A) ω ε,δ is a complete Kähler form on Y ε := X \ (Z ε ∪ D) for every δ > 0. (B) ω ε,δ ≥ ω on Y ε for every δ ≥ 0.
(C) For every point p in X, there exists a bounded function Ψ ε,δ on an open neighborhood B p such that ω ε,δ = i∂∂Ψ ε,δ on B p and Ψ ε,δ converges uniformly to a bounded function that is independent of ε.

For simplicity, we put H := he -ψ and H ε := h ε e -ψ . We consider a cohomology class β ∈ H q (X, K X ⊗ E ⊗ I(H)) such that β = 0 ∈ H q (X, K X ⊗ E ⊗ I(h)). By the De Rham-Weil isomorphism

H q (X, K X ⊗ E ⊗ I(H)) ∼ = Ker ∂ : L n,q (2) (E) H,ω → L n,q+1 (2)
(E) H,ω Im ∂ : L n,q-1 (2) (E) H,ω → L n,q (2) (E) H,ω
, the cohomology class β can be represented by a ∂-closed E-valued (n, q)-form u with u H,ω < ∞ (that is, β = {u}). Here L n,• (2) (E) H,ω is the L 2 -space of E-valued (n, •)-forms on X with respect to the L 2 -norm • H,ω defined by

• 2 H,ω := X | • | 2 H,ω dV ω ,
where dV ω := ω n /n! and n := dim X. For the L 2 -norm • Hε,ω ε,δ defined by

• 2 ε,δ := • 2 Hε,ω ε,δ := X | • | 2 Hε,ω ε,δ dV ω ε,δ , one can easily check that u ε,δ ≤ u H,ω ε,δ ≤ u H,ω < ∞. (29)
Indeed, the first inequality is obtained from property (b), and the second inequality is obtained from property (B) for ω ε,δ (for example see [START_REF] Fujino | Injectivity theorem for pseudo-effective line bundles and its applications[END_REF]Lemma 2,4]). In particular, we see that u belongs to the L 2 -space L n,q (2) (E) ε,δ := L n,q

(2) (Y ε , E) Hε,ω ε,δ of E-valued (n, q)-forms on Y ε (not X) with respect to • ε,δ . By the orthogonal decomposition (see for example [Mat16a, Proposition 5.8])

L n,q (2) (F ) ε,δ = Im ∂ ⊕ H n,q ε,δ (F ) ⊕ Im ∂ * ε,δ
, the E-valued form u can be decomposed as follows :

u = ∂w ε,δ + u ε,δ for some w ε,δ ∈ Dom ∂ ⊂ L n,q-1 (2) (E) ε,δ , and u ε,δ ∈ H n,q ε,δ (E). ( 30
)
Here ∂ * ε,δ is (the maximal extension of) the formal adjoint of the ∂-operator and H n,q ε,δ (E) is the space of harmonic forms on Y ε , that is,

H n,q ε,δ (E) := {w ∈ L n,q
(2) (E) ε,δ | ∂w = 0 and ∂ * ε,δ w = 0.}. Proposition 3.4 (resp. Proposition 3.5) can be proved by the same method as in [FM16, Proposition 5.4, 5.6, 5.7] (resp. [FM16, Proposition 5.9, 5.10]), so we omit the proofs here. for every relatively compact set K ⋐ X \ D, then the cohomology class β is zero in H q (X, K X ⊗ E ⊗ I(H)). Here • K,hε,ω ε,δ denotes the L 2 -norm on K with respect to h ε (not H ε ) and ω ε,δ .

Proposition 3.5. There exists v ε,δ ∈ L n,q-1

(2) (E) hε,ω ε,δ satisfying the following properties : ∂v ε,δ = u ε,δ and lim δ→0 v ε,δ ε,δ is bounded by a constant independent of ε. (31) Remark 3.6. In general, we have

L n,• (2) (E) ε,δ = L n,• (2) (E) Hε,ω ε,δ L n,• (2)
(E) hε,ω ε,δ , and thus v ε,δ may not be L 2 -integrable with respect to H ε .

For the above solution v ε,δ of the ∂-equation, by using the density lemma, we can take a family of smooth E-valued forms {v ε,δ,k } ∞ k=1 with the following properties :

v ε,δ,k → v ε,δ and ∂v ε,δ,k → ∂v ε,δ = u ε,δ in L n,• (2) (E) hε,ω ε,δ . (32)
Now we consider the level set

X c := {x ∈ X | -|s| 2 b < c} ⋐ X \ D for a negative number c. The set of the critical values of |s| 2
b is of Lebesgue measure zero from Sard's theorem. Hence, for a given relatively compact K ⋐ X \ D, we can choose -1 ≪ c < 0 such that

K ⋐ X c := {x ∈ X | -|s| 2
b < c} and d|s| 2 b = 0 at every point in ∂X c . Then, by [Mat16b, Proposition 2.5, Remark 2.6] (see also [START_REF] Folland | The Neumann problem for the Cauchy-Riemann complex[END_REF](1.3.2) Proposition]), we obtain and * denotes the Hodge star operator with respect to ω ε,δ . Note that dV ε,δ = dS ε,δ ∧ d|s| 2 b . One can easily see that lim

∂v ε,δ,k , u ε,δ X d ,hε,ω ε,δ = v ε,δ,k , ∂ * hε,ω ε,δ u ε,δ X d ,hε,ω ε,δ -((v ε,δ,k , (∂|s| 2 b ) * u ε,δ )) ∂X d ,hε,ω ε,δ (33
k→∞ ∂v ε,δ,k , u ε,δ X d ,hε,ω ε,δ = ∂v ε,δ , u ε,δ X d ,hε,ω ε,δ = u ε,δ , u ε,δ X d ,hε,ω ε,δ
by (32), and thus it is sufficient to show that the right hand side of equality (33) converges to zero. For this purpose, we first prove the following proposition.

η(Θ hε (E) + γ) -i∂∂η ≥ η -2 i∂η ∧ ∂η -εω ≥ η -2 i∂η ∧ ∂η -εω ε,δ .
Since u ε,δ is harmonic with respect to H ε and ω ε,δ , we have ∂ * ε,δ u ε,δ = 0 and ∂u ε,δ = 0. Further we have γ = i∂∂ψ on X \ D. Therefore we obtain 0

≥ - √ ηD ′ * u ε,δ 2 ε,δ = √ η∂u ε,δ 2 ε,δ + √ η∂ * ε,δ u ε,δ 2 ε,δ - √ ηD ′ * u ε,δ 2 ε,δ = ηΘ Hε -i∂∂η Λu ε,δ , u ε,δ ε,δ + 2Re ∂η ∧ ∂ * ε,δ u ε,δ , u ε,δ ε,δ = ηΘ Hε -i∂∂η Λu ε,δ , u ε,δ ε,δ ≥ η -2 i∂η ∧ ∂η)Λu ε,δ , u ε,δ ε,δ -εq u ε,δ 2 ε,δ .
from the twisted Bochner-Kodaira-Nakano identity (see [Ohs04, Lemma 2.1] or [Fuj13, Proposition 2.20. 2.21]). On the other hand, one can easily check that

η -2 i∂η ∧ ∂η)Λu ε,δ , u ε,δ ε,δ = η -1 (∂η) * u ε,δ 2 ε,δ = η -1 * ∂η * u ε,δ 2 ε,δ , ∂η = -χ ′ (σ)∂σ = - χ ′ (σ) (|s| 2 b + c) ∂|s| 2 b .
By the above arguments, we conclude that

εq u ε,δ 2 ε,δ ≥ χ ′ (σ) η(|s| 2 b + c) (∂|s| 2 b ) * u ε,δ 2 ε,δ .
It follows that the left hand side converges to zero from u H,ω ≥ u ε,δ ε,δ . Further the function χ

′ (σ)/η(|s| 2 b + c) is bounded below since we have 1 5 > |s| 2 b , C > η, χ ′ (σ) > 1
for some constant C. This completes the proof.

Finally we prove the following proposition by using Proposition 3.7. = v ε,δ K,hε,ω ε,δ ∂ * hε,ω ε,δ u ε,δ K,hε,ω ε,δ . Since lim δ→0 v ε,δ K,hε,ω ε,δ can be bounded by a constant that is independent of ε, it is sufficient to show that lim ε→0 lim δ→0 ∂ * hε,ω ε,δ u ε,δ K,hε,ω ε,δ = 0. We have e -ψ/2 = 1/|s| b < C K on K for some constant C K > 0, since K is a relatively compact set in X \ D. Hence we see that Remark 3.9. In the case of a holomorphically convex manifold, a proof based on injectivity theorems can be obtained by a slight modification of the above proof. The only problem is that an E-valued differential form u representing a given cohomology class is not necessarily L 2 -integrable but just locally L 2 -integrable. Since X admits a holomorphic map π : X → S to a Stein space S, the form u is L 2 -integrable with respect to the metric he -ψ e -Φ for a suitable psh exhaustion function Φ on X. Then it is not hard to check that our arguments still work by replacing h with he -Φ .

Remark 3.10. It would be interesting to know whether the hypothesis that ψ has analytic singularities is really needed. The main statement still makes sense when ψ has arbitrary analytic singularities, and one may thus guess that the result can be extended by performing a further regularization of ψ.

Proposition 3. 4 .

 4 If we have limε→0 lim δ→0 u ε,δ K,hε,ω ε,δ = 0,

  ) for almost all d ∈ ]c -a, c + a[, where a is a sufficiently small positive number. Here ∂ * hε,ω ε,δ is the formal adjoint of the ∂-operator in L n,• (2) (E) hε,ω ε,δ and ((•, •)) ∂X d ,hε,ω ε,δ is the inner product on the boundary ∂X d defined by ((a, b)) ∂X d ,hε,ω ε,δ := ∂X d a, b hε,ω ε,δ dS ε,δ , for smooth E-valued forms a, b, where dS ε,δ denotes the volume form on ∂X d defined by dS ε,δ := - * d|s| 2 b / d|s| 2 b hε,ω ε.δ

  b ) * u ε,δ ε,δ = 0.Proof of Proposition 3.7. By property (d) and property (B), we have

Proposition 3. 8 .

 8 (i) For a relatively compact set K ⋐ X \ D, we havelim ε→0 lim δ→0 lim k→0 v ε,δ,k , ∂ * hε,δ u ε,δ K,hε,ω ε,δ = 0. (ii) For almost all d ∈ ]c -a, c + a[, we have lim ε,δ,k , (∂|s| 2 b ) * u ε,δ )) ∂X d ,hε,ω ε,δ = 0.Proof of Proposition 3.8. In general, we have the formula D ′ g u = G -1 ∂(Gu) for a smooth hermitian metric g, where G is a local function representing g. Let G ε be a local function representing h ε . We remark that G ε e -ψ is a local function representing H ε . By the definition of∂ * ε,δ = ∂ * Hε,ω ε,δ , we have 0 = ∂ * ε,δ u ε,δ = - * D ′ Hε * u ε,δ = - * (G ε e -ψ ) -1 ∂(G ε e -ψ * u ε,δ ),and thus we obtain∂ * hε,ω ε,δ u ε,δ = - * (G ε ) -1 ∂(e ψ G ε e -ψ * u ε,δ ) = - * ∂e ψ * u ε,δ e -ψ = -(∂|s| 2 b ) * u ε,δ e -ψ . Now we have lim k→∞ v ε,δ,k , ∂ * hε,ω ε,δ u ε,δ K,hε,ω ε,δ ≤ lim k→∞ v ε,δ,k K,hε,ω ε,δ ∂ * hε,ω ε,δ u ε,δ K,hε,ω ε,δ

  ∂ * hε,ω ε,δ u ε,δ K,hε,ω ε,δ = -(∂|s| 2 b ) * u ε,δ e -ψ K,hε,ω ε,δ ≤ C K (∂|s| 2 b ) * u ε,δ K,ε,δ .We obtain the first statement (i) since the right hand side converges to zero by Proposition 3.7. Now we prove statement (ii). By the Cauchy-Schwarz inequality, we have((v ε,δ,k , (∂|s| 2 b ) * u ε,δ )) ∂X d ,hε,ω ε,δ 2 ≤ ((v ε,δ,k , v ε,δ,k )) ∂X d ,hε,ω ε,δ (((∂|s| 2 b ) * u ε,δ , (∂|s| 2 b ) * u ε,δ )) ∂X d ,hε,ω ε,δ . By Fubini's theorem, we obtain d∈]c-a,c+a[ ((v ε,δ,k , v ε,δ,k )) ∂X d ,hε,ω ε,δ dS ε,δ = c-a<-|s| 2 b <c+a |v ε,δ,k | 2 hε,ω ε,δ dV ε,δ ≤ v ε,δ 2 hε,δ .By Fatou's lemma, we see thatd∈]c-a,c+a[ ε,δ,k , v ε,δ,k )) ∂X d ,hε,ω ε,δ dS ε,δ ≤ limTherefore the integrand of the left hand side is finite for almost all d ∈ (c -a, c + a). On the other hand, by the same argument, we see thatd∈]c-a,c+a[ 2 b ) * u ε,δ , (∂|s| 2 b ) * u ε,δ )) ∂X d ,hε,ω ε,δ dS ε,δ ≤ lim ε→0 lim δ→0 (∂|s| 2 b ) * u ε,δ hε,ω ε,δ = 0.Therefore the integrand of the left hand side is zero for almost all d ∈ (c -a, c + a). This completes the proof.Theorem 1.1 is now a consequence of Proposition 3.4, Proposition 3.8, and equation (33).

This work was supported by the Agence Nationale de la Recherche grant "Convergence de Gromov-Hausdorff en géométrie kählérienne" (ANR-GRACK), the European Research Council project "Algebraic and Kähler Geometry" (ERC-ALKAGE, grant No. 670846 from September 2015), the Japan Society for the Promotion of Science Grant-in-Aid for Young Scientists (B) (JSPS, grant No. 25800051) and the JSPS Program for Advancing Strategic International Networks to Accelerate the Circulation of Talented Researchers.


Xin Fang

Jihong Wen
email: wenjihong@vip.sina.com

Bernard Bonello

Jianfei Yin

Dianlong Yu

Wave propagation in one-dimensional nonlinear acoustic metamaterials

Keywords: nonlinear acoustic metamaterial, wave propagation, band structure, nonlinear dynamics, chaos

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Acoustic metamaterials [1][2][3][4][5][6] (AMs) are artificial media that gain their properties from structure rather than composition. They generally feature built in resonators with subwavelength dimensions, rendering the concepts of effective mass density and effective elastic constants relevant for characterizing them. Because of their ability to manipulate the propagation of sound, AMs have become a hot topic during the last decade [6]. By now, most studies [8][9][10][11][12] are dealing with linear AMs (LAMs) which are based on a locally resonant (LR) mechanism [3] to yield negative indexes, super-lensing effects [12], or wave guiding. Different from the Bragg mechanism [13], LR paves the way to manipulate low-frequency waves with small "meta-atoms" [12]. However, the LR bandgap is generally narrow [6] and heavy "atoms" are necessary to enlarge this bandgap.

The concept of nonlinear metamaterial was first introduced in electromagnetism to investigate the photonic metamaterials in 2003 [14]. Nonlinear responses can be constructed on purpose, such as tenability, electromagnetically induced transparency, plasmonics, active media, etc [15,16]. In contrast, it is only recently that nonlinear acoustic metamaterials (NAMs) have emerged.

The presence of nonlinear media in linear phononic crystals can be utilized to realize interesting devices, such as acoustic diodes [17,18]. Compared with NAMs, studies on one-dimensional (1D) nonlinear periodic structures [19] (NPSs) have a longer history that can date back to 1955 when the studies on the Fermi-Pasta-Ulam (FPU) models [START_REF] Fermi | Studies of Nonlinear Problems[END_REF][START_REF] Gallavotti | The Fermi-Pasta-Ulam Problem: A Status Report[END_REF] that promoted the nonlinear physics. These nonlinear periodic chains are suitable for modeling a number of physical systems arising in different scientific contexts [22,23]. Based on the FPU model, Nesterenko [START_REF] Nesterenko | Dynamics of Heterogeneous Materials[END_REF] experimentally highlighted solitons in the granular crystal (a kind of NPS) interacting in a nonlinear way through Hertzian contact [START_REF] Daraio | [END_REF]26]. Moreover, a bistable lattice (another NPS) is studied recently [27]. By now, may new physical properties, that are different from the linear phononic crystals, are found in NPSs, such as unidirectional transition [27], discrete breathers [23,28], waves coupling [29], subharmonic frequencies [30,31], soliton waves [32,33] and surface waves [34]. Both simulations [22,23] and experiments [35,36] demonstrate that bifurcations in NPSs would be relevant to those properties. Acoustic devices such as diodes [36] and lenses [37] can be built upon them. The strong NPS has been proved to increase the velocity of sound and therefore the acoustic impedance [38]. Moreover, the bandgap properties in NPSs attract much attention [39][40][41][42][43][44]. The perturbation approach and the harmonic balance method are adopted to [40][41][42] studied the amplitude-dependent dispersions, stop band properties, and wave beaming in granular crystals; experimental works highlighted the role played by the critical amplitude in energy transmission [43] and bifurcation-induced bandgap reconfiguration [44] in NPSs. Actually, the granular crystals are suitable for ultrasonic applications; it is hard to consider them at low frequency regime because of the high contact stiffness they inherently feature.

Because of their promising applications, AMs with both low-frequency and broadband properties attract much attention. However, the mechanisms for both properties simultaneously happen are difficult to realize.

LAMs consist of linear "meta-atom", but when the linear "meta-atom" gets nonlinear in NAM, properties of the wave propagation get different patterns. In our recent works [45][46], the wave propagation in diatomic and tetratomic NAMs are analyzed using the homotopy analysis method and we found that the chaotic bands resulting from bifurcations can significantly enlarge the width of the wave-suppressing bands. This finding demonstrates that chaos is a novel and promising mechanism to achieve simultaneously low-frequency and broad band in both mono-bandgap NAMs and multi-bandgap NAMs, and still that a strong nonlinearity is beneficial to expand the bandwidth by several times. However, there are lots of phenomena arising in NAMs that have not yet been fully explained nor demonstrated. For example, why are the responses in the first pass band similar to ones observed with LAMs? Under which conditions can the elastic energy propagate in the bandgap? When will the wave be amplified by chaos? In the tetratomic system, why the nonlinearity has larger influence on the nonlinear LR bandgaps than it has for the Bragg gaps? Furthermore, some problems remain still unexplored, for example, the influence of the nonlinearity on the structure of the bands, the features of chaos and their differences, the jumps, the nonlinear resonances and their stabilities and the influence of the damping on chaos.

In this paper we attempt to answer these questions with the help of the frequency response analysis, the bifurcation theories, the Lyapunov exponents (LEs) and the fractal dimensions.

Models

Diatomic model

The 1D damped diatomic NAM model is illustrated in figure 1. The nonlinear "meta-atom" of this periodic model consists of a linear oscillator m and a damped nonlinear Duffing oscillator m r with cubic nonlinear stiffness k 1 Δ+k 2 Δ 3 , where Δ denotes the relative displacement and k 1 (k 2 ) symbolizes the linear (nonlinear) stiffness. In addition, a linear viscous damping term c  is taken into consideration.

3 0 1 1 1 2 3 12 ( 2 ) ( ) ( ) ( ) ( ) ( ) ( ) n n n n n n n n n n n n n n n n n x x x x y x y x y x y y x y x y x                                 (1)
where the dot denotes the derivation with respect to time t. The definitions of the generalized parameters are: mass ratio η=m r /m, stiffness ratio Applying the Bloch theorem for periodic structure, the displacement reads ( ) ( )exp( ) r a r i a  uu κ , where κ is the wave vector and a is the lattice constant. We can solve the dispersion relations with perturbation approach, harmonic balance method or homotopy analysis method. The algorithms for these approaches are described in [46]. The dispersion solutions are eigen waves in the infinite metamaterial without external input.

β 0 =k 0 /m= 2 s  , β 1 =k 1 /m r = 2 r  , β 2 =k 2 /m r .
We have considered a diatomic NAM with a finite length: 8 unit cells are involved in the model. Actually, further increasing the number of cells hardly changes the intrinsic characteristics of the NAM but would decrease the accuracy of the bifurcation diagram. Therefore, there are 16 degrees-of-freedom (DoF) and the state space has 32 dimensions (32D, nD symbolizes n-dimensional). In practice, applications would imply the dynamic of waves in a finite system with external excitations. The monochromatic sinusoidal wave u 0 (t) drives the system from the left end. The right end is free and the terminal linear oscillator is analyzed.

There is only one LR bandgap in the diatomic model. In practice, besides the LR bandgap, there are different Bragg bandgaps in the metamaterials. The nonlinearity would give rise to different behaviors in the multiple bandgaps. Though the tri-atomic model presents the case having a LR bandgap and a Bragg bandgap, here we discuss a more complex tetratomic model to figure out the influences on multiple bandgaps.

Tetratomic model

The model of the tetratomic NAM is illustrated in figure 2. A "molecule" of this model is composed of three identical linear "atoms" and one nonlinear "atom". Their parameters are labeled in the figure, where x n , y n , z n , r n denote the displacements of corresponding oscillators in the nth cell, respectively. The tetratomic NAM has more complicated dynamics. The generalized motion equation of the nth cell can be simplified as

m m m k 0 k 0 k 0 k 0 m m m k 0 k 0 k 0 n n+1
3 0 1 1 2 0 01 3 12 ( 2 ) ( ) ( ) ( ) ( 2 ) ( 2 ) ( ) ( ) (
)

n n n n n n n n n n n n n n n n n n n n n n n n n x y z x r x r x r x y x z y z x y z r r x r x r x                                            (2)
The definitions of other parameters in Eq. ( 2) are identical with those in Eq. (1) The algorithms to analyze the bifurcation and characteristics of chaos in the tetratomic model are also identical with the diatomic model. The difference is the chain we considered here has only 4 cells, so the state space is also 32D. The parameters of the model are: , and Ω=ω/ω s . ω s is 3 times of ω rs , which is similar with the continuous metamaterials which requires soft resonators.

m=1, m r =2, η=2, β 0 =225= 2 s  , β 1 =25= 2 sr  , β 2 =5×10
A larger but still weak damping ζ=0.02 is adopted. Assuming that the driving amplitude is A 0 =0.01 and the nonlinear strength is σ=0.6, so it is strongly nonlinear system in this case.

Theories

Method in frequency domain

(1) Some Algebraic laws The analysis of the nonlinear dynamics of high dimensional systems generally involves high order operations on matrices and vectors. For convenience, let us define at first some algebraic operations.

We define the "element product" of two nD columns vectors

The element product of two vectors can be abbreviated as xy . Therefore, the notation for the mth power of

( ) ( ) [[ ]], [[ ]] m m m     A x y A x A x A x yx (5)
These algebraic rules allow for a convenient analysis of the high-dimensional nonlinear equations and lead to compact formulas.

(2) Frequency responses It is not an easy task to achieving the analytical solutions of Eqs. ( 1) and ( 2) for a high-dimensional system. However, the frequency responses can be approximated with the help of the numerical integral method, the averaging method, or the harmonic balance method (HBM). The latter can be implemented to find the approximate steady frequency responses. To this end, we define the coordinate transformation ˆn

where M, C, K and N are mass, damping, stiffness and nonlinear coefficient matrices of the whole transformed system; f stands for the node force vector applied on every masses. Let us assume that the steady response has the form cos sin tt

  y = a b
, where a and b are constant vectors. The first-order HBM procedure leads to a system of algebraic equations:

    2 2 2 2 2 2 [ ] 3 ( ) 4 [ ] 3 ( ) 4                  K M a Cb N a b a f K M b Ca N a b b 0 (7)
The solutions of Eq. ( 7) would accurately describe the responses 22  Y = a b to the driving force. However, HBM method does not directly allow investigating the stability of the solutions. This is the reason why we adopted in this work the harmonic average approach (HAA) [48,49] instead, to study the frequency responses.

Within this approach, the solution is assumed to have the form:

( ) cos ( )sin ( )sin ( ) cos tt tt           y = u v y = u v (8)
where t   . The derivatives with respect to time t of the formula in (8) are

                   y = u v v u y = v u u v (9)
Comparing the expressions of y in Eqs. ( 8) and ( 9), we obtain ( cos sin )

    M u v 0 (10)
The further substitution of ( 8) and ( 9) into (6) gives another form of the equation of motion:

    2 3 2 3 3 2 3 2 2 cos sin [ ]( cos sin ) ( ) cos ( )sin {( 3 ) cos ( 3 )sin 3 sin 3 cos } cos                              M v u K M u v C u v v u N u uv v u v u v uv = f (11)
Then, assuming that u and v are constant, by calculating (11)×sinθ- (10)×cosθ and integrating the result from 0 to 2π, we obtain

  1 2 2 2 2 ( ) 3 ( ) 4            u M K M v Cu N v u v (12)
Similarly, calculating (11)×cosθ+ (10)×sinθ and integrating over the interval [0, 2π] gives:

  1 2 2 2 2 ( ) 3 ( ) 4              v M K M u Cv N v u u f (13)
The steady solutions correspond to the condition  u v 0 and their expressions are therefore:

    2 2 2 2 2 2 ( ) 3 ( ) 4 0 ( ) 3 ( ) 4                  K M v Cu N v u v K M u Cv N v u u f (14)
6

The amplitude of the response is

22  Y u v .
The derivation of these solutions using HAA turns the steady response problem into an equilibrium problem of differential equations. Actually, the solutions in Eqs. ( 14) and ( 7) are equal. However, HAA allows for analyzing the stability of the solutions of both the systems of differential equations ( 12) and ( 13) via the Jacobian matrix. Omitting the terms in 2ω≠0 in the left sides of ( 12) and ( 13) since they will not influence the properties of the solutions, the Jacobian matrix J is derived by performing the derivations of ( 12) and ( 13) with respect to the vectors u and v, namely

1 2 2 2 1 2 2 2 3 [[ ]] 2 3 [[3 ]] 4 3 [[3 ]] 4 3 [[ ]] 2                            M 0 C N uv K M N v u J 0 M K M N u v C N uv (15)
Therefore, with the solutions coming from Eq. ( 15), their stability can be determined: if the real part of an eigenvalue of J is positive, the corresponding steady solution is unstable; if a real eigenvalue goes from negative to positive, saddle-node (SN) bifurcation occurs, and if a pair of conjugate complex eigenvalues crosses the imaginary axis, Hopf bifurcation appears [START_REF] Yuri | Elements of Applied Bifurcation Theory[END_REF]. The succinct expressions derived above are universal for the NAMs and differential equations with cubic nonlinearity.

The pseudo-arclength continuation algorithm is further adopted to solve Eq. ( 14). Within this frame, Newton method (MATLAB fsolve) is used to find the solutions. Let us notice that if C=0, this formula reveals that v must read v=0.

Bifurcation analysis method in time domain

From now on, the periodic solutions will be mostly derived using time domain methods, which are parts of software such as AUTO [START_REF] Doedel | Continuation and bifurcation software for ordinary differential equations[END_REF]. Within this kind of approach, the integrations of the equations of motion are made in the state space. Furthermore, the method deriving LEs [START_REF] Wolf | [END_REF]53] is also based on numerical integrations. The incident elastic wave transforms the model in a non-autonomous system. However, when calculating the spectra of LEs and bifurcation diagrams, one must transform the non-autonomous system to an autonomous one. To calculate the former, the input wave is expressed as 00 sin , , (0) 0 uA

       (16)
i.e. it becomes a boundary value problem. Therefore, it becomes a 33D system in this situation.

When calculating the bifurcation diagrams, a differential system is added.

                          (17)
This system has a unique asymptotically stable solution, sin , cos tt

     .
The boundary value at a specified time t 0 is defined as 0

0 0 0 sin , cos tt      . Therefore, 00
uA   and the whole system becomes 34D.

We used the program AUTO [START_REF] Doedel | Continuation and bifurcation software for ordinary differential equations[END_REF] to analyze the bifurcations of periodic solutions. AUTO is based on a Newton iterative scheme to find the solutions. The continuation method is then used to find the branches of solutions. The bifurcations and stabilities of periodic solutions are identified by Floquet multipliers μ i . Let us recall here that, for a continuous differential system, there are three types of codimension-1 bifurcations:

saddle-node (SN), period doubling (PD) and torus (TR) bifurcations [START_REF] Strogatz | Nonlinear dynamics and chaos[END_REF]. The bifurcation associated with the appearance of μ i =1 (μ i =-1) is called a SN (PD) bifurcation; the bifurcation corresponding to the presence of conjugate multipliers μ 1,2 =exp(±iθ 0 ) , 0<θ 0 <π, is called a torus bifurcation. An invariant torus is a quasiperiodic solution. Both PD bifurcations and breakdowns of invariant tori may induce chaos [START_REF] Yuri | Elements of Applied Bifurcation Theory[END_REF]. At present, most studies on the bifurcation are focusing on the low-dimensional systems. In our 34D model, there are multiple pairs of conjugate multipliers.

Moreover, the QR decomposition algorithm proposed by wolf et al. [START_REF] Wolf | [END_REF] is adopted to compute the spectra of LEs of the 33D system. There are 33 LEs: λ 1 , λ 2 ,…λ 33 (in descending order). λ 1 is named as the largest LE (LLE). All non-autonomous systems have at least one zero LE that corresponds to the t-component. Therefore, if λ 1 =0, we discard it and make λ n * =λ n+1 be the new spectrum. If λ 1 >0, the motion is Lyapunov chaotic, and if there are more than two positive λ i , the motion is described as Lyapunov hyperchaotic.

Let the quantity

i s     
be the sum of all the positive LEs. The strength of chaos is determined by the average ratio of the exponential divergence of neighbor orbits in chaotic attractor. Therefore, s  may be used to describe the strength of chaos and the larger s  , the stronger the chaos is. Actually, noting that a larger λ 1 corresponds to a larger s  , it can be more conveniently used to characterize the chaos. Accordingly, if 1   →0, it is weak chaos which would behave as a quasiperiodic motion.

Finally, the Lyapunov dimension d LD is defined as

1 1 j LD i j i dj     
, where j is the maximum value of i that makes 0 and Ω c is moved to 2.277, as shown in figure 3(a) [45].

i    , i.e.
Actually, the pass bands of the finite LAM are composed of discrete resonances that become denser when increasing the periodic numbers. As shown in figure 3(a,b), although there is no damping, the appearance of nonlinearities significantly reduces and even suppresses resonances. However, there are some differences between the numerical integral results and HAA results obtained with Newton method, especially in the second passband (i.e. the optical band). There are several reasons for that: (1) there is an intrinsic randomness in the chaotic responses, so A max manifests the extreme case but this does not mean that the system gets a steady amplitude A max ; (2) the solution of HAA is an approximately averaged result at a specified frequency and a more accurate result should require high-order harmonic components. Despite all this, HAA well depicts the nonlinear responses and the role played by the nonlinearities as well.

Different from the linear resonance with infinite amplitudes, the nonlinear modes have finite phase volumes so that the amplitudes of the 16 resonances keep bounded [45]. Increasing the amplitude reduces the wave transmissibility, and meanwhile the nonlinear upshot of suppressing resonances extends to the first passband (i.e. the acoustic band). This result is further demonstrated by numerical integrations. resonances, the stable steady solutions become unstable when the frequency exceeds a threshold. Their amplitudes increase then along a branch until a stable solution appears at a certain frequency, where the amplitudes of the stable solutions switch to a lower branch. Subsequently, the amplitude increases along this branch and a "neighbor peak" with lower amplitude germinates.

When the driving amplitude increases to A 0 =0.005, the stable and unstable solutions appear alternately. As a consequence the high amplitudes are unstable and the total frequency width of the unstable regions is broader than the width of the stable regions. However, the "neighbor peak" disappears so a curve of nonlinear resonance becomes a single trajectory bending to the right. These results illustrate the fact that, near a resonance, a high dimensional NAM has a behavior different from that of a single DoF or even of a 2DoF system [START_REF] Nayfeh | Nonlinear oscillations[END_REF]. Although it is difficult to determine whether the amplitude may generate a "jump" at the catastrophe frequency, the broadband unstable solutions provide an opportunity for the chaos to appear. In contrast to the steady amplitude, the damping has a significant impact on the characteristics of the chaotic attractors. As shown in figure 4, in the chaotic region the LLE remarkably decreases upon damping as weak as ζ=0.01. Therefore, a strong damping effect could allow a chaotic solution to turn into a periodic solution.

Hereafter, we consider a weak damping ζ=0.01 in the diatomic system.

We have chosen three illustrative frequencies Ω=0.678, 1.124, Ω=2.2, whose locations are labeled in figure 3(a). Bifurcation diagrams and chaotic characteristics are illustrated in figure 5 for each of these frequencies.

For the conservative 34D system, there are dense TR points along the branches of periodic solutions because of the multiple pairs of conjugate multipliers, which indicates that there are extensive quasiperiodic solutions in the system. The damping does not alter the branches and stabilities of periodic solutions in the bifurcation diagram, because the linear damping x  does not change the properties of the equilibriums in the time-domain algorithm (the Jacobian matrix from state space). However, the phase volumes of the damped systems shrinks in the long-term responses, which results in that the Floquet multipliers on the unit circle enter the circle.

Therefore, the damping ζ=0.01 reduces the number of TR points along the branches, and the PD points found in the conservative system disappear. In some frequency range, the TR point disappears too, as illustrated by figure 5(c). In the main, the damping reduces the number of quasi-periodic solutions in the system.

The unstable periodic solution is an essential condition for the chaos to appear. The bifurcation diagram, mean amplitude that directly comes out from numerical integration, LEs and d LD in figure 5 are well consistent with each other.

In the case Ω=0.678 the first unstable solution appears at A 0 =1.867×10 -3 marked by the vertical black line in figure 5. Within the subsequent domain, the stable and unstable periodic solutions appear alternately. When A 0 >6.686×10 -3 , the periodic solution enters into a constantly unstable region, in which multiple branches are generated with many TR points. The terminal points of the new branches are SN points. In the interval 1×10 -3 <A 0 <1.867×10 -3 , 1   →0 and therefore the response is close to the quasiperiodic responses. Both behaviors of LEs and d LD attest that the periodic motion turns into chaos when the driving amplitude still increases. Actually, both 1   and d LD tend to increase as A 0 increases, meaning that the weak chaos becomes a strong chaos. However, LEs and d LD have local fluctuations. For the amplitude at which chaos appears, only one LE is positive, but there are quickly more than two positive LEs, i.e. hyperchaos occurs. Therefore, for the 33D NAM model, this is a general rule that the system is hyperchaotic when chaos occurs. Moreover, the variation laws of d LD are identical with λ 1 in the chaotic domain, and when strong chaos occurs d LD →N d =33. Combing the bifurcation diagram and LE spectrum, it is known that chaos at TR points would be induced by the breakdowns of invariant tori. A rigorous demonstration of this mechanism is based on the KAM theorem [START_REF] Yuri | Elements of Applied Bifurcation Theory[END_REF]. Such a demonstration falls out of the scope of this article, which rather focuses on the influence of chaos on the responses. At other unstable points, the chaos arises from period-doubling bifurcations. The waterfall plot of power spectra illustrates this process, as shown in figure 6. When A 0 →0, we are dealing with a 1:1 resonance. In this situation, the energy gets localized at the driving frequency Ω e . Further increasing A 0 leads to odd harmonic waves 3Ω e and then 5Ω e to appear, so the propagating waves become quasiperiodic. At the critical amplitude A 0c , the propagating wave cascades into chaos. A 0c has the same behavior as LEs. In the route toward chaos, energy localization switches to energy dispersion and therefore the energy is pumped [START_REF] Fang | [END_REF] and spread within a broad high frequency passband and even within the stop bands. If the NAM is excited by the wave at the bottom of the LR bandgap, the periodic solution exhibits multiple folded branches that induce multiple jumps, as shown by the case Ω=1.124. There are still some discrete TR points on the unstable branches. The SN bifurcation points act as the terminals of branches. The first SN bifurcation occurs at A 0 =3.807×10 -3 . When A 0 <3.807×10 -3 , the periodic solutions remain stable and their amplitude tends to zero; both LEs and d LD reveal that the long-term motion of the system is periodic. The amplitudes issued from numerical integration also tend to zero. This result means that the LR bandgap remains a complete stop band and that the bandgap is linearly stable in this interval. After this SN point, the periodic solutions become constantly unstable and a jump occurs at this point. Both the bifurcation diagram and the mean amplitude show that the motion jumps to a high-energy chaotic branch at this SN point. Then, as illustrated both by LEs and by d LD , the propagation through the NAM is along a chaotic way. Meanwhile, the hardened nonlinear stiffness shifts upwards the lower boundary of the LR bandgap. Along this chaotic branch, the mean amplitude has almost constant value with LEs and d LD fluctuating in a neighborhood of the corresponding constants. Another SN bifurcation point appears at A 0 =7.676×10 -3 , where the motion jumps to a higher chaotic orbit. However, this point does not fit well the long-term motion, as attested by the jumps of the mean amplitude to a higher-energy chaotic orbit at A 0 =0.013. Like the mean amplitude, LEs and d LD keep almost constant. These results show that the system has multiple approximately constant states in this frequency range, and that only discrete values of these states can be achieved: the system behaves as a "quantum object", as sketched in figure 7. Therefore, when manipulating waves in the bandgaps via amplitudes, the system can switch suddenly between the stop state (i.e. the bounded state) and the propagation states (i.e. the excited quantum states). These characteristics may be helpful to realize acoustic devices with small dimensions, as for example acoustic diodes or switches. In figure 5(b), it should be noted that, at some value of A 0 , the amplitude may jump from the highest obit (the third one) to a lower one. However, these jumps are unstable and a small perturbation may stimulate the jump back to the highest orbit.

The case Ω=2.2 is similar to the case Ω=1.124 except that the motion turns from periodic into chaotic for only one jump at the first unstable periodic solution. The trend in bifurcation diagram agrees well with the mean amplitude variations. After the amplitude at first slowly increases and then keeps constant as do both LEs and d LD . This suggests a mechanism to suppress elastic waves in strongly NAM: strengthening the nonlinearity decreases the wave transmissibility. More bifurcation diagrams at other frequencies are depicted in figure 8. According to the statements above, it can be anticipated that the chaos would occur at the first unstable periodic solutions. In the examined amplitude range, neither bifurcations nor unstable solutions occur on the low-frequency periodic orbits (Ω≤0.3) and the amplitude linearly increases with driving amplitude, i.e. the transmissibility is constant and the linear regime is preserved. At higher frequency the motion turns to weak chaos that behaves as a quasiperiodic wave, which explains why the responses of NAM in the low-frequency parts of the acoustic branch are similar to the ones of LAM. The transmission coefficient shows this feature (see above and figure 3). In the acoustic branch, the strong chaos can be observed at the frequencies near the LR bandgap.

Bifurcation diagrams in the optical branch have similar characteristics: firstly, the amplitudes increase along an orbit where both periodic and chaotic behaviors may occur. However, when the motions jump to a higher orbit at SN point, the amplitudes of the periodic solutions do not increase with the input wave. Instead, they remain constant or even they decrease.

In the bandgap, the periodic orbits at Ω=1.25 feature folding and multiple jumps, as it is the case for Ω=1.124. However, larger critical amplitude is required to force the wave into chaos. The reason is that the influences of the nonlinearity on the bandgap spread from low to high frequencies. This is further demonstrated by the bifurcation diagram at Ω=1.4: although there are unstable solutions, the amplitudes are so small that only subharmonic and superharmonic waves can propagate and this frequency interval keeps its ability to reflect the incident waves. In the high-frequency stop band, as shown by the bifurcation diagram at Ω=2.4, much stronger nonlinearity is needed to generate unstable solutions.

Based on all these results, the band structure of a strongly diatomic NAM is sketched in figure 9.

Tetratomic model

The dispersion relations of this model were elaborated in Ref. [46]. In the linear regime, there are three As it is the case with the diatomic model, the tetratomic NAM also proves that the broadband elastic waves are suppressed. The results from HAA establish that a strong nonlinearity causes the wave transmissibility to decrease. However, the elastic waves can still propagate in the nonlinear LR bandgap. The continuation method emphasizes the surprising nonlinear modes that compose the unstable branch and a neighbor peak. Overall, the behaviors are similar to that of diatomic NAM. However, in the amplitude range A 0 ≤0.01, the nonlinear modes with only one branch are not found. In the discrete models, in conjunction with the results of the diatomic NAM, we can conclude that the hardened NAM actually influences all the passbands that are higher than the "nonlinear LR bandgap", including a small region below and near this bandgap. We have further studied bifurcation and chaos to better understand the wave propagation properties in the system. To this end, we firstly have investigated the wave propagation in the pass bands. We have considered three frequencies Ω=0.2463, 0.672, 1.42 located in the 1 st , 2 nd , 3 rd pass bands respectively, (see figure 10). The properties of the fourth passband are similar with the third one. The bifurcation diagrams, mean amplitudes, LEs and LDs are shown in figure 11(a-c). Comparing these results, it can be seen that increasing the driving amplitude stimulates the motion into chaos. Moreover, there are differences between the chaotic critical points in bifurcation diagrams and LEs which are caused by: (1) the damping that greatly influences the long-term motion;

(2) the different principles implemented by the two algorithms that lead to different critical points.

In the case Ω=0.2463, the bifurcation diagram is consistent with the mean amplitude: between 0.006<A 0 <0.009 the general chaos (only one positive LE) arises, but the motion turns into hyperchaos when A 0 further increases. However, 1   <0.03 and d LD <<33, so that the wave undergoes a weak hyperchaos whose behavior is similar to the quasiperiodic one. Therefore, it can be deduced that the low-frequency waves in the first passband are quasiperiodic or weakly chaotic.

In fact, the propagation of the waves in the subsequent three passbands is similar except for some local differences. Along with the increasing driving amplitude, the variation law of the wave state in whole is "periodic → weakly chaotic → strongly hyper-chaotic". However, in these three passbands, the LLEs λ 1 for the strong chaos are 10 times larger than that calculated for the weak chaos in the first passband and d LD →33. The mechanisms that suppress resonances in the passbands are identical to the ones observed with the diatomic NAM: chaos induced by period-doubling bifurcation; the wave amplitudes along the chaotic branches slowly increase, remain constant or even decrease with A 0 . Moreover, as illustrated in figure 11(b,c), periodic windows in the chaotic region are observed, which means that LEs and LDs do not monotonously vary. Close to Ω=1.42, T A is equal to the transmission coefficient for the corresponding LAM. In fact, there are other specific frequency domains. The comparison with the results obtained with the 10-cells model [46] indicates that these domains relate to the density of the linear modes that is determined by the length of chain.

In the 4-cells model, this domain is a non-resonant region. The periodic solutions for Ω=1.42 have two branches: a low-energy one and a high-energy one. Along the low-energy branch, there is an unstable domain in between two stable domains. If the motion had been periodic for A 0 ≤0.01, the mean displacement of the last oscillator should have been linear against A 0 . Instead, both the mean amplitude and LEs indicate that the system enters a weak chaos regime in the range 0.006<A 0 <0.009. However, this weakly chaotic motion is different from that described above. Here, the oscillators jump to a high-energy unstable orbit. In the subsequent range 0.009<A 0 ≤0.011, the system jumps back to the low-energy weakly chaotic state, featuring a quasiperiodic behavior. Therefore, when the amplitude is A 0 =0.01, the linear and nonlinear transmissibilities near Ω=1.42 are equal. Further increasing the driving amplitude makes the system to jump to the strongly hyper-chaotic state, featuring constant amplitude and high-energy, except in a narrow window around A 0 =0.014 where the motion gets periodic. . However, in this small interval, amplitudes still are in the high-energy orbit.

The jumps between low-energy and high-energy orbits explains why the waves are amplified under non-resonant conditions. This amplification scheme can be applied to design a broadband wave amplifier. For finite LAM, the waves can be amplified at the discrete resonant frequencies in the passbands. In contrast, a four-cell NAM can allow for a homogeneously amplification in a broad passband. Another question we have addressed is: why is the nonlinear LR bandgap less efficient than other bandgaps to suppress the elastic waves in NAM? To answer this question, we have set the frequencies to Ω=0.35 and Ω=1. As shown in figure 12, as for the diatomic model, there are multiple bifurcations and branches. However, the stable periodic solutions disappear for the amplitude as low as A 0 =0.003 because of the strong nonlinearity and therefore this bandgap closes at a certain driving amplitude and is replaced by a high energy chaotic orbit, as it is the case in the diatomic model. In contrast, for the BG bandgaps away from this nonlinear bandgap, much higher amplitudes are needed to get the unstable periodic solution. In the example illustrated in figure 11(d), when Ω=1 in the first BG bandgap, a periodic motion is observed and the chaos is weak ( 1   <0.04) if the amplitude A 0 <0.01, as for the quasi-linear state in the interval A 0 =0.01~0.015. It is only as A 0 >0.015, that the motion is attracted to the high energy hyper-chaotic state. Combing the mechanisms in diatomic model, one concludes that the nonlinear LR bandgap has weaker ability to reflect the incident waves when it jumps to the high-energy orbits. Based on the results from this paper and Ref. [46], and combing with some predications, we sketch the band structure of the tetratomic NAM in figure 13. For the LAM here, the low-frequency LR bandgap is narrow. The strong nonlinearity shifts the nonlinear LR upwards, meanwhile, the responses in this gap jumps to the excited state. Below the nonlinear gap, the band structures are similar with the diatomic ones. However, the second and the third linear resonant passbands are replaced by the strong hyperchaotic passbands with lower transmission.

When the nonlinearity is weak, BG1, BG2 and the fourth resonant passband would remain the linear or quasi-linear state; but when the nonlinearity goes to much stronger, BG1 and BG2 would also get multi-state behaviors, meanwhile, the fourth passband become chaotic passband and its upper boundary spreads to the high-frequency stop band. The detailed statements about the band structures should base on the structure, the frequency and the amplitude.

Conclusions

The propagation of waves in the nonlinear acoustic metamaterials (NAMs) is fundamentally different from that in the conventional LAMs. However, those features are still not fully understood. In this work we investigate the elastic wave propagations in the 1D diatomic and tetratomic NAM models. We further demonstrate that nonlinear effects can greatly suppress elastic waves in a broad band. This band consists of bandgaps and chaotic pass-bands. The nonlinear wave behaviors, band structures, the bifurcations and the chaos are studied to demonstrate the novel mechanisms that can manipulate wave propagations.

HAA combined with the pseudo-arclength continuation are employed to calculate the frequency responses and nonlinear modes. Bifurcations and chaos are analyzed using both the continuation algorithm and the spectra of LEs (and d LD ). Our results show that the nonlinear resonances feature multiple branches with unstable peaks and a "neighbor peak". We further have demonstrated that the period doubling process leads the motion from being periodic to chaos and damping has a significant influence on the characteristics of chaotic attractors. Moreover, due to dispersion process in chaotic regime, the localized energy spreads in the broadband high frequency passbands and even in the stop bands. The trend for the positive LEs and d LD is to increase with the amplitude so the weak chaos (or hyperchaos) turns into strong hyperchaos.

Band structures of both diatomic and tetratomic NAMs are studied. The bifurcations and unstable solutions do not occur on the low-frequency periodic orbits, so the responses in these regions are similar to the ones of LAMs. Actually, the hardened NAM influences all the passbands that are higher than the nonlinear LR bandgap, including a small region below this bandgap.

In the nonlinear LR bandgap, we put into evidence "quantum" behaviors because of jumping bifurcations between low-energy and high-energy orbits, whose propagation states (excited states) and stop state (bounded state) have discrete characteristics and switch suddenly. This behavior also explained why the nonlinear LR gets weaker ability to reflect the incident waves. Moreover, jumps in passbands would amplify non-resonant waves.

This work provides a theoretical base for the understandings and applications of NAMs.

Figure 1 .

 1 Figure 1. Model of the diatomic NAM

Figure 2 .

 2 Figure 2. Model of the tetratomic NAM

. 4 .

 4 For a high-dimensional system (>2D), d LD can accurately describe the fractal nature of the chaotic attractor. For a chaotic attractor, 0<d LD <N d , where N d is the dimension of the system. For periodic and quasiperiodic motions, d LD =0. A larger d LD manifests a more complex chaos. Properties of wave propagation in the NAM models4.1. Diatomic modelBandgap properties of the diatomic NAM model have been laid out in Ref.[45]. It is shown that for this system, the generalized frequency range of the LR bandgap is [1.088, 1.5] , and the cutoff frequency of the pass band is Ω c =2.251. The transmission coefficient defined as T A =A max /A 0 , where A max denotes the maximum amplitude of the output wave, i.e., the maximum amplitude of the last linear oscillator in the chain. A max represents the extreme condition in the long-term chaotic responses. T A can be derived via three different approaches as shown in figure3. The hardened nonlinear stiffness leads the lower bound frequency of the LR bandgap to shift upward. On the other hand, when the driving amplitude is set to A 0 =0.005, the bandgap calculated by homotopy analysis method is [1.2, 1.758] 

Figure 3 .

 3 Figure 3. The frequency responses of diatomic LAM and NAMs. (a) maximum responses solved directly with numericalintegral method; here A 0 =0.005, ζ=0, "L" and "N" symbolize "linear" and "nonlinear" models, respectively; the chain has 4 or 8 cells, where "nc" means there are n cells in the chain; (b) the solutions found with Newton iterative method based on Eq. (14); both the model without damping (ζ=0) and with a weak damping ζ=0.01 are considered. (c) the response of the damped NAM in the optical branch; the solutions are solved with pseudo-arclength continuation algorithm, and the stabilities are determined with Jacobian matrix in Eq.(15); the stabilities are distinguished with different colors.

Figure 4 .

 4 Figure 4. Influence of damping on the largest Lyapunov exponents (Ω=0.678).

Figure 5 .

 5 Figure 5. Bifurcation diagrams and chaos of diatomic model. (a) Ω=0.678, (b) Ω=1.124, (c) Ω=2.2. ζ=0.01. In each figure,the four diagrams are corresponding to the bifurcation diagram, the mean amplitude of the last linear oscillator, the Lyapunov exponents (λ 1 , λ 3 , λ 5 ) and the Lyapunov dimension d LD , respectively. In bifurcation diagrams, the solid black (dashed green) lines represent the stable (unstable) periodic solutions, and the red dots (blue asterisks) denote the SN (TR)

Figure 6 .

 6 Figure 6. 3D waterfall plot of power spectrum density (PSD). The excitation frequency is Ω e =0.678. The different colors correspond to the value log 10 (PSD).

Figure 7 .

 7 Figure 7. The sketch plot of quantum states in NAM. The solid circles represent the excited quantum states and the dashed curves denote the fluctuations in the neighborhood of these states. The shade circle and the black point in the center represent the bounded quantum state.

Figure 8 .

 8 Figure 8. Bifurcation diagrams of the diatomic model under different frequencies. The frequencies are labeled in the corresponding plots. The locations of these frequencies are marked in figure 3(a). Other specifications are identical with figure 5.

Figure 9 .

 9 Figure 9. The sketch of the band structure of diatomic NAM model. The green lines represent the transmission of NAM.

  bandgaps for this metamaterial: a LR bandgap and two Bragg (BG) bandgaps. The frequency ranges of the linear bandgaps are: LR, [0.3065, 0.4195]; BG1, [1, 1.078]; BG2, [1.732, 1.777]. It is shown in Ref. [46] that, in the nonlinear regime, BG1 and BG2 are not altered if the amplitude A 0 =0.01, but the LR bandgap is shifted to [0.42, 0.5807]. The cutoff frequency of the passband is Ω=2.

Figure 10 .

 10 Figure 10. The frequency responses of tetratomic LAM and NAMs. ζ=0.02 (a) T A and H(ω) solved directly with numerical integral method; here A 0 =0.01, "L" and "N" symbolize "linear" and "nonlinear" models, respectively; (b) the solutions are found with the Newton iterative method; (c) the solutions are solved with pseudo-arclength continuation algorithm; the stabilities are distinguished with different colors.

Figure 11 .

 11 Figure 11. Bifurcation diagrams and chaos of tetratomic model. (a) Ω=0.2463, (b) Ω=0.672, (c) Ω=1.42, (d) Ω= 1. Other specifications of these figures are same with figure 5. z 4 denotes the displacement of the last linear oscillator.

Figure 12 .

 12 Figure 12. Bifurcation diagrams at frequency Ω=0.35.

Figure 13 .

 13 Figure 13. The sketch of the band structure of tetratomic NAM model. The blue lines represent the transmission of strongly nonlinear AM model. The dashed blue lines correspond to the much stronger nonlinearity.

  The generalized frequency is Ω=ω/ω s and the damping coefficient is ζ=c/m r . The strength factor of nonlinearity is defined as σ=3β 2 A 0 2 /β 1 , where A 0 stands for the amplitude of the incident wave u 0 =A 0 sin(ωt) at the left end. It is a weak nonlinearity if σ <<1. The values of the parameters used in this model are: m=1, m r =0.5, β 0 =10π, β 1 =15π, β 2 =10 5 . If A 0 =0.005, σ=0.16, it is moderate nonlinearity. To increase σ, one generally needs to increase the drive because it is difficult to

increase k 2 in the macroscopic mechanical systems. Just Recently, Huang et al

[47]

used chemical bonding interactions to generate a huge k 2 on the size scale about 1 μm.

  4

y z xz yz , and associative  xyz zyx .

							x
	is	x	m		[] m T i x	. It can be easily established that the element product is commutative  xy yx , distributive
	()  x Similarly, we define the element product of two matrices as being  	AB		[] ij ij ab	, where the dot "。" cannot
	be omitted. It is straightforward to show the element product of two matrices obeys the commutative,
	distributive and associative laws as well. However, the product between a matrix and an element product of two
	vectors does generally not observe the associative law, i.e. ( A x y	) ( )  Ax y , which could become an issue
	in practice. To overcome this obstacle and to ensure that this operation gets associative, we set the following
	transformation
						( A x y	) ( [[ ]])  A x y	(4)
	where the notation [[ ]] x stands for the square matrix built upon the vector x, namely
						[[ ]] [  x	x x	x	T ]
							n
	With this definition, together with the commutativity, it comes ( [[ ]]) A x y A y x . Moreover, it can be ( [[ ]]) 
	easily shown that the partial derivative of the element product is
							1

Acknowledgements

This research was funded by the National Nature Science Foundation of China (Project Nos. 51405502 and 51275519).