format: text

text:
# Charge-density wave formation in Sr14Cu24O41

\[

## Abstract

The electrodynamic response of the spin-ladder compound Sr14-xCaxCu24O41 ($`x=0,3,9`$) has been studied from radiofrequencies up to the infrared. At temperatures below 250 K a pronounced absorption peak appears around $`12`$ cm-1 in Sr14Cu24O41 for the radiation polarized along the chains/ladders ($`𝐄𝐜`$). In addition a strongly temperature dependent dielectric relaxation is observed in the kHz - MHz range. We explain this behavior by a charge density wave which develops in the ladders sub-system and produces a mode pinned at $`12`$ cm-1. With increasing Ca doping the mode shifts up in frequency and eventually disappears for $`x=9`$ because the dimensionality of the system crosses over from one to two dimensions, giving way to the superconducting ground state under pressure.

\] Low-dimensional quantum spin systems with spins and charges arranged in chains and/or ladders are under intensive study theoretically as well as experimentally . A number of intriguing magnetic and electronic properties are inherent to those systems because of their specific arrangement of spins and charges; also the interplay between spin and charge degrees of freedom produces a variety of unusual phenomena . Most important, it was predicted  that due to the presence of a spin gap and paired charge carriers, doped even-leg ladder compounds can produce superconductivity reminiscent to that in the underdoped high-Tc cuprates; alternatively, they can undergo a charge-density-wave (CDW) transition.

The spin-ladder system Sr14-xCaxCu24O41, where the superconductivity was discovered under pressure , has attracted much attention in the last years. Its structure contains mutually interpenetrating subsystems of CuO2 chains and Cu2O3 two-leg ladders . The parent compound Sr14Cu24O41 is inherently doped with holes and has a rather high dc conductivity along the $`c`$ axis (ladder/chain direction) of about $`300(\mathrm{\Omega }\mathrm{cm})^1`$ at 300 K, with an anisotropy $`\sigma _c`$:$`\sigma _a`$:$`\sigma _b`$ 1000:100:1  ($`a`$ denotes the rung direction and $`b`$ is the direction perpendicular to the ladders’ plane). When Sr is substituted by the isovalent Ca the conductivity along all three axes increases. It is under debate whether the electronic properties of Sr14-xCaxCu24O41 are determined by a single-particle response or can be influenced by collective effects, like CDWs . Even in the ‘simple’ Sr14Cu24O41 it is not clear what the mechanism of dc conductivity is, how the holes redistribute between chains and ladders when temperature and doping are changed; what is the nature and the interrelation of spin gaps and optical pseudogaps, the role of collective excitations, and the effect of the dimensionality on charge transport and superconductivity. In order to address these questions we have investigated the charge dynamics of Sr14-xCaxCu24O41 by means of optical spectroscopy. In this Letter we present clear evidence that a CDW forms in this quasi-one dimensional compound. Increasing the dimensionality by Ca doping and pressure suppresses the CDW ground state and makes the system superconducting.

The high-quality single crystals (cylinders of about 5 mm in diameter and 7 mm long) were grown by a traveling-solvent floating-zone method. For the dc resistivity measurements four electrical contacts to the sample were obtained with silver paste applied directly on the surface and heated to 600 K in oxygen. In the radiofrequency range (1 Hz to 1 MHz) the spectra of the real and imaginary parts of the dielectric constant ($`ϵ^{}+iϵ^{\prime \prime }`$) were obtained from the complex admittance measured at $`60\mathrm{K}<T<110`$ K by a technique described in . For frequencies 5 to 25 cm-1 a coherent source quasioptical spectrometer was utilized for direct measurements of $`ϵ^{}`$ and $`ϵ^{\prime \prime }`$ . From the far infrared (FIR) up to 10 000 cm-1 we determined the polarized reflection by a Fourier transform spectrometer. The combined data sets were analyzed by the Kramers-Kronig relations in order to obtain the spectra of $`ϵ^{}(\omega )`$, $`ϵ^{\prime \prime }(\omega )`$, and conductivity $`\sigma (\omega )`$. In addition, the non-linear transport of Sr14Cu24O41 was investigated at various temperatures.

The frequency dependences of the $`c`$-axis conductivity and the real and imaginary parts of the dielectric constant of Sr14Cu24O41 are displayed in Fig. 1. Our findings agree with previous infrared experiments . At 300 K the conductivity is Drude-like and coincides with the dc values. Weak phonon features are seen below 1000 cm-1, together with a bump around 2000 cm-1 which is caused by the holes in the ladders . Lowering the temperature, the overall conductivity from dc up to the FIR decreases making the phonons more pronounced; the bump looses its intensity. The new features we have observed are the intensive mode in the FIR range at $`1015\mathrm{cm}^1`$, and a strongly temperature dependent relaxation in the radiofrequency spectra of $`ϵ^{}(\omega )`$ and $`ϵ^{\prime \prime }(\omega )`$. No signs of the FIR mode are found in the $`a`$ and $`b`$ directions . Predictions of a mode existing below 50 cm-1 were made earlier on the basis of infrared measurements in pure  and Ca-doped  Sr14Cu24O41; indications of a peak around 14 cm-1 were also seen in the room-temperature Raman spectra . We observe the development of the mode below approximately 250 K. As the temperature decreases to 60 K (inset of Fig. 1) it gets more pronounced and softens slightly from 14 to 12 cm-1; basically no changes are found going down further to 5 K. While at high temperatures the shape of the mode can be fitted with a single Lorentzian (plus a Drude term for free carriers), it becomes more complicated at low temperatures .

The radiofrequency relaxation is clearly observed as a peak in the $`ϵ^{\prime \prime }`$ spectrum accompanied by a strong dispersion of $`ϵ^{}`$. These spectra can be described phenomenologically in terms of a generalized Debye expressions

$$ϵ^{}(\omega )+iϵ^{\prime \prime }(\omega )=\mathrm{\Delta }ϵ/[1+(i\omega \tau _1)^{1\alpha }].$$
(1)
Here $`\mathrm{\Delta }ϵ`$ is the strength of the relaxation, $`\tau _1`$ is a mean relaxation time and $`1\alpha `$ describes the relaxation time distribution. The solid lines in Fig. 1 result from fits of the spectra by Eq. (1). The so-obtained temperature dependences of the relaxation parameters are shown in Fig. 2 together with the dc conductivity of Sr14Cu24O41.

We ascribe the 12 cm-1 mode to the CDW which develops on the ladders, where the majority of delocalized charges sit , and is pinned by imperfections; correspondingly, we assign the low frequency mode to the CDW relaxation broadened due to normal carriers, which flow around the pinned CDW and produce screening currents. Before we come to a detailed analysis, let us summarize the arguments for our interpretation: (a) It is unlikely that the discovered FIR mode is of simple phonon origin. While practically no change is seen in the phonon spectrum of Sr14-xCaxCu24O41 (above 100 cm-1) when going from $`x=0`$ to $`x=3`$ , this mode shows a strong shift from 12 to 23 cm-1 (Fig.1) and disappears with further doping ($`x=9`$). We also rule out that the resonance is produced by pairs which form a correlated Luttinger liquid – a possibility proposed recently  – because the mode appears at elevated temperatures and cannot be fitted by the corresponding expressions. (b) The CDW is one of the possible ground states in a spin-ladder system like Sr14Cu24O41  where the electronic bands are essentially of one-dimensional character . (c) We observe the CDW-related absorptions only along the highly conducting $`c`$ axis. (d) As expected for the formation of a CDW , experimental indications for the charge order  and lattice superstructure  are seen in the ladders of Sr14-xCaxCu24O41. The phonon spectrum shows signatures of a Brillouin zone folding below 200 K  which gives evidence for a lattice superstructure. (e) The resonance at 12 cm-1 is characterized by an enhanced effective mass, typical for collective excitations like the CDW. (f) We do not observe any resonance for $`x=9`$ when the system looses its one-dimensionality . (g) The absorption features we observe in Sr14Cu24O41 are extremely reminiscent of the CDW response in most well-known one-dimensional compounds like K0.3MoO3, NbSe3, TaS3, and (TaSe4)2I . As we will discuss next, this holds qualitatively (Fig. 1) as well as quantitatively.

Utilizing a standard theory for one-dimensional CDW systems we now will describe the dynamical response of Sr14Cu24O41. Following Littlewood’s model  we consider a deformable CDW which interacts with the uncondensed carriers in conditions of non-uniform pinning. It accounts for the existence of two modes: transverse and longitudinal. The former couples to electromagnetic radiation and yields a high frequency pinned CDW mode, whereas the latter couples to an electrostatic potential and results in an overdamped low frequency relaxation around $`1/(2\pi \tau _1)`$ due to screening effects. In the real systems, absorptions arise from the mixing of a longitudinal plasma-like oscillations of the CDW (with the screened plasma frequency $`\mathrm{\Omega }_p^2=\rho _c^2/ϵ_zm^{}`$) with a transverse pinned mode at $`\mathrm{\Omega }_0=\sqrt{V_0/m^{}}`$ due to symmetry breaking (non-uniform pinning); this mixing can account for the unusual spectral shape of the CDW resonance we observe in Sr14Cu24O41. The pinning potential is given by $`V_0`$; the charge density by $`\rho _c=ne`$; $`ϵ_z`$ is the high-frequency dielectric constant, and $`m^{}`$ is the effective CDW mass. First we note that the position of the relaxation peak in $`ϵ^{\prime \prime }(\omega )`$ is thermally activated (Fig. 2) with the same activation energy as the dc conductivity, as predicted by the model . Further, estimating $`n10^{21}\mathrm{cm}^3`$ as the carrier concentration , taking $`\sigma _z=10^4(\mathrm{\Omega }\mathrm{cm})^1`$ as the plateau conductivity between the two excitations and $`1/(2\pi c\tau _1)=1.7\times 10^5\mathrm{cm}^1`$, $`\mathrm{\Omega }_0/(2\pi c)=12\mathrm{cm}^1`$ (Fig. 1), $`ϵ_z=5`$ , we get an effective mass $`m^{}200m_0`$, a plasma frequency $`\mathrm{\Omega }_p/(2\pi c)190\mathrm{cm}^1`$, and $`n^{}=m^{}\mathrm{\Omega }_p^2/(4\pi e^2)2\times 10^{19}\mathrm{cm}^3`$ carriers condensed, i.e., only 2 % of the charge carriers contribute to the CDW. The CDW plasma frequency can be independently estimated from the oscillator strength of the mode $`\mathrm{\Omega }_p/(2\pi c)=\mathrm{\Omega }_0/(2\pi c\sqrt{ϵ_z})=120\mathrm{cm}^1`$ , which is in good agreement with the value of $`190\mathrm{cm}^1`$ obtained above.

We assign the low-temperature decrease of the spectral weight in the FIR (below 1000 cm-1, see Fig. 1 and ) to the opening of a CDW pseudogap $`E_g/hc2000\mathrm{cm}^1`$ (i.e., $`E_g0.27`$ eV) and determine the CDW mean-field temperature $`T_{\mathrm{CDW}}^{MF}=E_g/3.5k_B900`$ K; this implies that CDW fluctuations should be observed well above the room temperature. Similar rather high mean-field temperatures were estimated in other quasi one-dimensional conductors like TTF-TCNQ ($`T_{\mathrm{CDW}}^{MF}450`$ K) or K0.3MoO3 (330 K) . The pseudogap should govern the activated behavior of the dc conductivity (Fig. 2, $`170\mathrm{K}<T<250`$ K). The fact that $`\sigma _{dc}(T)`$ is not simply activated for $`250\mathrm{K}<T<300`$ K can be explained with CDW fluctuations which produce a collective contribution to the dc transport. This sort of paraconductivity was also detected in other one-dimensional CDW conductors, like K0.3MoO3 or TTF-TCNQ .

Below approximately 170 K the activation energy of $`\sigma _{dc}(T)`$ changes from $`E_g=0.27`$ eV to 0.12 eV. A number of magnetic, electronic and structural properties display an anomaly around 170 K  which may indicate some type of phase transition. We suggest that this transition is due to a second CDW which develops in the chain sub-system of Sr14Cu24O41. Mobile holes in the chains, essential for the CDW formation, have been suggested on the basis of ESR measurements . The additional peak observed in the microwave range at $`1.7\mathrm{cm}^1`$ for $`T<170`$ K , shown in Fig. 1, thus corresponds to the second pinned CDW mode; its spectral weight is much smaller compared to the FIR mode, due to the lower charge concentration in the chains. Similar to the FIR mode reported here, this microwave resonance is seen in pure Sr14Cu24O41 and for small Ca doping of $`x=3`$ while is disappears for larger $`x`$. We note that signs of this resonance lying below 8 cm-1 are also clearly seen in our low temperature reflectivity spectra  for $`x=3`$.

