text:
# 1 Introduction

## 1 Introduction

### 1.1 The model system

Gap solitons (GS) is a common name for solitary waves in nonlinear models which feature one or more  gaps in their linear spectrum, see review papers  and , respectively, for these two cases. A soliton may exist if its frequency belongs to the gap, as then it does not decay into linear waves.

Gaps in the linear spectrum are a generic phenomenon in two- or multi-component systems, as intersection of dispersion curves belonging to different components is generically prevented by a linear coupling between the components. Excluding cases when the system’s linear spectrum is unstable (which are possible in a fluid dynamics application ), the intersection avoidance alters the spectrum so that a gap opens in place of the intersection. Approximating the two dispersion curves, that would intersect in the absence of the coupling, by straight lines, and assuming a generic cubic \[$`\chi ^{(3)}`$\] nonlinearity, one arrives at a generalized massive Thirring model (GMTM), which has a family of exact GS solutions that completely fill the gap . The model has a direct application to nonlinear optics, describing co-propagation of forward- and backward-traveling electromagnetic waves in a fiber with a resonant Bragg grating (BG). Gap solitons, first predicted theoretically, were observed in experiments with light pulses launched into a short piece of the BG-equipped fiber  (in fact, optical solitons that were first observed in the BG fiber  were, strictly speaking, not of the GS type, but more general ones, whose central frequency did not belong to the fiber’s bandgap).

GS are known not only in optics but also in other physical settings, for instance, in density-stratified fluid flows, where dispersion curves pertaining to two different internal-wave modes often exhibit near-intersections. Again taking into regard weak nonlinearity, one can predict the occurrence of GS in density-stratified fluids .

In this work, we aim to consider GS that may exist in a generic situation of the next type, when the underlying system contains three wave components, and the corresponding dispersion curves intersect at three nearly-coincident points, unless linear coupling terms are taken into account. Situations of this type can readily occur in the above-mentioned density-stratified fluid flows, since tuning of two suitable external parameters can often lead to a near-coincidence in the linear phase speeds of three independent internal wave modes, for certain wavenumbers (see ). Indeed, similar considerations can be applied to many other physical systems. As for the usual GS systems, the generic wave-envelope model can be expected to contain cubic nonlinear terms.

In optics, a $`\chi ^{(3)}`$-nonlinear model with three linearly coupled waves is possible too, in terms of either temporary evolution of fields in a planar nonlinear waveguide equipped with a triple BG in the form of three systems of parallel scores, or spatial evolution of stationary fields in a bulk waveguide with a similar triple BG consisting of three systems of parallel interfaces between layers. The latter realization seems natural enough, as it strongly resembles photonic-crystal fibers, which have recently attracted a great deal of interest . Note that the former version of the model is a generalization of a three-wave model for a $`\chi ^{(2)}`$-nonlinear planar waveguide with an ordinary BG, which was introduced in . Both versions of the proposed model are illustrated in Fig. 1, where the periodic lattice shows the triple BG. In the case of the temporal evolution, Fig. 1 displays the planar waveguide, while in the case of the spatial evolution, it is a transverse cross-section of the bulk waveguide.

We stress that the lattice in Fig. 1 is not completely symmetric; although the triangular cells of the lattice are equilateral ones, the two diagonal sub-gratings are assumed, in the general case, to have the strength (contrast of the refractive index) smaller than the horizontal one. The bold triangles inscribed into the two triangular cells illustrate the resonant Bragg reflections that give rise to linear couplings between the waves. Then, neglecting intrinsic dispersion or diffraction of the waves in comparison with the strong artificial dispersion/diffraction induced by the Bragg reflections, normalized equations governing the spatial evolution of the three fields whose Poynting vectors are shown by three bold arrows in Fig. 1, are

$$i\left(\frac{u_1}{t}\frac{u_1}{x}\frac{1}{\sqrt{3}}\frac{u_1}{y}\right)+u_2+\kappa u_3+\left(|u_1|^2+2|u_2|^2+2|u_3|^2\right)u_1=0,$$
(1)
$$i\left(\frac{u_2}{t}+\frac{u_2}{x}\frac{1}{\sqrt{3}}\frac{u_2}{y}\right)+u_1+\kappa ^{}u_3+\left(|u_2|^2+2|u_1|^2+2|u_3|^2\right)u_2=0,$$
(2)
$$i\left(\frac{u_3}{t}+\frac{2}{\sqrt{3}}\frac{u_3}{y}\right)+\kappa ^{}u_1+\kappa u_2+\left(|u_3|^2+2|u_1|^2+2|u_2|^2\right)u_3=\omega _0u_3.$$
(3)
Here, the evolution variable $`t`$ is the proper time in the case of the temporal evolution in the planar waveguide, or the coordinate $`z`$ in the direction perpendicular to the plane of the figure in the case of the spatial evolution in the bulk waveguide. In the latter case, the beam enters the medium through the plane $`z=0`$ and evolves along the coordinate $`z`$, that is why it plays the role of the evolution variable in the equations (the initial conditions necessary to launch a soliton will be discussed in more detail below). The relative coefficients in front of the $`x`$\- and $`y`$\- derivative terms correspond to the geometry in Fig. 1, the coefficient of the walk-off in the $`x`$-direction in the first two equations being normalized to be $`\pm 1`$. The coefficient of the BG-induced linear conversion between the waves $`u_1`$ and $`u_2`$ is normalized to be $`1`$, while the parameter $`\kappa `$ (which is complex, in the most general case, but see the discussion below) accounts for the linear conversion between these waves and the third wave $`u_3`$, and the usual ratio $`1:2`$ between the coefficients of the self-phase modulation (SPM) and cross-phase modulation (XPM) is adopted. Lastly, $`\omega _0`$ is a frequency/wavenumber mismatch between the third and the first two waves, which is caused by the above-mentioned asymmetry between the diagonal and horizontal subgratings, as well as by other reasons.

As we mentioned above, the Bragg constant $`\kappa `$ in Eqs. (1) - (3), which couples the field $`u_3`$ to the pair $`u_{1,2}`$, is complex in the general case (note that the constant of the Bragg coupling between the fields $`u_1`$ and $`u_2`$ might also be complex in its primary form, and making it equal to $`1`$ in Eqs. (1) and (2) involves opposite constant phase shifts of the fields $`u_1`$ and $`u_2`$, which is why $`\kappa `$ and $`\kappa ^{}`$ appear exactly as in Eqs. (1) - (3)). However, assuming that each score, the families of which constitute the triple grating (in the case of the temporal-domain evolution in the planar waveguide), gives rise to simple reflection described by the classical Fresnel formulas, it is easy to conclude that all the coupling constants are real and positive, provided that either the light is polarized orthogonally to the waveguide’s plane, and the reflection takes place from a less optically dense material (i.e., the “score” is, literally, a shallow trough on the surface of the planar waveguide), or the light is polarized parallel to the plane of the waveguide, and the reflection is from a more optically dense material. Similarly, in the case when the same equations describe the spatial evolution of the optical fields in the layered bulk medium, one may assume that either the light is polarized in the $`z`$ -direction, and seams between the layers are filled with a material (for instance, air) which is optically less dense than the bulk medium, or the polarization is orthogonal to the $`z`$-axis (i.e., it is parallel to the plane of Fig. 1), and the material filling the inter-layer seams is optically denser than the host medium.

In the present paper, we focus on this case, which was described above in detail for the realizations of the model in terms of both planar and bulk optical waveguides, and which corresponds to $`\kappa `$ real and positive in Eqs. (1) - (3). Note, incidentally, that the case when $`\kappa `$ is real and negative can be reduced to the same case simply by reversing the sign in the definition of $`u_3`$.

The model displayed in Fig. 1 may be further generalized by introducing an additional asymmetry, which will remove the equality between the horizontal side of the lattice’s triangular cell and its diagonal sides. Then, the simultaneous fulfillment of the Bragg-reflection conditions for the waves $`u_{1,2}`$ and $`u_3`$ can be secured by making the waveguide anisotropic. However, such a generalization goes beyond the scope of this paper.

For the physical realization of the model, Eqs. (1) - (3) must be supplemented by initial conditions at $`t=0`$ in the case of the temporal evolution in the planar waveguide, or boundary conditions at $`z=0`$ in the case of the spatial evolution in the bulk medium. It is sufficient to assume that, at $`t=0`$, a single wave component (for instance, $`u_3`$) is launched into the waveguide. The linear-coupling terms in the equations will then start to generate the other components, and, if solitons that might exist in this model are stable (see below), they may self-trap from such an initial beam.

Bearing in mind also the above-mentioned application to internal waves in stratified fluids, as well as similar realizations in other physical media, Eqs. (1) - (3) may be naturally extended by introducing more general SPM and XPM coefficients, as in applications other than nonlinear optics, the ratios between the XPM and SPM coefficients may be different from those adopted above. Thus the generalized system of equations takes the following form, in which we confine consideration to $`y`$-independent solutions (the consideration of possible three-dimensional solitons in the case of $`y`$-dependent fields is not an objective of this work),

$$i(\frac{u_1}{t}\frac{u_1}{x})+u_2+\kappa u_3+\alpha \left(\alpha \sigma _1|u_1|^2+\alpha |u_2|^2+|u_3|^2\right)u_1=0,$$
(4)
$$i(\frac{u_2}{t}+\frac{u_2}{x})+u_1+\kappa u_3+\alpha \left(\alpha \sigma _1|u_2|^2+\alpha |u_1|^2+|u_3|^2\right)u_2=0,$$
(5)
$$i\frac{u_3}{t}+\kappa \left(u_1+u_2\right)+\left(\sigma _3|u_3|^2+\alpha |u_1|^2+\alpha |u_2|^2\right)u_3=\omega _0u_3.$$
(6)
where, in accord with the discussion above, we set $`\kappa `$ to be real and positive.

The coefficients $`\sigma _{1,3}`$ and $`\alpha `$ in Eqs. (4) - (6) are the generalized SPM and XPM coefficients, respectively. In particular, $`\alpha `$ is defined as a relative XPM coefficient between the first two waves and the third wave. In fact, the coefficients $`\sigma _1`$ and $`\sigma _3`$ both may be normalized to be $`\pm 1`$, unless they are equal to zero; however, it will be convenient to keep them as free parameters, see below (note that the SPM coefficients are always positive in the optical models, but in those describing density-stratified fluids they may have either sign). In optical models, all the coefficients $`\alpha `$ and $`\sigma _{1,3}`$ are positive. However, in the models describing the internal waves in stratified fluids, there is no inherent restriction on their signs, and some of them may indeed be negative.

The symmetry between the walk-off terms in Eqs. (4) and (5) is not really essential, and we will comment later on the more general case when these terms are generalized as follows:

$$\frac{u_1}{x}c_1\frac{u_1}{x},+\frac{u_2}{x}+c_2\frac{u_1}{x},$$
(7)
where $`c_1`$ and $`c_2`$ are different, but have the same sign. As for Eq. (ref3), it is obvious that the walk-off term in this equation, if any, can always be eliminated by means of a straightforward transformation.

We have kept only the most natural nonlinear SPM and XPM terms in Eqs. (4) - (6), i.e., the terms of the same types as in the standard GMT model. Additional terms, including nonlinear corrections to the linear couplings (e.g., a term $`|u_1|^2u_2`$ in Eq. (4)) may appear in more general models, such as a model of a deep (strong) BG .

Equations (4) - (6) conserve the norm, which has the physical meaning of energy in optics,

$$N\underset{n=1,2,3}{}_{\mathrm{}}^+\mathrm{}\left|u_n(x)\right|^2𝑑x,$$
(8)
the Hamiltonian,

$`H`$ $``$ $`H_{\mathrm{grad}}+H_{\mathrm{coupl}}+H_{\mathrm{focus}},`$ (9)
$`H_{\mathrm{grad}}`$ $``$ $`{\displaystyle \frac{i}{2}}{\displaystyle _{\mathrm{}}^+\mathrm{}}\left(u_1^{}{\displaystyle \frac{u_1}{x}}u_2^{}{\displaystyle \frac{u_2}{x}}\right)𝑑x+\mathrm{c}.\mathrm{c}.,`$ (10)
$`H_{\mathrm{coupl}}`$ $``$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}\left[u_1^{}u_2+\kappa u_3^{}\left(u_1+u_2\right)\right]𝑑x+\mathrm{c}.\mathrm{c}.,`$ (11)
$`H_{\mathrm{focus}}`$ $``$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}\left[{\displaystyle \frac{1}{2}}\alpha ^2\sigma _1\left(\left|u_1\right|^4+\left|u_2\right|^4\right)+{\displaystyle \frac{1}{2}}\sigma _3\left|u_3\right|^4+\alpha ^2\left|u_1\right|^2\left|u_2\right|^2+\alpha \left|u_3\right|^2\left(\left|u_1\right|^2+\left|u_2\right|^2\right)\right]𝑑x,`$ (12)
and the momentum, which will not be used here. In these expressions, the asterisk and $`\mathrm{c}.\mathrm{c}.`$ both stand for complex conjugation, $`H_{\mathrm{grad}}`$, $`H_{\mathrm{coupl}}`$, and $`H_{\mathrm{focus}}`$ being the gradient, linear-coupling, and self-focusing parts of the Hamiltonian. To obtain the Eqs. (4) - (6)from the Hamiltonian, the conjugate pairs of the variables are defined, in a standard fashion, as $`u_n,u_n^{}`$.

