text:
# Specific heat and thermal conductivity in the vortex state of the two-gap superconductor MgB2

## I INTRODUCTION

Although superconductivity in MgB2  has been established as conventional phonon-mediated s-wave superconductivity, it is often difficult to analyze the results of measurements of various physical quantities because there exist two gaps of different magnitude associated with different bands. Their ratio is estimated as $`r=\mathrm{\Delta }_0^S/\mathrm{\Delta }_0^L0.30.4`$ where the larger gap $`\mathrm{\Delta }_0^L`$ is associated with the two-dimensional $`\sigma `$-bands and the smaller gap $`\mathrm{\Delta }_0^S`$ with the three-dimensional $`\pi `$-bands. Evidence for two gaps is provided by the rapid rise of the specific heat coefficient $`\gamma _s(H)`$ at very low fields. Recent measurements of the in-plane thermal conductivity $`\kappa _s(H)`$ for fields both parallel and perpendicular to the c-axis provide more evidence for the existence of two gaps because they show a very unusual field dependence at low temperatures.  For increasing field, $`\kappa _s(H)`$ drops rapidly and reaches a minimum at a relatively low field, then exhibits an S-shape behavior where a steep rise in the low-field region is followed by a downward curvature up to about $`H_{c2}/2`$ and finally an upward curvature near $`H_{c2}`$ for $`𝐇𝐜`$. The initial drop is attributed to the decrease of the phononic contribution $`\kappa _{ph}(H)`$ which is caused by the rapid enhancement of scattering of phonons by the quasiparticles in the vortex cores associated with the $`\pi `$-band. The steep increase in the low field region is attributed to the rapid rise of the electronic contribution $`\kappa _e(H)`$ which arises from the the rapid release of mobile quasiparticles in the vortex lattice asociated with the $`\pi `$-band. The final upward curvature of $`\kappa (H)`$ near $`H_{c2}`$ for $`𝐇𝐜`$ is presumed to be due to the mobile carriers of the vortex lattice associated with the $`\sigma `$-band.

The field dependence of $`\gamma _s(H)`$ and $`\kappa _{es}(H)`$ can be explained qualitatively in terms of a two-band model with two energy gaps of different magnitude in the two bands where the band with strong pairing (the L, or $`\sigma `$-band) is responsible for the superconductivity and the superconductivity in the second band (the S, or $`\pi `$-band) is induced by Cooper pair tunneling.  The vortices in the S-band have large vortex cores of radius $`\xi _0^S`$ and already start to overlap in weak fields. This leads to a strong suppression of the small gap at a ”virtual” upper critical field $`H_{c2}^S34`$ kOe which is much smaller than $`H_{c2}^{(c)}35`$ kOe. The large radius of the vortex core and the field-induced suppression of the smaller gap is consistent with recent scanning tunneling spectroscopy (STS) measurements.

Theories of the electronic thermal conductivity $`\kappa _e`$  and the phononic thermal conductivity $`\kappa _{ph}`$  in the mixed state have been developed on the basis of the Brandt-Pesch-Tewordt (BPT)  and the equivalent Pesch (P)  approximation schemes. Linear response equations for the vortex state  yield another expression for $`\kappa _e`$ which has been applied to d-wave pairing superconductivity in the cuprates.  In the BPT-approximation the normal and anomalous Green’s functions $`G`$ and $`F`$ are derived from Gorkov’s integral equations with kernel $`\mathrm{\Delta }(r_1)G_\omega ^0(r_1r_2)\mathrm{\Delta }^{}(r_2)\mathrm{exp}(i_1^22eA𝑑s)`$ where $`\mathrm{\Delta }(r)`$ is Abrikosov’s vortex lattice order parameter and $`G_\omega ^0`$ is the quasihole propagator. Thus these Green’s functions take into account the effects of supercurrent flow and Andreev scattering in the vortex cores as well as quasiparticle transfer between the vortices. The expression for $`\kappa _e`$ , which has been derived from the Kubo formula in analogy to the procedure in Ref. , consists essentially of the $`\omega `$-integral with the factor $`(\omega /T)^2\text{sech}^2(\omega /2T)`$ determining the temperature dependence, the ”destructive” coherence factor arising from the combination $`(GG^{}FF^{})`$, and the total scattering rate in the denominator consisting of the sum of the quasiparticle scattering rates $`\gamma `$ and $`\gamma _A(\omega ,H)`$ due to impurity and Andreev scattering. In the present paper we shall also calculate the scattering rate $`1/\tau _{es}`$ due to thermal scattering in the vortex state, but this becomes comparable to $`\gamma `$ only at higher temperatures. The expression for $`\kappa _{ph}`$ contains, in addition to the sum of scattering rates of phonons by different defects, the scattering rate $`1/\tau _{ps}`$ for phonons scattered by quasiparticles in the vortex state. The latter expression was obtained by replacing in the BCS-expression  the density of states (DOS) and ”coherence” functions by the expressions resulting from the P-approximation for the vortex state.

The results for $`\gamma _s(H)/\gamma _n`$, $`\kappa _{es}(H)/\kappa _{en}`$, and $`\tau _{pn}/\tau _{ps}`$, which were obtained in Ref.  for a single cylindrical or spherical band and an isotropic s-wave pairing gap, turned out to depend sensitively on the reduced impurity scattering rate $`\delta =\mathrm{\Gamma }/\mathrm{\Delta }_0=(\pi /2)\xi _0/\mathrm{}`$ where $`\xi _0`$ is the zero-field coherence length and $`\mathrm{}`$ is the mean free path. If one takes for the $`\pi `$-band in MgB$`_2`$ $`\xi _\pi 500\text{Å}`$ from the STS measurements  and a mean free path $`\mathrm{}_\pi 500800\text{Å}`$, one arrives at a reduced sacttering rate $`\delta _\pi 1`$ which is much larger than the largest value of $`\delta =0.5`$ taken in Ref. . For such a large impurity scattering rate it is necessary to renormalize the impurity scattering self energy and to take into account the corresponding vertex corrections. The vertex corrections yield a renormalization of the gap function  which, near $`H_{c2}`$, tends to the renormalization function which has been derived in the calculation of $`H_{c2}`$.  It turns out that this gap renormalization substantially reduces the dependence of the functions $`\gamma _s(H)/\gamma _n`$, $`\kappa _{es}(H)/\kappa _{en}`$, and $`\tau _{pn}/\tau _{ps}(H)`$ on impurity scattering.

Recently it has been shown that numerical solutions of the Eilenberger equations yield a spatial average of the density of states (DOS) which is approximated very well by the DOS obtained with the BPT- and P-approximations over the whole field range between $`H_{c2}`$ and $`H_{c1}`$.  This result motivates us to calculate the field dependence of the central parameter of the theory, $`\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }\mathrm{\Lambda }/v`$, from the gap equation in the P-approximation  by including the impurity scattering vertex corrections in analogy to those in the gap equation of the BPT-approximation.  Here $`\mathrm{\Delta }^2`$ is the spatial average of $`|\mathrm{\Delta }(r)|^2`$ where $`\mathrm{\Delta }(r)`$ is Abrikosov’s vortex lattice order parameter, $`\mathrm{\Lambda }=(2eH)^{1/2}=(\pi /2)^{1/2}a`$ where $`a`$ is the lattice constant of the vortex lattice, and $`v`$ is the Fermi velocity perpendicular to the uniform field H. It turns out that one can fit the exact results of these calculations very well with the relation between $`\stackrel{~}{\mathrm{\Delta }}^2`$ and $`H/H_{c2}`$ given by the Ginzburg-Landau relationship if one uses an appropriate value of the Abrikosov parameter $`\beta _A=<|\mathrm{\Delta }|^4>/(<|\mathrm{\Delta }|^2>)^2`$.