Although the presented results clearly point towards the collective CDW ground state in the spin-ladder system Sr14Cu24O41 and can even be well described by Littlewood’s model, a few open questions remain. A signature of an absorption line at 14 cm-1 was already observed in the Raman spectra at room temperature . While the CDW phase fluctuations should be infrared active and the CDW amplitude fluctuations only Raman active , it is likely that both excitations mix due to symmetry breaking; hence the resultant excitation is observable by both techniques (but still only along the conducting $`c`$ direction). Interestingly, the Raman mode seems to be pinned even at 300 K, while we do not see any sign of the optical resonance at this temperature probably due to the large screening by free carriers. The CDW can be depinned by an external electric field leading to non-linear transport . We do observe a field-dependent conductivity in our system, however, the non-linearity is very small and the threshold current cannot be clearly determined. An estimate of the threshold field $`E_T0.1`$ V/cm (Fig. 2) cannot be related to the pinning frequency $`\mathrm{\Omega }_0`$ by the standard expression  $`m^{}\mathrm{\Omega }_0^2=2eE_T/\lambda `$, with the wavelength $`\lambda `$ of the CDW taken of the order of the distance between Cu atoms on the ladder ($``$4 Å); a similar discrepancy was noticed for the microwave mode . The effective mass $`m^{}200m_0`$ is small compared to standard CDW materials where up to the order of ten thousand has been reported . The mentioned points can be ascribed to an unusual CDW state formed in Sr14Cu24O41, as compared to typical one-dimensional conductors. This is most likely due to the two interacting sublattices on the chains and on the ladders. Since each is characterized by magnetic, structural, and charge ordering, the ground state of the system is determined by a mixture of correspondent order parameters leading to excitations which may be rather unconventional. Detailed low-energy studies of the interplay between the spin and charge subsystems in (Sr,Ca,La)14Cu24O41 for different doping and in magnetic fields are required. In addition, the internal deformations of the CDW and the pinning mechanism may be more complicated in this system of interacting chains and ladders compared to strictly one-dimensional compounds. It is not surprising that the observed CDW phenomenon is strongly modified compared to the conventional CDW compounds.

In summary, we presented evidence of a charge-density wave formation in the ladder subsystem in Sr14-xCaxCu24O41. For $`x=0`$ we observed a pinned CDW mode at 12 cm-1 below 250 K and a strongly temperature dependent relaxation in the radiofrequency range due to single-particle screening currents. The pinned mode shifts up to 23 cm-1 for $`x=3`$ and disappears when the Ca doping reaches $`x=9`$. From the two competing ground states (density wave and superconductivity) expected for spin-ladder compounds, the CDW instability dominates when Sr14-xCaxCu24O41 is highly one-dimensional ($`0x3`$) but shows a tendency to disappear when stronger Ca doping ($`x=9`$) makes it more two-dimensional. When external pressure increases the coupling (and thus the dimensionality) even further, the system finally becomes superconducting.

We thank P. Lemmens and A. Maeda for discussions and exchange of results prior to publication. The work was supported by the Deutsche Forschungsgemeinschaft (DFG) and Ministry of Science of RH.

Fig.1. (a) Frequency dependent conductivity of Sr14Cu24O41 measured at different temperatures along the highly conducting $`c`$ direction. The arrows indicate the dc conductivity; the solid lines represent fits by the generalized Debye relaxation, the dashed lines are guides to the eye. The peak around 1 cm-1 corresponds to the data from . The inset shows how the shape of the pinned CDW resonance varies with temperature. (b) Real and imaginary parts of the dielectric constant as a function of frequency; the solid lines are fits by the generalized Debye expression. The inset exhibits the evolution of the spectral weight in the FIR range with decreasing temperature: the dotted line shows the pinned CDW mode of Sr11Ca3Cu24O41.

Fig.2. Temperature dependences of the dc conductivity of Sr14Cu24O41 along the $`c`$ axis and of the parameters which characterize the dielectric relaxation at radiofrequencies \[Eq. (1)\]: the inverse time $`1/(2\pi \tau _1)`$ and the strength $`\mathrm{\Delta }ϵ`$ of the dielectric relaxation, and the parameter 1-$`\alpha `$ which describes the symmetrical broadening of the relaxation time distribution. The dc conductivity is thermally activated $`\sigma _{dc}\mathrm{exp}\{\mathrm{\Delta }/k_BT\}`$ between 250 K and 170 K with an energy gap $`\mathrm{\Delta }0.27`$ eV ($`\mathrm{\Delta }/hc2200`$ cm-1) in agreement with the onset of the low-temperature conductivity spectrum in Fig. 1 (inset of the lower panel) and is identified as a pseudogap due to CDW in the ladders. Below 170 K a change of the activation energy to $`\mathrm{\Delta }0.12`$ eV is observed due to the onset of a CDW in the chains. The solid lines correspond to an activated behaviors of $`\sigma _{dc}(T)`$ and $`1/\tau _1(T)`$. The inset shows an electric field dependent resistance of Sr14Cu24O41 measured along the $`c`$ direction at 87 K.


text:
# On the existence of self-similar spherically symmetric wave maps coupled to gravity

## 1 Introduction

In this paper we continue our investigations, started in  (referred to as I), of wave maps coupled to gravity, that is solutions of Einstein’s equations with an $`SU(2)`$ sigma field as matter. We found numerically in I that for $`\alpha <1/2`$ ($`\alpha `$ is the dimensionless coupling constant) the model admits a countable family of continuously self-similar (CSS) solutions, labeled by an integer nodal index $`n=0,1,\mathrm{}`$, that are analytic inside the past light cone of the singularity. We also provided evidence that the $`n`$th CSS solution can be continued up to the future light cone of the singularity if $`\alpha <\alpha _n`$, where $`\{\alpha _n\}`$ is an increasing sequence bounded above by $`1/2`$. The purpose of this paper is to make the results of I into theorem–proof rigorous mathematics. This is accomplished by applying a shooting argument to the resulting dynamical system. We note that the case $`\alpha =0`$ was previously analyzed in .

The physical importance of the CSS solutions considered here was discussed in I, in particular we conjectured that in a certain parameter range ($`\alpha _0<\alpha <\alpha _1`$) the $`n=1`$ solution is a critical solution at the threshold of black hole formation. This conjecture has been recently confirmed in numerical studies of the critical behaviour  and in the linear stability analysis . As far as we know, this is the only case where the existence of a self-similar solution, which was numerically found as the critical solution in gravitational collapse, has been established rigorously.

## 2 Setup

For the reader’s convenience we repeat from I the basic setting for the problem. Let $`X:MN`$ be a map from a spacetime $`(M,g_{ab})`$ into a Riemannian manifold $`(N,G_{AB})`$. Wave maps coupled to gravity are defined as extrema of the action

$$S=_M\left(\frac{R}{16\pi G}+L_{WM}\right)𝑑v_g$$
(1)
with the Lagrangian density

$$L_{WM}=\frac{f_\pi ^2}{2}g^{ab}_aX^A_bX^BG_{AB}.$$
(2)
Here $`G`$ is Newton’s constant and $`f_\pi ^2`$ is the wave map coupling constant. The product $`\alpha =4\pi Gf_\pi ^2`$ is dimensionless. The field equations derived from (1) are the wave map equation

$$\mathrm{}_gX^A+\mathrm{\Gamma }_{BC}^A(X)_aX^B_bX^Cg^{ab}=0,$$
(3)
where $`\mathrm{\Gamma }_{BC}^A(X)`$ are the Christoffel symbols of the target metric $`G_{AB}`$ and $`\mathrm{}_g`$ is the d’Alembertian associated with the metric $`g_{ab}`$, and the Einstein equations $`R_{ab}\frac{1}{2}g_{ab}R=8\pi GT_{ab}`$ with the stress-energy tensor

$$T_{ab}=f_\pi ^2\left(_aX^A_bX^B\frac{1}{2}g_{ab}(g^{cd}_cX^A_dX^B)\right)G_{AB}.$$
(4)
As a target manifold we take the three-sphere $`S^3`$ with the standard metric in polar coordinates $`X^A=(F,\mathrm{\Theta },\mathrm{\Phi })`$

$$G_{AB}dX^AdX^B=dF^2+\mathrm{sin}^2F(d\mathrm{\Theta }^2+\mathrm{sin}^2\mathrm{\Theta }d\mathrm{\Phi }^2).$$
(5)
For the domain manifold we assume spherical symmetry and use Schwarzschild coordinates

$$g_{ab}dx^adx^b=e^{2\delta }Adt^2+A^1dr^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2),$$
(6)
where $`\delta `$ and $`A`$ are functions of $`(t,r)`$. Next, we assume that the wave maps are corotational, that is

$$F=F(t,r),\mathrm{\Theta }=\theta ,\mathrm{\Phi }=\varphi .$$
(7)
Equation (3) reduces then to the single semilinear wave equation

$$\mathrm{}_gF\frac{\mathrm{sin}(2F)}{r^2}=0,$$
(8)
where

$$\mathrm{}_g=e^\delta _t(e^\delta A^1_t)+\frac{e^\delta }{r^2}_r(r^2e^\delta A_r),$$
(9)
and the Einstein equations become

$`_tA`$ $`=`$ $`2\alpha rA(_tF)(_rF),`$ (10)
$`_r\delta `$ $`=`$ $`\alpha r\left((_rF)^2+A^2e^{2\delta }(_tF)^2\right),`$ (11)
$`_rA`$ $`=`$ $`{\displaystyle \frac{1A}{r}}\alpha r\left(A(_rF)^2+A^1e^{2\delta }(_tF)^2+2{\displaystyle \frac{\mathrm{sin}^2F}{r^2}}\right).`$ (12)
These equations are invariant under dilations $`(t,r)(\lambda t,\lambda r)`$ so it is natural to look for continuously self-similar (CSS) solutions, that is solutions which are left invariant by the action of the homothetic Killing vector $`K=t_t+r_r`$. To study such solutions it is convenient to use similarity variables $`\rho =r/(t)`$ and $`\tau =\mathrm{ln}(t)`$. Then $`K=_\tau `$, so CSS solutions do not depend on $`\tau `$. Assuming this and using an auxiliary function $`Z=e^\delta \rho /A`$, we reduce equations (8-12) to the system of ordinary differential equations (where prime is $`d/d\rho `$)

$`F^{\prime \prime }`$ $`+`$ $`{\displaystyle \frac{2}{\rho }}F^{}\alpha (1+Z^2)\rho F_{}^{}{}_{}{}^{3}{\displaystyle \frac{\mathrm{sin}(2F)}{A\rho ^2(1Z^2)}}=0,`$ (13)
$`A^{}`$ $`=`$ $`2\alpha \rho AF_{}^{}{}_{}{}^{2},`$ (14)
$`\rho Z^{}`$ $`=`$ $`Z(1+\alpha (1Z^2)\rho ^2F_{}^{}{}_{}{}^{2}),`$ (15)
$`\rho A^{}`$ $`=`$ $`1A\alpha \left(\rho ^2A(1+Z^2)F_{}^{}{}_{}{}^{2}+2\mathrm{sin}^2F\right).`$ (16)
The combination of (14) and (16) yields the constraint

$$1A2\alpha \mathrm{sin}^2F+\alpha A\rho ^2F_{}^{}{}_{}{}^{2}(1Z^2)=0.$$
(17)
This system of equations has a fixed singularity at the center $`\rho =0`$ and moving singularities at points where $`Z(\rho )=\pm 1`$ and/or $`A(\rho )=0`$. In terms of the similarity coordinate $`\rho `$, the metric (6) takes the form

$$ds^2=A^1(1Z^2)\rho ^2dt^2+2A^1t\rho dtd\rho +A^1t^2d\rho ^2+t^2\rho ^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2),$$
(18)
hence the hypersurfaces $`Z=\pm 1`$ are null (provided that $`A>0`$). The first $`\rho _1`$ where $`Z(\rho _1)=1`$ is the locus of the past light cone of the singularity at the origin $`(t=0,r=0)`$ (in what follows we shall refer to the past and future light cones of the singularity as to the past and future self-similarity horizons (SSH)). By rescaling, $`\rho \rho /\rho _1`$, one can always locate the past self-similarity horizon at $`\rho _1=1`$, that is $`Z(1)=1`$. To ensure regularity of solutions in the interval $`0\rho 1`$, the equations (13-17) must be supplemented by the boundary conditions at both endpoints

$`F(0)`$ $`=`$ $`0,F^{}(0)=c,Z(0)=0,A(0)=1,`$ (19)
$`F(1)`$ $`=`$ $`{\displaystyle \frac{\pi }{2}},F^{}(1)=b,Z(1)=1,A(1)=12\alpha ,`$ (20)
where $`c`$ and $`b`$ are free parameters.

Our main result is the following theorem:

###### Theorem 1.

For any $`0\alpha <1/2`$ and any nonegative integer $`n`$, the equations (13-17) have an analytic solution $`(F_n,A_n,Z_n)`$ which satisfies the boundary conditions (19-20) and has precisely $`n`$ oscillations of $`F_n(\rho )`$ around $`\pi /2`$.

In the next section we shall prove this theorem using a shooting technique. The numerical evidence for Theorem 1 was given in I . The case $`\alpha =0`$ was proved previously in  so hereafter we assume that $`0<\alpha <1/2`$.

## 3 Proof of Theorem 1

For convenience we rewrite equations (13-15) in terms of $`H=F\pi /2`$:

$`H^{\prime \prime }`$ $`+`$ $`{\displaystyle \frac{2}{\rho }}H^{}\alpha (1+Z^2)\rho H_{}^{}{}_{}{}^{3}+{\displaystyle \frac{\mathrm{sin}(2H)}{A\rho ^2(1Z^2)}}=0,`$ (21)
$`A^{}`$ $`=`$ $`2\alpha \rho AH_{}^{}{}_{}{}^{2},`$ (22)
$`\rho Z^{}`$ $`=`$ $`Z(1+\alpha (1Z^2)\rho ^2H_{}^{}{}_{}{}^{2}),`$ (23)
The constraint becomes

$$12\alpha A+2\alpha \mathrm{sin}^2H+\alpha A\rho ^2H_{}^{}{}_{}{}^{2}(1Z^2)=0.$$
(24)
The initial conditions at $`\rho =0`$ are

$$H(0)=\frac{\pi }{2},H^{}(0)=c,A(0)=1,Z(0)=0,Z^{}(0)=1.$$
(25)
Note that the above equations have a residual scaling symmetry $`\rho \lambda \rho `$. The initial condition $`Z^{}(0)=1`$ is imposed temporarily in order to fix the scale. We shall refer to solutions of equations (21-24) satisfying the initial conditions (25) as to $`c`$-orbits. In the appendix we show that $`c`$-orbits exist locally and are analytic in $`\rho `$ and $`c`$. Now we shall show that $`c`$-orbits can be extended up to a point $`\rho _1`$ at which $`Z(\rho _1)=1`$.

###### Proposition 2.

For any $`0<\alpha <1/2`$ and $`c>0`$ there is a $`\rho _1(c)(\sqrt{12\alpha },1)`$, such that the $`c`$-orbit is defined for all $`\rho <\rho _1`$ and $`lim_{\rho \rho _1}Z(\rho )=1`$. Furthermore, the following limits exist

$$\frac{\pi }{2}<\overline{H}\stackrel{def}{=}\underset{\rho \rho _1}{lim}H(\rho )<\frac{\pi }{2},\overline{A}\stackrel{def}{=}\underset{\rho \rho _1}{lim}A(\rho )=12\alpha \mathrm{cos}^2\overline{H},\underset{\rho \rho _1}{lim}(1Z^2)H_{}^{}{}_{}{}^{2}=0.$$

*Proof:* Let the maximum domain of definition of the $`c`$-orbit be $`0\rho <\rho _1`$ and assume that $`Z(\rho )<1`$ in this interval. Then, from constraint (24) we have $`A12\alpha >0`$ and hence $`\overline{A}=lim_{\rho \rho _1}A(\rho )>0`$ ($`\overline{A}`$ exists since $`A(\rho )`$ is monotone decreasing). By (23) $`Z^{}0`$, hence $`\overline{Z}=lim_{\rho \rho _1}Z(\rho )`$ exists. If $`\overline{Z}<1`$, then from constraint (24) $`H_{}^{}{}_{}{}^{2}`$ is bounded so $`\overline{H}=lim_{\rho \rho _1}H(\rho )`$ exists, which in turn implies, again by (24), that $`lim_{\rho \rho _1}H^{}`$ exits. Thus, $`H,H^{},A`$, and $`Z`$ all have finite limits at $`\rho _1`$ and therefore the $`c`$-orbit may be continued beyond $`\rho _1`$ contradicting the maximality of $`\rho _1`$. We conclude that $`\overline{Z}=1`$.

Now, we must show that $`\overline{H}(\pi /2,\pi /2)`$ exists. Since $`\overline{Z}=1`$, we may no longer conclude that $`H_{}^{}{}_{}{}^{2}`$ is bounded but from equation (22) we get $`(\mathrm{ln}A)^{}=2\alpha \rho H_{}^{}{}_{}{}^{2}`$, so $`H_{}^{}{}_{}{}^{2}`$ is integrable near $`\rho _1`$ which implies that $`H^{}`$ is absolutely integrable ($`|H^{}|<1+H_{}^{}{}_{}{}^{2}`$) and thus $`\overline{H}`$ exists. From constraint (24), $`H(\rho )=\pm \pi /2`$ for some $`0<\rho <\rho _1`$ is not possible since $`1A>0`$. Thus, $`\pi /2<H(\rho )<\pi /2`$ and so $`\pi /2\overline{H}\pi /2`$. In fact, for $`\rho \rho _1/2`$ we have $`1A\sigma >0`$, so $`2\alpha \mathrm{cos}^2H\sigma >0`$ (remember that we assume $`\alpha >0`$), hence $`H`$ is uniformly bounded away from $`\pm \pi /2`$, and thus $`\pi /2<\overline{H}<\pi /2`$.

To prove $`\overline{A}=12\alpha \mathrm{cos}^2\overline{H}`$, note that by (24) $`d=lim_{\rho \rho _1}H_{}^{}{}_{}{}^{2}(1Z^2)`$ exits and is finite. Hence, by (23) $`lim_{\rho \rho _1}Z^{}`$ exists and is finite so $`1Z^2=O(\rho \rho _1)`$ near $`\rho _1`$. If $`d0`$, then $`H^{}(\rho )d/(\rho _1\rho )`$ would not be integrable near $`\rho _1`$, thus $`d`$ must be zero. Inserting this into (24) we get $`\overline{A}=12\alpha \mathrm{cos}^2\overline{H}`$.

Next, $`(Z/\rho )^{}>0`$ by (23) and $`lim_{\rho 0}(Z/\rho )=1`$ by L’Hôpital’s rule, hence $`Z\rho `$ for all $`\rho >0`$, and thus $`\rho _11`$. Finally, from (22) and (23)

$$\left(\frac{AZ^2}{\rho ^2}\right)^{}=\frac{2Z^4A\alpha H_{}^{}{}_{}{}^{2}}{\rho }<0,$$
(26)
and since $`lim_{\rho 0}(AZ^2/\rho ^2)=1`$, we have $`(AZ^2/\rho ^2)1`$ and hence $`\rho _1>\sqrt{A}>\sqrt{12\alpha }`$.

If $`Z(\rho _2)=1`$ for some $`\rho _2<\rho _1`$, we replace $`\rho _1`$ by $`\rho _2`$ in the above arguments.

###### Corollary 3.

The function $`\rho _1(c)`$ is continuous.

*Proof:* Let $`\stackrel{~}{c}`$ be given and let $`ϵ>0`$. By Proposition 2, $`\rho _1(\stackrel{~}{c})`$ is defined. The function $`Z(\rho )`$ is monotone increasing for $`\rho <\rho _1(\stackrel{~}{c})`$, so $`Z(\rho _1(\stackrel{~}{c})ϵ,\stackrel{~}{c})<1`$, hence for all $`c`$ sufficiently close to $`\stackrel{~}{c}`$, $`Z(\rho _1(\stackrel{~}{c})ϵ,c)<1`$, and thus $`\rho _1(c)>\rho _1(\stackrel{~}{c})ϵ`$. To show that $`\rho _1(c)<\rho _1(\stackrel{~}{c})+ϵ`$ for all $`c`$ sufficiently close to $`\stackrel{~}{c}`$, we assume otherwise and get a contradiction. By the mean-value theorem $`Z(\rho _1(\stackrel{~}{c})+ϵ,c)Z(\rho ,c)=Z^{}(\xi ,c)(\rho _1(\stackrel{~}{c})+ϵ\rho )`$. By continuity $`Z(\rho ,c)`$ is close to $`Z(\rho ,\stackrel{~}{c})`$ and $`Z(\rho ,\stackrel{~}{c})`$ is close to $`1`$ if $`\rho `$ is close to $`\rho _1(\stackrel{~}{c})`$, hence $`Z(\rho ,c)`$ is arbitrarily close to $`1`$. But $`Z^{}(\rho ,c)>Z(\rho ,c)/\rho >1`$, so $`Z(\rho _1(\stackrel{~}{c})+ϵ,c)>Z(\rho ,c)+ϵ>1`$, which is a contradiction. Thus, $`\rho _1(c)<\rho _1(\stackrel{~}{c})+ϵ`$.

###### Lemma 4.

$`H^{}(\rho )`$ is bounded near $`\rho _1`$ if and only if $`\overline{H}=0`$.

*Proof:* Suppose that $`\overline{H}0`$ and $`H^{}(\rho )`$ is bounded. Then, in (21) we have

$$H^{\prime \prime }=\text{bounded terms}\frac{\mathrm{sin}2H}{A\rho ^2(1Z^2)}\frac{d}{\rho _1\rho },$$
(27)
where $`d0`$. This contradicts that $`H^{}(\rho )`$ is bounded near $`\rho _1`$ and concludes the ”only if” part of Lemma 4.

Suppose now that $`H(\rho _1)=0`$ and $`H^{}(\rho )`$ is unbounded. Without loss of generality we consider the case that $`H(\rho )<0`$ and $`H^{}(\rho )>0`$ near $`\rho _1`$. Dividing equation (21) by $`H^{}`$ and integrating from $`\rho `$ to $`\rho _1`$ we obtain

$$_\rho ^{\rho _1}\left(\frac{H^{\prime \prime }}{H^{}}+\frac{2}{\rho }\alpha (1+Z^2)\rho H_{}^{}{}_{}{}^{2}+\frac{\mathrm{sin}(2H)}{H^{}A\rho ^2(1Z^2)}\right)𝑑\rho =0.$$
(28)
The first integral is divergent because $`lim_{\rho \rho _1}\mathrm{ln}H^{}=\mathrm{}`$. The second and the third terms are integrable (remember that $`H_{}^{}{}_{}{}^{2}`$ is integrable). Thus, to get a contradiction it suffices to show that the last term is integrable. We write this term as

$$\frac{\mathrm{sin}(2H)}{H^{}A\rho ^2(1Z^2)}=\frac{\mathrm{sin}(2H)}{HA\rho ^2}\frac{H}{(1Z^2)H^{}}.$$
(29)
The first factor is continuous and we now show that the second factor is also continuous. Applying L’Hôpital’s rule we get

$$\underset{\rho \rho _1}{lim}\frac{H}{(1Z^2)H^{}}=\underset{\rho \rho _1}{lim}\frac{H^{}}{2ZZ^{}H^{}+(1Z^2)H^{\prime \prime }}=\underset{\rho \rho _1}{lim}\frac{1}{2ZZ^{}+(1Z^2)H^{\prime \prime }/H^{}}.$$
(30)
Next, using (21) we get

$$(1Z^2)\frac{H^{\prime \prime }}{H^{}}=\frac{2(1Z^2)}{\rho }+\alpha \rho (1+Z^2)(1Z^2)H_{}^{}{}_{}{}^{2}\frac{\mathrm{sin}(2H)}{A\rho ^2H^{}}.$$
(31)
In the limit $`\rho \rho _1`$, the first term on the rhs of (31) obviously goes to zero, the second does by Proposition 2, and the third does by the assumption that $`H^{}\mathrm{}`$. Thus, the limit (30) is finite and consequently so is (29). This contradicts (28) and thus concludes the proof of the ”if” part of Lemma 4.

###### Corollary 5.

A $`c`$-orbit which has $`\overline{H}(c)`$=0 is analytic on the whole interval $`0\rho \rho _1`$.

*Proof:* The boundedness of $`H^{}(\rho )`$ implies by (21) that $`H^{\prime \prime }>2H^{}/\rho `$ is bounded below (remember that $`H(\rho )<0`$ and $`H^{}(\rho )>0`$ near $`\rho _1`$), hence $`lim_{\rho \rho _1}H^{}(\rho )`$ exists. Having that, it is easy to show by applying L’Hôpital’s rule to $`lim_{\rho \rho _1}(H^{(k)}(\rho _1)H^{(k)}(\rho ))/(\rho _1\rho )`$ for $`k=0,1`$ that the solution $`(H,A,Z)`$ is $`C^2`$ near $`\rho _1`$. By a routine contraction mapping argument one can show that $`C^2`$ solutions are unique, hence a $`c`$-orbit must belong to the one-parameter family of analytic solutions from Proposition 14 (see the appendix). Next, we describe the behaviour of $`c`$-orbits for small and large values of the shooting parameter $`c`$. We define a nodal number of a $`c`$-orbit $`N(c)=`$ number of zeros of the function $`H(\rho )`$ on the interval $`0\rho <\rho _1`$. We first show that $`c`$-orbits with small $`c`$ have no nodes.

###### Proposition 6.

If $`c`$ is sufficiently small then $`N(c)=0`$.

*Proof:* For $`c=0`$ we have $`H(\rho )\pi /2`$ and $`Z(\rho )=\rho `$ so $`\rho _1(c=0)=1`$. By continuity, for any $`ϵ>0`$ and sufficiently small $`c`$ we can find $`\rho _0`$ such that $`1ϵ<\rho _0<\rho _1(c)<1`$ and $`H(\rho _0)<\pi /2+ϵ`$. We know from the proof of Proposition 2 that $`lim_{\rho \rho _1}\sqrt{\rho _1\rho }H^{}=0`$, hence