Our objective is to find various types of solitons existing in the generic three-wave system (4) - (6) and investigate their stability. Focusing first on the case (suggested by the analogy with GMTM) when the SPM term in Eq. (6) may be neglected (i.e., $`\sigma _3=0`$), in section 3 we find a general family of zero-velocity solitons in an exact analytical form. We will demonstrate that they are of two drastically different types: regular GS, and cuspons, i.e., solitons with a cusp singularity at the center, in which the soliton amplitude is finite, but the derivative is infinite; further, the energy of the cuspons is finite. Cuspons are known to exist in degenerate models without linear terms (except for the evolution term such as $`u/t`$), i.e., without a linear spectrum, a well-known example being the exactly integrable Camassa-Holm (CH) equation  (see also ). Our model resembles the CH one in the sense that both give rise to coexisting solutions in the form of regular solitons and cuspons. The cause for the existence of these singular solitons in our model is the fact that, looking for a zero-velocity soliton solution, one may eliminate the field $`u_3`$ by means of an algebraic relation following, in this case, from Eq. (6). The subsequent substitution of that result into the first two equations (4) and (5) produces a rational nonlinearity in them, the corresponding rational functions featuring a singularity at some (critical) value of the soliton’s amplitude. If the amplitude of a regular-soliton solution is going to exceed the critical value, it actually cannot exist, and in the case when $`\sigma _3=0`$ it is replaced by a cuspon, whose amplitude is exactly equal to the critical value.

In the limit $`\kappa 0`$, which corresponds to the vanishing linear coupling between the first two and third waves, the cuspon resembles a peakon, which is a finite-amplitude solitary wave with a jump of its first derivative at the center. Note that peakon solutions, coexisting with regular solitons (this property is shared by our model), are known in a slightly different (also integrable) version of the CH equation, see, e.g., Refs. . We also note that soliton solutions with a discontinuity in the first derivative have been found in the BG model (which does contain a linear part) in the case where the grating parameter changes abruptly .

Then, we show that, when the SPM term is restored in Eq. (6) (i.e., $`\sigma _30`$; the presence or absence of the SPM terms $`\sigma _1`$ in Eqs. (4) and (5) is not crucially important), the system supports a different set of soliton solutions. These are regular GS and, depending on the sign of certain parameters, a family of peakons, which, this time, appear as generic solutions, unlike the case $`\sigma _3=0`$, when they only exist as a limiting form of the cuspon solutions corresponding to $`\kappa 0`$. As far as we know, the model formulated in the present work is the first spatially uniform non-degenerate one (i.e., a model with a non-vanishing linear part) which yields both cuspons and peakons.

### 1.2 Stability of the solitons and spatiotemporal collapse

As concerns the dynamical stability of the various solitons in the model (4) -(6), in this work we limit ourselves to direct simulations, as a more rigorous approach, based on numerical analysis of the corresponding linear stability-eigenvalue problem , is technically difficult in the case of cuspons and peakons (results of such an analysis, based on the Evans-function technique, will be presented elsewhere). In fact, direct simulations of perturbed cuspons and peakons is a hard problem too, but we have concluded that identical results concerning the stability are produced (see section 3 below) by high-accuracy finite-difference and pseudo-spectral methods (each being implemented in more than one particular form), which lends the results credibility. A general conclusion is that the regular solitons are always stable. As for the cuspons and peakons, they may be either stable or unstable.

If the cusp is strong enough, the numerical results presented below demonstrate that the instability of a cuspon initiates formation of a genuine singularity, i.e., onset of a spatiotemporal collapse  in the present one-dimensional model. Before proceeding to the consideration of solitons in the following sections, it is relevant to discuss collapse phenomenon in some detail.

A simple virial-type estimate for the possibility of the collapse can be done, assuming that the field focuses itself in a narrow spot with a size $`(t)`$, amplitude $`\mathrm{}(t)`$, and a characteristic value $`K(t)`$ of the field’s wavenumber . The conservation of the norm (8) imposes a restriction $`\mathrm{}^2LN`$, i.e., $`LN/\mathrm{}^2`$. Next, the self-focusing part (9) of the Hamiltonian (9), which drives the collapse, can be estimated as

$$H_{\mathrm{focus}}\mathrm{}^4LN\mathrm{}^2.$$
(13)
On the other hand, the collapse can be checked by the gradient term (10) in the full Hamiltonian, that, in the same approximation, can be estimated as $`H_{\mathrm{grad}}\mathrm{}^2KLNK`$. Further, Eqs. (4) - (6) suggest an estimate $`K\mathrm{}^2`$ for a characteristic wavenumber of the wave field \[the same estimate for $`K`$ follows from an expression (21) for the exact stationary-soliton solution given below\], thus we have $`H_{\mathrm{grad}}N\mathrm{}^2`$. Comparing this with the expression (13), one concludes that the parts of the Hamiltonian promoting and inhibiting the collapse scale the same way as $`\mathrm{}\mathrm{}`$ (or $`L0`$), hence a weak collapse  may be possible (but does not necessarily take place) in systems of the present type. We stress that, in one-dimensional models of GS studied thus far and based on GMTM, collapse has never been reported. The real existence of the collapse in the present one-dimensional three-wave GS model, which will be shown in detail below as a result of numerical simulations, is therefore a novel dynamical feature, and it seems quite natural that cuspons and peakons, in the case when they are unstable, play the role of catalysts stimulating the onset of the collapse. The possibility of a real collapse in a 1D system is quite interesting by itself, and also because experimental observation of spatio-temporal self-focusing in nonlinear optical media is a subject of considerable interest, see, e.g., Ref. .

## 2 Analytical solutions

### 2.1 The dispersion relation

The first step in the investigation of the system is to understand its linear spectrum. Substituting $`u_{1,2,3}\mathrm{exp}(ikxi\omega t)`$ into Eqs. (4 -6), and omitting nonlinear terms, we arrive at a dispersion equation,

$$(\omega ^2k^21)(\omega \omega _0)=\mathrm{\hspace{0.17em}2}\kappa ^2(\omega 1).$$
(14)
If $`\kappa =0`$, the third wave decouples, and the coupling between the first two waves produces a commonly known gap, so that the solutions to Eq. (14) are $`\omega _{1,2}=\pm \sqrt{1+k^2}`$ and $`\omega _3=\omega _0`$. If $`\kappa 0`$, the spectrum can be easily understood by treating $`\kappa `$ as a small parameter. However, the following analysis is valid for all values of $`\kappa `$ in the range $`0<\kappa ^2<1`$.

First, consider the situation when $`k=0`$. Three solutions of Eq. (14) are then

$$\omega =1,\omega =\omega _\pm (\omega _01)/2\pm \sqrt{(\omega _0+1)^2/4+2\kappa ^2}.$$
(15)
It can be easily shown that $`\omega _{}<\mathrm{min}\{\omega _0,1\}\mathrm{max}\{\omega _0,1\}<\omega _+`$, so that one always has $`\omega _{}<1`$, while $`\omega _+<1`$ if $`\omega _0<1\kappa ^2`$, and $`\omega _+>1`$ if $`\omega _0>1\kappa ^2`$. Next, it is readily seen that, as $`k^2\mathrm{}`$, either $`\omega ^2k^2`$, or $`\omega \omega _0`$. Each branch of the dispersion relation generated by Eq. (14) is a monotonic function of $`k^2`$. Generic examples of the spectrum are shown in Fig. 2, where the panels (a) and (b) pertain, respectively, to the cases $`\omega _0<1\kappa ^2`$ with $`\omega _+<1`$, and $`\omega _0>1`$ with $`\omega _+>1`$. The intermediate case, $`1\kappa ^2<\omega _0<1`$, is similar to that shown in panel (a), but with the points $`\omega _+`$ and $`1`$ at $`k=0`$ interchanged. When $`\omega _0<1`$, the upper gap in the spectrum is $`\mathrm{min}\{\omega _+,1\}<\omega <`$ $`\mathrm{max}\{\omega _+,1\}`$, while the lower gap is $`\omega _{}<\omega <\omega _0`$. When $`\omega _0>1`$, the upper gap is $`\omega _0<\omega <\omega _+`$, and the lower one is $`\omega _{}<\omega <1`$.

### 2.2 Gap solitons

The next step is to search for GS solutions to the full nonlinear system. In this work, we confine ourselves to the case of zero-velocity GS, substituting into Eqs. (4) - (6)

$$u_n(x,t)=U_n(x)\mathrm{exp}(i\omega t),n=1,2,3,$$
(16)
where it is assumed that the soliton’s frequency $`\omega `$ belongs to one of the gaps. In fact, even the description of zero-velocity solitons is quite complicated. Note, however, that if one sets $`\kappa =0`$ in Eqs. (4) - (6), keeping nonlinear XPM couplings between the first two and third waves, the gap which exists in the two-wave GMT model remains unchanged, and the corresponding family of GS solutions does not essentially alter, in accord with the principle that nonlinear couplings cannot alter gaps or open a new one if the linear coupling is absent ; nevertheless, the situation is essentially different if $`\kappa `$ is vanishingly small, but not exactly equal to zero, see below.

The substitution of (16) into Eqs. (4) and (5) leads to a system of two ordinary differential equations for $`U_1(x)`$ and $`U_2(x)`$, and an algebraic relation for $`U_3(x)`$,

$$iU_1^{}=\omega U_1+U_2+\kappa U_3+\alpha \left(\alpha \sigma _1|U_1|^2+\alpha |U_2|^2+|U_3|^2\right)U_1,$$
(17)
$$iU_2^{}=\omega U_2+U_1+\kappa U_3+\alpha \left(\alpha \sigma _1|U_2|^2+\alpha |U_1|^2+|U_3|^2\right)U_2,$$
(18)
$$(\omega _0\omega +\sigma _3|U_3|^2+\alpha |U_1|^2+\alpha |U_2|^2)U_3=\kappa (U_1+U_2),$$
(19)
where the prime stands for $`d/dx`$. To solve these equations, we substitute $`U_{1,2}=A_{1,2}(x)\mathrm{exp}\left(i\varphi _{1,2}(x)\right)`$ with real $`A_n`$ and $`\varphi _n`$. After substituting the expression (19) into equations (17,18), and some simple manipulations, it can be found that $`\left(A_1^2A_2^2\right)^{}=0`$ and $`\left(\varphi _1+\varphi _2\right)^{}=0`$. Using the condition that the soliton fields vanish at infinity, we immediately conclude that

$$A_1^2(x)=A_2^2(x)S(x);$$
(20)
as for the constant value of $`\varphi _1+\varphi _2`$, it may be set equal to zero without loss of generality, so that $`\varphi _1(x)=\varphi _2(x)\varphi (x)/2`$, where $`\varphi (x)`$ is the relative phase of the two fields. After this, we obtain two equations for $`S(x)`$ and $`\varphi (x)`$ from Eqs. (17) and (18),

$$\varphi ^{}=2\omega 2\mathrm{cos}\varphi 2\alpha ^2\left(1+\sigma _1\right)SS^1U_3^2\left(\omega _0\omega \sigma _3U_3^2\right),$$
(21)
$$S^{}=2S\mathrm{sin}\varphi 2\kappa \sqrt{S}U_3\mathrm{sin}\left(\varphi /2\right),$$
(22)
and Eq. (19) for the third wave $`U_3`$ takes the form of a cubic algebraic equation

$$U_3\left(\omega _0\omega 2\alpha S\sigma _3|U_3|^2\right)=\mathrm{\hspace{0.17em}2}\kappa \sqrt{S}\mathrm{cos}\left(\varphi /2\right),$$
(23)
from which it follows that $`U_3`$ is a real-valued function.

This analytical consideration can be readily extended for more general equations (4) and (5) that do not assume the symmetry between the waves $`u_1`$ and $`u_2`$, i.e., with the group-velocity terms in the equations altered as in Eq. (7). In particular, the relation (20) is then replaced by $`c_1A_1^2(x)=c_2A_2^2(x)S(x)`$. The subsequent analysis is similar to that above, and leads to results for the asymmetric model that are qualitatively similar to those presented below for the symmetric case.