We have also investigated the temperature dependence of the phonon relaxation time ratio $`\tau _{pn}/\tau _{ps}(H)`$ and find that this ratio increases rapidly to one with increasing $`H`$ in proportion to $`[\gamma _s(H)/\gamma _n]^2`$ for temperatures $`T<(1/2)\mathrm{\Delta }(T,H)`$. This yields a rapid decrease of $`\kappa _{ph}`$ for small fields and relatively low temperatures if one takes for $`\mathrm{\Delta }`$ the small gap $`\mathrm{\Delta }^S(T,H)`$ of the $`\pi `$-band. In this way one can explain the observed drop of $`\kappa `$ for small fields and temperatures near $`T_c/6`$ which is attributed to the drop of $`\kappa _{ph}`$.

Since the interband impurity scattering between the $`\sigma `$\- and $`\pi `$-band is small,  one can add, in a crude approximation, our results obtained for a single isotropic gap in a cylindrical or spherical Fermi surface (FS) band by weighting them with the corresponding density of states. If one takes into account that the ”virtual” upper critical field for the three-dimensional $`\pi `$-band is much smaller than the upper critical field of the two-dimensional $`\sigma `$-band, one can explain qualitatively the measured field dependence of the specific heat and thermal conductivity in MgB$`_2`$.

## II Theory of the specific heat and electronic and phononic thermal conductivity in Abrikosov’s vortex lattice

The basis of our theory is the approximation of Brandt, Pesch, and Tewordt (BPT)  for the normal and anomalous Green’s functions $`G`$ and $`F`$ in the mixed state. A simplified version of this theory has been derived by Pesch (P)  from the quasiclassical Eilenberger equations. In this version the $`\epsilon `$-integrals of the spatial averages of the spectral functions, $`\text{Im}G/\pi `$ and $`\text{Im}F/\pi `$, denoted by $`A`$ and $`B`$, take on the following form:

$$A(\mathrm{\Omega };\stackrel{~}{\mathrm{\Delta }},\theta )N(\omega ,\theta )/N_0=\text{Re}\left[1+\frac{8\stackrel{~}{\mathrm{\Delta }}_r^2}{\mathrm{sin}^2\theta }\left[1+i\sqrt{\pi }zw(z)\right]\right]^{1/2},$$
(1)
$$B(\mathrm{\Omega };\stackrel{~}{\mathrm{\Delta }},\theta )=\text{Re}\left[\frac{i\sqrt{\pi }\mathrm{\hspace{0.17em}2}(\stackrel{~}{\mathrm{\Delta }}_r/\mathrm{sin}\theta )w(z)}{\left\{1+(8\stackrel{~}{\mathrm{\Delta }}_r^2/\mathrm{sin}^2\theta )[1+i\sqrt{\pi }zw(z)]\right\}^{1/2}}\right],$$
(2)
$$z=2[\mathrm{\Omega }\stackrel{~}{\mathrm{\Delta }}+i(\mathrm{\Lambda }/v)\gamma ]/\mathrm{sin}\theta ;\theta =\mathrm{}(𝐩,𝐇),$$
(3)
$$\mathrm{\Lambda }=(2eH)^{1/2},\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }\mathrm{\Lambda }/v,\mathrm{\Omega }=\omega /\mathrm{\Delta },$$
(4)
$$\stackrel{~}{\mathrm{\Delta }}^2=(H_{c2}H)/6\beta _AH,\mathrm{\Delta }^2=\mathrm{\Delta }_0^2(T)\left[1(H/H_{c2})\right],$$
(5)
$$w(z)=\mathrm{exp}(z^2)\mathrm{erfc}(iz).$$
(6)
$$\stackrel{~}{\mathrm{\Delta }}_r=\stackrel{~}{\mathrm{\Delta }}/D,D(\mathrm{\Omega })=1_0^{\pi /2}𝑑\theta \mathrm{\hspace{0.17em}2}(\mathrm{\Lambda }\gamma /v)\sqrt{\pi }w(z).$$
(7)
Here $`\mathrm{\Delta }_0(T)`$ is the s-wave gap in zero field, $`H`$ is the spatial average of the magnetic field, $`\mathrm{\Delta }^2=<|\mathrm{\Delta }|^2>`$, where $`\mathrm{\Delta }(r)`$ is Abrikosov’s vortex lattice order parameter, and $`\beta _A`$ is the Abrikosov parameter. For the total impurity scattering rate $`\gamma `$ we employ the Born approximation, $`\gamma =\mathrm{\Gamma }\overline{A}`$ where $`\mathrm{\Gamma }=1/2\tau _n`$ and $`\overline{A}`$ is the average of $`A`$ over the Fermi surface which is calculated self-consistently. For a spherical Fermi surface (FS), $`v\mathrm{sin}\theta `$ is the component of the Fermi velocity $`𝐯(𝐩)`$ perpendicular to the field $`𝐇`$. For a cylindrical FS parallel to $`𝐜`$ and field $`𝐇||𝐜`$, one has to take $`\theta =\pi /2`$ and $`v`$ equal to the Fermi velocity in the ab-plane. For a general FS, $`v\mathrm{sin}\theta `$ is replaced by the component $`v_{}(𝐩)`$ perpendicular to the field. The gap renormalization function $`D`$ in Eq.(7) was first calculated by U. Brandt  from the ladder summation of impurity interaction lines bridging the vertices at $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }^{}`$ occuring in the self energy part $`G_\omega ^0\mathrm{\Delta }G_\omega ^0\mathrm{\Delta }^{}G_\omega ^0`$ in Gorkov’s integral equation for $`G`$. In the P-approxiamtion scheme $`D`$ takes the form given in Eq.(7) which, near $`H_{c2}`$, tends to the function introduced in Ref. .

The field dependence of the DOS, A, and the ”coherence factor” function, B, are governed by the quantity $`\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }\mathrm{\Lambda }/v=(2\pi )^{1/2}a/\xi `$ where $`a=(2/\pi )^{1/2}\mathrm{\Lambda }`$ is the lattice constant of the Abrikosov vortex lattice and $`\xi =v/\mathrm{\Delta }\pi `$ is the effective coherence length in the plane perpendicular to H . For $`HH_{c2}`$ and $`H0`$ one obtains the limiting values $`\stackrel{~}{\mathrm{\Delta }}0`$ and $`\stackrel{~}{\mathrm{\Delta }}\mathrm{}`$, respectively. The relationships between $`\stackrel{~}{\mathrm{\Delta }}`$ and $`H/H_{c2}`$ and between $`\mathrm{\Delta }`$ and $`H/H_{c2}`$ given in Eq.(5) have been derived from the GL-theory which is valid near $`H_{c2}`$. To get the general relationship at lower fields we make use of the gap equation in the P-approximation. At $`T=0`$ this equation can be rewritten in the following form:

$$_0^{\mathrm{}}𝑑\mathrm{\Omega }[B(\mathrm{\Omega })C(\mathrm{\Omega })]=0.$$
(8)
Here $`B(\mathrm{\Omega })`$ is the spectral function of $`F`$ given in Eq. (2) and $`C(\mathrm{\Omega })`$ is obtained by taking the limit of $`B`$ for $`HH_{c2}`$ where $`\mathrm{\Delta }0`$:

$$C(\mathrm{\Omega })=\text{Re}\left[i\sqrt{\pi }2(H/H_{c2})^{1/2}(\stackrel{~}{\mathrm{\Delta }}_r/\mathrm{sin}\theta )w\left[(H/H_{c2})^{1/2}z\right]\right]|_{\mathrm{\Lambda }(H_{c2})}.$$
(9)
The gap equation Eq. (8) in the P-approximation corresponds to the gap equation in the BPT-approximation  with the function $`B`$ replaced by the spectral function of $`F`$ derived for the theory of NMR.  Thus Eq. (8) is an implicit equation for $`H/H_{c2}`$ as a function of $`\stackrel{~}{\mathrm{\Delta }}`$ for given impurity scattering rate $`\delta =\mathrm{\Gamma }/\mathrm{\Delta }_0`$. This calculation is rather complicated because the impurity scattering self energy $`\gamma =\mathrm{\Gamma }\overline{A}(\mathrm{\Omega })`$ occuring in the expression for $`B(\mathrm{\Omega })`$ in Eq. (2) has to be calculated self-consistently for each value of $`\mathrm{\Omega }`$ from the expression for $`A(\mathrm{\Omega })`$ (see Eq. (1)). For $`𝐇||𝐜`$ and a cylindrical FS one has to set $`\theta =\pi /2`$ in Eq.(8), and for a spherical FS Eq.(8) is averaged over the polar angle $`\theta `$. In Fig. 1 we show our results for $`H/H_{c2}`$ versus $`\stackrel{~}{\mathrm{\Delta }}`$ obtained from Eq.(8) for $`\theta =\pi /2`$ and for $`\delta `$ = 0.1, 0.5, and 1.0. These results can be fitted very well by the relationship $`H/H_{c2}=(1+6\beta _A\stackrel{~}{\mathrm{\Delta }}^2)^1`$ \[see Eq.(5)\] over the whole field range if one takes $`\beta _A`$ = 1.1, 1.4, and 1.7 for $`\delta `$ = 0.1, 0.5, and 1.0. Similar results are obtained for a spherical FS if one averages Eq.(8) over the polar angle $`\theta `$ (see the dashed curves in Fig. 1 for $`\delta `$ = 0.1, 0.5, and 0.8 with $`\beta _A`$ = 2.3, 2.8, and 3.2).

In Fig. 2 we have plotted our numerical results for the normalized DOS at zero energy, $`A(\mathrm{\Omega }=0)=\gamma _s(H)/\gamma _n`$, versus $`H/H_{c2}`$ where $`\gamma `$ is the specific heat coefficient. Here we have approximated $`\stackrel{~}{\mathrm{\Delta }}`$ occuring in Eq. (1) by the analytical expression in Eq. (5) with $`\beta _A`$ = 1.1, 1.4, and 1.7 for bare impurity scattering rates $`\delta =\mathrm{\Gamma }/\mathrm{\Delta }_0=`$0.1, 0.5, and 1.0 (solid curves, from top to bottom). Here we have employed the Born approximation for the total impurity scattering self energy, $`\gamma =\mathrm{\Gamma }\overline{A(0)}`$, where $`\overline{A(0)}`$ is calculated self-consistently. This means that $`\delta `$ is replaced by $`\delta \overline{A(0)}`$. For $`\delta =0.1`$ and $`\theta =\pi /2`$ the function $`A(0)`$ is close to the clean limit result shown in Ref.  for the s-wave pairing state. The latter function has been shown to be very close to the numerical solution of the Eilenberger equations over the whole field range from $`H_{c2}`$ to $`H=0`$ corresponding to $`H_{c1}`$. One can see that all curves in Fig. 2 tend to zero in the limit $`H/H_{c2}0`$, which is due to the gap renormalization. For a spherical FS one has to average Eq.(1) over the polar angle $`\theta `$. For the plot of $`A`$ versus $`H/H_{c2}`$ one needs Eq.(5) where $`\beta _A`$ is calculated from the $`\theta `$-angle average of Eq.(8). The results are shown by the dashed curves in Fig. 2 for $`\delta `$ = 0.1, 0.5, and 0.8 with the $`\beta _A`$-values used in Fig. 1. It is interesting that the curve for $`\delta `$ = 0.8 is almost linear over the whole field region.

We turn now to the theory of the electronic thermal conductivity $`\kappa _e`$ in the vortex state which has been developed in Ref.  by replacing the zero-field Green’s functions $`G`$ and $`F`$ in the theory of Ref.  by the Green’s functions of the BPT-approximation. This theory has been applied to superconducting states with nodes in the gap such as Sr2RuO$`_4`$, the cuprates, and UPt$`_3`$. To save space we omit here the expression for $`\kappa _e`$. The main features of this expression are the following. Most important is the term Im$`\epsilon _0`$ in the denominator of the $`\omega `$-integral with the well-known factor $`\omega ^2\text{sech}^2(\omega /2T)`$. This term Im$`\epsilon _0`$ corresponds to the scattering rate of quasiparticles due to impurity and Andreev scattering. The other term corresponds to the coherence factor of BCS theory. In the zero field limit both terms in the expression of Ref.  tend correctly to those occuring in the expression of Ref. . The physical meaning of Im$`\epsilon _0`$ is the following. The equation for the position, $`\epsilon _0`$, of the pole of the Green’s function $`G`$ in the BPT-approximation  yields Im$`\epsilon _0=\gamma +\gamma _A`$ where $`\gamma `$ is the total impurity scattering rate and $`\gamma _A`$ is the imaginary part of the quasiparticle self energy $`\mathrm{\Sigma }_\omega `$ at $`\epsilon _0`$. From the kernel of Gorkov’s integral equation for $`G`$ one sees that, in the spatial representation, $`\mathrm{\Sigma }_\omega (r_1,r_2)=V(r_1,r_2)G_\omega ^0(r_1r_2)`$, where $`V(r_1,r_2)=\mathrm{\Delta }(r_1)\mathrm{\Delta }^{}(r_2)\mathrm{exp}(ı_1^22eA𝑑s)`$ and $`\mathrm{\Delta }(r)`$ is Abrikosov’s vortex lattice order parameter. From this expression one recognizes that $`\gamma _A=\text{Im}\mathrm{\Sigma }_\omega `$ is the scattering rate for converting a quasiparticle into a quasihole by Andreev reflection at $`\mathrm{\Delta }^{}(r_2)`$ and then back into a quasiparticle at $`\mathrm{\Delta }(r_1)`$. The Fourier transform of $`\mathrm{\Sigma }_\omega `$ with respect to the difference coordinate $`(r_1r_2)`$ is equal to

$$\mathrm{\Sigma }_\omega (p)=d^3p^{}V(p^{})G_\omega ^0(pp^{}),$$
(10)
with

$$V(p^{})=8\pi ^2\mathrm{\Delta }^2\mathrm{\Lambda }^2\delta (p_z^{})\mathrm{exp}\left[\mathrm{\Lambda }^2\left(p_x^{\mathrm{\hspace{0.17em}2}}+p_y^{\mathrm{\hspace{0.17em}2}}\right)\right](\widehat{𝐳}𝐇).$$
(11)
For a sperical Fermi surface the result is

$$\mathrm{\Sigma }_\omega (p)=i\sqrt{\pi }\mathrm{\Delta }^2(\mathrm{\Lambda }/v\mathrm{sin}\theta )w[(\omega +i\gamma +\epsilon )\mathrm{\Lambda }/v\mathrm{sin}\theta ]$$
(12)
where $`\epsilon `$ is the normal state energy measured from the Fermi energy, $`\theta =\mathrm{}(𝐩,𝐇)`$, and $`\gamma `$ is the total impurity scattering self energy.

In the $`\omega 0`$ limit corresponding to $`T0`$ one obtains the following explicit expression for Im$`\epsilon _0`$:

$`\text{Im}\epsilon _0`$ $`=`$ $`\gamma +\gamma _A;\gamma _A=(v/\mathrm{\Lambda })\beta `$ (13)
$`\beta `$ $`=`$ $`\stackrel{~}{\mathrm{\Delta }}_r^2\left\{\left[(\gamma \mathrm{\Lambda }/v)^2+\stackrel{~}{\mathrm{\Delta }}_r^2+\eta \mathrm{sin}^2\theta \right]^{1/2}+(\gamma \mathrm{\Lambda }/v)\right\}^1;(\pi ^1\eta 1/2).`$ (14)
All the other quantities $`\stackrel{~}{\mathrm{\Delta }}`$ and $`\beta `$ occuring in the expression for $`\kappa _e`$ (see Eq.(13) of Ref. )) are also renormalized by the replacement $`\stackrel{~}{\mathrm{\Delta }}\stackrel{~}{\mathrm{\Delta }}_r`$ (see Eq.(7) for $`\mathrm{\Omega }=0`$). Note that $`\stackrel{~}{\mathrm{\Delta }}a/\xi `$ and $`\mathrm{\Gamma }\mathrm{\Lambda }/va/\mathrm{}`$ . In Fig. 3 we have plotted our results for $`\kappa _{es}/\kappa _{en}`$ versus $`H/H_{c2}`$ in the limit $`T0`$. The bare impurity scattering rates $`\delta =\mathrm{\Gamma }/\mathrm{\Delta }_0`$ are 0.1, 0.5, and 1.0 (from bottom to top). For a spherical FS one has to evaluate the average of $`\kappa _e`$ over the polar angle $`\theta `$ which yields nearly the same results as shown in Fig. 3 for the plots of $`\kappa _{es}/\kappa _{en}`$ vs $`H/H_{c2}`$ if one uses Eq.(5) with the $`\beta _A`$ calculated from the average of Eq.(8) (see the dashed curves in Fig. 3 for $`\delta `$ = 0.1, 0.5, and 0.8 with the $`\beta _A`$-values given in Fig. 1). All curves tend to zero in the limit $`H/H_{c2}0`$. One sees that, for low impurity scattering rates, $`\kappa _{es}/\kappa _{en}`$ exhibits a strong upward curvature near $`H_{c2}`$ which is caused by the rapid reduction of the Andreev scattering rate $`\gamma _A`$ with decreasing $`\stackrel{~}{\mathrm{\Delta }}`$ or increasing field as can be seen from Eq. (14). For increasing impurity scattering rate $`\gamma `$ this effect is diminished because $`\gamma _A`$ becomes smaller in comparison to $`\gamma `$ as can be seen from the expression in Eq. (14).

The main problem in analyzing the measured thermal conductivity $`\kappa `$ in the mixed state is to separate the electronic and phononic contributions $`\kappa _e`$ and $`\kappa _{ph}`$.  The general expression for $`\kappa _{ph}`$ (see Eq. (4) of Ref. ) contains, in addition to other phonon relaxation times due to defect scattering, the relaxation time $`\tau _p(\omega )`$ due to scattering of phonons of frequency $`\omega `$ by quasiparticles which changes drastically in the BCS state. The expression for $`1/\tau _{ps}`$ in the BCS state  has been extended to the vortex lattice state by replacing the BCS DOS $`E/\epsilon `$, where $`\epsilon =(E^2\mathrm{\Delta }^2)^{1/2}`$, and the function $`\mathrm{\Delta }/\epsilon `$ occuring in the BCS coherence factor, with the spectral functions $`A`$ and $`B`$ of Eqs. (1) and (2). In the zero field limit the latter expressions tend correctly to the BCS functions with $`\omega E`$. The resulting ratio of phonon relaxation times in the normal and mixed states is given by (see Eq. (3) of Ref. ):

$`\tau _{pn}/\tau _{ps}=g(\mathrm{\Omega }_0)`$ $`=`$ $`[1\mathrm{exp}(\mathrm{\Omega }_0/y)](2/\mathrm{\Omega }_0){\displaystyle _0^{\mathrm{}}}𝑑\mathrm{\Omega }f[(\mathrm{\Omega }\mathrm{\Omega }_0/2)/y]f[(\mathrm{\Omega }+\mathrm{\Omega }_0/2)/y]`$ (16)
$`\times \left[A(\mathrm{\Omega }\mathrm{\Omega }_0/2)A(\mathrm{\Omega }+\mathrm{\Omega }_0/2)B(\mathrm{\Omega }\mathrm{\Omega }_0/2)B(\mathrm{\Omega }+\mathrm{\Omega }_0/2)\right].`$
Here, $`\mathrm{\Omega }_0=\omega _0/\mathrm{\Delta }`$ is the phonon energy $`\omega _0`$ divided by $`\mathrm{\Delta }=\mathrm{\Delta }(T,H)`$, f is the Fermi function, and $`y=T/\mathrm{\Delta }(T,H)`$. In Ref.  we calculated $`g`$ in the limit $`T0`$. It was found that in the limit $`\stackrel{~}{\mathrm{\Delta }}\mathrm{}`$, or $`H/H_{c2}0`$, g tends correctly to the BCS step functon with a step at $`\mathrm{\Omega }_0=2`$, or $`\omega _0=2\mathrm{\Delta }`$. This step function is washed out as $`\stackrel{~}{\mathrm{\Delta }}0`$, or $`H/H_{c2}1`$, and it tends to the constant, 1. Since it is important to know the temperature dependence of $`\kappa _{ph}`$, and thus of $`\tau _{pn}/\tau _{ps}=g`$, for analyzing the measured $`\kappa `$ as a function of $`T`$ at fixed field,  we have now calculated the $`T`$-dependence of the expression in Eq.(16). For this purpose we make use of the fact that, according to the general expression for $`\kappa _{ph}`$ (see Eq.(4) of Ref. ), the maximum of the integrand in that x-integral occurs at about $`x_m=\omega _{0m}/T4`$. Therefore it is a good approximation to make the variable transformations $`\mathrm{\Omega }_0=x_my`$ and $`\mathrm{\Omega }=ty`$ in Eq.(16). This yields $`\tau _{pn}/\tau _{ps}=g(x_my)`$ as a function of the $`T`$-dependent quantity $`y`$ for fixed $`\stackrel{~}{\mathrm{\Delta }}`$, or field $`H/H_{c2}`$. In Fig. 4 we have plotted $`g(x_my)`$ versus $`y`$ for constant $`x_m=4`$ and constant $`\stackrel{~}{\mathrm{\Delta }}=`$ 0.2, 0.3, and 0.6, or $`H/H_{c2}=`$ 0.78, 0.62, and 0.29 (see Fig. 1) for $`\delta =0.1`$. In the limit $`H/H_{c2}0`$, one obtains a step function with a step at $`y=1/2`$, that is, $`\mathrm{\Omega }_0=4y=2`$ or $`\omega _0=2\mathrm{\Delta }`$, as it should be. One sees, that for intermediate values of the field, $`\tau _{pn}/\tau _{ps}=g`$ considered as a function of $`T`$ is smeared out and tends to one for high fields. In the limit $`y0`$, or $`\mathrm{\Omega }_00`$ Eq.(16) yields $`g(\mathrm{\Omega }_0=0)=[A(\mathrm{\Omega }_0=0,\stackrel{~}{\mathrm{\Delta }},\theta )]^2`$. If the scattering rate $`1/\tau _{ps}=(1/\tau _{pn})g`$ occuring in the denominator of the $`\omega `$-integral in the expression for $`\kappa _{ph}`$ (see Eq.(4) in Ref. ) is comparable to or larger than the other scattering rates due to phonon-defect scattering, we expect that $`\kappa _{ph}`$ rapidly decreases for increasing field for temperatures smaller or close to $`T(1/2)\mathrm{\Delta }(T,H)`$. In the measurements of $`\kappa `$ one observes a fast drop of $`\kappa `$ in small fields at temperatures around $`T_c/6`$  which can be attributed to the fast drop of $`\kappa _{ph}`$ due to scattering by quasiparticles in the $`\pi `$-band vortex lattice.

It is interesting to compare the temperature dependence of the phonon lifetime $`\tau _p`$ with the lifetime $`\tau _e`$ of a quasiparticle due to scattering by acoustic phonons. This scattering rate has been derived for the BCS state in Ref.  (see Eq.(3.9)). We extend this expression for the scattering rate $`1/\tau _{es}`$ in anology to our procedure for $`1/\tau _p`$ from the BCS to the vortex state by replacing the density of states by the DOS $`A(\mathrm{\Omega })`$ in Eq.(1) and the coherence factor terms by $`B(\mathrm{\Omega })`$ in Eq.(2). Then we obtain:

$`1/\tau _{es}(\omega )`$ $`=`$ $`C[f(\beta \omega )]^1A^1(\mathrm{\Omega }){\displaystyle _0^{\mathrm{}}}d\omega _0\omega _0^2\{A(\mathrm{\Omega })A(\mathrm{\Omega }\mathrm{\Omega }_0)B(\mathrm{\Omega })B(\mathrm{\Omega }\mathrm{\Omega }_0)`$ (19)
$`\times f[\beta (\omega \omega _0)][1+b(\beta \omega _0)]`$
$`+[A(\mathrm{\Omega })A(\mathrm{\Omega }+\mathrm{\Omega }_0)B(\mathrm{\Omega })B(\mathrm{\Omega }+\mathrm{\Omega }_0)]f[\beta (\omega +\omega _0)]b(\beta \omega _0)\}.`$
Here $`C`$ is a constant arising from the electron-phonon matrix elements, $`\beta =1/T`$, $`b(x)=(e^x1)^1`$, and $`\omega _0`$ is the phonon frequency. We use this expression now in the calculation of the electronic thermal conductivity $`\kappa _e`$. We mentioned above that the integrand of the $`\omega `$-integral in the expression for $`\kappa _e`$ given in Ref.  contains the term Im$`\epsilon _0=\gamma +\gamma _A`$ in the denominator which corresponds to the sum of relaxation rates of a quasiparticle due to impurity and Andreev scattering. The effect of thermal scattering by phonons can now be taken into account by adding $`1/\tau _{es}`$ to Im$`\epsilon _0`$. The temperature dependence of $`\kappa _e`$ is governed by the factor $`(\omega /T)^2\text{sech}^2(\omega /2T)`$ in the integrand of the $`\omega `$-integral which yields the main contribution of the integrand from the vicinity of $`x=\omega /T`$ at about $`x_m=2`$. Therefore we can estimate the temperature dependence of $`1/\tau _{es}`$ by carrying out the following variable transformations in Eq.(19): $`x=\omega /T`$, $`t=\omega _0/T`$, $`\mathrm{\Omega }=xy`$, and $`\mathrm{\Omega }_0=ty`$ with $`y=T/\mathrm{\Delta }(T,H)`$. This yields $`1/\tau _{es}`$ as a function of $`\mathrm{\Omega }=xy`$ and $`x`$. The normal state relaxation rate $`1/\tau _{en}`$ is obtained from Eq.(19) by setting $`A=1`$ and $`B=0`$. In Fig. 5 we have plotted $`\tau _{en}/\tau _{es}`$ as a function of $`y`$ for constant $`x=x_m=2`$ and $`\stackrel{~}{\mathrm{\Delta }}=`$ 0.2, 0.3, and 0.6 ($`H/H_{c2}`$= 0.78, 0.62, and 0.29, for $`\delta =0.1`$, see Fig. 1). One sees that for increasing $`\stackrel{~}{\mathrm{\Delta }}`$, or decreasing field, this function tends to a step function where the step occurs at $`y=1/2`$, that is, at a quasiparticle energy $`\mathrm{\Omega }=\omega /\mathrm{\Delta }=1`$. This result agrees with the previous result for the BCS state.  One recognizes that the variations of $`\tau _{en}/\tau _{es}`$ as a function of $`y`$ are much smaller than those of the phonon lifetime ratio $`\tau _{pn}/\tau _{ps}`$ shown in Fig. 4 where, for increasing $`\stackrel{~}{\mathrm{\Delta }}`$, a step occurs at about $`y=\mathrm{\Omega }/4=1/2`$, i.e., at a phonon energy $`\mathrm{\Omega }_0=\omega _0/\mathrm{\Delta }=2`$. The contribution of $`1/\tau _{es}`$ to the relaxation rate Im$`\epsilon _0`$ occuring in the electronic thermal conductivity $`\kappa _e`$ increases with temperature like $`T^3`$ because $`1/\tau _{en}`$ is proportional to $`T^3`$. In the presence of moderate impurity scattering this becomes comparable with $`\gamma `$ only at temperatures such that $`y=T/\mathrm{\Delta }(T,H)>1`$ where, according to Fig. 5, the ratio $`\tau _{en}/\tau _{es}`$ is nearly one. The expression for $`1/\tau _{es}`$ in Eq.(19) has been calculated in analogy to Ref.  for quasiparticle scattering by acoustic phonons. Instead of this spectrum, $`C\omega _0^2`$, one should employ in the expression for $`1/\tau _{es}`$ the transport Eliashberg functions $`\alpha _{tr}^2F_{nn}(\omega _0)`$, where $`n=\sigma ,\pi `$, which have been calculated for MgB$`_2`$.  Since the coupling constant $`\lambda _{\sigma \sigma }=0.81>\lambda _{\pi \pi }=0.41`$, we conclude that for the $`\sigma `$-band the field dependence of $`1/\tau _{es}`$ should be taken into account while, for the $`\pi `$-band, $`1/\tau _{es}`$ can safely be approximated by $`1/\tau _{en}`$.

## III Conclusions

We have calculated the density of states $`A(\omega )`$, the electronic thermal conductivity $`\kappa _e`$, and the scattering rates $`1/\tau _p(\omega )`$ of a phonon and $`1/\tau _e(\omega )`$ of a quasiparticle for a single isotropic gap over the entire field range of the mixed state. That the BPT  and P  methods for these calculations are valid approximation schemes has been proved by comparison of the results for $`A`$ with the solutions of the Eilenberger equations for Abrikosov’s vortex lattice order parameter.  By solving the gap equation including impurity self energy and vertex corrections (see Fig. 1), we have shown that the relation between the central parameter of the theory, $`\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }\mathrm{\Lambda }/v`$, and $`H/H_{c2}`$ can be very well approximated by the Ginzburg-Landau (GL) relation \[see Eq.(5)\] over the entire field range if one uses an appropriate value of the Abrikosov parameter $`\beta _A`$. This result is connected with the fact that the Abrikosov order parameter is a solution of the linearized GL equations which allows extension of the BPT method over the entire region of linear magnetization.

The second main result of the present paper is that the vertex corrections due to impurity scattering  substantially reduce the dependence of the functions $`\gamma _s(H)/\gamma _n`$, $`\kappa _{es}(H)/\kappa _{en}`$, and $`\tau _{pn}/\tau _{ps}(H)`$ on the amount of impurity scattering. The impurity scattering rate $`\delta =\mathrm{\Gamma }/\mathrm{\Delta }_0=(\pi /2)(\xi _0/\mathrm{})`$ is renormalized here in the Born approximation $`\mathrm{\Gamma }\gamma =\mathrm{\Gamma }\overline{A}(\omega )`$ where $`\overline{A}(\omega )`$ is calculated self-consistently from Eq.(1). The vertex corrections yield a renormalization of the gap $`\stackrel{~}{\mathrm{\Delta }}_r`$ due to the function $`D`$ \[see Eq.(7)\]. As can be seen from Fig. 2, the functions $`\gamma _s(H)/\gamma _n`$ for $`\delta =`$0.1, 0.5, and 1.0 are not very different from each other. The results are however quite different if we set the renormalization function $`D`$ equal to one. Then the slope of $`\gamma _s(H)/\gamma _n`$ versus $`H/H_{c2}`$ at $`H=0`$ rises rapidly with $`\delta `$ and becomes very large for $`\delta =1`$. It should be noticed that all curves in Fig. 2 tend to zero in the zero field limit as they should to be in accordance with Anderson’s theorem for ”dirty” superconductors. Our result that the effect of impurity scattering on the field dependence of $`\gamma _s(H)/\gamma _n`$ for a single isotropic gap is rather small and does not yield a rapid rise for small fields, gives support for the two-gap model in interpreting the measurements of $`\gamma _s(H)/\gamma _n`$ on single crystals of MgB$`_2`$.  According to this model, the first steep increase of $`\gamma _s(H)/\gamma _n`$, which is nearly the same for applied fields parallel and perpendicular to the c-axis, arises from the behavoir of the small gap associated with the three dimensional $`\pi `$-band. This in effect closes for increasing field due to the overlap of the vortex cores with large radius $`\xi _\pi 500\text{Å}`$  at a relatively small ”virtual” upper critical field $`H_{c2}^S`$ 3 - 4 kOe. The remaining superconducting contribution above $`H_{c2}^S`$ arises from the large gap associated with the two dimensional $`\sigma `$-band which closes at an upper critical field $`H_{c2}^{(c)}`$ 32 kOe and $`H_{c2}^{(ab)}`$ 160 kOe. The field dependence of the $`\pi `$-band and the $`\sigma `$-band contributions to $`\gamma _s(H)/\gamma _n`$ is presumably given by the almost linear curves in Fig. 2 for $`\delta =0,8`$ and $`\delta =0.5`$ with $`H_{c2}`$ equal to $`H_{c2}^S`$ and $`H_{c2}^{(c)}`$ or $`H_{c2}^{(ab)}`$, respectively.

The results for $`\kappa _{es}(H)/\kappa _{en}`$ shown in Fig. 3 depend on the impurity scattering rate $`\delta `$ primarily because the effect of Andreev scattering on the upward curvature becomes relatively smaller as $`\delta `$ increases. Although the slope at $`H=0`$ is relatively small even for $`\delta =1`$, it becomes very large at $`\delta =1`$ if the renormalization of the gap $`\stackrel{~}{\mathrm{\Delta }}_r=\stackrel{~}{\mathrm{\Delta }}/D`$ due to the vertex correction $`D`$ in $`\kappa _e`$ \[see Eq.(14)\] is neglected. This new result has the consequence that the observed steep rise of $`\kappa _e(H)`$ for small fields parallel or perpendicular to the c-axis in single-crystalline MgB$`_2`$  can again be explained only in terms of the two-gap model where the smaller gap in the three dimensional $`\pi `$-band closes in effect at a ”virtual” upper critical field $`H_{c2}^S`$3 - 4 kOe due to the overlap of the vortex cores with large radius $`\xi _\pi `$.  At higher fields $`\kappa _e(H)`$ tends to saturate until, in the vicinity of $`H_{c2}^{(c)}`$32 kOe, the contribution to $`\kappa _e`$ from the larger gap associated with the cylindrical $`\sigma `$-band exhibits an upward curvature. The contributions to $`\kappa _e`$ arising from the $`\pi `$\- and $`\sigma `$-bands are presumably given by the curves in Fig. 3 for $`\delta =0.8`$ and $`\delta =0.5`$ with the corresponding $`H_{c2}`$ equal to $`H_{c2}^S`$ and $`H_{c2}^{(c)}`$ where the curve for $`\delta =0.5`$ yields the measured upward curvature near $`H_{c2}`$.

The rapid decrease of $`\kappa `$ observed at low fields, which is attributed to the decrease of the phononic thermal conductivity $`\kappa _{ph}`$,  can be explained by the field and relevant temperature dependence of the ratio of phonon lifetimes $`\tau _{pn}/\tau _{ps}(H)=g(H)`$ shown in Fig. 4. One sees that $`g(H)`$ tends rapidly to the constant one for increasing field and values of $`y<1/2`$ corresponding to temperatures $`T<(1/2)\mathrm{\Delta }(T,H)`$. If we assume that $`\mathrm{\Delta }`$ is the smaller gap associated with the $`\pi `$-band, then the $`T`$ satisfying this condition are much smaller than $`T_c`$. However, at the lowest temperatures, this effect on $`\kappa _{ph}`$ is strongly reduced because the scattering rate $`1/\tau _{ps}`$ occuring in the denominator of the $`\omega `$-integral for $`\kappa _{ph}`$ is multiplied by $`T/T_c`$ in comparison to the scattering rate due to sample boundary scattering (see Eq.(4) of Ref. ).

We have calculated $`\kappa _e(H)`$ in the limit $`T0`$, or $`\omega 0`$. For finite $`T`$ one has to solve the transcendental equation determining the pole $`\epsilon _0`$ of the BPT-Green’s function  for each finite value of $`\omega `$ occuring in the $`\omega `$-integral for $`\kappa _e`$. It should be pointed out that this calculation is avoided in the expression for $`\kappa _e`$ which was derived from linear response theory.  For higher temperatures one should add to the sum of scattering rates $`\gamma `$ and $`\gamma _A`$ due to impurity and Andreev scattering \[see Eq.(14)\], the scattering rate $`1/\tau _{es}(H)`$ due to scattering of quasiparticles by phonons. The field and relevant temperature dependence of the ratio $`\tau _{en}/\tau _{es}(H)`$ is shown in Fig. 5 as a function of $`y=T/\mathrm{\Delta }(T,H)`$. One sees that this ratio is nearly equal to one for values of $`y>1`$ corresponding to temperatures $`T>\mathrm{\Delta }(T,H)`$. For the small gap associated with the $`\pi `$-band this condition is probably satisfied at those temperatures where $`1/\tau _{en}T^3`$ becomes comparable to the intraband impurity scattering $`\gamma `$ for the $`\pi `$-band. However, for the large gap associated with the $`\sigma `$-band, this condition might not be satisfied. We conclude from these estimates that, for the $`\pi `$-band, $`1/\tau _{es}`$ can be approximated by $`1/\tau _{en}`$ while, for the $`\sigma `$-band, the field dependence of $`1/\tau _{es}(H)`$ should be taken into account. These considerations are supported by detailed calculations which yield a larger intraband impurity scattering rate $`\gamma `$ for the $`\pi `$-band than for the $`\sigma `$-band,  and a larger scattering rate $`1/\tau _{en}`$ for the $`\sigma `$-band than for the $`\pi `$-band because the intraband electron-phonon coupling constant satisfies $`\lambda _{\sigma \sigma }>\lambda _{\pi \pi }`$.

Our assumption that we can simply add our results obtained for two different isotropic gaps and different upper critical fields in the $`\pi `$\- and $`\sigma `$-bands is a crude approximation in view of the actual situation in MgB$`_2`$. It is true that the interaction due to interband impurity scattering between the $`\sigma `$\- and $`\pi `$-bands is small.  However, we have neglected the pairing interaction between the two bands which presumably leads to Cooper pair tunneling from the $`\sigma `$-band with strong pairing to the $`\pi `$-band where it leads to giant vortex cores and where superconductivity is maintained above the virtual $`H_{c2}^S`$ up to $`H_{c2}`$.  Nevertheless we believe that our clear-cut results for the mixed state of a single isotropic gap are helpful in analyzing the complicated behavior of different physical quantities in MgB$`_2`$.

Another approximation is that we have neglected higher harmonic oscillator components N of Abrivosov’s vortex lattice order parameter which lead to a distortion of the lattice at lower fields.  For an s-wave superconductor the BPT method breaks down at low fields where properties are determined by the states bound to the vortex cores. However, STS measurements in MgB2 show an absence of localized states in the cores.

###### Acknowledgements.

We would like to thank T. Dahm and K. Scharnberg for valuable discussions.


format: text

text:
# Extended Recurrence Plot Analysis and its Application to ERP Data

## 1 Introduction

Neurons are known to be nonlinear devices because they become activated when their somatic membrane potential crosses a certain threshold \[Kandel et al., 1995\]. This nonlinearity is one of the essentials in neural modelling which leads to the sigmoidal activation functions of neural networks \[Amit, 1989\]. The activity of large formations of neurons is macroscopically measurable as the electroencephalogram (EEG) at the human scalp which results from a spatial integration of postsynaptic potentials \[Nunez, 1981\]. However, it is an unsolved problem whether the EEG should be treated as a time series stemming from a linear or a nonlinear dynamical system. Applying nonlinear techniques of data analysis to EEG measurements has a long tradition. Most of these efforts have been done by computing the correlation dimension of spontaneous EEG \[e. g. Babloyantz et al., 1985; Rapp et al., 1986; Gallez and Babloyantz, 1991; Lutzenberger et al., 1992; Pritchard and Duke, 1992\]. Theiler et al.  applied the technique of surrogate data to correlation dimensions of EEG and reported that there is no evidence of low-dimensional chaos but of significance for nonlinearity in the data.

While correlation dimensions are only well defined for stationary time series generated by a low-dimensional dynamical system moving around an attractor, these measures fail in investigating event-related brain potentials \[ERPs, Sutton et al., 1965\] since they are nonstationary by definition. Traditionally, ERP waveforms are determined by computing an ensemble average over a collection of stimulus time locked EEG trials. This is based on the following assumptions: (1) the presentation of stimuli of the same kind is followed by the same sequence of processing steps, (2) these processing steps always lead to activation of the same brain structures, (3) this activation always elicits the same pattern of electrophysiological activity, which can be measured at the scalp \[Rösler, 1982\] and (4) spontaneous activity is stationary and ergodic \[beim Graben et al., 2000\].

By averaging the data-points time-locked to the stimulus presentation (cf. Oddball experiment) it is possible to filter out the signal (ERP) of the noise (spontaneous activity). In the next step the functional significance of a component is assessed. Antecedent conditions of the occurrence of a component and variables, which influence its parameters are defined. Now the commonalities of these factors are identified. The generalization of all empirically found influencing factors leads to a more abstract cognitive theory of the functional meaning of a event-related potential component and makes it usable for the validation of models of cognitive processes.

The disadvantage of the averaging method is the high number of trials needed to reduce the signal-to-noise-ratio \[Kutas and Petten, 1994\]. This is crucial for example for clinical studies, for studies with children and for studies, in which repeating a task would influence the performance. So it is desirable to find new ways of analyzing event-related activity on a single-trial basis. Applying nonlinear methods to electrophysiological data could be one way of dealing with this problem.

To compute dimensions of ERPs, Molnár et al.  used the pointwise dimensions and reported a drop of the pointwise dimension as a function of time corresponding to the P300 component observed in the Oddball experiment. Recently, concepts of information theory have been introduced to analyse ERPs. On one hand this is the wavelet entropy of Quiroga et al.  and on the other hand symbolic dynamics of EEG and ERP \[beim Graben et al., 2000; Frisch et al., 2002; Steuer, 2002; Schack, 2002\].

A further promising approach is the recurrence quantification analysis (RQA), which is based on the quantification of the diagonal oriented structures in recurrence plots \[RPs, Webber Jr. and Zbilut, 1994; Zbilut and Webber Jr., 1992\]. The RQA was broadly applied in a wide field of the analysis of physiological data \[e. g. Casdagli, 1997; Faure and Korn, 1998; Thomasson et al., 2001; Marwan et al., 2002\]. The important advantage of methods based on the quantification of RPs is that the required data length can be relatively short. However, the measures of the classical RQA are only able to recognize transitions between periods and chaos and vice versa \[Trulla et al., 1996\]. In this work, we will use recently introduced additional measures based on RPs in order to find chaos-chaos transitions in physiological data. These new measures use the vertical structures in the RP and are able to identify laminar states \[Marwan et al., 2002\].

In the first section we will give a short introduction into RPs and their quantification analysis. In the next section we will introduce the new measures and finally we will apply them to event related potential data gained from the Oddball experiment.

## 2 Recurrence Plots and Their Quantification

The method of recurrence plots (RP) was introduced to visualize the time dependent behavior of the dynamics of systems, which can be pictured as a trajectory in the phase space \[Eckmann et al., 1987\]. It represents the recurrence of the $`m`$-dimensional phase space trajectory $`\stackrel{}{x}_i^m`$ ($`i=1,\mathrm{},N`$, time discrete) to a certain state. The main step of this visualization is the calculation of the $`N\times N`$-matrix

$$𝐑_{i,j}:=\mathrm{\Theta }(\epsilon _i\stackrel{}{x}_i\stackrel{}{x}_j),i,j=1\mathrm{}N,$$
(1)
where $`\epsilon _i`$ is a state dependent cut-off distance, $``$ is the norm of vectors, $`\mathrm{\Theta }`$ is the Heaviside function and $`N`$ is the number of states. The phase space vectors for one-dimensional time series $`u_i`$ from observations can be reconstructed with the Taken’s time delay method $`\stackrel{}{x}_i=(u_i,u_{i+\tau },\mathrm{},u_{i+(m1)\tau })`$ with dimension $`m`$ and delay $`\tau `$ \[Kantz and Schreiber, 1997\]. The recurrence plot exhibits characteristic large-scale and small-scale patterns which are caused by typical dynamical behavior \[Eckmann et al., 1987; Webber Jr. and Zbilut, 1994\], e. g. diagonals (similar local time evolution of different parts of the trajectory) or horizontal and vertical black lines (state does not change for some time).

Zbilut and Webber have developed the recurrence quantification analysis (RQA) to quantify an RP \[Webber Jr. and Zbilut, 1994; Zbilut and Webber Jr., 1992\]. They defined measures using the recurrence point density and diagonal structures in the recurrence plot, the recurrence rate $`RR`$ (density of recurrence points), the determinism $`DET`$ (ratio of recurrence points forming diagonal structures to all recurrence points), the maximal length of diagonal structures $`L_{max}`$ (or their averaged length $`L`$), the Shannon entropy $`ENT`$ of the distribution of the diagonal lengths and the trend $`TREND`$ (paling in the RP). The computation of these measures in shifted windows along the main diagonal of the RP enables one to find characteristic excursions of the trajectory in the phase space of the considered systems.

Trulla et al. have applied these measures in order to find transitions in dynamical systems \[Trulla et al., 1996\]. They have showed, that the RQA is able to find transitions between chaos and order (periodical states). But they could not find the chaos-chaos transitions.

## 3 Laminarity and Trapping Time

We have recently introduced two additional measures which are based on the vertical structures in the RP \[Marwan et al., 2002\]. We define these measures analogous to the definition of $`DET`$ and $`L`$, but we consider the distribution $`P(v)`$ of the length of the vertical structures in the RP.

First, the laminarity $`LAM`$

$$LAM:=\frac{_{v=2}^NvP(v)}{_{v=1}^NvP(v)},$$
(2)
is the ratio of recurrence points forming vertical structures to all recurrence points and represents the probability of occurrence of laminar states in the system, but it does not describe the length of these laminar phases. It will decrease if the RP consists of more single recurrence points than vertical structures.

Next, the trapping time $`TT`$

$$TT:=\frac{_{v=2}^NvP(v)}{_{v=2}^NP(v)},$$
(3)
is the averaged length of the vertical structures. The measure $`TT`$ contains information about the amount and the length of the laminar phases.

The difference between these measures and the traditional RQA measures is their ability to find transitions between chaos and chaos \[Marwan et al., 2002\]. For example, such transitions can be found in the logistic map $`x_{n+1}=ax_n\left(1x_n\right)`$ with increasing control parameter $`a[0,4]`$ and $`x_n[0,1]`$. For such trajectories $`x(a)`$ which contain laminar states (e. g. $`a=3.678,3.791,3.927`$), $`LAM`$ and $`TT`$ show pronounced maxima (Fig. 1). The application of these measures to heart rate variability data, has shown, that they are able to detect and quantify laminar phases before a life-threatening cardiac arrhythmia and, thus, to enable a prediction of such an event \[Marwan et al., 2002\]. These findings can be of importance for the therapy of malignant cardiac arrhythmias.

In the next section we will apply this extended RQA to physiological data.

## 4 Event Related Potentials

### 4.1 The Oddball experiment

As mentioned in the Introduction, the Oddball experiment studies brain potentials during a stimulus presentation.

The measurement of the EEG was done with 31 electrodes/ channels (Tab. 1). The first 25 electrodes were localized as shown in Figure 2; the others were reference electrodes. The sample interval for the measurements was 4 ms.

Probands were seated in a dimly lit room in front of a monitor and were instructed to count tones of high pitch. Each subject was tested in nine blocks. The blocks varied in the probability of occurrence of the higher tones from 10 to 90 %. Each block contained at least 30 target tones. Response was given in a three alternative choice (using cursor keys of the keyboard). During the test, the EEG was recorded. The stimuli were computer-generated beeps of 100 ms length. Tones were either high (1400 Hz) or low (1000 Hz). They were presented with an interstimulus interval of 1000 ms.

After computing event-related voltage averages for the experimental manipulations (10 % up to 90 % target probability) one can observe a P300 ERP component whose amplitude is anti-correlated to the probability of the stimuli (surprise ERP, Fig. 3).

The P300 component of the ERP was the first potential discovered to vary in dependence on subject-internal factors like attention and expectation instead of physical characteristics \[Sutton et al., 1965\]. The amplitude of the P300 component is highly sensitive to novelty of an event and its relevance. So this component is assumed to reflect the updating of the environmental model of the information processing system \[context updating, Donchin, 1981; Donchin and Coles, 1988\].

### 4.2 Data analysis

Our focus will be directed to the ERP data of two extreme event probabilities. Henceforth, the time (measured in ms) is denoted as $`t`$, the trial number as $`i`$ and the electrode as $`e`$ (the allocation of the electrode numbers with their notion, see Fig. 2).

The first set ERP90 contains 40 trials of ERP data for an event frequency of 90 % and the second set ERP10 contains 31 trials for an event frequency of 10 %. Both data sets can be rather well discriminated in the N100 and P300 components by the average over all trials (Fig. 3). As expected, both components have increased for lower event probabilities (ERP10). The maxima of the P300 are located around the central and central-parietal electrodes. However, the single trials do not obtain such a clear result. The P300 component is only well pronounced in 15 trials. When the single trials are observed, then extreme values can also occur in the ERP90 data and vanish in the ERP10 data (Fig. 4). We applied also a statistical variance-based T-test to the single trial ERP data. However, this method could also not clearly distinguish the single trials.

The recurrence quantification analysis (RQA) is based on the structures obtained by recurrence plots (RPs). The RPs were firstly computed for the means of ERP90 and ERP10 over all trials and then for the single trials. This was done with the embedding parameters $`m=3`$, $`\tau =3`$ and $`\epsilon =10\%`$ (fixed amount of nearest neighbours). The embedding parameters were estimated by using the standard methods false nearest neighbours (dimension) and mutual information (delay) \[Kantz and Schreiber, 1997\]. Due to the N100 and the P300 components in the data, the RPs show varying structures changing in time (Fig. 5). Diagonal structures and clusters of black points occur. The nonstationarity of the data around the N100 and P300 causes extended white bands along these times in the RPs. However, the clustered black points around 300 ms occur in almost all RPs of the ERP10 data set.

The RQA was computed from the RPs of ERP90 and ERP10 for the single trials, in sliding windows over the RPs (which have the dimension $`m=3`$) with a length of 240 ms and with a shifting step of 8 ms. This window length corresponds with a data length of 60 values.

The mean of all RQA variables of ERP10 reveal typical structures in the data (Fig. 6, right column). They indicate the transitions corresponding to the N100 and P300 components around the central electrodes. The RQA variables for the ERP90 do not reveal these transitions (Fig. 6, left column). The onset of the increasing of the parameter is about 120 ms before the event. This is due to the windowed analysis of the RPs (240 ms windows). We have chosen the middle of the RP window for the time, what results in a 120 ms earlier onset of the RQA variables.

The four RQA variables are quite different, especially in their amplitude. For ERP10, $`LAM`$ and $`TT`$ are the best pronounced parameter and have two distinct maxima at some electrodes; $`DET`$ and $`L`$ reveal these maxima at these electrodes too, but are lesser pronounced (Fig. 6). These maxima occur at the transition around 100 ms and 300 ms after the event and occur at the electrodes F3, F4, FZ, C3, FCZ, PZ, POZ and PO3. Differences between the various transitions found by these measures also occur in time and brain locations (electrodes). But, the study was not detailed enough in order to give reliable results.

The analysis of the single trials achieves similar results (Figs. 7 and 8 show the results for selected trials). The $`LAM`$ clearly found the N100 and P300 components for ERP10 in 26 trials (of 31), but not for the ERP90 trials. The other measures have lesser maxima and, thus, are not suitable for such recognition.

This result indicates that our introduced measures of complexity (especially $`LAM`$) are able to recognize transitions in brain potentials, which are caused by e. g. stimulative events. These transitions can be found in the single trials, which is an improvement to the classical method of averaging all observations.

## 5 Summary

We have applied an extended recurrence quantification analysis (RQA) to physiological event related potential data (ERP). The classical RQA consists of measures which are mainly based on diagonal structures in the recurrence plots (RPs), e. g. the determinism ($`DET`$), which is the ratio of recurrence points located on connected diagonal structures in the RP, and the averaged diagonal line length ($`L`$). We have extended the RQA with two recently introduced measures, the laminarity ($`LAM`$) and the trapping time ($`TT`$). These measures are analogously defined as $`DET`$ and $`L`$, but provided by the vertical structures in recurrence plots. Whereas the classical RQA enables the identification of period-chaos transitions, the new measures make the identification of chaos-chaos transitions and laminar states possible.

The classical method to study ERP data is to average them over many trials. Our aim was to study the single trials in order to find transitions in the data.

The application of the extended RQA to ERP data has discriminated the single trials with a distinct P300 component due to a high surprise moment (less frequent events) against such trials with a low surprise moment (high frequent events). Considering the raw ERP10 data, the P300 component can only be found in the half of all trials. Also a statistical variance test fails to distinguish cleary the trials. The $`LAM`$ is the most pronounced parameter in this analysis. It measures the ratio of recurrence points located on connected vertical structures in an RP. This structures correspond with laminarity within the underlying process. In the ERP data, the $`LAM`$ reveals transitions from less laminar states to higher laminar states after the occurrence of the event and a transition from higher laminar states to less laminar states after about 400 ms. These transitions occur around bounded brain areas (parietal to frontal along the central axis). The comparable measures $`DET`$/ $`LAM`$ and $`L`$/ $`TT`$ are quite different in their amplitude. There should also be differences in time and brain location of the found transitions.

These results show that the measures based on vertical RP structures make the identification of transitions possible, which are not found by the classical RQA measures. They indicate transitions in the brain processes into laminar states due to the surprising moment of observed events.

A future work will be concerned with the development of a statistical evaluation of these results. Furthermore, this investigation has to be extended to ERP data gained from other frequent events and a detailed study of the comparable measures $`DET`$/ $`LAM`$ and $`L`$/ $`TT`$ should give hints about the different transitions in the brain processes.

## 6 Acknowledgments

This work was partly supported by the priority programme SPP 1097 of the German Science Foundation (DFG). We gratefully acknowledge the colleagues of the Nonlinear Dynamics Group Potsdam for useful discussions and support of this work.