$$H(\rho _1)H(\rho _0)=_{\rho _0}^{\rho _1}H^{}(\rho )𝑑\rho <const_{\rho _0}^{\rho _1}\frac{d\rho }{\sqrt{\rho _1\rho }}<const\sqrt{ϵ}.$$
(32)
Thus, $`H(\rho )`$ stays arbitrarily close to $`\pi /2`$ all the way up to $`\rho _1`$ if $`c`$ is sufficiently small and therefore $`N(c)=0`$. We remark that using a scaling argument one can derive the precise asymptotic behaviour of $`c`$-orbits for small $`c`$. We omit this argument since it is not needed for the proof. We show next that $`c`$-orbits with large $`c`$ have arbitrarily many nodes.

###### Proposition 7.

$`N(c)\mathrm{}`$ for $`c\mathrm{}`$.

*Proof:* We rescale the variables, setting

$$x=c\rho ,\stackrel{~}{H}(x)=H(\rho ),\stackrel{~}{A}(x)=A(\rho ),\stackrel{~}{Z}(x)=cZ(\rho ).$$
(33)
Then, equations (21-24) become

$`\stackrel{~}{H}^{\prime \prime }`$ $`+`$ $`{\displaystyle \frac{2}{x}}\stackrel{~}{H}^{}\alpha (1+{\displaystyle \frac{\stackrel{~}{Z}^2}{c^2}})xH_{}^{}{}_{}{}^{3}+{\displaystyle \frac{\mathrm{sin}(2\stackrel{~}{H})}{\stackrel{~}{A}x^2(1\frac{\stackrel{~}{Z}^2}{c^2})}}=0,`$ (34)
$`\stackrel{~}{A}^{}`$ $`=`$ $`2\alpha x\stackrel{~}{A}\stackrel{~}{H^{}}^2,`$ (35)
$`x\stackrel{~}{Z}^{}`$ $`=`$ $`\stackrel{~}{Z}(1+\alpha (1{\displaystyle \frac{\stackrel{~}{Z}^2}{c^2}})x^2\stackrel{~}{H^{}}^2),`$ (36)
with the constraint

$$12\alpha \stackrel{~}{A}+2\alpha \mathrm{sin}^2\stackrel{~}{H}+\alpha \stackrel{~}{A}x^2\stackrel{~}{H^{}}^2(1\frac{\stackrel{~}{Z^2}}{c^2})=0,$$
(37)
and the initial conditions at $`x=0`$

$$\stackrel{~}{H}(0)=\frac{\pi }{2},\stackrel{~}{H}^{}(0)=1,\stackrel{~}{A}(0)=1,\stackrel{~}{Z}(0)=0,\stackrel{~}{Z}^{}(0)=1.$$
(38)
As $`c\mathrm{}`$, the solutions of equations (34)-(38) tend uniformly on compact intervals to solutions of the limiting equations

$`h^{\prime \prime }`$ $`+`$ $`{\displaystyle \frac{2}{x}}h^{}\alpha xh_{}^{}{}_{}{}^{3}+{\displaystyle \frac{\mathrm{sin}(2h)}{ax^2}}=0,`$ (39)
$`a^{}`$ $`=`$ $`2\alpha xah_{}^{}{}_{}{}^{2},`$ (40)
$`xz^{}`$ $`=`$ $`z(1+\alpha x^2h_{}^{}{}_{}{}^{2}),`$ (41)
with the constraint

$$12\alpha a+2\alpha \mathrm{sin}^2h+\alpha ax^2h_{}^{}{}_{}{}^{2}=0,$$
(42)
and the same initial conditions at $`x=0`$

$$h(0)=\frac{\pi }{2},h^{}(0)=1,a(0)=1,z(0)=0,z^{}(0)=1.$$
(43)
We observe first that the function $`a(x)`$ is monotone decreasing by (40) and bounded below, $`a>12\alpha `$, by (42). Thus, no singularity can develop due to $`a`$ going to zero. Also, by (42) no singularity can develop due to $`h^{}`$ becoming unbounded. Thus, solutions exist for all $`x>0`$ (assuming the existence of a solution for small $`x`$). In order to complete the proof it is sufficient to show that the function $`h(x)`$ has an infinite number of zeros for $`x>0`$. Since $`a<1`$, it follows from (42) that $`\pi /2<h(x)<\pi /2`$ for all $`x>0`$. To show that $`h(x)`$ oscillates around zero we consider three cases:
(i) Assume that $`lim_x\mathrm{}h(x)`$ does not exist. Then, there must be a sequence $`\mathrm{}x_k<y_k<x_{k+1}<y_{k+1}<\mathrm{}`$ such that $`h`$ has a local minimum at $`x_k`$ and a local maximum at $`y_k`$. By (39), $`h^{}(x_k)=0,h^{\prime \prime }(x_k)0`$ imply that $`\mathrm{sin}(2h(x_k))0`$, hence $`h(x_k)0`$. By a similar argument, $`h(y_k)0`$. Thus, $`h(x)`$ has a zero in each interval $`x_k<x<y_k`$.
(ii) Assume that a nonzero $`lim_x\mathrm{}h(x)`$ exists. Then, from (42) $`lim_x\mathrm{}x^2h_{}^{}{}_{}{}^{2}`$ exists and, in fact, equals zero because $`lim_x\mathrm{}h(x)`$ exists. This implies by (39) that $`lim_x\mathrm{}x^2h^{\prime \prime }(x)=\mathrm{sin}(2h(\mathrm{}))/A(\mathrm{})0`$, hence $`lim_x\mathrm{}x^2h_{}^{}{}_{}{}^{2}(x)0`$. Thus the case (ii) does no arise.
(iii) Assume that $`lim_x\mathrm{}h(x)=0`$. We define the rotation function $`\theta (x)`$ by

$$\mathrm{tan}\theta (x)=\frac{xh^{}(x)}{h(x)},\theta (0)=0.$$
(44)
*Remark 1.* The rotation function $`\theta (x)`$ is well defined because the case $`h(x)=h^{}(x)=0`$ is impossible for solutions satisfying the initial conditions (43). To see this, assume that $`h(x_0)=h^{}(x_0)=0`$ for some $`x_0>0`$. Then, by (42) $`a(x_0)=12\alpha `$ and the unique solution with these initial conditions at $`x_0`$ is $`h(x)=0,h^{}(x)=0,a(x)=12\alpha `$ for all $`x`$, contradicting the initial conditions (43). We want to show that $`lim_x\mathrm{}\theta (x)=\mathrm{}`$. Using (39) we obtain

$$x\theta ^{}(x)=\mathrm{sin}^2\theta \frac{\mathrm{sin}2h}{2h}\frac{2\mathrm{cos}^2\theta }{a}\frac{(12\alpha \mathrm{cos}^2h)\mathrm{sin}\theta \mathrm{cos}\theta }{a}.$$
(45)
Under the assumption $`lim_x\mathrm{}h(x)=0`$, it follows from (42) that $`lim_x\mathrm{}a(x)=12\alpha `$, hence for sufficiently large $`x`$

$$\theta ^{}(x)\frac{1}{x}\left(\mathrm{sin}^2\theta +\mathrm{sin}\theta \mathrm{cos}\theta +\frac{2\mathrm{cos}^2\theta }{12\alpha }\right)<\frac{3}{4x},$$
(46)
so $`lim_x\mathrm{}\theta (x)=\mathrm{}`$. Thus, given any integer $`k`$ there exists an $`x_k`$ such that $`h(x)`$ has at least $`k`$ zeroes for $`x<x_k`$. By continuous dependence on initial conditions, we may choose $`c>x_k/\sqrt{12\alpha }`$ so that the $`c`$-solution has $`k`$ zeroes also for $`x<x_k`$. In terms of the variable $`\rho =x/c`$ the $`c`$-solution has $`k`$ zeroes for $`\rho <\sqrt{12\alpha }<\rho _1(c)`$. This completes the proof of Proposition 7. Next, we need two lemmas which tell us how the number of nodes $`N(c)`$ changes under small variations of $`c`$.

###### Lemma 8.

If $`\overline{H}(\stackrel{~}{c})=0`$, then $`N(c)=N(\stackrel{~}{c})`$ or $`N(c)=N(\stackrel{~}{c})+1`$ for $`c`$ sufficiently close to $`\stackrel{~}{c}`$.

*Proof:* First note that if $`H(\rho ,\stackrel{~}{c})`$ has a zero at some $`\rho _0<\rho _1(\stackrel{~}{c})`$, then $`H^{}(\rho _0,\stackrel{~}{c})0`$ (see Remark 1) so by continuity of $`H(\rho ,c)`$ with respect to $`c`$, $`H(\rho ,c)`$ also has a zero if $`c`$ is sufficiently close to $`\stackrel{~}{c}`$. Thus $`N(c)N(\stackrel{~}{c})`$ and it suffices to show that $`N(c)N(\stackrel{~}{c})+1`$. Let $`\stackrel{~}{a}<\rho _1(\stackrel{~}{c})`$ be the last node of the $`\stackrel{~}{c}`$-orbit, that is $`H(\stackrel{~}{a},\stackrel{~}{c})=0`$ and, for concreteness, $`H(\rho ,\stackrel{~}{c})<0`$ for $`\stackrel{~}{a}<\rho <\rho _1`$. By continuity with respect to $`c`$, $`H(\rho ,c)`$ will also have a zero at $`a`$ near $`\stackrel{~}{a}`$ if $`c`$ is near $`\stackrel{~}{c}`$. In order to prove that $`H(\rho ,c)`$ cannot have more than one zero in the interval $`a<\rho <\rho _1(c)`$, we now show that if $`H(\rho ,c)`$ becomes positive for some $`\rho >a`$, then it would not have time to change the sign again before going singular. Assume for contradiction that there is a segment $`a<\rho _N\rho \rho _D`$ of the $`c`$-orbit in which the function $`H(\rho )`$ is monotone decreasing from a local maximum $`H(\rho _N)>0`$ to $`H(\rho _D)=0`$.

We define

$$W=\frac{1}{2}\rho ^2AH_{}^{}{}_{}{}^{2}(1Z^2)+\mathrm{sin}^2H.$$
(47)
From (24) $`W=(A1+2\alpha )/(2\alpha )`$, hence by (22) $`W^{}<0`$. We have

$$\frac{H_{}^{}{}_{}{}^{2}}{W\mathrm{sin}^2H}=\frac{2}{\rho ^2A(1Z^2)},\text{so}\frac{H^{}}{\sqrt{W\mathrm{sin}^2H}}=\frac{\sqrt{2}}{\rho \sqrt{A(1Z^2)}}.$$
(48)
Integrating the left-hand side from $`\rho _N`$ to $`\rho _D`$, we get (using $`H_N=H(\rho _N)`$)

$$_{\rho _N}^{\rho _D}\frac{H^{}d\rho }{\sqrt{W\mathrm{sin}^2H}}=_0^{H_N}\frac{dH}{\sqrt{W\mathrm{sin}^2H}}_0^{H_N}\frac{dH}{\sqrt{\mathrm{sin}^2H_N\mathrm{sin}^2H}}>\frac{\pi }{2},$$
(49)
where the first inequality follows from $`W(\rho )W(\rho _N)=\mathrm{sin}^2H_N`$ (since $`W^{}`$ decreases) and the second inequality is a simple calculation using a substitution $`\mathrm{sin}H=u\mathrm{sin}H_N`$ (remember that $`H_N<\pi /2`$).

Next, we derive an upper bound for the integral of the right-hand side of (48). We have

$$_{\rho _N}^{\rho _D}\frac{d\rho }{\rho \sqrt{A(1Z^2)}}\frac{1}{\rho _N\sqrt{12\alpha }}_{\rho _N}^{\rho _D}\frac{d\rho }{\sqrt{1Z^2}}\frac{1}{\rho _N\sqrt{12\alpha }}_{\rho _N}^{\rho _D}\frac{d\rho }{\sqrt{1Z}}.$$
(50)
We showed above that $`Z^{}>1`$, hence $`1Z\rho _1\rho `$. Therefore

$$_{\rho _N}^{\rho _D}\frac{d\rho }{\sqrt{1Z}}_{\rho _N}^{\rho _D}\frac{d\rho }{\sqrt{\rho _1\rho }}=2(\sqrt{\rho _1\rho _N}\sqrt{\rho _1\rho _D})<2\sqrt{\rho _1\rho _N}.$$
(51)
By continuity of solutions with respect to $`c`$ and by Corollary 3, $`\rho _N`$ is arbitrarily close to $`\rho _1(c)`$ if $`c`$ is sufficiently close to $`\stackrel{~}{c}`$, hence it follows from (51) that the integral of the right-hand side of (48) is arbitrarily small. This is in contradiction with (49), hence $`H(\rho ,c)`$ cannot have a second additional zero, which completes the proof of Lemma 8.

###### Lemma 9.

If $`\overline{H}(\stackrel{~}{c})0`$, then $`N(c)=N(\stackrel{~}{c})`$ for $`c`$ sufficiently close to $`\stackrel{~}{c}`$.