Equations (21) and (22) have a Hamiltonian structure, as they can be represented in the form

$$\frac{dS}{dx}=\frac{H}{\varphi },\frac{d\varphi }{dx}=\frac{H}{S},$$
(24)
with the Hamiltonian

$$H=\mathrm{\hspace{0.17em}2}S\mathrm{cos}\varphi +\alpha ^2\left(1+\sigma _1\right)S^2+2\omega S+U_3^2\left(\omega _0\omega 2\alpha S\right)\frac{3}{2}\sigma _3U_3^4,$$
(25)
which is precisely a reduction of the Hamiltonian (9) of the original system (4) - (6) for the solutions of the present type. Note that $`H`$ is here regarded as a function of $`S`$ and $`\varphi `$, and the relation (23) is regarded as determining $`U_3`$ in terms of $`S`$ and $`\varphi `$. We stress that the dependence $`U_3(S,\varphi )`$ was taken into account when deriving the Hamiltonian representation (24).

For soliton solutions, the boundary conditions at $`x=\pm \mathrm{}`$ yield $`H=0`$ so that the solutions can be obtained in an implicit form,

$$2S\mathrm{cos}\varphi +\alpha ^2\left(1+\sigma _1\right)S^2+2\omega S+U_3^2\left(\omega _0\omega 2\alpha S\right)\left(3/2\right)\sigma _3U_3^4=\mathrm{\hspace{0.17em}0}.$$
(26)
In principle, one can use the relations (23) and (26) to eliminate $`U_3`$ and $`\varphi `$ and so obtain a single equation for $`S`$. However, this is not easily done unless $`\sigma _3=0`$ (no SPM term in Eq. (ref3)), and so we proceed to examine this special, but important, case first. Note that the no-SPM case also plays an important role for GMTM, which is exactly integrable by means of the inverse scattering transform just in this case .

### 2.3 Cuspons, the case $`\sigma _3=0`$

Setting $`\sigma _3=0`$ makes it possible to solve Eq. (23) for <sub>3</sub> explicitly in terms of $`S`$ and $`\varphi `$,

$$U_3=\frac{2\kappa \sqrt{S}\mathrm{cos}\left(\varphi /2\right)}{\omega _0\omega 2\alpha S}.$$
(27)
For simplicity, we also set $`\sigma _1=0`$ in Eqs. (4) and (5) and subsequent equations, although the latter assumption is not crucially important for the analysis developed below. Indeed, the analysis is based on the fact that the field $`U_3`$ can be explicitly eliminated by means of Eq. (27), which is not affected by $`\sigma _1`$. If $`\sigma _1`$ is kept in the system, it merely renormalizes some coefficients in the formulas derived below.

At the next step, one can also eliminate $`\varphi `$, using Eqs. (26 ) and (27), to derive a single equation for $`S`$,

$$\left(dS/dx\right)^2=\mathrm{\hspace{0.17em}4}S^2F(S),$$
(28)
$$F(S)(1\omega \frac{1}{2}\alpha ^2S)\left[2\left(1+\frac{\kappa ^2}{\omega _0\omega 2\alpha S}\right)(1\omega \frac{1}{2}\alpha ^2S)\right].$$
(29)
The function $`F(S)`$ has either one or three real zeros $`S_0`$. One is

$$S_{01}=\mathrm{\hspace{0.17em}2}\left(1\omega \right)/\alpha ^2,$$
(30)
and the remaining two, if they exist, are real roots of the quadratic equation,

$$(2+2\omega +\alpha ^2S_0)(\omega _0\omega 2\alpha S_0)+4\kappa ^2=0.$$
(31)
Only the smallest positive real root of Eq. (31), to be denoted $`S_{02}`$ (if such exists), will be relevant below. Note, incidentally, that $`F(S)`$ cannot have double roots. It is easy to see that a consequence of this fact is that Eq. (28) cannot generate kink solutions, which have different limits as $`x\pm \mathrm{}`$, for both of which the right-hand-side of (28) must have a double zero.

For a bright-soliton solution of (28), we need first that $`F(0)>0`$ (in this work, we do not consider dark solitons, nor “anti-dark” solitons, i.e., solitons on top of a finite-amplitude flat background, a reason being that there is little chance that the flat background would be modulationally stable). Comparing the condition $`F(0)>0`$ with the expressions given in Section 2.1 for the gaps in the linear spectrum, it is readily shown that this condition is exactly equivalent to requiring that $`\omega `$ belongs to either the upper or the lower gap of the linear spectrum. We note that the coupling to the third wave gives rise to nonlinearity of the rational type in the expression (29), despite the fact that the underlying system (4) - (6) contains only cubic polynomial nonlinear terms. Even if the coupling constant $`\kappa `$ is small, it is clear that the rational nonlinearity may produce a strong effect in a vicinity of a critical value of the squared amplitude at which the denominator in the expression (29) vanishes,

$$S_{\mathrm{cr}}=\left(\omega _0\omega \right)/2\alpha ,$$
(32)
where one must have $`\alpha (\omega _0\omega )>0`$ (otherwise, this critical value is not relevant).

If $`S_{\mathrm{cr}}>0`$, the structure of the soliton crucially depends on whether, with an increase of $`S`$, the function $`F(S)`$ defined by Eq. (29 ) first reaches zero at $`S=S_0>0`$ (i.e. either $`S=S_{01}`$ or $`S=S_{02}`$, whichever is the smaller positive value), or, instead, it first reaches the singularity at $`S=S_{\mathrm{cr}}`$, i.e., whether $`0<S_0<S_{\mathrm{cr}}`$, or $`0<S_{\mathrm{cr}}<S_0`$. In the former case, the existence of $`S_{\mathrm{cr}}`$ plays no role, and the soliton is a regular one, having the amplitude $`\sqrt{S_0}`$. This soliton may be regarded as obtained by a smooth deformation from the usual GS known in GMTM at $`\kappa =0`$.

In the case $`0<S_{\mathrm{cr}}<S_0`$, as the soliton cannot have an amplitude larger than $`\sqrt{S_{\mathrm{cr}}}`$, the amplitude takes this critical value. The soliton is singular in this case, being a cuspon (see details below), but, nevertheless, it is an absolutely relevant solution. The remaining possibilities are that either $`S_{\mathrm{cr}}<0`$ and $`S_0>0`$, or vice versa; then the soliton may only be, respectively, regular or singular. Of course no soliton exists if both $`S_0`$ and $`S_{\mathrm{cr}}`$ are negative. Further, using the symmetries of the equations, it is readily shown that for all these soliton solutions, $`S(x)`$ is symmetric about its center, which may be set at $`x=0`$, that is, $`S(x)`$ is an even function of $`x`$. For the cuspon solutions, and for those regular solutions whose squared amplitude is $`S_{01}`$, it can also be shown that the phase variable $`\psi (x)=\varphi (x)\pi `$ and $`U_3(x)`$ are odd functions of $`x`$, while for those regular solutions whose squared amplitude is $`S_{02}`$ the phase variable $`\varphi (x)`$ and $`U_3(x)`$ are, respectively, odd and even functions of $`x`$.

It is now necessary to determine which parameter combinations in the set ($`\omega ,\omega _0,\alpha `$) permit the options described above. The most interesting case occurs when $`\omega _0>\omega `$ (so that $`\omega `$ belongs to the lower gap, see Fig. 2) and $`\alpha >0`$ (the latter condition always holds in the applications to nonlinear optics). In this case, it can be shown that the root $`S_{02}`$ of Eq. (31) is not relevant, and the options are determined by the competition between $`S_{01}`$ and $`S_{\mathrm{cr}}`$. The soliton is a cuspon ($`0<S_{\mathrm{cr}}<S_{01}`$) if

$$\alpha (\omega _0\omega )<\mathrm{\hspace{0.17em}4}(1\omega ).$$
(33)
In effect, the condition (33) sets an upper bound on $`\alpha `$ for given $`\omega _0`$ and $`\omega `$. In particular, this condition is always satisfied if $`0<\alpha <4`$.