*Proof:* Without loss of generality we assume that $`\overline{H}(\stackrel{~}{c})<0`$. As above let $`\stackrel{~}{a}<\rho _1(\stackrel{~}{c})`$ be the last node of the $`\stackrel{~}{c}`$-orbit, that is $`H(\stackrel{~}{a},\stackrel{~}{c})=0`$ and $`H(\rho ,\stackrel{~}{c})<0`$ for $`\stackrel{~}{a}<\rho \rho _1`$. Let $`a`$ be the corresponding zero of $`H(\rho ,c)`$ for $`c`$ near $`\stackrel{~}{c}`$. We want to show that $`H(\rho ,c)`$ cannot have an extra zero for $`\rho >a`$. Suppose for contradiction that $`H(b,c)=0`$ for some $`b>a`$. Then, there must be a $`\delta <b`$ such that $`H(\delta ,c)=\overline{H}(\stackrel{~}{c})`$. Let us integrate the identity

$$\frac{H^{}}{\sqrt{W\mathrm{sin}^2H}}=\frac{\sqrt{2}}{\rho \sqrt{A(1Z^2)}}$$
(52)
from $`\delta `$ to $`b`$. For the left hand side we get

$$_\delta ^b\frac{H^{}d\rho }{\sqrt{W\mathrm{sin}^2H}}=_0^{\overline{H}}\frac{dH}{\sqrt{W\mathrm{sin}^2H}}.$$
(53)
From Proposition 2 we know that $`lim_{\rho \rho _1}(1Z^2)H_{}^{}{}_{}{}^{2}=0`$ so $`W(\rho ,\stackrel{~}{c})<(1+ϵ/2)\mathrm{sin}^2\overline{H}`$ for $`\rho `$ near $`\rho _1`$ and hence $`W(\rho ,c)<(1+ϵ)\mathrm{sin}^2\overline{H}`$ for $`c`$ near $`\stackrel{~}{c}`$. Since $`W`$ is decreasing, $`W(\delta ,c)<W(\rho ,c)<(1+ϵ)\mathrm{sin}^2\overline{H}`$. Thus

$$_0^{\overline{H}}\frac{dH}{\sqrt{W\mathrm{sin}^2H}}_0^{\overline{H}}\frac{dH}{\sqrt{(1+ϵ)\mathrm{sin}^2\overline{H}\mathrm{sin}^2H}}\mathrm{arcsin}\left(\frac{1}{\sqrt{1+ϵ}}\right)>\frac{\pi }{2}$$
(54)
for sufficiently small $`ϵ`$, where the last but one inequality can be seen by substituting $`\mathrm{sin}H=u\mathrm{sin}\overline{H}`$ into the integral. By the same argument as in the proof of Lemma 8, the integral of the right hand side of (52) is $`O(\sqrt{\rho _1\rho })`$. By continuity of solutions with respect to $`c`$ and by Corollary 3, $`\delta `$ is arbitrarily close to $`\rho _1(c)`$ if $`c`$ is sufficiently close to $`\stackrel{~}{c}`$, hence the integral of the left hand side of equation (52) is arbitrarily small. This contradicts (54) and completes the proof of Lemma 9. Now we are ready to make a shooting argument. We define a set

$$C_0=\{c|N(c)=0\}$$
(55)
and let $`c_0=supC_0`$. The set $`C_0`$ is nonempty (by Proposition 6) and bounded above (by Proposition 7) so $`c_0`$ exists. We claim that the $`c_0`$-orbit has no nodes and satisfies the boundary condition $`\overline{H}(c_0)=0`$. To see this, note that the $`c_0`$-orbit cannot have a node because then by Lemmas 8 and 9 all nearby $`c`$-orbits would have a node so there would be an interval around $`c_0`$ without any elements of $`C_0`$ in it, contradicting the assumption that $`c_0`$ is the least upper bound. Thus, $`N(c_0)=0`$. Now, if $`\overline{H}(c_0)<0`$, then by Lemma 9 all nearby $`c`$-orbits would have no nodes, so there would be an interval around $`c_0`$ consisting of elements of $`C_0`$, contradicting the assumption that $`c_0`$ is an upper bound of $`C_0`$. Thus $`\overline{H}(c_0)=0`$.

Next, we define $`C_1=\{c>c_0|N(c)=1\}`$. This set is nonempty by the previous step and Lemma 8 and bounded above by Proposition 7, hence $`c_1=supC_1`$ exists. By the same argument as above, the $`c_1`$-orbit has exactly one node and satisfies $`\overline{H}(c_1)=0`$. The construction of subsequent $`c_n`$-orbits proceeds by induction. *Conclusion of the proof of Theorem 1:*

Returning to the original variable $`F(\rho )`$ and rescaling $`\rho \rho /\rho _1(c_n)`$ we get the solution of equations (13-17) which satisfies the boundary conditions (19) and (20) and has exactly $`n`$ intersections with the line $`F=\pi /2`$. By Corollary 5 this solution is analytic in the whole interval $`0\rho 1`$.

## 4 Beyond the past self-similarity horizon

In this section we consider the behaviour of the CSS solutions of Theorem 1 outside the past SSH, in particular we ask the question: do these solutions possess a regular future self-similarity horizon? Note that $`\rho =\mathrm{}`$ corresponds to the hypersurface $`(t=0,r>0)`$ so in order to analyze the global behaviour of solutions (for $`t>0`$) we need to go ”beyond $`\rho =\mathrm{}`$”. To this end we define, after I, a new coordinate $`x`$ by

$$\frac{d}{dx}=\rho Z\frac{d}{d\rho },x(\rho =1)=0.$$
(56)
We also define an auxiliary function $`w(x)=1/Z(\rho )`$. In these new variables, the past SSH where $`w=1`$ is at $`x=0`$, while the future SSH (if it exists) is at some $`x_A>0`$ where $`w(x_A)=1`$.

In terms of $`x`$ and $`w`$, the equations (21)-(23) become autonomous (where now prime is $`d/dx`$)

$`H^{\prime \prime }`$ $``$ $`2\alpha wH_{}^{}{}_{}{}^{3}+{\displaystyle \frac{\mathrm{sin}(2H)}{A(w^21)}}=0,`$ (57)
$`A^{}`$ $`=`$ $`2\alpha AwH_{}^{}{}_{}{}^{2},`$ (58)
$`w^{}`$ $`=`$ $`1+\alpha (1w^2)H_{}^{}{}_{}{}^{2}.`$ (59)
The constraint (24) becomes

$$12\alpha A+2\alpha \mathrm{sin}^2H+\alpha AH_{}^{}{}_{}{}^{2}(w^21)=0.$$
(60)
From (20) the initial conditions at $`x=0`$ are

$$H(x)bx,w(x)1x,A(x)12\alpha 2\alpha (12\alpha )b^2x.$$
(61)
We know from Theorem 1 that for each $`\alpha <1/2`$ there is an infinite sequence $`\{b_n(\alpha )\}`$ determining solutions which are regular inside the past SSH, that is for all $`x0`$ (note that $`\rho =0`$ corresponds to $`x=\mathrm{}`$). In I we showed that for $`x>0`$ the solutions starting from the past SSH with the initial conditions (61) tend in finite ”time” to $`w=1`$ if $`b`$ is small, or to $`w=+1`$ if $`b`$ is large. After I we shall refer to these two kinds of behaviour as to type A and type B solutions, respectively. Now we want to show that the solutions of Theorem 1 are of type A (and therefore possess the future SSH) provided that $`\alpha `$ is sufficiently small. Unfortunately, the shooting argument gives us insufficient information about the parameters $`b_n`$ so we cannot apply the above mentioned result of I to determine the character of solutions outside the past SSH. Instead, we shall make use of the obvious fact that for $`\alpha =0`$ all solutions are of type A.

###### Lemma 10.

For sufficiently small $`\alpha `$ the $`c_n`$-orbits of Theorem 1 (rescaled so that $`\rho _1(c)=1`$) have $`|b_n|`$ uniformly bounded above for all $`n`$.

*Proof:* It was shown in  (see Lemma 4 in that reference) that for $`\alpha =0`$ the solution to equations (57)-(61) for $`x<0`$ must exit the strip $`|H|\pi /2`$ if $`|b|`$ is too large, say $`|b|>B`$. By continuous dependence, the same is true for sufficiently small $`\alpha `$. But from Proposition 2 the $`c`$-orbits must stay in the strip $`|H|\pi /2`$ for all $`x<0`$. Thus, $`|b_n|B`$ for small $`\alpha `$.

###### Lemma 11.

If a solution to equations (57)-(60) has $`w(x_0)<0`$ and $`A(x_0)>2/3`$ for some $`x_0`$, then there is $`x_A>x_0`$ such that $`lim_{xx_A}w(x)=1`$, i.e., the solution is of type A.

*Proof:* By (58) $`A`$ is increasing for $`w<0`$. Thus, using equation (59) and the constraint (60) we get for $`x>x_0`$

$$w^{}=1+\alpha (1w^2)H_{}^{}{}_{}{}^{2}=1+\frac{1A2\alpha \mathrm{cos}^2H}{A}<2+\frac{1}{A}\frac{1}{2},$$
(62)
hence $`w`$ must hit $`1`$ for some finite $`x_A>x_0`$.

###### Proposition 12.

The $`c_n(\alpha )`$-orbits are of type A if $`\alpha `$ is sufficiently small.

*Proof:* For $`\alpha =0`$ and any $`b`$ we have $`w(x)=1x`$ and $`A(x)1`$; in particular $`A(3/2)=1>2/3`$ and $`w(3/2)=1/2<0`$. By continuous dependence on initial conditions there exists a $`\delta (b)`$ such that if $`\alpha <\delta (b)`$ and $`|bb^{}|<\delta (b)`$, then $`A(3/2,b^{})>2/3`$ and $`w(3/2,b^{})<0`$. This implies by Lemma 11 that the solutions corresponding to such values of $`\alpha `$ and $`b^{}`$ are of type A. By a standard theorem of advanced calculus there is a $`\delta ^{}>0`$ (independent of $`b`$) such that the solutions with $`\alpha <\delta ^{}`$ and $`|b|B`$ are of type A. By Lemma 10 any $`c_n`$-orbit has $`|b|B`$, so for $`\alpha <\delta ^{}`$ the $`c_n`$-orbits are of type A. By a similar argument as in the proof of Proposition 2 one can easily show that the type A solutions are generically only $`C^0`$ at the future SSH (for isolated values of $`\alpha `$ there are solutions that go smoothly through the future SSH). In I we showed that the leading order asymptotic behaviour at the future SSH is (using $`y=x_Ax`$)

$$w1+y,AA_02\alpha A_0C^2y\mathrm{ln}^2(y),HH_0Cy\mathrm{ln}(y),$$
(63)
where $`A_0=12\alpha \mathrm{cos}^2H_0`$, $`C=\mathrm{sin}(2H_0)/2A_0`$, and $`H_0`$ is a free parameter. Using this expansion one can check that the curvature is finite as $`y0`$ which means that the type A solutions are examples of naked singularities.

## Acknowledgments

We acknowledge the hospitality of the Erwin Schrödinger Institute for Mathematical Physics in Vienna, where part of this paper was produced. The research of PB was supported in part by the KBN grant 2 P03B 010 16.

## Appendix (local existence theorems)

In  (Proposition 1) Breitelohner, Forgács, and Maison have derived the following result concerning the behaviour of solutions of a system of ordinary differential equations near a singular point (see also  for a similar result):

###### Theorem (BFM).

Consider a system of first order differential equations for $`n+m`$ functions $`u=(u_1,\mathrm{},u_n)`$ and $`v=(v_1,\mathrm{},v_m)`$

$$t\frac{du_i}{dt}=t^{\mu _i}f_i(t,u,v),t\frac{dv_i}{dt}=\lambda _iv_i+t^{\nu _i}g_i(t,u,v),$$
(64)
where constants $`\lambda _i>0`$ and integers $`\mu _i,\nu _i1`$ and let $`C`$ be an open subset of $`R^n`$ such that the functions $`f`$ and $`g`$ are analytic in the neighbourhood of $`t=0,u=c,v=0`$ for all $`cC`$. Then there exists an $`n`$-parameter family of solutions of the system (64) such that

$$u_i(t)=c_i+O(t^{\mu _i}),v_i(t)=O(t^{\nu _i}),$$
(65)
where $`u_i(t)`$ and $`v_i(t)`$ are defined for all $`cC,|t|<t_0(c)`$ and are analytic in $`t`$ and $`c`$.

We shall use this theorem to prove existence of local solutions of equations (21)-(23) near the singular points $`\rho =0`$ and $`\rho =1`$.

###### Proposition 13.

The equations (21)-(23) admit a two-parameter family of local solutions near $`\rho =0`$

$`H(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi }{2}}+c\rho +O(\rho ^3),`$ (66)
$`A(\rho )`$ $`=`$ $`1\alpha c^2\rho ^2+O(\rho ^4),`$ (67)
$`Z(\rho )`$ $`=`$ $`d\rho +O(\rho ^3),`$ (68)
which are analytic in $`c,d`$ and $`\rho `$.

*Proof:* Using the variables

$$w_1=\frac{H+\pi /2}{\rho },w_2=H^{},w_3=\frac{1A}{\rho ^2},w_4=\frac{Z}{\rho }$$
(69)
we rewrite the equations (21)-(23) as the first order system