If, on the other hand, the condition (33) does not hold (i.e., $`0<S_{01}<S_{\mathrm{cr}}`$), we obtain a regular soliton. In a less physically relevant case, when again $`\omega _0>\omega `$ but $`\alpha <0`$, cuspons do not occur (as this time $`S_{\mathrm{cr}}<0`$, see Eq. (32), and only regular solitons may exist.

Next we proceed to the case $`\omega _0<\omega `$, so that $`\omega `$ is located in the upper gap of the linear spectrum. For $`\alpha >0`$, we have $`S_{\mathrm{cr}}<0`$, hence only regular solitons may occur, and indeed in this case there is always at least one positive root $`S_0`$, so a regular soliton does exist. If $`\alpha <0`$, then we have $`S_{\mathrm{cr}}>0`$, but if $`\omega _0<1\kappa ^2`$ (when also $`\omega <1`$), there is at least one positive root $`S_0<S_{\mathrm{cr}}`$; thus, only a regular soliton can exist in this case too. On the other hand, if $`\alpha <0`$ and $`\omega _0>1\kappa ^2`$ (and then $`\omega >1`$), there are no positive roots $`S_0`$, and so only cuspons occur.

Let us now turn to a detailed description of the cuspon’s local structure near its center, when $`S`$ is close to $`S_{\mathrm{cr}}`$. From the above analysis, one sees that cuspons occur whenever $`\omega `$ lies in the lower gap, with $`\omega _0>\omega `$ and $`\alpha >0`$, so that the criterion (33) is satisfied, or when $`\omega `$ lies in the upper gap with $`1\kappa ^2<\omega _0<\omega `$ and $`\alpha <0`$. To analyze the structure of the cuspon, we first note that, as it follows from Eq. (26), one has $`\mathrm{cos}\varphi =1`$ (i.e., $`\varphi =\pi `$) when $`S=S_{\mathrm{cr}}`$, which suggest to set

$$S_{\mathrm{cr}}S\delta \kappa ^2R,1+\mathrm{cos}\varphi \delta \rho ,$$
(34)
where $`\delta `$ is a small positive parameter, and the stretched variables $`R`$ and $`\rho `$ are positive. At the leading order in $`\delta `$, it then follows from Eq. (26) that $`\rho =\rho _0R`$, where

$$\rho _0\alpha ^3(S_{01}S_{\mathrm{cr}}).$$
(35)
As it follows from the above analysis, $`\rho _0`$ is always positive for a cuspon. We also stretch the spatial coordinate, defining $`x\delta ^{3/2}\kappa ^2y`$, the soliton center being at $`x=0`$. Since $`S(x)`$ is an even function of $`x`$, it is sufficient to set $`x>0`$ in this analysis. Then, on substituting the first relation from Eq. (34) into Eq. (28), we get, to the leading order in $`\delta `$, an equation

$$R\left(dR/dy\right)^2=\rho _0S_{\mathrm{cr}}^2/\alpha ^2K^2,$$
(36)
so that

$$R=\left(3Ky/2\right)^{2/3}.$$
(37)
Note that in the original unstretched variables, the relation (37 ) shows that, near the cusp,

$$S_{\mathrm{cr}}S(x)(3K\kappa x/2)^{2/3},$$
(38)
$$dS/dx(2/3)^{1/3}\left(K\kappa \right)^{2/3}x^{1/3},$$
(39)
and it follows from Eq. (27) that $`U_3`$ is unbounded near the cusp,

$$U_3(S_{\mathrm{cr}}/\alpha )(2\alpha \rho _0K^2/3\kappa x)^{1/3}.$$
(40)
In particular, Eq. (39) implies that, as $`K\kappa `$ decreases, the cusp gets localized in a narrow region where $`|x|_{}^<K^2\kappa ^2`$ (outside this region, $`\left|dS/dx\right|`$ is bounded and shows no cusp). Note that this limit can be obtained either as $`\kappa ^20`$, or as $`\rho _00`$ (recall that $`\rho _0`$ is defined in Eq. (35)).

An example of the cuspon is shown in Fig. 3. Although the first derivative in the cuspon is singular at its center, as it follows from Eq. (39) \[see also Fig. 3(a)\], and its $`U_3`$ component diverges at $`x0`$ as per Eq. (40), it is easily verified that the value of the Hamiltonian (9) \[and, obviously, the value of the norm (8) too\] is finite for the cuspon solution. These solitons are similar to cuspons found as exact solutions to the Camassa-Holm (CH) equation , which have a singularity of the type $`|x|^{1/3}`$ or $`|x|^{2/3}`$ as $`|x|0`$, cf. Eqs. (38) and (39). The CH equation is integrable, and it is degenerate in the sense that it has no linear terms except for $`u/t`$ (which makes the existence of the solution with a cusp singularity possible). Our three-wave system (4) - (6) is not degenerate in that sense; nevertheless, the cuspon solitons are possible in it because of the model’s multicomponent structure: the elimination of the third component generates the non-polynomial nonlinearity in Eqs. (17), (18), and, finally, in Eqs. (22) and (28), which gives rise to the cusp. It is noteworthy that, as well as the CH model, ours gives rise to two different coexisting families of solitons, viz., regular ones and cuspons. It will be shown below that the solitons of both types may be stable.

Of course, the presence of the singularities in $`U_3(x)`$ and $`dS/dx`$ at $`x0`$ suggests that higher-order terms, such as the higher-order dispersion, should be taken into regard in this case. The fact that the cuspon’s Hamiltonian converges despite these singularities, as well as a direct analysis, suggest that such higher-order terms will smooth the shape of the cuspon in a very narrow layer for small $`x`$, allowing for large but not diverging values of the fields. However, the small higher-order terms will not essentially alter the global shape of the cuspons. In the next section we will show that, in fact, the genuine generic singular solitons are (in the presence of the SPM terms) peakons, for which the singularities are much weaker, hence the latter issue is still less significant. Besides that, it appears to be an issue of principal interest to understand what types of solitons the system may generate without intrinsic dispersion (cf. the situation for the traditional GMTM, in which the spectrum of soliton solutions is completely altered by the addition of intrinsic dispersion ).

In the special case $`\kappa 1`$, when the third component is weakly coupled to the first two ones in the linear approximation (in terms of the optical model represented by Fig. 1, it is the case when the subgratings shown by the dashed lines are very weak), straightforward inspection of the above results shows that the cuspons look like peakons; that is, except for the above-mentioned narrow region of the width $`|x|\kappa ^2`$, where the cusp is located, they have the shape of a soliton with a discontinuity in the first derivative of $`S(x)`$ and a jump in the phase $`\varphi (x)`$, which are the defining features of peakons (). A principal difference of true peakons from cuspons is that the first derivative does not diverge inside a peakon, but is of course, discontinuous.

An important result of our analysis is that the family of solitons obtained in the limit $`\kappa 0`$ is drastically different from that in the model where one sets $`\kappa =0`$ from the very beginning. In particular, in the most relevant case, with $`\omega _0>\omega `$ and $`\alpha >0`$, the family corresponding to $`\kappa 0`$ contains regular solitons whose amplitude is smaller than $`\sqrt{S_{\mathrm{cr}}}`$; however, the solitons whose amplitude at $`\kappa =0`$ is larger than$`\sqrt{S_{\mathrm{cr}}}`$, i.e., the ones whose frequencies belong to the range (33) \[note that the definition of $`S_{\mathrm{cr}}`$ does not depend on $`\kappa `$ at all, see Eq. (32)\], are replaced by the peakons which are constructed in a very simple way: drop the part of the usual soliton above the critical level $`S=S_{\mathrm{cr}}`$, and bring together the two symmetric parts which remain below the critical level, see Fig. 3(b).

It is interesting that peakons are known as exact solutions to a version of the integrable CH equation slightly different from that which gives rise to the cuspons. As well as in the present system, in that equation the peakons coexist with regular solitons . In the next subsection, we demonstrate that the peakons, which are found only as limit-form solutions in the no-SPM case $`\sigma _3=0`$, become generic solutions in the case $`\sigma _30`$.

### 2.4 Peakons, the case $`\sigma _30`$

A natural question is whether the cuspon solutions are structurally stable, i.e., if they will persist on inclusion of terms that were absent in the analysis presented above (the other type of the stability, viz., dynamical stability against small initial perturbations, will be considered in the next section). Here, we address this issue by restoring the SPM term in Eq. (6), that is, we now set $`\sigma _30`$, but assume that it is a small parameter. Note that, in the application to nonlinear optics, one should expect that $`\sigma _3>0`$, but there is no such a restriction on the sign of $`\sigma _3`$ in the application to the flow of a density-stratified fluid. We still keep $`\sigma _1=0`$, as the inclusion of the corresponding SPM terms in Eqs. (4) and (5) amounts to trivial changes both in the above analysis, and in that presented below. On the other hand, we show below that the inclusion of the SPM term in Eq. (6) is a structural perturbation which drastically changes the character of the soliton solutions.

In view of the above results concerning the cuspons, we restrict our discussion here to the most interesting case when $`S(x)`$ is an even function of $`x`$, while $`\psi (x)=\varphi (x)\pi `$ and $`U_3(x)`$ are odd functions. In principle, one can use the relations (23) and (26) to eliminate $`\varphi `$ and $`U_3`$ and so obtain a single equation for $`S`$ (a counterpart to Eq. (28)), as it was done above when $`\sigma _3=0`$. However, when $`\sigma _30`$, it is not possible to do this explicitly. Instead, we shall develop an asymptotic analysis valid for $`x0`$, which will be combined with results obtained by direct numerical integration of Eqs. (21) and (22), subject of course to the constraints (23) and (26). Since singularities only arise at the center of the soliton (i.e., at $`x=0`$) when $`\sigma _3=0`$, it is clear that the introduction of a small $`\sigma _30`$ will produce only a small deformation of the soliton solution in the region where $`x`$ is bounded away from zero.

First, we consider regular solitons. Because the left-hand side of Eq. (23) is not singular at any $`x`$, including the point $`x=0`$, when $`\sigma _3=0`$, we expect that regular solitons survive a perturbation induced by $`\sigma _30`$. Indeed, if there exists a regular soliton, with $`S_0S(x=0)`$, and $`\varphi (x=0)=\pi `$ and $`U_3(x=0)=0`$, it follows from Eq. (26) that the soliton’s amplitude remains exactly the same as it was for $`\sigma _3=0`$, due to the fact that the regular soliton has $`U_3(x=0)=0`$.

Next, we turn to the possibility of singular solutions, that is, cuspons or peakons. Since we are assuming that $`S_0=S(x=0)`$ is finite, and that $`\varphi (x=0)=\pi `$, it immediately follows from Eq. (23) that when $`\sigma _30`$, $`U_3`$ must remain finite for all $`x`$, taking some value $`U_00`$, say, as $`x+0`$. Since $`U_3`$ is an odd function of $`x`$, and $`U_00`$, there must be a discontinuity in $`U_3`$ at $`x=0`$, i.e., a jump from $`U_0`$ to $`U_0`$. This feature is in marked contrast to the cuspons for which $`U_3`$ is infinite at the center, see Eq. (40). Further, it then follows from Eq. (22) that, as $`x0`$, there is also a discontinuity in $`dS/dx`$, with a jump from $`2\kappa U_0\sqrt{S_0}`$ to $`2\kappa U_0\sqrt{S_0}`$. Hence, if we can find soliton solutions of this type, with $`U_00`$, they are necessarily peakons, and we infer that cuspons do not survive the structural perturbation induced by $`\sigma _30`$.

Further, if we assume that $`U_00`$, then Eq. (23), taken in the limit $`x0`$, immediately shows that

$$2\alpha (S_{\mathrm{cr}}S_0)=\sigma _3U_0^2$$
(41)
(recall that $`S_{\mathrm{cr}}`$ is defined by Eq. (32)). Next, the Hamiltonian relation (26), also taken in the limit $`x0`$, shows that

$$\frac{\rho _0}{\alpha }S_0\alpha ^2S_0(S_{\mathrm{cr}}S_0)=\frac{1}{2}\sigma _3U_0^4,$$
(42)
where we have used Eq. (41) \[recall that $`\rho _0`$ is defined by Eq. (35)\]. Elimination of $`U_0`$ from Eqs. (41) and (42) yields a quadratic equation for $`S_0`$, whose positive roots represent the possible values of the peakon’s amplitude.

We recall that for a cuspon which exists at $`\sigma _3=0`$ one has $`\rho _0>0`$, i.e., the amplitude of the corresponding formal regular soliton exceeds the critical value of the amplitude, see Eq. (35). Then, if we retain the condition $`\rho _0>0`$, it immediately follows from Eqs. (41) and (42) that no peakons may exist if the SPM coefficient in Eq. (6) is positive, $`\sigma _3>0`$. Indeed, Eq. (41) shows that $`S_{\mathrm{cr}}S_0>0`$ if $`\sigma _3>0`$, which, along with $`\rho _0>0`$, leads to a contradiction in the relation (42).

Further, it is easy to see that a general condition for the existence of peakons following from Eqs. (41) and (42) is

$$\sigma _3\rho _0<0,$$
(43)
hence peakons are possible if $`\sigma _3<0`$, or if we keep $`\sigma _3>0`$ but allow $`\rho _0<0`$. In the remainder of this section, we will show that peakons may exist only if $`\rho _0>0`$. Hence, it follows from the necessary condition (43) that peakons may indeed be possible solely in the case $`\sigma _3<0`$. On the other hand, regular solitons do exist in the case $`\sigma _3>0`$ (i.e., in particular, in nonlinear-optics models), as they have $`U_0=0`$, hence neither Eq. (41 ) nor its consequence in the form of the inequality (43) apply to regular solitons. The existence of (stable) peakons for $`\sigma _3<0`$, and of (also stable) regular solitons for $`\sigma _3>0`$ will be confirmed by direct numerical results presented in the next section.

To obtain a necessary condition (which will take the form of $`\rho _0>0`$) for the existence of the peakons, we notice that the existence of any solitary wave implies the presence of closed dynamical trajectories in the phase plane of the corresponding dynamical system, which is based on the ordinary differential equations (21) and (22), supplemented by the constraint (23). Further, at least one stable fixed point (FP) must exist inside such closed trajectories, therefore the existence of such a stable FP is a necessary condition for the existence of any solitary wave.

The FPs are found by equating to zero the right-hand sides of Eq. (21 ) and (22), which together with Eq. (23) give three equations for the three coordinates $`\varphi ,S`$ and $`U_3`$ of the FP. First of all, one can find a trivial unstable FP of the dynamical system,

$$\mathrm{cos}\varphi =\frac{\omega +\kappa ^2/(\omega _0\omega )}{1+\kappa ^2/(\omega _0\omega )},S=\mathrm{\hspace{0.17em}0},$$

which does not depend on $`\sigma _3`$. Then, three nontrivial FPs can be found, with their coordinates $`\varphi _{}`$, $`S_{}`$ and $`U_3`$ given by the following expressions:

$$\varphi _{}^{(1)}=\pi ,S_{}^{(1)}=\frac{1\omega }{\alpha ^2}=\frac{1}{2}S_{01},U_3^{(1)}=\mathrm{\hspace{0.17em}0},$$
(44)
$$\varphi _{}^{(2)}=\pi ,(2\sigma _3)S_{}^{(2)}=2S_{\mathrm{cr}}\frac{\sigma _3}{2}S_{01},(2\sigma _3)\left[\alpha U_3^{(2)}\right]^2=\rho _0\alpha ^3S_{\mathrm{cr}},$$
(45)
$$(2\sigma _3)S_{}^{(3)}=2S_{\mathrm{cr}}\frac{1}{2}\sigma _3S_{01}+\frac{\kappa ^2}{\alpha },(2\sigma _3)\left[\alpha U_3^{(3)}\right]^2=\rho _0\alpha ^3S_{\mathrm{cr}}\alpha ^2\kappa ^2,$$

$$\mathrm{cos}\left(\varphi _{}^{(3)}/2\right)=\frac{1}{2}\kappa U_3^{(3)}/\sqrt{S_{}^{(3)}},$$
(46)
where the superscript is a number label for the FP. To be specific, we now consider the case of most interest, when both $`S_{01}>0`$ and $`S_{\mathrm{cr}}>0`$ . In this case, the FP given by Eqs. (44 ) exists for all $`\sigma _3`$ and all $`\rho _0`$. However, for small $`\sigma _3`$ (in fact $`\sigma _3<2`$ is enough) and small $`\kappa `$, the FPs given by Eqs. (45) and (46) exist only when $`\rho _0>0`$. Indeed, they exist only for $`\rho _0>\alpha ^3S_{01}`$ and $`\rho _0>\alpha ^3S_{01}+\kappa ^2`$, respectively, or, on using the definition (35 ) of $`\rho _0`$, when $`S_{01}>2S_{\mathrm{cr}}`$ and $`S_{01}>2S_{\mathrm{cr}}+\kappa ^2/\alpha `$, respectively.

Let us first suppose that $`\rho _0<0`$. Then there is only the single non-trivial FP, namely the one given by Eqs. (44). This FP is clearly associated with the regular solitons, whose amplitude at the crest is $`S_{01}`$. Hence, we infer that for $`\rho _0<0`$ there are no other solitary-wave solutions, and in particular, no peakons (and no cuspons when $`\sigma _3=0`$ either, in accordance with what we have already found in subsection 2.3 above). When combined with the necessary condition (43) for the existence of peakons, we infer that there are no peakons when $`\sigma _3>0`$, thus excluding peakons from applications to the nonlinear-optics models, where this SPM coefficient is positive. However, peakons may occur in density-stratified fluid flows, where there is no inherent restriction on the sign of $`\sigma _3`$. This case is considered below, but first we note that in the case $`\rho _0<0`$ and $`\sigma _3>0`$ (which includes the applications to nonlinear optics), the same arguments suggest that there may be periodic solutions with a peakon-type discontinuity at the crests; indeed, our numerical solutions of the system (21,22) show that this is the case.

Next, we suppose that $`\rho _0>0`$. First, if $`S_{01}<2S_{\mathrm{cr}}`$ then there is again the single non-trivial FP given by (44). But now, by analogy with the existence of cuspons when $`\rho _0>0`$ and $`\sigma _3=0`$ , , we infer that the solitary wave solution which is associated with this fixed point is a peakon, whose squared amplitude $`S_0`$ for small $`\sigma _3`$ is close to $`S_{\mathrm{cr}}`$, while the FP has $`S_{}^{(1)}=S_{01}/2<S_{\mathrm{cr}}`$.

If, on the other hand, $`S_{01}>2S_{\mathrm{cr}}`$, the FPs given by Eqs.(45) and (46) become available as well. We now infer that the peakon solitary-wave solution continues to exist, and for sufficiently small $`\sigma _3`$ and $`\kappa `$ it is associated with the FP given by Eq. (45). Although Eq. (45) implies that $`S_{}^{(2)}S_{\mathrm{cr}}`$, and the peakon’s squared amplitude $`S_0`$, determined by Eqs. ( 41) and (42), is also approximately equal to $`S_{\mathrm{cr}}`$, we nevertheless have $`S_0>S_{}^{(2)}`$ as required. Note that, in the present case, the FPs given by Eqs. (44) and (46) lie outside the peakon’s homoclinic orbit. In Fig. 4, we show a plot of a typical peakon obtained, in this case, by numerical solution Eqs. (21) and (22).

## 3 Numerical results

### 3.1 Simulation techniques

The objectives of direct numerical simulations of the underlying equations (4) - (6) were to check the dynamical stability of regular solitons, cuspons, and peakons in the case $`\sigma _3=0`$, and the existence and stability of peakons in the more general case, $`\sigma _30`$. Both finite-difference and pseudo-spectral numerical methods have been used, in order to check that the same results are obtained by methods of both types. We used semi-implicit Crank-Nicholson schemes, in which the nonlinear terms were treated by means of the Adams-Bashforth method.

The presence of singularities required a careful treatment of cuspon and peakon solutions. To avoid numerical instabilities due to discontinuities, we found it, sometimes, beneficial to add an artificial weak high-wavenumber viscosity to the pseudospectral code. This was done by adding linear damping terms to the right-hand side of Eqs. (4), (5) and (6), which have the form $`i\nu (k)k^2\widehat{u_n}`$ in the Fourier representation, where $`\widehat{u_n}`$ is the Fourier transform of $`u_n`$ ($`n=1,2,3`$). The high-pass filter viscosity $`\nu (k)`$ suppresses only high wavenumbers and does not act on others. In particular, we chose

$$\nu (k)=\{\begin{array}{cc}0,\hfill & \mathrm{if}|k|<\frac{5}{16}K\hfill \\ \eta (\frac{16|k|}{K}5),\hfill & \mathrm{if}\frac{5}{16}K<|k|<\frac{3}{8}K\hfill \\ \eta ,\hfill & \mathrm{if}|k|>\frac{3}{8}K\hfill \end{array},$$

where $`K`$ is the largest wavenumber in the actual numerical scheme, and $`\eta `$ is a small viscosity coefficient. We have found that $`\eta 10^5`$ was sufficient to avoid Gibbs’ phenomenon in long-time simulations.

When instabilities occur at a singular point (cusp or peak), it is hard to determine whether the instability is a real one, or a numerical artifact. Therefore, we checked the results by means of a finite-difference code which used an adaptive staggered grid; motivated by the analysis of the vicinity of the point $`x=0`$ reported above, we introduced the variable $`\xi x^{2/3}`$ to define an adaptive grid, and also redefined $`U_3\sqrt{\xi }\stackrel{~}{U}_3`$. In these variables, the cusp becomes a regular point. This approach was solely used to check the possible occurrence of numerical instabilities.

In the following subsections we present typical examples of the numerical results for both cases considered above, viz., $`\sigma _3=0`$ and $`\sigma _3<0`$, when, respectively, cuspons and peakons are expected.

### 3.2 The case $`\sigma _3=0`$

First, we report results obtained for the stability of regular solitary waves in the case $`\sigma _3=0`$. As initial configurations, we used the corresponding stationary solutions to Eqs. (21) and (22). To test the stability of the regular solitary waves, we added small perturbations to them. As could be anticipated, the regular solitary wave sheds off a small dispersive wave train and relaxes to a stationary soliton, see Fig. 5 (for a more detailed illustration of the generation of small radiated waves by a soliton, see also Fig. 8 below). If, however, a regular soliton is taken as an initial condition for parameter values inside, but close to the border of the cuspon region, it does not become unstable in this slightly modified section of parameter space (which only supports cuspons), but instead this soliton exhibits persistent internal vibrations, see an example in Fig. 6. These and many other simulations clearly show that the regular soliton is always stable, and, close to the parameter border with cuspons, it has a persistent internal mode.

It was shown analytically above that Eqs. (17) and (18) ( with $`\sigma _3=0`$) support peakons when $`\rho _0>0`$ and $`\rho _0\kappa ^2`$ is very small. Direct simulations show that peakons do exist in this case, and they may be either unstable or stable. In the case when they are unstable, a high-wavenumber instability develops around the central peak. In Fig. 7, we display the time evolution of a typical stable peakon.

Next, we look at what happens if we take a regular soliton as an initial condition in a section of the parameter space which supports only stable peakons. This enables us to study the competition of the structural stability of regular solitons (as confirmed in Fig. 6) and the stability of peakons. The initial condition is taken as a stationary regular soliton in the parameter region (close to the boundary of the peakon region) with $`\rho _0<0`$, whereas the simulations are run for values of the parameters corresponding to $`\rho _0>0`$, which only admits peakons and excludes regular solitons. Unlike the case shown in Fig. 6, the time evolution now does not exhibit internal vibrations. Instead, the pulse slowly decays into radiation. This outcome can be explained by the fact that the peakon’s norm (see Eq.(8)) turns out to be larger than that of the initial pulse in this case, hence its rearrangement into a stable peakon is not possible. An essential result revealed by the simulations is that cuspons may also be stable, a typical example being displayed in Fig. 8. In this figure, one can see a small shock wave which is initially generated at the cuspon’s crest. It seems plausible that this shock wave is generated by some initial perturbation which could be a result of the finite mesh size in the finite-difference numerical scheme employed for the simulations. In fact, the emission of a small-amplitude shock wave is quite a typical way of the relaxation of both cuspons and peakons to a final stable state. To make it sure that the shock wave is not an artifact generated by the numerical scheme, we have checked that its shape does not change with the increase of the numerical accuracy.

To further test the stability of the cuspons and peakons, in many cases we allowed the initially generated shock wave to re-enter the integration domain (due to periodic boundary conditions used in the simulations) and interact again with the cuspon or peakon. As a result, the stability of the solitons of these types has been additionally confirmed. An example of the spatial profile of the cuspon established after a long evolution is shown in Fig. 9. Both the stability of the cuspon, and the presence of a tiny shock wave are evident in the figure.

However, unlike the regular solitons, which were found to be always stable, the cuspons are sometimes unstable. Typically, their instability triggers the onset of spatiotemporal collapse, i.e., formation of a singularity in a finite time (see a discussion of the feasible collapse in systems of the present type, given in the Introduction). Simulations of the collapse were possible with the use of an adaptive grid. A typical example of the collapse is shown in Fig. 10, where the inset shows that (within the numerical accuracy available) the amplitude of the collapsing pulse indeed diverges in a finite time.

However, the collapse is not the only possible outcome of the instability. In some other cases, which are not displayed here, the instability of peakons could be quite weak, giving rise to their rearrangement into regular solitons by the shedding of a small amount of radiation.

### 3.3 The case $`\sigma _30`$

The predictions of the analysis developed above for the most general case, when the SPM terms are present in the model ($`\sigma _30`$), were also checked against direct simulations. As a result, we have found, in accord with the predictions, that only regular solitons exist in the case $`\sigma _3>0`$, while in the case $`\sigma _3<0`$ both regular solitons and peakons have been found as generic solutions. Further simulations, details of which are not shown here, demonstrate that both regular solitons and peakons are stable in this case.

## 4 Conclusion

In this work, we have introduced a generic model of three waves coupled by linear and nonlinear terms, which describes a situation when three dispersion curves are close to an intersection at one point. The model was cast into the form of a system of two waves with opposite group velocities that, by itself, gives rise to the usual gap solitons, which is further coupled to a third wave with zero group velocity (in the laboratory reference frame). Situations of this type are quite generic, being realizable in various models of nonlinear optics, density-stratified fluid flows, and in other physical contexts. In particular, two versions (temporal and spatial) of a nonlinear-optical model, which is based on a waveguide carrying the triple spatial Bragg grating, have been elaborated in the Introduction. Our consideration was focussed on zero-velocity solitons. In a special case when the self-phase modulation (SPM) is absent in the equation for the third wave, soliton solutions were found in an exact form. It was shown that there are two coexisting generic families of solitons: regular solitons and cuspons. In the special case when the coefficient of the linear coupling between the first two waves and the third one vanishes, cuspons are replaced by peakons. Direct simulations have demonstrated that the regular solitons are stable (in the case when the regular soliton is close to the border of the cuspon region, it has a persistent internal mode). The cuspons and peakons may be both stable and unstable. If they are unstable, they either shed off some radiation and rearrange themselves into regular solitons, or, in most typical cases, the development of the cuspon’s instability initiates onset of spatiotemporal collapse. Actually, the present system gives the first explicit example of collapse in one-dimensional gap-soliton models.

The most general version of the model, which includes the SPM term in the equation for the third wave, has also been considered. Analysis shows that cuspons cannot exist in this case, i.e., cuspons, although being possibly dynamically stable, are structurally unstable. However, depending on the signs of the SPM coefficient, and some combination of the system’s parameters, it was shown that a generic family of peakon solutions may exist instead. In accord with this prediction, the peakons have been found in direct simulations. The peakons, as well as the regular solitons, are stable in the system including the SPM term. We stress that peakons are physical solutions, as they have all their field components and their first derivatives finite.

The next step in the study of this system should be consideration of moving solitons, which is suggested by the well-known fact that the usual two-wave model gives rise to moving gap solitons too . However, in contrast to the two-wave system, one may expect a drastic difference between the zero-velocity and moving solitons in the present three-wave model. This is due to the reappearance of a derivative term in Eq. (6), when it is written for a moving soliton, hence solitons which assume a singularity or jump in the $`U_3`$ component, i.e., both cuspons and peakons, cannot then exist. Nevertheless, one may expect that slowly moving solitons will have approximately the same form as the cuspons and peakons, with the singularity at the central point replaced by a narrow transient layer with a large gradient of the $`U_3`$ field. Detailed analysis of the moving solitons is, however, beyond the scope of this work.

## Acknowledgements

We would like to thank Tom Bridges, Gianne Derks and Sebastian Reich for valuable discussions. B.A.M. and G.A.G. appreciate hospitality of the University of Loughborough (UK). The work of G.A.G. is supported by a European Commission Grant, contract number HPRN-CT-2000-00113, for the Research Training Network Mechanics and Symmetry in Europe (MASIE)

Fig. 1. A schematic representation of the optical model that gives rise to the three linearly coupled waves in a $`\chi ^{(3)}`$ waveguide. The figure shows either the planar waveguide with the triple Bragg grating, in the case of the temporal evolution of the fields, or the transverse cross section of the bulk medium of the layered (photonic-crystal-fiber) type, in the case of the spatial evolution along the coordinate $`z`$ perpendicular to the plane of the figure. The triangles formed by bold arrows illustrate how the linear couplings between the three waves, whose Poynting vectors are represented by the arrows, are induced by the Bragg reflections on the three gratings. The difference between the gratings represented by the continuous and dashed lines is in their strength (refractive-index contrast).

Fig. 2. Dispersion curves produced by Eq. (14) in the case $`\kappa =0.5`$: (a) $`\omega _0<1\kappa ^2`$; (b) $`\omega _0>1`$. The dashed line in each panel is $`\omega =\omega _0`$. The case with $`1\kappa ^2<\omega _0<1`$ is similar to the case (a) but with the points $`\omega _+`$ and $`1`$ at $`k=0`$ interchanged.

Fig. 3. The shape of the cuspon for $`\alpha =2.0`$, $`\omega _0=0.1`$, $`\omega =0.5`$, and (a) $`\kappa =0.5`$, i.e., in the general case, and (b) $`\kappa =0.1`$, i.e., for small $`\kappa `$. In the case (b) we also show the usual gap soliton (by the dashed line), the part of which above the critical value $`S=S_{\mathrm{cr}}`$ (shown by the dotted line) should be removed and the remaining parts brought together to form the peakon corresponding to $`\rho _0\kappa ^20`$.

Fig. 4. The shape of the peakon for the case when $`\sigma _3<0`$, where we plot $`|U_1|^2`$. The parameters are $`\sigma _3=0.01`$, $`\kappa =0.1`$, $`\alpha =2.0`$, $`\omega _0=0.1`$, and $`\omega =0.5`$. In this case, $`\rho _0=4.8`$.

Fig. 5. The shape of an initially perturbed regular soliton in the case $`\sigma _3=0`$ at $`t=5`$, which illustrates the stabilization of the soliton through the shedding of small-amplitude radiated waves. The plot displayed is the field $`\mathrm{Re}U_1(x)`$. The parameters are $`\kappa =0.01,\alpha =1.0,\omega _0=0.2`$, and $`\omega =0.9`$.

Fig. 6. Internal vibrations of an initially-perturbed regular soliton, which was taken close to the parameter boundary of the cuspon region. The plot shows the squared amplitude $`a\left|U_1(x=0)\right|^2`$ of the $`U_1(x)`$ field versus time. The parameters are $`\kappa =0.01,\alpha =1.9,\omega _0=1.5`$, and $`\omega =0.5`$, with $`\rho _0=0.095`$ (see Eq. (35).

Fig. 7. An example of a stable peakon. The plot shows the field $`\left|U_1\right|^2`$ versus $`x`$ and $`t`$. The parameters are $`\kappa =1.0,\alpha =1.95,\omega _0=1.5`$, and $`\omega =0.5`$, with $`\rho _0=0.04875`$.

Fig. 8. An example of a stable cuspon. The plot shows the field $`\left|U_1\right|^2`$ versus $`x`$ and $`t`$. The parameters are $`\kappa =1.0,\alpha =1.0,\omega _0=1.5`$, and $`\omega =0.5`$, with $`\rho _0=0.5`$.

Fig. 9. The spatial profile of the stable cuspon at $`t=20`$. The parameters are the same as in Fig. 8.

Fig. 10. The spatial profile is shown for an unstable cuspon in terms of $`\mathrm{Im}U_1`$ at $`t=10^3`$. The inset depicts the time-evolution of the maximum value of $`\left|U_1\right|^2`$. The transition to collapse is clearly seen as an explosive temporal behavior of the amplitude. The parameters are $`\kappa =0.01,\alpha =1.1,\omega _0=0.1`$, and $`\omega =0.3`$, with $`\rho _0=2.618`$.


format: text

text:
# Effects of pressure on the superconducting properties of magnesium diboride

## I INTRODUCTION

The recent discovery of superconductivity in MgB2 (Ref. naga, ) has attracted considerable interest in the study of this material, both to understand the mechanism of superconductivity and to explore other properties of MgB2 and related materials. The high transition temperature $`T_c`$ 40 K in this material offers other possibility for finding high-$`T_c`$ superconductivity in some binary intermetallic compounds besides cuprates and $`C_{60}`$-based compounds. Meanwhile, the high critical currents observed in MgB2 thin films eom  and wiresjin  reveal that MgB2 belongs to a new class of low-cost, high-performance superconducting materials for magnets and electronic applications.

Measurements of the isotope effect and of the influence of pressure on the transition temperature and critical field of superconductors yield information on the interaction causing superconductivity. Indeed, the pressure (or volume) and the mass number would seem to be the only variables whose effect might be capable of immediate theoretical interpretation. By observing how pressure changes the parameters of the lattice in the normal state and in the superconducting state, and comparing the measurements with the theoretical predictions, one can test the validity of some theoretical models. Olsen $`et`$ $`al.`$ olse1  have shown that the volume (V) dependence of effective interaction $`N(E_F)v`$, $`d\mathrm{ln}N(E_F)v/d\mathrm{ln}V`$, can be scaled well with the deviation $`\xi `$ from the full isotope effect where $`\xi `$ is defined by $`T_cM^{0.5(1\xi )}`$ in superconducting metals. Bud’ko $`et`$ $`al.`$ bud  and Hinks $`et`$ $`al.`$ hink  reported a sizeable isotope effect for B ($`\alpha _B`$=0.26(3) or 0.30(1)) in newly discovered superconductor MgB2. Although the total isotope coefficient $`\alpha `$=0.32(1) (Ref. hink, ) is smaller than the canonical BCS value of 0.5, it is the same as that in Cd (Ref. palm, ). The isotope effect along with other measurements such as inelastic neutron scattering,sst ; osbo  tunneling,kara  NMR (Ref. kote, ), and specific heatkrem ; walt ; bouq  confirmed that MgB2 is an electron-phonon mediated $`s`$-wave superconductor.

Soon after the discovery of superconductivity in MgB2, the effect of pressure on $`T_c`$ was studied by two groupsmont ; lore1  by resistivity or $`ac`$ susceptibility measurements. Both groups observed a decrease of $`T_c`$ with increasing pressure, with initial pressure derivative $`dT_c/dP`$ of –0.8 K/GPa (Ref. mont, ) or –1.6 K/GPa (Ref. lore1, ), respectively. Moreover, Monteverde $`et`$ $`al.`$ mont  found that the superconductivity is not destroyed by applying high pressure up to 25 GPa, at which point $`T_c`$ is as high as 21 K. A somewhat larger $`dT_c/dP`$ of –2.0 K/GPa was late reported by Saito $`et`$ $`al.`$sait  from high-pressure resistivity measurements. Using a He-gas apparatus, Tomita $`et`$ $`al.`$tomi  determined a $`dT_c/dP`$ of –1.11 K/GPa under pure hydrostatic pressure conditions. In order to find the reason why the reported values of $`dT_c/dP`$ are different among different groups, Lorenz $`et`$ $`al.`$ lore2  carried out high pressure experiments on MgB2 samples with different $`T_c`$’s at ambient pressure and different pressure media. $`T_c`$ was found to decrease linearly over the whole pressure range (0–1 GPa). In the He environment, the two samples with the initial $`T_c`$=39.2 and 37.4 K yield the pressure derivatives of –1.07 and –1.45 K/GPa, respectively. The former is obviously very close to that of Tomita $`et`$ $`al.`$tomi  The latter approaches their previous data,lore1  which was obtained by using the Fluorinert FC77 as pressure medium. They therefore concluded that the variation in the value of $`dT_c/dP`$ by various groups results from the differences in sample preparation conditions. The value of $`dT_c/dP`$–1.1 K/GPa is then confirmed to give the true hydrostatic pressure dependence of $`T_c`$ in MgB2.

Two theoretical models have been tried to describe the systematics of the behavior of $`T_c`$ under pressure in MgB2. Based on the theory of hole superconductivity, Hirsch hirs  predicted an increased $`T_c`$ with the decrease of B-B intraplane distance under the application of in-plane biaxial pressure. However, this prediction has not been confirmed experimentally yet. No uniaxial pressure measurement was reported due to the extreme difficulty in growing MgB2 single crystal. The experiments of hydrostatic pressure effect on $`T_c`$ do not particularly support this theory provided that no charge transfer between the Mg and B layers occurs. Alternatively, the experimental results have been analyzed lore1 ; tomi ; loa ; vogt  by using the McMillan formula mcmi  derived from Eliashberg theory,elia  supporting electron-phonon mediated superconductivity. Interestingly, Loa and Syassen loa  analyzed the pressure effect on $`T_c`$ from McMillan formula on the basis of their calculated elastic and electronic structure data. Assuming that the electron-ion matrix element $`I`$ is pressure independent, they found that the pressure effect on $`T_c`$ is in good agreement with experimental data by using a lattice Grüneisen parameter $`\gamma _G`$=1. These assumptions, however, deserve some refinements. Recent band structure calculations suggest that MgB2 is a traditional $`sp`$ metal superconductor.kort ; pick ; kong  The pressure dependence of $`I`$ has long been an interesting issue of the research of pressure effect in simple $`sp`$ metals.hodd ; trof ; coul ; rott ; daco  Ziman’s calculation of the electron-phonon interaction leads to $`<I^2>N(E_F)^2`$ at least in the limit of long wavelengths.zima  This then indicated that the consideration of the pressure dependence of $`I`$ would be important for better understanding the superconducting properties of MgB2 under pressure. On the other hand, it has been foundtomi  that the choice of lattice Grüneisen parameter $`\gamma _G`$ is crucially important in explaining both the magnitude and the sign of the pressure derivative of $`T_c`$ when using McMillan formula. The value of $`\gamma _G`$=1 in the calculation of Loa and Syassen is obviously lower than those reported recently.roun ; gonc  The pressure dependence of the effective electron-electron Coulomb repulsion $`\mu ^{}`$ appearing in the McMillan formula is usually neglected in previous studies due to the assumption of the small change of $`\mu ^{}`$ compared with that of the electron-phonon coupling parameter $`\lambda `$ (Ref. seid, ). However, the magnitude of $`\mu ^{}`$ is also of interest in connection with the possibility that superconductivity may be destroyed by pressure.seid ; bo  It was argued that the pressure dependence of $`\mu ^{}`$ makes a significant contribution to the behavior of $`T_c`$ under very high pressures and must be handled carefully.daco ; gubs ; papa

In this paper we discuss the pressure dependences of some interested superconducting properties in MgB2. The outline of this paper is as follows: In Sec. II we presented a theoretical approach for pressure effects on the superconducting properties in the simple $`sp`$ metals superconductors. Section III contains the theoretical results obtained and a comparison with experiments for MgB2. We draw conclusions in Sec. IV.

## II THEORETICAL FORMULATION

For our purposes, the relation between $`T_c`$ and microscopic parameters is given adequately by the McMillan equationmcmi

$$T_c=\frac{\mathrm{\Theta }_D}{1.45}\mathrm{exp}\left[\frac{1.04(1+\lambda )}{\lambda \mu ^{}(1+0.62\lambda )}\right],$$
(1)
which relates $`T_c`$ to the electron-phonon coupling parameter $`\lambda `$, the Coulomb repulsion strength $`\mu ^{}`$, and a temperature $`\mathrm{\Theta }_D`$ characteristic of the phonons.

Considering the variations of $`\mathrm{\Theta }_D`$, $`\lambda `$, and $`\mu ^{}`$ with pressure or volume and introducing parameters $`\phi =\mathrm{ln}\lambda /\mathrm{ln}V`$ and $`\varphi =\mathrm{ln}\mu ^{}/\mathrm{ln}V`$, we can get the pressure coefficient of $`T_c`$

$`{\displaystyle \frac{d\mathrm{ln}T_c}{dP}}={\displaystyle \frac{\gamma _G}{B_0}}{\displaystyle \frac{1.04\lambda (1+0.38\mu ^{})}{[\lambda \mu ^{}(1+0.62\lambda )]^2}}{\displaystyle \frac{\phi }{B_0}}`$
$`+{\displaystyle \frac{1.04\mu ^{}(1+\lambda )(1+0.62\lambda )}{[\lambda \mu ^{}(1+0.62\lambda )]^2}}{\displaystyle \frac{\varphi }{B_0}},`$ (2)
where $`B_01/\kappa _V=P/\mathrm{ln}V`$ is the bulk modulus and $`\mathrm{\Theta }_D`$ is assumed to be proportional to $`<\omega ^2>^{1/2}`$ and $`\gamma _G=\mathrm{ln}<\omega ^2>^{1/2}/\mathrm{ln}V`$ being the effective Grüneisen parameter.

It is well known that the usual BCS result for the energy gap can be expressed bybcs

$$\mathrm{\Delta }_0=2\mathrm{\Theta }_D\mathrm{exp}\left[\frac{1}{N(E_F)v}\right].$$
(3)
Where $`N(E_F)`$ is the electronic density of states at the Fermi energy $`E_F`$ and $`v`$ is the pairing potential arising from the electron-phonon interaction. If we renormalize the Morel-Anderson result moan  by introducing the renormalization parameter $`Z_n(0)1+\lambda `$ into their analysis, the effective interaction strength $`N(E_F)v`$ can be rewritten asleav

$$N(E_F)v=\frac{\lambda \mu ^{}}{1+\lambda }.$$
(4)
The logarithmic volume derivative of $`N(E_F)v`$ is then given by

$$\frac{d\mathrm{ln}N(E_F)v}{d\mathrm{ln}V}=\frac{\lambda (1+\mu ^{})}{(\lambda \mu ^{})(1+\lambda )}\phi \frac{\mu ^{}}{\lambda \mu ^{}}\varphi .$$
(5)
Considering the experimental observations of the pressure dependence of the energy gap of superconductor,jpf  we differentiate Eq. (3) with respect to pressure

$$\frac{d\mathrm{ln}\mathrm{\Delta }_0}{dP}=\frac{\gamma _G}{B_0}\frac{1}{B_0}\left[\frac{\lambda (1+\mu ^{})}{(\lambda \mu ^{})^2}\phi \frac{\mu ^{}(1+\lambda )}{(\lambda \mu ^{})^2}\varphi \right].$$
(6)
BCS expression for the critical field $`H_c`$ at absolute zero temperature isbcs

$$\frac{H_c(0)^2}{8\pi }=2N(E_F)\omega ^2\mathrm{exp}\left[\frac{2}{N(E_F)v}\right].$$
(7)
Differentiating Eq. (7) with respect to the pressure, one obtains an expression of the pressure coefficient of $`H_c(0)`$,

$$\frac{d\mathrm{ln}H_c(0)}{dP}=\frac{d\mathrm{ln}\mathrm{\Delta }_0}{dP}\frac{\gamma _N}{2B_0}+\frac{1}{2B_0},$$
(8)
where $`\gamma _N=\mathrm{ln}N(E_F)/\mathrm{ln}V`$.

The expressions for $`\gamma _G`$, $`\phi `$, and $`\varphi `$ can be integrated to give

$`\mathrm{\Theta }_D(V)`$ $`=`$ $`\mathrm{\Theta }_D(0)\left[V/V_0\right]^{\gamma _G}`$ (9)
$`\lambda (V)`$ $`=`$ $`\lambda (0)\left[V/V_0\right]^\phi `$
$`\mu ^{}(V)`$ $`=`$ $`\mu ^{}(0)\left[V/V_0\right]^\varphi .`$
Here $`V`$ and $`V_0`$ are the unit cell volumes under the applied pressure and at ambient pressure, respectively. These two volumes can be related according to the first-order Murnaghan equation of state $`V(P)=V(0)(1+B_0^{}P/B_0)^{1/B_0^{}}`$. The Eq. (9) is then rewritten as

$`\mathrm{\Theta }_D(P)`$ $`=`$ $`\mathrm{\Theta }_D(0)\left[1+{\displaystyle \frac{B_0^{}P}{B_0}}\right]^{\gamma _G/B_0^{}}`$ (10)
$`\lambda (P)`$ $`=`$ $`\lambda (0)\left[1+{\displaystyle \frac{B_0^{}P}{B_0}}\right]^{\phi /B_0^{}}`$
$`\mu ^{}(P)`$ $`=`$ $`\mu ^{}(0)\left[1+{\displaystyle \frac{B_0^{}P}{B_0}}\right]^{\varphi /B_0^{}}.`$
From Eqs. (1) and (10) we arrive at the expression for the pressure dependence of $`T_c`$

$$T_c(P)=T_c[\mathrm{\Theta }_D(P),\lambda (P),\mu ^{}(P)].$$
(11)
Knowing $`B_0`$, $`B_0^{}`$, $`\gamma _G`$, $`\gamma _N`$, $`\varphi `$, and $`\phi `$, one can evaluate the pressure effects on the superconducting properties, especially the behavior of $`T_c`$ under pressure. $`B_0`$ and $`B_0^{}`$ can be obtained from the compressibility data determined by neutron or synchrotron x-ray diffractions. A direct experimental determination of $`\gamma _G`$ can be made by measuring electron tunneling keel ; yama ; zava  or inelastic neutron scatteringlech  under high pressure. In general, for metals in which different techniques yield similar Grüneisen constants, a good approximation to $`\gamma _G`$ is provided by the room temperature value determined from the Grüneisen equation

$$\gamma _G=\frac{\alpha _VV_m}{\kappa _VC_p},$$
(12)
where $`\alpha _V`$ is the volume coefficient of thermal expansion, $`V_m`$ is the molar volume, and $`C_p`$ is the molar heat capacity at constant pressure. The approximation for $`\gamma _G`$ of Slater is derived from the pressure derivative of the bulk modulus slat

$$\gamma _G^S\frac{B_0^{}}{2}\frac{1}{6}=\frac{2}{3}\frac{1}{2}\frac{V^2P/V^2}{P/V}.$$
(13)
The formula for $`\mu ^{}`$ due to Morel and Anderson moan  used here is

$$\mu ^{}=\frac{\mu }{1+\mu \mathrm{ln}(E_F/\omega _{ph})},$$
(14)
with $`\mu =0.5\mathrm{ln}[(1+a^2)/a^2]`$ and $`a^2=\pi e^2N(E_F)/k_F^2`$, from which we evaluate the volume dependence of $`\mu ^{}`$ as

$$\varphi =\mu ^{}\left[\frac{2}{3}\gamma _G\frac{1e^{2\mu }}{2\mu ^2}(\gamma _N+\frac{2}{3})\right].$$
(15)
Here the variation of $`k_F`$ with volume has been calculated from the fundamental definition $`k_F=(3\pi ^2Z/V)^{1/3}`$ with $`Z`$ the valency. Unfortunately, to the best of our knowledge $`\gamma _N`$ has never measured directly for any superconductor in the case of free electron gas it would have a value of $`2/3`$. Using the expression given by Migdal migdal  for the electronic specific heat coefficient $`\gamma `$ one obtains for the electronic Grüneisen parameter

$$\gamma _e=\frac{\mathrm{ln}\gamma }{\mathrm{ln}V}=\gamma _N+\frac{\lambda }{1+\lambda }\phi .$$
(16)
The electronic Grüneisen parameter $`\gamma _e`$ is usually deduced from measurements through the simple relation barr

$$\gamma _e=\frac{\alpha _eV_m}{\kappa _VC_e}.$$
(17)
Here $`\alpha _e`$ is the contribution to the expansion coefficient from the electrons at lower temperatures and $`C_e`$ is the electronic heat capacity. A theoretical estimate of $`\gamma _e`$ can also be given from the measurement of the volume dependence of the orbital susceptibility fawc  or from band-structure considerations.flet

The electron-phonon coupling parameter $`\lambda `$ can be expressed as

$$\lambda =\frac{N(E_F)}{M<\omega ^2>}\frac{\eta }{M<\omega ^2>},$$
(18)
where $`<I^2>`$ is the mean-square electron-ion matrix element and $`M`$ the ionic mass. The McMillan-Hopfield parameter $`\eta `$ (or $`N(E_F)<I^2>`$) has been regarded as a local “chemical” property of an atom in a crystal. Allen and Dynes alle  pointed out that $`\eta `$ is the most significant single parameter in understanding the origin of the high $`T_c`$ of conventional superconductors. For strong coupling systems, variation in $`\eta `$ is more important than variation of $`<\omega ^2>`$ in causing $`T_c`$ to change. Softening of $`<\omega ^2>`$ often does enhance $`T_c`$, but very high $`T_c`$ should be caused more by large $`\eta `$ than by small $`<\omega ^2>`$.

The logarithmic volume derivative of $`\lambda `$, $`\phi `$, is then obtained

$$\phi =\frac{\mathrm{ln}\eta }{\mathrm{ln}V}+2\gamma _GS+2\gamma _G.$$
(19)
In order to understand how the electronic contribution $`\eta =N(E_F)<I^2>`$ varies with volume, we use the Gaspari-Gyoriffy theory gasp  for $`\eta `$, i.e.,

$$\eta =\frac{k_F^2}{\pi ^2N(E_F)}\underset{l}{}\frac{2(l+1)\mathrm{sin}^2(\delta _{l+1}\delta _l)N_lN_{l+1}}{N_l^1N_{l+1}^1},$$
(20)
where $`N_l`$ is the $`l`$th angular momentum component of the density of states, $`N_l^1`$ is the $`l`$th component of the single scatterer density of states evaluated at $`E_F`$, and $`\delta _l`$ the phase shift.

For simple metals, the scatterers are assumed to be weak. We can take $`N_l=N_l^1`$ and approximate $`\mathrm{sin}^2(\delta _{l+1}\delta _l)`$ by $`(\delta _{l+1}\delta _l)^2`$. Eq. (20) is rewritten as

$$\eta =\frac{k_F^2}{\pi ^2N(E_F)}\underset{l}{}2(l+1)(\delta _{l+1}\delta _l)^2.$$
(21)
This expression is identical to the pseudopotential formula of McMillan.mcmi ; evan  Assuming that the phase shift $`\delta _l`$ does not vary very much under pressure for simple $`sp`$ superconductors, we then obtain

$$S=\gamma _N\frac{2}{3}.$$
(22)
The form in Eq. (22) is the same as that of Baryakhtar and Makarov,bary  who used the constant of the electron-phonon interaction of Fröhlich and Mitra.floh  The expression is an improvement over the expressions of $`S`$=0 and $`S=`$4/3 obtained by Olsen $`et`$ $`al.`$ olsen1  and Seiden,seid  respectively. It is interesting to notice that substitution of $`\gamma _N`$=2/3 into Eq. (22) yields $`S=4/3`$. Eq. (22) reduces to the expression of Seiden,seid  who modified McMillan’s expression for $`\lambda `$ somewhat by considering the effects of a real lattice spectrum as opposed to the Jellium model. Since the electronic Grüneisen parameter $`\gamma _e`$ usually varies among different metals even in the simple non-transition elements,palm ; barr ; flet  we believe that Eq. (22) should provide a more reasonable value of $`S`$ compared with Seiden’s formula.

## III RESULTS AND DISCUSSION

Using the experimental value of $`T_c`$=39.25 K (Refs. tomi, ; lore2, ), and the theoretical estimates of $`\lambda `$=0.87 and $`\mu ^{}`$=0.10 (Ref. kong, ), we got $`\mathrm{\Theta }_D`$=860 K from Eq. (1) for MgB2. We believe that all these parameters, which will enter our calculations, are reliable. For example, the inelastic neutron scattering measurements osbo  provide an estimate of $`\lambda `$0.9, which is close to that we used. The calculated value of $`\mathrm{\Theta }_D`$=860 K is in the range from 746 to 1050 K determined from the specific heat measurements.krem ; walt ; bouq

We took the structural parameters $`B_0`$=147.2 GPa and $`B_0^{}`$=4 from the measurements under the pure hydrostatic pressures up to 0.62 GPa (Ref. jorg, ) and under high pressures up to 15 GPa (Refs. vogt, ; gonc, ), respectively. To our knowledge no inelastic neutron scattering or tunneling data exist for MgB2 under hydrostatic pressure. We have to use Eq. (12) or (13) for estimating the lattice Grüneisen parameter $`\gamma _G`$. The measurements of heat capacityswif  give a $`C_p`$ of 47.80 J/(K mol) at $`T`$=298.16 K. $`V_m`$=1.75$`\times 10^5`$ m3/mol, $`\kappa _V`$=6.79$`\times 10^{12}`$ Pa-1, and $`\alpha _V`$=2.22$`\times 10^5`$ K-1 can be drawn from the neutron diffraction data. jorg  We therefore obtained $`\gamma _G`$=1.2 by using Eq. (12). Based on the first-order Murnaghan equation for $`V(P)`$ and the Slater expression of Eq. (13), we got a somewhat larger $`\gamma _G`$ of 1.83 compared to that from Eq. (12). For most simple metals, there is no much difference between the room temperature lattice Grüneisen parameter given through Eq. (12) and the Slater relation.seid ; gsch  It was found that the Slater expression usually can yield the reasonable values of $`\gamma _G`$ for most metals.moru  The only uncertainty entering Eq. (12) in our calculation comes from the indirect measurements of the linear coefficients of thermal expansion.jorg  Roundy $`et`$ $`al.`$ roun  reported a value of $`\gamma _G`$2.3 from $`ab`$ $`initio`$ calculations, which is close to our calculated $`\gamma _G^S`$ according to Eq. (13). Meanwhile, Goncharov $`et`$ $`al.`$ gonc  determined a large $`E_{2g}`$ mode Grüneisen parameter of 2.9$`\pm `$0.3 from the measurements of Raman spectra under pressure. This value is obviously larger than those derived from Eqs. (12) and (13).

In the calculations of electronic density of states, Loa and Syassanloa  found that $`N(E_F)`$ decreases with pressure at a rate of $`d\mathrm{ln}N(E_F)/dP=`$–3.1$`\times `$10-3 GPa-1. Combining this calculated value and experimental value of $`B_0`$, we got $`\gamma _N`$=0.46. The volume dependence of $`\mu ^{}`$ is then derived from (15) once having the values of $`\gamma _G`$ and $`\mu `$. For simple $`sp`$ metals, $`a^2`$ has a typical value of 0.4 (Ref. moan, ), which yields $`\mu =`$0.63. The volume dependence of $`\lambda `$, $`\phi `$, is therefore determined from Eqs. (19) and (22). Using $`\gamma _N`$=0.46, we have $`S`$=–1.13 for MgB2, which is smaller in magnitude than –4/3 in Seiden’s formula for simple metals.seid  For the transition metals, Hopfield hopf  commented that $`S`$ is a relatively constant quantity with a value of about –3.5. The values of $`S`$=–3.5$``$–3.1 obtained by inverting the measured $`dT_c/dP`$ for YNi2-xMxB2C (M=Co and Cu) (Ref. loon, ) are comparable to that of the transition metals, but are larger in magnitude compared to that of MgB2.

With the parameters determined above, we have calculated the pressure derivatives of $`T_c`$ for MgB2 by using Eq. (2). In Fig. 1 we plotted $`dT_c/dP`$ as a function of $`\gamma _G`$ in the interested range. It is interesting to note that $`\gamma _G`$ plays a predominant role for the pressure effect of $`T_c`$. For the four different $`\gamma _G`$’s considered here $`dT_c/dP`$ are negative. The values obtained from $`\gamma _G`$=1.83 and 2.3 are –0.78 and –1.12 K/GPa, respectively. These are close to the hydrostatic pressure value of –1.1 K/GPa (Refs. tomi, ; lore2, ). Thus the hydrostatic pressure results can be reproduced in terms of our present model by using the values of $`\gamma _G`$ obtained from either the Slater relation or $`ab`$ $`initio`$ calculation. It is difficult to obtain the measured results by using $`\gamma _G`$=1 as suggested by Loa and Syassen.loa  We noticed that a $`\gamma _G`$ of 2.27 is necessary so as to account for the pressure effect on $`T_c`$ for MgB2. As emphasized above, all quantities entering Eq. (12) are experimental values and only $`\alpha _V`$ was taken from indirect measurements. Thus it is highly expected to operate the thermal expansion measurement to yield a direct $`\alpha _V`$. The present results indicate that the range from $`\gamma _G`$=1.83 to 2.3 should cover the reasonable choices for the lattice Grüneisen parameters.

To verify these results, and also to study the behavior of $`T_c`$ as a function of pressure, we have performed explicit calculation based on Eq. (11). The theoretical results in the pressure range from 0 to 1.0 GPa are shown in Fig. 2. The experimental data points of Tomita $`et`$ $`al.`$ tomi  and Lorenz $`et`$ $`al.`$ lore2  measured under hydrostatic pressure conditions are also plotted for comparison. It is clearly seen that our calculations agree well with the experiments.

In Fig. 3 we presented the calculated results as well as the experimental data points of Monteverde $`et`$ $`al.`$ mont  and Deemyad $`et`$ $`al.`$deem  measured in the relatively high pressure region. Here we assume that phase transitions of all kinds do not occur under pressure range that we consider. We noticed that the experimental data points of Deemyad $`et`$ $`al.`$deem  and the sample 4 of Monteverde $`et`$ $`al.`$mont  are situated well between the two theoretical curves calculated by using $`\gamma _G`$=1.83 and 2.3, respectively. Interestingly, the agreement between our theoretical curve calculated by using $`\gamma _G`$=1.83 and the experimental data points of other samples of Monteverde $`et`$ $`al.`$mont  is seen to be reasonable, although there are some scatters among different samples and the reason is not clear. Furthermore, although the pressure measurements are limited to the region below 33 GPa, it is seen from the inset of Fig. 3 that Eq. (11) continues to describe the pressure dependence of $`T_c`$ as high as 100 GPa. Even at this point, the superconductivity is not destroyed by pressure in newly discovered superconductor MgB2. There was a discrepancy on whether pressure can destroy superconductivity.seid ; bo ; olsen2 ; smithc  However, our results support the conclusion of Olsen and collaboratorsbo ; olsen2  that the possibility of destruction of superconductivity by the application of sufficiently high pressure most likely does not exist. It follows from the comparison of Figs. 2 and 3 that the the pressure effect on $`T_c`$ indeed provides a support to the electron-phonon mediated superconductivity in MgB2.

In Fig. 4 we presented the normalized $`\lambda `$ and $`\mu ^{}`$ as a function of pressure up to 30 GPa, calculated from Eq. (10) by using $`\gamma _G`$=2.3. The Coulomb pseudopotential $`\mu ^{}`$ increases slightly with pressure. Whereas $`\lambda `$ changes significantly with pressure. The contribution from $`\mu ^{}`$(P) to the variation of $`T_c`$ with pressure is much less important than that of $`\lambda `$(P). Thus in the range from 0 to 30 GPa the pressure effect of $`T_c`$ for MgB2 is dominated by the competition of $`\lambda `$ and $`\mathrm{\Theta }_D`$ (or $`<\omega ^2>^{1/2}`$).

Table I contains the calculated values of pressure dependences of superconducting parameters for MgB2 from Eqs. (5), (6), (8) and (16) by using $`\gamma _G=1.83`$ and $`\gamma _G=2.3`$, respectively. The reliable values of $`\gamma _e`$ is readily determined using Eqs. (16), (19), and (22). We obtained $`\gamma _e`$=1.64 and 2.07 for MgB2, which are close to $`\gamma _e`$=1.7 for Pb and for Sn, $`\gamma _e`$=2.0 for Al (Ref. smithc, ). The negative sign for $`d\mathrm{ln}H_c(0)/dP`$ predicted for MgB2 is in agreement with the measurements for all simple elements superconductors with the exception of thallium.levyo  For the simple $`sp`$ metals superconductors, Rohrerrohrer  has demonstrated that $`d\mathrm{ln}N(E_F)v/d\mathrm{ln}V`$ must have approximately a value of 2.0. However, it was realized levyo  that the transition metals fail to show such simple behavior. Our estimated $`d\mathrm{ln}N(E_F)v/d\mathrm{ln}V`$=1.71 and 2.34 for MgB2 are comparable to those obtained for simple $`sp`$ metals superconductors.levyo ; rohrer  Early measurements for most simple metals trof ; jpf ; zava  show that there is difference between the quantities in $`d\mathrm{ln}\mathrm{\Delta }_0/dP`$ and $`d\mathrm{ln}T_c/dP`$. This can be understood with the aid of the results by Geǐlikman and Kresin,beili  that is, $`2\mathrm{\Delta }_0/k_BT_c=3.52[1+5.3(T_c/\omega _{ph})^2\mathrm{ln}\omega _{ph}/T_c]`$. The calculated data of MgB2 listed in Table I make it possible to support this theory. Since the phonon spectrum shifts under pressure, it follows that for all superconductors with $`2\mathrm{\Delta }_0/k_BT_c>3.52`$ a change of $`2\mathrm{\Delta }_0/k_BT_c`$ under pressure can be expected. It is interesting from the viewpoint of experiment to investigate the tunnel characteristics of MgB2 under hydrostatic pressure.

Figure 5 is a plot of $`d\mathrm{ln}N(E_F)v/d\mathrm{ln}V`$ versus the deviation $`\xi `$ from the full isotope effect for nine simple $`sp`$ metals superconductors as well as MgB2. The experimental values of $`d\mathrm{ln}N(E_F)v/d\mathrm{ln}V`$ for simple metals are chosen from the work of Olsen, Andres, and Geballe.olsen1  The experimental results for isotope effect exponent $`\alpha `$ and its deviation $`\xi `$ are taken from the works in Refs. \[ palm, ; maxw, ; hake, ; fass1, ; fass2, ; meser, \]. There are no experimental data available for Al and In now, we took the calculated results from Leavens and Carbotte.leav  We summarized these results for simple $`sp`$ metals in Table II. Olsen $`et`$ $`al.`$ olse1  suggested that $`d\mathrm{ln}N(E_F)v/d\mathrm{ln}V`$ is related to the isotope effect exponent $`\alpha `$ in metals superconductors. As seen from Fig. 5 the relation between $`d\mathrm{ln}N(E_F)v/d\mathrm{ln}V`$ and $`\xi `$ is not very clear when more data included. An isotope effect $`\alpha `$=0.32(1) (Ref. hink, ) in MgB2 is consistent with what appears to be a systematic variation of $`\alpha `$ across the non-transition elements. It is well known that deviations of the isotope effect exponent from $`1/2`$ are a measure of the relative strengths of the Coulomb and phonon-mediated electron-electron interactions. It is indicated, from the compared values of $`\xi `$ for MgB2 (Ref. hink, ) with Zn (Ref. fass1, ), Cd (Ref. palm, ), and Al (Ref. leav, ), that MgB2 should be a medium coupling superconductor.

## IV SUMMARY AND CONCLUSIONS

The major conclusions given by present investigation can be summarized as follows:

(i) A simple expression was derived for the pressure dependences of superconducting properties in simple $`sp`$ superconductor on the basis of McMillan equation. The logarithmic volume derivatives of $`\lambda `$, $`\mu ^{}`$, and $`\mathrm{\Theta }_D`$ can be self-consistently determined from experiments and theories. We gave an expression for $`\phi `$ from the theory of Gaspari and Gyorffy. gasp  The theory of Morel and Anderson moan  was used to obtain $`\varphi `$, which makes it possible to investigate the pressure dependence of $`\mu ^{}`$. Neglecting the pressure dependence of $`\mu ^{}`$, the present theoretical model can be reduced to the two popular models of Seiden seid  and Baryakhtar and Makarov bary  when taking $`\gamma _N=`$ 2/3 and neglecting the direct electron-electron interaction, respectively. Furthermore, we obtained an explicit expression for the change of $`T_c`$ as a function of pressure with the help of Murnaghan equation. The present model enables us to study the pressure behaviors of some interested superconducting parameters such as the zero temperature energy gap $`\mathrm{\Delta }_0`$, the critical field at absolute zero temperature $`H_c(0)`$, the effective interaction strength $`N(E_F)v`$, and the electronic specific heat coefficient $`\gamma `$.

(ii) We investigated the pressure effects on superconducting properties in the newly discovered superconductor MgB2 using our simple approach. It was found that the hydrostatic pressure derivative of $`T_c`$ can be reproduced by using the values of $`\gamma _G`$ obtained from either the Slater relation or $`ab`$ $`initio`$ calculation. The calculated $`d\mathrm{ln}N(E_F)v/d\mathrm{ln}V`$2.0 in MgB2 is close to those obtained in simple $`sp`$ superconductors. The quantitative agreement for the variation of $`T_c`$ with pressure in the low pressure region as well as high pressure region is very good when comparing our theoretical results with experimental data measured by three groups. The predicted values of $`d\mathrm{ln}H_c(0)/dP`$, $`d\mathrm{ln}\mathrm{\Delta }_0/dP`$, and $`\gamma _e`$ are also comparable to those in simple $`sp`$ metals superconductors. All these characteristic pressure behaviors allow us to conclude that MgB2 should be a simple electron-phonon mediated $`sp`$ superconductor and the mechanism in simple $`sp`$ metals superconductors is also responsible for the superconductivity in MgB2.

###### Acknowledgements.

The authors acknowledge useful discussions with O. K. Andersen, O. Jepsen, Y. Kong, R. K. Kremer, and K. Syassen. We are indebted to J. S. Schilling and J. D. Jorgensen for allowing us to use their experimental data prior to publication as well as their valuable comments on the manuscript. XJC thanks the MPG for financial support.

Note added.– After submission of this manuscript, the authors have learned that the superconductivity is not destroyed up to 44 GPa where $`T_c`$ is still as high as 12 K. struz  The intrinsic $`dT_c/dP1.1`$ K/GPa under hydrostatic pressure conditions has recently been reported by other three groups.struz ; schla ; tang  We also have learned a possible explanation given by Tissen $`et`$ $`al.`$tissen  for the large $`dT_c/dP`$ observed in the low $`T_c=37.4\pm 0.1`$ samples.