$`\rho w_1^{}`$ $`=w_1+w_2,`$ $`\rho w_2^{}`$ $`=2w_12w_2+\rho ^2h_1,`$
$`\rho w_3^{}`$ $`=2w_3+2\alpha w_2^2+\rho ^2h_2,`$ $`\rho w_4^{}`$ $`=\rho ^2h_3,`$ (70)
where the functions $`h_i`$ are analytic near $`\rho =0`$. Next, substituting

$`w_1`$ $`=u_1v_1,`$ $`w_2`$ $`=u_1+2v_1,`$
$`w_3`$ $`=v_2+\alpha (u_1^22v_1^28u_1v_1),`$ $`w_4`$ $`=u_2`$ (71)
we put (Appendix (local existence theorems)) into the form (64)

$`\rho u_1^{}`$ $`=\rho ^2f_1,`$ $`\rho u_2^{}`$ $`=\rho ^2f_2,`$
$`\rho v_1^{}`$ $`=3v_1+\rho ^2g_1,`$ $`\rho v_2^{}`$ $`=2v_2+\rho ^2g_2,`$ (72)
where the functions $`f_i,g_i`$ are analytic in an open neighbourhood of $`\rho =0,u_1=c,u_2=d,v_i=0`$ for any $`c`$ and $`d`$. Thus, according to the BFM theorem there exists a two-parameter family of solutions such that

$`u_1`$ $`=c+O(\rho ^2),`$ $`u_2`$ $`=d+O(\rho ^2),`$ (73)
$`v_1`$ $`=O(\rho ^2),`$ $`v_2`$ $`=O(\rho ^2),`$ (74)
which is equivalent to (66)-(68).

###### Proposition 14.

The equations (21)-(23) admit a one-parameter family of local solutions near $`\rho =1`$

$`H(\rho )`$ $`=`$ $`b(\rho 1)+O((\rho 1)^2),`$ (75)
$`A(\rho )`$ $`=`$ $`12\alpha 2\alpha (12\alpha )b^2(\rho 1)+O((\rho 1)^2),`$ (76)
$`Z(\rho )`$ $`=`$ $`\rho +O((\rho 1)^2)`$ (77)
which are analytic in $`b`$ and $`\rho `$.

*Proof:* We define the variables

$`u`$ $`=H^{},`$ $`v_1`$ $`={\displaystyle \frac{H}{\rho 1}}H^{},`$ (78)
$`v_2`$ $`={\displaystyle \frac{(12\alpha )A}{\rho 1}}2\alpha (12\alpha )H_{}^{}{}_{}{}^{2},`$ $`v_3`$ $`={\displaystyle \frac{Z1}{\rho 1}}1.`$ (79)
Then, the equations (21)-(23) take the form (using $`t=\rho 1`$)

$$tu^{}=tf,tv_i^{}=v_i+tg_i,$$
(80)
where the functions $`f`$ and $`g_i`$ are analytic in an open neighourhood of $`t=0,u=b,v_i=0`$ for any $`b>0`$. Thus, according to the BFM theorem there exists a one-parameter family of solutions such that

$$u(t)=b+O(t),v_i(t)=O(t),$$
(81)
which is equivalent to (75)-(77).


text:
# A simple model for 1/𝑓^𝛼 noise

\[

## Abstract

We present a simple stochastic mechanism which generates pulse trains exhibiting a power law distribution of the pulse intervals and a $`1/f^\alpha `$ power spectrum over several decades at low frequencies with $`\alpha `$ close to one. The essential ingredient of our model is a fluctuating threshold which performs a Brownian motion. Whenever an increasing potential $`V(t)`$ hits the threshold, $`V(t)`$ is reset to the origin and a pulse is emitted. We show that if $`V(t)`$ increases linearly in time, the pulse intervals can be approximated by a random walk with multiplicative noise. Our model agrees with recent experiments in neurobiology and explains the high interpulse interval variability and the occurrence of $`1/f^\alpha `$ noise observed in cortical neurons and earthquake data.

\]

The omnipresence of $`1/f^\alpha `$ noise in nature is one of the oldest puzzles in contemporary physics still lacking a generally accepted explanation. The phenomenon is characterized by a certain behavior of the respective time signal: The power spectrum $`S(f)`$ is proportional to $`1/f^\alpha `$ at low frequencies $`f`$ with $`\alpha 1`$. Examples include the light of quasars , electrical measurements , music and speech , human cognition  and coordination , the current through ion channels , network traffic , burst errors in communication systems , freeway traffic , granular flow , *etc*. The time signals of a large number of these systems  resemble a pulse train consisting of individual, largely identical events which occur at discrete times. This is especially true for spike trains of single nerve cells for which $`1/f^\alpha `$ noise has been observed in various brain structures . The reported exponents, which depend both on the presence or absence of a sensory stimulus  and on the state of the animal (REM sleep vs awake state) , vary from 0.68 to 1.38. The power-law behavior for spike train power spectra lies within the range 0.01 to 10 Hz, extending typically over 2 decades. In almost all cases the upper limit of the observed time over which fractal correlations exist is imposed by the duration of the recording.

In this Letter, we propose a simple mechanism for generating pulse trains with $`1/f^\alpha `$ behavior in systems with a threshold-controlled dynamics like, e.g., neurons and earthquake faults. Our model is based on an integrate-and-fire (IaF) mechanism and consists of a single unit characterized by two variables (see Fig. 1): The voltage $`V(t)`$ and the threshold $`C(t)`$. Initially, the voltage is below the threshold. Then, the voltage increases monotonically in time — in the simplest case just linearly — and the threshold evolves according to a Brownian motion with diffusion constant $`D`$ within reflecting boundaries $`V_0<C_l<C(t)<C_u`$. As soon as $`V(t)`$ has reached the threshold, the voltage is reset to $`V_0`$ and a pulse of unit height is emitted. In this way, a pulse train is generated. Note that the threshold is *not* reset to its initial value.

Such a model is often used to describe single neurons: $`V(t)`$ is the membrane potential and the emitted pulse is the generated action potential. However, it is usually assumed that the threshold is constant in time. Recent investigation have shown that this is not true for cortical neurons *in vivo*  and *in vitro* . Additionally, there is evidence that the spike trains of auditory neurons  and of neurons in the mesencephalic reticular formation  are not renewal, i.e., successive time intervals between spikes are correlated. These facts are incorporated in our model in the simplest possible way. Moreover, the effects of dead time or absolute refractoriness, which limits the rate at which a neuron can fire, are automatically included in our model via the lower bound $`C_l`$. The upper bound $`C_u`$ prevents an infinite time difference between two pulses.

As one can see from Fig. 2, the power spectrum of the pulse train generated by our model with linear increasing voltage shows a $`1/f^\alpha `$ decay over several decades . The frequency below which white noise behavior is observed is determined by $`C_uC_l`$ and goes to zero for $`C_uC_l\mathrm{}`$. The exponent $`\alpha `$ is not universal and increases if the ratio $`(C_lV_0)/\sqrt{D}`$ increases. In the limit $`(C_lV_0)\sqrt{D}`$, we find $`\alpha =0.5`$. This result is explained by the fact that this limit corresponds to the case in which the waiting times between pulses are purely determined by the diffusive dynamics of the threshold. Hence, each waiting time is merely the first return time of a Brownian motion which obeys a power-law distribution with exponent $`1.5`$. This implies $`\alpha =0.5`$ .

To compare our results with real neurons, consider the parameters of curve (b): The dead time of a neuron is typically of order of milliseconds and the maximal time difference between two spikes of the order of seconds. This is exactly the ratio between $`C_l`$ and $`C_u`$. Moreover, we can now identify one unit of time in our model with 5 milliseconds real time. Hence, the $`1/f^\alpha `$ behavior in (b) is found for $`f<`$ 11 Hz. This and the fact that $`\alpha 1.02`$ reproduces the experimental results very well.

Extensive numerical investigations have shown that our model is very stable with respect to variations of the dynamics. Different forms of voltage increase, e.g., a linear, a squared, or a square-root increase, give similar results. The assumption of a monotonical increase of the voltage can also be dropped. A small amount of noise can be added to the voltage signal without altering our findings. Substituting the reflecting boundaries by a confining potential does not change our results, either. This points towards a generic behavior.

To obtain the power spectrum of our model analytically, consider the signal $`X(t)=_k\delta (tt_k)`$ where $`t_k`$ denotes the time of occurrence of the $`k`$th pulse. It follows as shown in

$`S(f)`$ $`=`$ $`\underset{T\mathrm{}}{lim}{\displaystyle \frac{1}{2T}}\left|{\displaystyle \underset{T}{\overset{T}{}}}𝑑tX(t)\mathrm{exp}^{i2\pi ft}\right|^2,`$ (1)
$`=`$ $`\underset{T\mathrm{}}{lim}{\displaystyle \frac{1}{2T}}{\displaystyle \underset{k}{}}{\displaystyle \underset{q}{}}\mathrm{exp}^{i2\pi f(t_{k+q}t_k)}.`$ (2)
With $`\overline{I}=lim_T\mathrm{}\frac{1}{2T}(k_{\mathrm{max}}k_{\mathrm{min}}+1)`$, this leads to

$`S(f)`$ $``$ $`\overline{I}{\displaystyle \underset{q}{}}\mathrm{exp}^{i2\pi f(t_{k+q}t_k)},`$ (3)
where $`\mathrm{}`$ denotes the average over the ensemble and over $`k`$. Hence, we need to know the probability distributions $`\mathrm{\Psi }_q(\tau )d\tau `$ of the time differences between pulses $`\tau =t_kt_{k+q}`$ for all integers $`q`$. $`\mathrm{\Psi }_q`$ is merely the $`q`$th passage times density function $`g_q`$ averaged over the stationary probability distribution of “initial” states $`V(0)`$

$`\mathrm{\Psi }_q(\tau )=g_q(\tau |V(0))_{V(0)}.`$ (4)
Since $`g_q`$ can be computed from the first passage times density function (FPTDF), we will first focus on the latter.

The FPDTF of our model with *linear* voltage increase can be obtained by mapping the model to an IaF model with constant threshold $`C=(C_u+C_l)/2`$: The voltage $`\stackrel{~}{V}(t)`$ is defined as the sum of $`V(t)`$ and $`CC(t)`$. This means that $`\stackrel{~}{V}(t)`$ fluctuates around $`V(t)`$. It also implies that $`\stackrel{~}{V}(t)`$ is reset to $`V_0+CC(t_k)`$ after the $`k`$th threshold crossing. Hence, the correlations are now encoded in the fluctuating reset. In conclusion, $`\stackrel{~}{V}(t)`$ behaves almost as a Brownian motion with drift to an absorbing barrier with a reset following each barrier crossing. The only difference is that the stochastic process is restricted to the interval $`[V(t)+CC_u,V(t)+CC_l]`$ at time $`t`$ which makes an exact mathematical treatment difficult. Neglecting the restriction for a moment and setting $`C=0`$, we obtain from the associated Fokker-Planck equation the FPTDF

$`g_1(\tau |\stackrel{~}{V}(0))={\displaystyle \frac{\stackrel{~}{V}(0)}{\sqrt{2\pi D}\tau ^{3/2}}}\mathrm{exp}{\displaystyle \frac{(\stackrel{~}{V}(0)+\tau )^2}{2D\tau }},`$ (5)
and $`g_10`$ for $`\tau <0`$. Here, $`\stackrel{~}{V}(0)`$ is the initial distance. Eq. (5) is, of course, just an approximation to the FPTDF of our model. This can already be seen from the fact that the FPTDF of our model is identically zero for $`\tau <(C_lV_0)`$ and $`\tau >(C_uV_0)`$. However, as will be shown below this approximation proves to be very useful and Eq. (5) can even be simplified to a Gaussian with the same mean and variance:

$`g_1(\tau |\stackrel{~}{V}(0))\mathrm{exp}{\displaystyle \frac{(\stackrel{~}{V}(0)+\tau )^2}{2D|\stackrel{~}{V}(0)|}}.`$ (6)
From the FPTDF one could now compute higher passage times density functions $`g_n`$ in principle by convolution

$`g_n(\tau |\stackrel{~}{V}(0))={\displaystyle _0^\tau }g_1(t|\stackrel{~}{V}(0))g_{n1}(\tau t|t)𝑑t.`$ (7)
However, this convolution resists an analytical treatment due to the Markovian character of our process and the form of Eqs. (6,7). The former manifests itself in the fact that the first passage time directly gives the initial distance between the threshold and the voltage for the next passage problem. The reason for this is that the norm of the $`k`$th reset of the voltage $`\stackrel{~}{V}`$ equals the time difference between the $`k`$th and the $`(k1)`$th pulse due to the linear increase of the voltage $`\stackrel{~}{V}`$ with slope 1. Hence the reset depends on the last passage time only and the whole stochastic point process is totally described by the FPTDF plus its Markovian property.

These properties enable us to interpret the stochastic process generated by the Gaussian approximation (6) as a random walk of the inter-spike-intervals (ISI) $`\tau _k`$:

$`\tau _{k+1}=\tau _k+\sqrt{D\tau _k}\eta _k.`$ (8)
Here, $`\eta _k`$ denotes the white noise source. Note the special kind of multiplicative noise which distinguishes our model from the one in . The variance of the step length is proportional to the current “position”. Hence, the origin is a fixed point of the random walk. To omit this difficulty and in spirit of our original model, we consider the random walk $`\tau _k`$ to be confined by two reflecting boundaries, i.e., $`0<(C_uV_0)<\tau _k<(C_lV_0)`$. The power spectrum of the pulse-train generated by such a random walk with $`t_{k+1}=t_k+\tau _k`$ shows a clear $`1/f^\alpha `$ behavior with the same exponent as for our model (see Fig. 3). This is also true if the new inter-spike-interval is chosen from an inverse Gaussian distribution given in Eq. (5) with an initial condition depending on the last interval. Hence, these approximations seem to be justified. This is further confirmed by the stationary ISI distribution. For the random walk, the ISI distribution function $`P(\tau )`$ is proportional to $`\tau ^1`$ . Simulations show that this is in excellent agreement with our model and the inverse Gaussian approximation.

We have to point out that the behavior of $`P(\tau )`$ depends on the specific kind of voltage increase. For a voltage increase with $`(tt_{last})^\beta `$ with $`0<\beta <2`$, we find $`P(\tau )\tau ^{\beta 2}`$ for our model. Substituting the reflecting boundaries by a potential of the form $`a/x+bx^2`$ for example, we still find power law tails in the ISI distribution. This is exactly what measurements show for cortical neurons , especially in the mesencephalic reticular formation . Hence, our model seems to be especially well-suited to explain the occurrence of $`1/f^\alpha `$ noise in that formation. In general, IaF models can not explain the high interspike interval variability exhibited by cortical neurons . A fluctuating threshold as described by our model, however, can solve this long-standing issue.

For renewal processes the ISI distribution function and the FPTDF are identical and completely describe the process. This is the case for a standard IaF neuron with reset of the voltage to the origin after each barrier crossing which is described by the FPTDF in Eq. (5). Such a model can neither explain a $`1/f^\alpha `$ signal nor a power-law decay of the ISI distribution. In contrast to our model, Usher and co-workers showed that fractal behavior might be a consequence of the global activity dynamics of a network of IaF neurons . Due to experimental limitations, however, the link between $`1/f^\alpha `$ single-unit power spectra and macroscopic activity dynamics remains, as of now, a conjecture. Another attempt to explain the phenomenon of $`1/f^\alpha `$ noise is based on fractal and fractal-rate stochastic point processes . Certain types of these processes properly characterize the statistical properties found in different experiments. However, in many cases it is not clear *a priori* why the real system should generate such a process. This is especially true for the clustering Poisson process which was applied to explain $`1/f^\alpha `$ noise in the mesencephalic reticular formation . In contrast to that, our model provides a simple explanation of $`1/f^\alpha `$ noise in integrate-and-fire systems in general and, hence, applies to neurons in particular. The crucial assumption is a fluctuating threshold. Such a behavior is related to models presented in  and was already considered in  to explain the power law decay of neuronal ISI distribution functions. However, the stochastic point process was still assumed to be renewal and could not explain a $`1/f^\alpha `$ behavior of the power spectral density function. Consequently, the second crucial assumption of our model is the Markovian character of the process in agreement with Refs. . To directly verify our model for neurons in the mesencephalic reticular formation, one should measure the evolution of the reset from the threshold potential.

All three main characteristics of our model, i.e., a $`1/f^\alpha `$ spectrum, a power-law decay of the ISI distribution, and a Markovian dependence of the ISI, can not only be observed in nerve cells but also in other systems. For example, the pulse signal, defined by associating the occurrence of pulses of unit height with earthquakes in the Mojave region , shows $`1/f^\alpha `$ fluctuations with $`\alpha =1.3`$ and has an ISI distribution with exponent $`1.1`$ (Fig. 4). Moreover, there is evidence that the ISI do not obey a renewal process. Rather a Markovian dynamics can be found . These findings imply that the effective evolution of the ISI is similar to the one generated by our model and caricatured by Eq. (8).

To summarize, we have presented a simple stochastic model which is able to transform integrated white noise with a $`1/f^2`$ power spectrum into noise with $`1/f`$ tail. The basic mechanism is similar to what is expected to describe the dynamics of single neurons and, thus, our model can explain the behavior of neurons in the central-nervous-system . The analysis of earthquake data even suggests that the effective dynamics of the ISI is the same for many systems exhibiting $`1/f^\alpha `$ noise.

We thank A. Aertsen and C. Goltz for useful discussions. J. D. would like to thank the Land Schleswig-Holstein, Germany, for financial support.


text:
# LMA parameters and non-zero 𝑈_{𝑒3} effects on atmospheric 𝜈 data?

## 1 Introduction

Recent results on atmospheric neutrinos  as well as results from the long base-line experiment K2K  further confirmed the interpretation of the atmospheric neutrino anomaly in terms of $`\nu _\mu \nu _\tau `$ oscillations with maximal or close to maximal mixing and mass squared difference in the interval, $`\mathrm{\Delta }m_{atm}^2=(1.54)\times 10^3\mathrm{eV}^2,\mathrm{sin}^22\theta _{atm}>0.88,`$ at 90 % C.L. .

A sub-dominant oscillation of electron neutrinos is not excluded yet. It seems that there is an excess of the $`e`$like events in the low energy part of the sub-GeV sample ($`p<0.4`$ GeV, where $`p`$ is the momentum of lepton). In comparison with predictions based on the atmospheric neutrino flux from Ref. the excess is about (12 - 15)%. For higher energies, the excess is much smaller.

Can these results be related to the $`\nu _e`$oscillations? What could be the implications of the positive answer? We have some preliminary results that we will discuss in next sections.

## 2 Formalism

In the three neutrino schemes which explain the atmospheric and solar neutrino data, there are two possible channels of the $`\nu _e`$ oscillations:

1. $`\nu _e`$oscillations driven by $`\mathrm{\Delta }m_{atm}^2`$ responsible for dominant mode of the atmospheric neutrino oscillations . These oscillations require non-zero value of $`U_{e3}`$. The effects are restricted by the CHOOZ result .

2. $`\nu _e`$oscillations driven by the solar mass splitting $`\mathrm{\Delta }m_{}^2`$ .

The detailed study of the effect have been done in our previous paper  where we have shown that neutrino oscillations with parameters in the LMA MSW allowed region  $`\mathrm{\Delta }m_{}^2=(230)10^5`$ eV2, $`\mathrm{sin}^22\theta _{}>0.65`$, favored by analyzes of solar neutrino data from SNO  and Super-Kamiokande  data, can lead to an observable excess of the e-like events in the sub-GeV atmospheric neutrino sample.

It was shown that the excess is determined by the two neutrino transition probability $`P_2`$ and the “screening” factor:

$$\frac{F_e}{F_e^0}1=P_2(rc_{23}^21),$$
(1)
where $`F_e`$ and $`F_e^0`$ are the electron neutrino fluxes with and without oscillations and $`r`$ is the ratio of the original muon and electron neutrino fluxes. In the sub-GeV region $`r2`$, so that the screening factor is zero when the $`\nu _\mu \nu _\tau `$ mixing is maximal. We show in Fig. 1 our previous results compared with the latest data on Super-Kamiokande .

In previous studies the effects of oscillations driven by the solar and atmospheric $`\mathrm{\Delta }m^2`$ have been considered separately: The studies of the $`\mathrm{\Delta }m_{atm}^2`$driven oscillations where performed in the framework of the so called “one level dominating scheme” when the effect of solar mass splitting between two lightest states, $`\mathrm{\Delta }m_{21}^2`$, is neglected. In studies of the solar $`\mathrm{\Delta }m_{21}^2`$ driven oscillations it was assumed that $`U_{e3}`$ is negligible.

In this paper we study the effects of the interplay of oscillations with the LMA parameters and non-zero $`U_{e3}`$.

## 3 $`U_{e3}`$ and induced interference

We consider the three-flavor neutrino system with hierarchical mass squared differences: $`\mathrm{\Delta }m_{21}^2=\mathrm{\Delta }m_{}^2<<\mathrm{\Delta }m_{31}^2=\mathrm{\Delta }m_{atm}^2`$. The evolution of the neutrino vector of state $`\nu _f(\nu _e,\nu _\mu ,\nu _\tau )^T`$ is described by the equation

$$i\frac{d\nu _f}{dt}=\left(\frac{UM^2U^{}}{2E}+V\right)\nu _f,$$
(2)
where $`E`$ is the neutrino energy and $`M^2=diag(0,\mathrm{\Delta }m_{21}^2,\mathrm{\Delta }m_{31}^2)`$ is the diagonal matrix of neutrino mass squared eigenvalues. $`V=diag(V_e,0,0)`$ is the matrix of matter-induced neutrino potentials with $`V_e=\sqrt{2}G_FN_e`$, $`G_F`$ and $`N_e`$ being the Fermi constant and the electron number density, respectively. The mixing matrix $`U`$ is defined through $`\nu _f=U\nu _{mass}`$, where $`\nu _{mass}(\nu _1,\nu _2,\nu _3)^T`$ is the vector of neutrino mass eigenstates. It can be parameterized as $`U=U_{23}U_{13}U_{12}.`$ The matrix $`U_{ij}=U_{ij}(\theta _{ij})`$ performs the rotation in the $`ij`$\- plane by the angle $`\theta _{ij}`$. Here we have neglected possible CP-violation effects in the lepton sector.

### 3.1 Propagation basis

The dynamics of oscillations is simplified in the “propagation” basis $`\stackrel{~}{\nu }=(\stackrel{~}{\nu }_e,\stackrel{~}{\nu }_2,\stackrel{~}{\nu }_3)^T`$, which is related with the flavor basis by $`\nu _f=\stackrel{~}{U}\stackrel{~}{\nu }`$. We define the propagation basis in such a way that projection matrix $`\stackrel{~}{U}`$ equals: $`\stackrel{~}{U}=U_{23}U_{13}`$ . The propagation basis can be introduced in the following way. First, let us perform the rotation $`\nu _f=U_{23}U_{13}\nu ^{}`$. Using Eq. (2) we find that in the basis $`\nu ^{}`$ the Hamiltonian takes the form,

$$H^{}\left(\begin{array}{cc}H_2& 0\\ 0& \mathrm{\Delta }m_{31}^2/2E+V_es_{13}^2\end{array}\right),$$
(3)
where $`H_2=U_{12}M_2U_{12}^{}/2E+V_ec_{13}^2`$, and $`M_2=diag(0,\mathrm{\Delta }m_{21}^2)`$. We neglect off-diagonal terms in the evolution equation, Eq. (2).

The evolution matrix S in the propagation basis $`(\stackrel{~}{\nu }_e,\stackrel{~}{\nu }_\mu ,\stackrel{~}{\nu }_\tau )`$ has the following form:

$$\stackrel{~}{S}\left(\begin{array}{ccc}\stackrel{~}{A}_{ee}& \stackrel{~}{A}_{e\mu }& 0\\ \stackrel{~}{A}_{\mu e}& \stackrel{~}{A}_{\mu \mu }& 0\\ 0& 0& \stackrel{~}{A}_{\tau \tau }\end{array}\right),$$
(4)
where $`A_{\tau \tau }\mathrm{exp}(i\mathrm{\Delta }m_{31}^2L/2E),`$ and $`L`$ is the total distance traveled by the neutrinos.

### 3.2 Flavor transitions and interference

Let us find the probabilities of ($`\nu _\mu \nu _e`$) oscillations, $`P_{\mu e}`$, and ($`\nu _e\nu _e`$) oscillations, $`P_{ee}`$, relevant for our problem. The $`S`$matrix in the flavor basis equals: $`S=\stackrel{~}{U}\stackrel{~}{S}\stackrel{~}{U}^{}`$, and we find

$`P_{\mu e}=\left|s_{13}c_{13}s_{23}\stackrel{~}{A}_{ee}+c_{13}c_{23}\stackrel{~}{A}_{\mu e}\right|^2+s_{13}^2c_{13}^2s_{23}^2,`$
and $`P_{ee}=c_{13}^4(1\stackrel{~}{P}_{\mu e})+s_{13}^4`$. For sub-GeV sample oscillations driven by $`\mathrm{\Delta }m_{31}^2`$ are averaged out, so that there is no interference effect due to state $`\stackrel{~}{\nu }_\tau `$. At the same time, according to (3.2) the amplitudes $`\stackrel{~}{A}_{ee}`$ and $`\stackrel{~}{A}_{\mu e}`$ interfere. It this interference which produces effect we are interested in this paper. Notice that amplitudes $`\stackrel{~}{A}_{ee}`$ and $`\stackrel{~}{A}_{\mu e}`$ are both due to solar oscillation parameters. However their interference appears due to presence of the third neutrino (non-zero $`s_{13}`$). In what follows we will call the interference of the amplitudes (with solar oscillation parameters) due to non-zero $`U_{e3}s_{13}`$ as induced interference.

Combining $`P_{\mu e}`$ and $`P_{ee}`$, the excess of the $`\nu _e`$flux equals:

$`{\displaystyle \frac{F_e}{F_e^0}}1=(rc_{23}^21)\stackrel{~}{P}_{\mu e}rs_{13}c_{13}^2\mathrm{sin}2\theta _{23}Q`$
$`s_{13}^2\left[2W_{23}+\stackrel{~}{P}_{\mu e}(r2)\right]+s_{13}^4W_{23}(2\stackrel{~}{P}_{\mu e})`$
and $`QRe(\stackrel{~}{A}_{ee}^{}\stackrel{~}{A}_{\mu e})`$ and $`W_{23}(1rs_{23}^2)`$. The first term on the left hand side (zero order in $`s_{13}^2`$) corresponds to the contribution we have discussed in . The second term is the effect of the induced interference. Let us stress its properties: 1). The interference term depends on $`s_{13}`$ linearly. So its effect may not be strongly suppressed even for small $`s_{13}`$. The interference depends on the sign of $`s_{13}`$, also does not have screening factor, and its smallness is mainly due to smallness of $`s_{13}`$. 3). Beside this term has opposite signs for neutrinos and anti-neutrinos.

We have calculated dependences of the excess of the $`e`$like events on the zenith angle of electron, $`\mathrm{\Theta }_e`$. The procedure was described before in Ref. . In Fig. 2 we show the zenith angle dependences of the excess of the e-like events for different values of oscillation parameter.

Concluding, we show that if the LMA solution is the correct one and for $`\theta _{23}=45^{}`$ ( in this case the effects due the oscillations driven by $`\mathrm{\Delta }m_{21}^2`$ only are suppressed) we can have a direct way to determine $`U_{e3}`$, from the electron neutrino zenith distribution as is shown in Figure 2.


text:
# Energy bandstructures of MgB2 and the related compounds

## 1 Introduction

From the discovery of the superconductivity in MgB2, a lot of theoretical and experimental researches have been performed in order to reveal the mechanism of the really high superconducting transition temperature ($`T_{\mathrm{sc}}`$) in metals and find a route to a new superconductor with much higher $`T_{\mathrm{sc}}`$. It has been indicated that MgB2 is a phonon-mediated BCS s-wave superconductor by theoretical calculations, the NMR study and the high-resolution photoemission study, then detailed phonon dispersion is calculated. On the other hand, multiple superconducting gaps are also indicated by the specific heat measurement, the tunneling spectroscopy study and the high-resolution photoemission study. In the electron doped system Mg1-xAlxB2, the structural instabilityhas been observed around $`x`$ 10%, resulting in the loss of superconductivity.

In this paper, the Fermi surface property in MgB2 is theoretically obtained for the further microscopic measurements of single crystals. The bandstructures for iso-electronic compounds LiBC and MgB2C2 are calculated for the comparison. Total energy calculations for Mg1-xAlxB2 are performed to clarify the structural instability.

## 2 Method of Calculation

Band structure calculations are carried out by using a full potential APW (FLAPW) method with the local density approximation (LDA) for the exchange correlation potential. For the LDA, the formula proposed by Gunnarsson and Lundqvist is used. For the band structure calculation, we used the program codes; TSPACE and KANSAI-99.

We used in the calculations the experimentally observed lattice constants: a= 3.083Å and c=3.527Å for MgB2. Muffin-tin (MT) radii are set as 0.4979 a and 0.2742 a for Mg and B, respectively. Core electrons (Ne-core for Mg, He-core for B) are calculated inside the MT spheres in each self-consistent step. The LAPW basis functions are truncated at $`|𝐤+𝐆_i|3.33\times 2\pi /𝐚`$, corresponding to 225 LAPW functions at the $`\mathrm{\Gamma }`$ point. The sampling points are uniformly distributed in the irreducible 1/24 th of the Brillouin zone, 333 k-points (divided by 12, 12, 12) for both the potential convergence and the final band structure. The similar parameters are used in calculations for LiBC and Mg1-xAlxB2.

MgB2C2 crystallizes in a base centered orthorhombic structure containing 8 formula units in the primitive unit cell. In this case, calculations are performed by using about 1,600 LAPW basis functions at the uniformly distributed 39 k-points for the potential convergence and 369 k-points (divided by 8, 8, 8) for the final band structure in the irreducible 1/8 th of the Brillouin zone.

## 3 Bandstructure and Fermi surface of MgB2

The calculated energy bandstructure and the density of states are shown in Fig. 1 and Fig. 2, respectively. They are very similar to the previous calculations, though Fig. 1 shows bandstructures along all the symmetry axes; the $`\mathrm{\Gamma }`$-K-M-$`\mathrm{\Gamma }`$ lines is parallel to the A-H-L-A lines on the basal plane and the $`\mathrm{\Lambda }`$, U and P axes are along the c-axis (see Fig. 3 (a)).

Let us see the dispersion along the P axis; the dispersive band crossing the Fermi level. The band is originated from B-$`p_z`$($`\pi `$) band and doubly degenerated. These B-$`p_z`$($`\pi `$) states are non-bonding along the P axis within the honeycomb lattice, and anti-bonding (bonding) at the K (H) point for the inter-layer (along the c-axis) mixing. While B-$`p_{xy}`$($`\sigma `$) states are split to bonding and anti-bonding sates within the honeycomb lattice, then the top of the bonding states are located just above the Fermi level along the $`\mathrm{\Lambda }`$ axis. While the bottom of the anti-bonding bands are located around 1.1 Ry. The energy and its components foe the states near the Fermi level are listed in Table 1. Eventually, there are four Fermi surfaces; two cylindrical hole Fermi surfaces along the $`\mathrm{\Lambda }`$ axis originating in the $`\sigma `$ band, and from the $`\pi `$ bands one hole around the T’ axis one electron along S’ axis, as shown in Fig. 3. The total number of holes and electrons are compensated.

The unique characteristics in the bandstructure of MgB2 is the existence of the cylindrical $`\sigma `$-holes (Fig. 3 (a) and (b)). The situation has been discussed as the relative shift of the $`\sigma `$ and $`\pi `$ bands. In fact, the totally non-bonding $`\sigma `$ states (about 0.91 Ry) is much higher in energy than the non-bonding $`\pi `$ state (0.65 Ry; the mid-point on the P axis). It should be ascribed to the ionic Mg2+ strongly affecting the crystalline field for the B-$`p`$ states. Moreover, here it is emphasized that the 3D network in the $`\pi `$ band is suggested by the dispersive band along the P axis.

As was already discussed, graphite (primitive graphite C2) shows a 2D bandstructure. In graphite, the inter-layer distance ($`d_{layer}=`$3.348Å) is smaller than in MgB2 ($`d_{layer}=`$3.528Å), but the planar bond length ($`d_{bond}=`$1.418 Å) is much smaller than in MgB2 ($`d_{bond}=`$1.780 Å). The planar bond length determines the total band width in each case, because of the much smaller distance and more coordination in plane, though it depends on how extended the wave functions are. As a result, the valence band width is narrower in MgB2 than graphite. However, what should be noticed here is the ratio of $`d_{layer}`$ and $`d_{bond}`$, which is 2.36 in graphite and 1.98 in MgB2. The much smaller ratio in MgB2 causes more 3D feature, resulting in the dispersive $`\pi `$ band along the P axis, which obtains both the hole and electron Fermi surfaces. The $`\sigma `$ bands keep the 2D character, because the $`p_{xy}`$ orbitals are extended in a layer. Eventually, the strong 3D feature in low-lying $`\pi `$ bands brings about 2D Fermi surface of the $`\sigma `$ band in MgB2, as shown in Fig. 3. The 3D feature and the low-lying in $`\pi `$ bands are ascribed to the existence of Mg2+ ion.

The calculated specific heat coefficient $`\gamma _{band}`$ is 1.73 mJ/mol K2, while the experimental value $`\gamma _{exp}`$ is 2.6 mJ/mol K2. Therefore, the averaged electron-phonon mass enhancement factor $`\lambda `$ is derived as 0.50 from $`\gamma _{exp}=(1+\lambda )\gamma _{band}`$.

The calculated angular dependence of the de Haas-van Alphen (dHvA) frequencies are shown in Fig. 4. The dHvA frequencies and the cyclotron masses are summarized in Table 2. Several theoretical studies have pointed out the importance of electron-phonon coupling in the $`\sigma `$ band for the superconductivity, then the superconducting gap in MgB2 can be expected to be larger on the $`\sigma `$ Fermi surfaces than on the $`\pi `$ Fermi surfaces. It is also expected to be larger mass enhancement in the cyclotron masses in the $`\sigma `$ Fermi surfaces, though they have not yet observed so far, unfortunately. Recently MgB2 single crystals are grown by a several groups, then the dHvA oscillations is expected to be detected soon. The experimentally observed cyclotron masses will be the direct evidence for the two superconducting gaps in MgB2.

## 4 Bandstructures of LiBC and MgB2C2

It is worth to see bandstructures of the related compounds. LiBC and MgB2C2 are iso-electric compounds with MgB2, though the crystal structure is slightly different. LiBC can be obtained, if Mg and one of B in MgB2 are replaced by Li and C. The stacking along the c-direction has not been experimentally determined, so we assume B-B-B…(C-C-C…) stacking. MgB2C2 crystallizes in a orthorhombic structure, as shown in Fig. 5, in which there are a distorted B-C honeycomb lattice layer and a Mg layer. In LiBC, $`d_{layer}=`$3.53 Å, $`d_{bond}=`$1.59 Å and $`d_{layer}/d_{bond}=`$2.22. While, in MgB2C2, $`d_{layer}=`$3.74 Å, $`d_{bond}=`$1.58 Å in average and $`d_{layer}/d_{bond}=`$2.37. Therefore, both compounds are expected to be less 3D than MgB2.

Figure 6 shows the bandstructure for LiBC. The space group of LiBC is P6̄$`m`$2, then the degeneracy along the P-axis is lifted. Li1+ ion has less attractive for the $`\pi `$ bands, then the bonding $`\sigma `$ bands get under the Fermi level. Eventually, there are hole and electron Fermi surfaces both from the $`\pi `$ bands, as shown in Fig. 7.

More interesting bandstructure is found in MgB2C2. As shown in Fig. 8, MgB2C2 becomes a semiconductor in consistent with the measured conductivity. The opening the gap may be due to the less anisotropic hybridyzation in the lower symmetry. However, the dispersion of top of the valence bands are very flat along the $`\mathrm{\Lambda }`$ axis, which is similar to the $`\sigma `$ bands in MgB2. Cylindrical Fermi surface, which is relevant in the superconductivity in MgB2, can be obtained with hole doping.

## 5 Lattice Constants and Structural Instability in Mg1-xAlxB2

With substituting Al for Mg in MgB2, the superconducting transition temperature gradually decreases and structural instability appears with Al content from 10 % to 20%. After the instability, the superconducting transition temperature rapidly decreases. Partial Al substitution for Mg means electron doping to the MgB2 bandstructure, leading to vanishing $`\sigma `$ hole Fermi surfaces. One can expect the structural phase instability is related to the disappearance of the $`\sigma `$ hole.

To reveal the structural instability, the total energy calculations with a several concentration and the lattice constant are performed. In the calculations, the lattice constant a is assumed to be constant, as its small Al concentration dependence (a=3.083Å in MgB2 to 3.006 Å in AlB2) and the virtual lattice approximation is adopted, in which the atom with the atomic number of $`12.x`$ is used for the case of Al concentration $`x`$.

The total energies are obtained at the points in Fig. 9, then the optimized lattice constant c for each concentration $`x`$ are determined by using parabolic fitting. The total energy as a function of c is well reproduced as a parabolic function, so that the present calculations cannot predict the observed structural instability for any Al concentration. There is another structural instability reported around $`x=`$70%, where the superconductivity disappears. Nevertheless the optimized c decreases monotonically with increasing $`x`$, as in Fig. 9. Concerning that Aluminum ordering is observed for $`x=0.5`$, the virtual lattice approximation may not be applied to the system Mg1-xAlxB2.

For AlB2, with using the measured a=3.006Å, the optimized c is successfully obtained as 3.245Å for the experimental value 3.251Å. With the value a=3.083Å for MgB2, which is used in Fig. 9, c=3.241Å for AlB2. In each case, the error is less than about 0.3% for AlB2. However, for MgB2, the optimized c is 3.457Å which is 2% smaller than the measured value 3.527Å. The lattice constants decrease in low temperatures, but c is reported as 3.51504(4)Å at 15K, which is only 0.3 % smaller than in room temperature. This 2% error should make some sense and corresponds to the jumps of two discontinuities (structural instabilities) around 20% and 70% Al concentration. The larger inter-layer distance (c) could decrease the number of the $`\sigma `$ holes. Therefore there may exist some mechanism beyond the LDA that favours the small number of the $`\sigma `$ holes and is related to the superconductivity. This should be solved in the future.

## Acknowledgments

The author should thank to J. Akimitsu, K. Ishida, Y. Kitaoka, T. Takahashi, H. Kohno, K. Miyake, M. Imada, Y. Kuramoto, M. Shirai for fruitful, enlightening and stimulating discussions. This work was partly supported by a Grant-in-Aid for Science Research from the Ministry of Education, Culture, Sports, Science and Technology, Japan.