text:
# Fermi surface topology and the upper critical field in two-band superconductors - application to MgB2

## Abstract

Recent measurements of the anisotropy of the upper critical field $`B_{c2}`$ on MgB2 single crystals have shown a puzzling strong temperature dependence. Here, we present a calculation of the upper critical field based on a detailed modeling of bandstructure calculations that takes into account both the unusual Fermi surface topology and the two gap nature of the superconducting order parameter. Our results show that the strong temperature dependence of the $`B_{c2}`$ anisotropy can be understood as an interplay of the dominating gap on the $`\sigma `$-band, which possesses a small $`c`$-axis component of the Fermi velocity, with the induced superconductivity on the $`\pi `$-band possessing a large $`c`$-axis component of the Fermi velocity. We provide analytic formulas for the anisotropy ratio at $`T=0`$ and $`T=T_c`$ and quantitatively predict the distortion of the vortex lattice based on our calculations.

Our understanding of the physical properties of the recently discovered superconductivity in MgB2 has made rapid progress since its discovery Akimitsu . Its high critical temperature $`T_c=39`$ K can be understood as arising from strong conventional electron-phonon coupling to a high frequency phonon mode Dolgov ; Liu ; Heid . Its pairing symmetry seems to be of conventional $`s`$-wave type Manzano ; Yang . However, in contrast to conventional superconductors, a number of recent experiments indicate that there exist two gaps of different size in this compound Szabo ; Giubileo ; Bouquet ; Sologubenko ; Bouquet2 ; Iavarone . This possibility is supported by band structure calculations, which have shown that the Fermi surface of this compound consists of four bands: two $`\sigma `$-type two-dimensional cylindrical hole sheets and two $`\pi `$-type three dimensional tubular networks Shulga ; Kortus  in good overall agreement with recent de Haas-van Alphen experiments Carrington . Microscopic calculations of the superconducting gap based on band structure calculations have shown recently that indeed one should expect a big superconducting gap living on the $`\sigma `$ bands and a smaller one, induced by interband pairing interaction, living on the $`\pi `$ bands Liu ; Louie . Impurity scattering, which in conventional superconductors tends to average out strongly differing gap values, in this case becomes ineffective, because the $`\sigma `$ and $`\pi `$ bands possess different symmetries, making interband impurity scattering much weaker than intraband impurity scattering Mazin .

Recent measurements of the upper critical field $`B_{c2}`$, particularly its anisotropy, on single crystal MgB2 have shown a puzzling strong temperature dependence of the anisotropy ratio $`B_{c2}^{ab}/B_{c2}^c`$ between the ab-plane and the c-axis upper critical field Angst ; Eltsev ; Lyard . In conventional systems this ratio rarely changes by more than 10 to 20 percent as a function of temperature. In MgB2 changes by more than a factor of 2 have been observed. It has been shown that such a strong temperature dependence of the anisotropy ratio can be obtained within a strongly anisotropic single gap model, possessing a small gap in c-axis direction and a big gap in ab-plane direction Posa1 . However, the gap anisotropy would have to be a factor of 10, which is too big as compared with experimental values. In addition, this scenario would be inconsistent with penetration depth studies which clearly indicate the presence of a small gap within the ab-plane Jin ; Posa2 .

Here, we want to present first calculations of the upper critical field $`B_{c2}`$, which take into account the multi-band Fermi surface structure seriously. We show that the strong temperature dependence of the anisotropy of $`B_{c2}`$ can be traced back to the influence of the two topologically very different Fermi surface types. While the cylindrical Fermi surface sheets are dominating the behavior of $`B_{c2}`$ at low temperatures, leading to a large anisotropy, at temperatures approaching $`T_c`$ the $`\pi `$ bands due to their much larger c-axis Fermi velocity start to play a more important role, strongly reducing the $`B_{c2}`$ anisotropy. Thus, the strong temperature dependence of the $`B_{c2}`$ anisotropy appears as a cross-over from a low temperature $`\sigma `$ band dominated regime to a higher temperature mixed $`\sigma `$ and $`\pi `$ band regime.

In order to include the multi-band Fermi surface structure in the calculation of the upper critical field, we start our investigation from the fully momentum dependent multi-band formulation of the quasiclassical (Eilenberger) theory of the upper critical field Rieck . For that purpose we have to solve the linearized multi-band gap equation in the presence of an external magnetic field, which reads

$$\mathrm{\Delta }_\alpha (\stackrel{}{r})=\pi T\underset{\alpha ^{}}{}\underset{|\omega _n^{}|<\omega _c}{}\lambda ^{\alpha \alpha ^{}}f_\alpha ^{}(\stackrel{}{r},\widehat{k}^{};\omega _n^{})_\alpha ^{}$$
(1)
Here, $`f_\alpha `$ is the anomalous Eilenberger propagator on the Fermi surface sheet denoted by $`\alpha `$. $`\mathrm{\Delta }_\alpha `$ is the gap function on that Fermi surface sheet and $`\lambda ^{\alpha \alpha ^{}}`$ is the pairing interaction, which becomes a matrix in the band indices. The brackets $`\mathrm{}_\alpha ^{}`$ denote a Fermi surface average over momentum $`\widehat{k}^{}`$ of the Fermi surface sheet $`\alpha ^{}`$. In Eq. (1) we have already assumed that the gaps are isotropic $`s`$-wave on the Fermi surfaces, as indicated by experiment, but may have different values on different sheets. At $`B_{c2}`$ the anomalous Eilenberger propagator has to be determined from the linearized Eilenberger equation (for $`\omega _n>0`$)

$$\left\{\omega _n+\stackrel{}{v}_{F,\alpha }\left[\frac{\mathrm{}}{2}\stackrel{}{}i\frac{e}{c}\stackrel{}{A}\left(\stackrel{}{r}\right)\right]\right\}f_\alpha (\stackrel{}{r},\widehat{k};\omega _n)=\mathrm{\Delta }_\alpha (\stackrel{}{r})$$
(2)
Here, $`\stackrel{}{v}_{F,\alpha }`$ is the (momentum dependent) Fermi velocity on Fermi surface sheet $`\alpha `$ and $`\stackrel{}{A}`$ the vector potential due to the magnetic field $`\stackrel{}{B}=\stackrel{}{}\times \stackrel{}{A}`$. Eqs. (1) and (2) constitute an eigenvalue problem for $`\mathrm{\Delta }_\alpha (\stackrel{}{r})`$ . For a given temperature $`T`$ the solution $`\mathrm{\Delta }_\alpha (\stackrel{}{r})`$ which solves Eqs. (1) and (2) for the highest value of $`B`$ determines $`B_{c2}`$.

Usually this eigenvalue problem is solved by a Landau level expansion of $`\mathrm{\Delta }_\alpha (\stackrel{}{r})`$ above the Abrikosov ground state of the vortex lattice. It has been shown recently, however, that a variational ansatz for $`\mathrm{\Delta }_\alpha (\stackrel{}{r})`$ corresponding to a distorted Abrikosov lattice leads to much better results for strongly anisotropic systems Posa1  and we will adopt that method here. In this method the ansatz reads $`\mathrm{\Delta }_\alpha (\stackrel{}{r})=\mathrm{\Delta }_\alpha \psi _\mathrm{\Lambda }^\tau (\stackrel{}{r})`$ with $`\psi _\mathrm{\Lambda }^\tau (x,y)=\psi _\mathrm{\Lambda }(e^\tau x,e^\tau y)`$. Here, $`\psi _\mathrm{\Lambda }`$ is the usual Abrikosov groundstate and $`\tau `$ is a variational parameter describing the distortion of the vortex lattice. The undistorted lattice corresponds to $`\tau =0`$ and $`\tau `$ has to be determined by maximizing $`B_{c2}`$. Introducing this ansatz into Eqs. (1) and (2) and using standard operator techniques (see e.g. Refs. Rieck ; Schopohl2 ; Graser ; SunMaki ) we are led to the following eigenvalue problem for $`\mathrm{\Delta }_\alpha `$ in band space:

$$\mathrm{\Delta }_\alpha =\underset{\alpha ^{}}{}\lambda ^{\alpha \alpha ^{}}\left[\frac{1}{\lambda _+}\mathrm{ln}\frac{T}{T_c}l_\alpha ^{}(\tau ,\frac{B_{c2}}{T^2})\right]\mathrm{\Delta }_\alpha ^{}$$
(3)
Here, $`\lambda _+`$ is the highest eigenvalue of $`\lambda ^{\alpha \alpha ^{}}`$, which determines $`T_c`$. The function $`l_\alpha `$ is given by the expression

$$l_\alpha =_0^{\mathrm{}}\frac{du}{\mathrm{sinh}u}1e^{u^2\frac{eB_{c2}}{8\pi ^2T^2}\left(e^{2\tau }v_{F1,\alpha }^2(\widehat{k})+e^{2\tau }v_{F2,\alpha }^2(\widehat{k})\right)}_\alpha $$
(4)
Here, $`v_{F1,\alpha }`$ and $`v_{F2,\alpha }`$ are the components of the Fermi velocity perpendicular to the magnetic field $`\stackrel{}{B}`$ on Fermi surface $`\alpha `$. According to Eq. (3) the criterion for $`B_{c2}`$ is that the highest eigenvalue of the matrix $`\lambda ^{\alpha \alpha ^{}}\left[\frac{1}{\lambda _+}\mathrm{ln}\frac{T}{T_c}l_\alpha ^{}(\tau ,\frac{B_{c2}}{T^2})\right]`$ becomes 1. Apparently, at $`T=T_c`$ this is fulfilled for $`B_{c2}=0`$, because $`l_\alpha 0`$.

Eqs. (3) and (4) allow to determine the temperature dependence and angular dependence of $`B_{c2}`$ from microscopic grounds. The material parameters we need for the solution are the Fermi velocities $`\stackrel{}{v}_{F,\alpha }(\widehat{k})`$ and the coupling matrix $`\lambda ^{\alpha \alpha ^{}}`$, which can be taken from bandstructure calculations. In order to simplify the analysis we restrict ourselves to two relevant bands, because the two $`\sigma `$-bands and the two $`\pi `$-bands are very similar Liu . The $`\sigma `$-band can be described to a good approximation by a cylindrical Fermi surface with a small $`c`$-axis hopping parameter. The $`\pi `$-band can be modeled by a half-torus as shown in Fig. 1 (for comparison see the $`\mathrm{\Gamma }`$-point centered Fermi surfaces obtained from bandstructure calculations in Fig. 1 in Ref. Shulga ). For the Fermi velocities on these two Fermi surfaces we thus write $`\stackrel{}{v}_{F,\pi }(\theta ,\varphi )=v_{F\pi }(\mathrm{cos}\theta \mathrm{cos}\varphi ,\mathrm{cos}\theta \mathrm{sin}\varphi ,\mathrm{sin}\theta )`$ and $`\stackrel{}{v}_{F,\sigma }(k_c,\varphi )=v_{F\sigma }(\mathrm{cos}\varphi ,\mathrm{sin}\varphi ,ϵ_c\mathrm{sin}ck_c)`$. Here, $`\varphi `$ is the azimuthal angle within the $`ab`$-plane, $`\theta [\frac{\pi }{2},\frac{3\pi }{2}]`$ the polar angle of the torus, $`k_c`$ the $`c`$-axis component of the momentum, and $`c`$ the lattice constant in $`c`$-direction. The dimensionless parameter $`ϵ_c`$ describes the small $`c`$-axis dispersion of the cylinder. The Fermi surface averages over the cylinder and the torus are then given by $`\mathrm{}_\sigma =\frac{c}{4\pi ^2}_{\pi /c}^{\pi /c}𝑑k_c_0^{2\pi }𝑑\varphi \mathrm{}`$ and $`\mathrm{}_\pi =\frac{1}{2\pi ^2}_{\pi /2}^{3\pi /2}𝑑\theta _0^{2\pi }𝑑\varphi \frac{1+\kappa \mathrm{cos}\theta }{12\kappa /\pi }\mathrm{}`$, where $`\kappa `$ is the ratio of the two radii of the torus. The parameters of the Fermi velocities can be found from bandstructure calculations. Taking the values given in Ref. Brinkman , we determine $`v_{F\pi }=8.210^5`$m/s, $`v_{F\sigma }=4.410^5`$m/s, and $`ϵ_c=0.23`$. The ratio of the two radii of the torus can be estimated from the Fermi surface crossings in Ref. Kortus  to be about $`\kappa =0.25`$.

For a $`2\times 2`$ matrix $`\lambda ^{\alpha \alpha ^{}}`$ the criterion that the biggest eigenvalue of Eq. (3) becomes 1 leads to the equation $`(1\eta )l_\sigma +\eta l_\pi +\mathrm{ln}t=\mathrm{\Lambda }_\pm \left(l_\sigma +\mathrm{ln}t\right)\left(l_\pi +\mathrm{ln}t\right)`$. Here, $`t=T/T_c`$, $`\lambda _{}`$ is the smaller eigenvalue of $`\lambda ^{\alpha \alpha ^{}}`$ and $`\eta =\frac{\lambda ^{\pi \pi }\lambda _{}}{\lambda _+\lambda _{}}`$ is a dimensionless parameter describing the interband coupling strength of the two bands ($`\eta =0`$ corresponds to no coupling, $`\eta =0.5`$ to maximum coupling). From bandstructure calculations in Ref. Liu  the following effective matrix elements can be obtained: $`\lambda ^{\sigma \sigma }=0.959`$, $`\lambda ^{\sigma \pi }=0.222`$, $`\lambda ^{\pi \sigma }=0.163`$, and $`\lambda ^{\pi \pi }=0.278`$. From these we find $`\lambda _+=1.008`$, $`\lambda _{}=0.228`$, and $`\eta =0.064`$. As it turns out the remaining parameter $`\mathrm{\Lambda }_\pm =\frac{\lambda _+\lambda _{}}{\lambda _+\lambda _{}}`$ only weakly affects the results as soon as $`\lambda _{}`$ is sufficiently smaller than $`\lambda _+`$, as is the case here.

We have calculated $`B_{c2}`$ numerically for the material parameters above optimizing the parameter $`\tau `$ such that $`B_{c2}`$ is maximized. In Fig. 2 we show our result for the anisotropy ratio $`\mathrm{\Gamma }=B_{c2}^{ab}/B_{c2}^c`$ as a function of temperature for different values of the interband coupling strength $`\eta `$. For $`\eta =0`$, when there is no coupling between the two bands, the temperature dependence of $`\mathrm{\Gamma }`$ is determined by the cylindrical $`\sigma `$-band because of its stronger pairing interaction. Here, $`\mathrm{\Gamma }`$ changes only by 20 percent, as one expects for an isotropic single gap superconductor. When $`\eta `$ is increased, however, the temperature dependence of $`\mathrm{\Gamma }`$ becomes more pronounced. Our result for the parameters given above is shown as the dashed line. Note, that our calculation is parameter free, relying only on the parameters given by bandstructure calculations. For comparison, also the experimental data by Lyard et al. Lyard  are shown (solid circles). As will become clear below, the most important parameters determining the temperature dependence of $`\mathrm{\Gamma }`$ are the interband coupling strength $`\eta `$ and the $`c`$-axis dispersion parameter $`ϵ_c`$. If we allow these two parameters to vary somewhat, we can obtain an excellent fit of the experimental data. The dotted line shows our result for $`\eta =0.121`$ and $`ϵ_c=0.182`$, where the other parameters have been kept constant. We mention here that this set of parameters also gives a correspondingly good fit of the temperature dependences of $`B_{c2}^c`$ and $`B_{c2}^{ab}`$ separately including the upward curvature of $`B_{c2}^{ab}`$ that has been noted in the experiments Graser2 . This is just an immediate consequence of the strong temperature dependence of the anisotropy ratio and $`B_{c2}^c`$ varying linearly near $`T_c`$ because of an absence of vortex lattice distortion in this field direction.

Apparently, the anisotropy ratio $`\mathrm{\Gamma }`$ turns out to be more sensitive to interband coupling $`\eta `$ close to $`T_c`$ than at $`T=0`$. In order to get a better physical understanding for this behavior, we want to discuss some limits, in which analytical solutions for $`B_{c2}`$ can be obtained. At first we consider the limit $`TT_c`$. In this case Eqs. (3) and (4) can be solved exactly and we obtain for the angular dependence of $`B_{c2}`$ as a function of the angle $`\beta `$ the magnetic field makes with the $`c`$-axis:

$$\frac{B_{c2}(\beta )}{B_{c2}(\beta =0)}=(\mathrm{cos}^2\beta +A\mathrm{sin}^2\beta )^{1/2}$$
(5)
where $`A=2\frac{\left(1\eta \right)v_{c,\sigma }^2_\sigma +\eta v_{c,\pi }^2_\pi }{\left(1\eta \right)v_{ab,\sigma }^2_\sigma +\eta v_{ab,\pi }^2_\pi }`$. Here, $`v_{c,\alpha }^2_\alpha `$ is the average of the squared $`c`$-axis component of the Fermi velocity over Fermi surface $`\alpha `$, while $`v_{ab,\alpha }^2_\alpha `$ is the corresponding in-plane quantity. The distortion parameter $`\tau `$ is given by $`e^{2\tau }=\frac{B_{c2}(\beta )}{B_{c2}(\beta =0)}`$. When we introduce the numerical values for $`v_{F\pi }`$, $`v_{F\sigma }`$ and $`\kappa `$ given above we find for the anisotropy ratio

$$\frac{B_{c2}^{ab}}{B_{c2}^c}=\frac{1}{\sqrt{A}}=\frac{1}{ϵ_c}\sqrt{\frac{1+0.63\eta }{1+\eta \left(\frac{3.69}{ϵ_c^2}1\right)}}$$
(6)
For $`ϵ_c=0.23`$ this result is shown in Fig. 3 as a function of $`\eta `$ (solid line). From Eq. (6) we see that at $`\eta =0`$ the anisotropy is given by the $`c`$-axis dispersion of the cylindrical Fermi surface and actually diverges for $`ϵ_c0`$. (Note, that this divergence could not have been obtained from a Landau level expansion above the Abrikosov ground state and our variational ansatz above is crucial for this result). Once the interband coupling $`\eta `$ is increased, the anisotropy quickly reduces. This reduction becomes sizeable already, when $`\eta \frac{v_{c,\sigma }^2_\sigma }{v_{c,\pi }^2_\pi }=\frac{ϵ_c^2}{3.69}`$. Thus, it is the small $`c`$-axis Fermi velocity of the cylindrical Fermi surface as compared to the one of the $`\pi `$-band which leads to a high sensitivity of $`\mathrm{\Gamma }`$ near $`T_c`$ and a small interband coupling is already sufficient to make the influence of the $`\pi `$-band visible.

In the limit $`T0`$ the integration in Eq. (4) can be performed and we find $`l_\alpha (\tau ,\frac{B_{c2}}{T^2})=\frac{1}{2}\mathrm{ln}\left(\frac{\gamma eB_{c2}}{2\pi ^2T^2}\right)+\frac{1}{2}\mathrm{ln}\left(e^{2\tau }v_{F1,\alpha }^2+e^{2\tau }v_{F2,\alpha }^2\right)_\alpha `$ where $`\mathrm{ln}\gamma =0.577`$ is Euler’s constant. Since the pairing interaction is dominant in the $`\sigma `$-band, we can employ two approximations: at first we set $`\lambda _{}=0`$. This corresponds to assuming that superconductivity in the $`\pi `$-band is completely induced by interband coupling. As a second approximation we assume that the vortex lattice distortion $`\tau `$ is also dominated by the $`\sigma `$-band. Then we can find $`\tau `$ by just optimizing it at $`\eta =0`$ and then use it for evaluation at $`\eta >0`$. With these two approximations the logarithmic averages in $`l_\alpha `$ can be performed and we finally get

$$\frac{B_{c2}^{ab}\left(T=0\right)}{B_{c2}^c\left(T=0\right)}=\frac{1.246}{ϵ_c}e^{\eta \left(0.4822\mathrm{ln}ϵ_c\right)}$$
(7)
The parameter $`\tau `$ is found to be $`\tau =\frac{1}{2}\mathrm{ln}ϵ_c`$, which shows that at low temperatures there is no simple relation anymore between $`\tau `$ and the anisotropy ratio. Also we mention here that at low temperatures Eq. (5) does not hold anymore. In Fig. 3 Eq. (7) is shown for $`ϵ_c=0.23`$ as the dotted line. In order to demonstrate the quality of these approximations also the fully numerical result without these approximations is shown as the dashed line. Apparently, now the anisotropy $`\mathrm{\Gamma }`$ is much less sensitive to interband coupling and only varies on a scale given by $`\eta 1/\left(0.4822\mathrm{ln}ϵ_c\right)`$. The reason for this is that $`ϵ_c`$ now only comes in logarithmically.

Eq. (7) shows that at $`T=0`$ the divergence for $`ϵ_c0`$ appears at all $`\eta <0.5`$, while at $`T=T_c`$ this divergence only shows up at $`\eta =0`$. Numerically we observe that at finite $`\eta `$ and $`ϵ_c=0`$ there is a certain temperature at which $`B_{c2}^{ab}`$ diverges. This temperature becomes smaller when $`\eta `$ is increased. In Fig. 4 we show the distortion of the vortex lattice at high magnetic field and low temperature that we expect from our calculation. When the magnetic field is directed along the $`c`$-axis of the crystal a regular Abrikosov lattice is expected as shown in Fig. 4 (left panel). However, when the field is directed within the $`ab`$-plane we expect a distortion of $`e^\tau =\sqrt{ϵ_c}=0.48`$ as shown in Fig. 4 (right panel). This prediction can be checked by neutron scattering or STM tunneling Eskildsen .

To summarize, we have calculated the anisotropy of the upper critical field for the two band Fermi surface topology shown in Fig. 1. Using parameters from bandstructure calculations for MgB2 we find a strong temperature dependence of the upper critical field in agreement with recent measurements on MgB2 single crystals. Fine tuning of the parameters can yield a very good fit of the experimental data. We observe that the small $`c`$-axis dispersion of the $`\sigma `$-band leads to a high sensitivity of the anisotropy ratio on the interband coupling strength near $`T_c`$, but not at $`T=0`$. This suggests that the interplay of these two quantities leads to the strong temperature dependence of the upper critical field in MgB2.

We would like to thank F. Bouquet, O.V. Dolgov, M.R. Eskildsen, K. Maki, A.I. Posazhennikova and L. Tewordt for valuable discussions. Thanks are also due to S. Graser for his help.


format: text

text:
# Why helicity injection causes coronal flux tubes to develop an axially invariant cross-section

## 1
Introduction

Remarkably detailed images of coronal loops provided by the TRACE spacecraft \[Aschwanden et al (2000)\] indicate that these loops have cross-sectional area varying by only 10-20% over their entire length; this surprising behavior cannot be explained by potential or force-free field models \[Klimchuk (2000\]. We present here a model explaining this behavior. Since the loops are typically twisted by a fraction of a turn over their length, they contain magnetic helicity and so our model should be intrinsic to helicity injection processes.

A coronal loop will be represented here by an axisymmetric flux tube with straight axis (cf. top, Fig. 1) and cylindrical coordinates $`(r,\varphi ,z)`$ will be used. The loop is characterized by an initially potential poloidal flux function

$$\psi (r,z)=_0^rB_z(r^{},z)2\pi r^{}𝑑r^{}$$
(1)
with associated poloidal field

$$𝐁_{pol}=\frac{1}{2\pi }\psi \times \varphi .$$
(2)
Axial non-uniformity corresponds to $`\psi `$ being $`z`$-dependent and bulging corresponds to $`\psi ^1^2\psi /z^2>0.`$ Similarly, the poloidal current is written as

$$I(r,z)=_0^rJ_z(r^{},z)2\pi r^{}𝑑r^{}$$
(3)
with associated poloidal current density

$$𝐉_{pol}=\frac{1}{2\pi }I\times \varphi .$$
(4)
We show that establishing a steady-state current $`I`$ involves three sequential stages having distinct physics. The first stage, “ramp-up”, has physics akin to a linear Alfvén wave, but it is assumed that the ramp-up rate is sufficiently slow that the effective Alfvén wavelength is infinite. This means that retarded time effects due to wave propagation issues are negligible and the current ramps up everywhere simultaneously as in an ordinary electrical circuit. The second stage, “axial flow”, has $`I/t=0`$, but is not in MHD equilibrium because unbalanced, non-conservative $`𝐉\times 𝐁`$ forces exist which drive plasma flows. The third stage, “stagnation”, involves convection of magnetic flux by the flows, plasma heating as a result of flow stagnation, and straightening of the $`\psi `$ profile until MHD equilibrium is established.

## 2 First stage (Ramp-up)

We represent the current ramp-up by the time-dependence

$$I(t)=I_0\frac{1+\mathrm{tanh}(t/\tau )}{2}$$
(5)
where the ramp-up time is assumed to be much longer than the time it takes for an Alfvén wave to propagate the length $`h`$ of the flux tube, i.e., $`\tau >>h/v_A.`$ Alfvén wave propagation effects are therefore unimportant in which case the system behaves like an electric circuit. From Ampere’s law the toroidal magnetic field is

$$B_\varphi (r,z,t)=\frac{\mu _0I(r,z,t)}{2\pi r}.$$
(6)
The toroidal component of Faraday’s law is

$$\frac{E_r}{z}\frac{E_z}{r}=\frac{B_\varphi }{t}.$$
(7)
We note that $`B_\varphi `$ has minimal $`z`$ dependence and that Ohm’s law implies $`E_z0`$. Thus, integration of Eq.(7) with respect to $`z`$ gives

$$E_r\frac{\mu z}{2\pi r}\frac{I}{t}$$
(8)
where on the basis of symmetry the location $`z=0`$ is set to be at the axial midpoint of the flux tube. Since there is no axial force in this stage, $`U_z`$remains zero and so the radial component of the ideal Ohm’s law gives

$$U_\varphi =\frac{E_r}{B_z}\frac{\mu z}{2\pi rB_z}\frac{I}{t},$$
(9)
showing that $`U_\varphi `$is finite only when $`I`$ is changing. The change in $`U_\varphi `$ (toroidal acceleration) implies the existence of a radial current determined from the toroidal component of the equation of motion

$$\rho \frac{U_\varphi }{t}=J_rB_z.$$
(10)
This current is just the polarization current \[Chen(1984), p.40\]

$$J_r=\frac{\rho }{B_z}\frac{U_\varphi }{t}=\frac{\rho }{B_z^2}\frac{E_r}{t}=\frac{\mu z\rho }{2\pi rB_z^2}\frac{^2I}{t^2}.$$
(11)
The transient toroidal velocity given by Eq.(9) results in an azimuthal displacement of the plasma,

$$r\mathrm{\Delta }\varphi =_0^tU_\varphi 𝑑t=\frac{\mu _0zI}{2\pi rB_z}=\frac{zB_\varphi }{B_z}$$
(12)
showing that the plasma motion in this stage follows the twisting of the magnetic field (the field line can be thought of as being frozen to the plasma so that when the field line twists, so does the plasma). Thus $`J_r`$ is finite only when $`I`$ is changing and $`J_r`$ is first negative and then positive, corresponding to toroidal acceleration followed by toroidal deceleration. The $`r`$ direction here is really a proxy for the $`\psi `$ direction, since the polarization current is in the $`\psi `$ direction (because the flux tube is long and slender, the $`r`$ direction is nearly the same as the $`\psi `$ direction).

Once $`I`$ has been established, both $`U_\varphi `$ and the polarization current $`J_r`$ remain zero. The poloidal flux function $`\psi `$ at this stage has not changed from its original value — all that has happened is that a toroidal field has been added so that the total magnetic field is now

$$𝐁=\frac{1}{2\pi }\left(\psi \times \varphi +\mu _0I\varphi \right).$$
(13)
The original potential flux tube has become twisted as shown in the bottom sketch of Fig. 1. Since $`\psi `$ is unchanged from its potential value, the bottom sketch in Fig.1 has the same poloidal profile (envelope) as the top sketch. This means that

$$\mu _0J_\varphi =r\varphi \times \frac{\left(\psi \times \varphi \right)}{2\pi }=\frac{r}{2\pi }\left(\frac{1}{r^2}\psi \right)=0.$$
(14)
## 3 Second stage (flow acceleration)

The magnetic force is

$`𝐉\times 𝐁`$ $`=\left(𝐉_{pol}+𝐉_{tor}\right)\times \left(𝐁_{pol}+𝐁_{tor}\right)`$
$`=𝐉_{pol}\times 𝐁_{pol}+𝐉_{pol}\times 𝐁_{tor}+𝐉_{tor}\times 𝐁_{pol}.`$ (15)
The term $`𝐉_{pol}\times 𝐁_{pol}`$ provides acceleration in the toroidal direction and involves the component of $`𝐉_{pol}`$ which is perpendicular to $`𝐁_{pol}.`$ However, in the previous section we showed that the component of $`𝐉_{pol}`$ which is perpendicular to $`𝐁_{pol}`$ is the polarization current and that this current scales as $`^2I/t^2.`$ Thus, when $`I`$ is in steady state there is no polarization current and no toroidal acceleration, and so $`U_\varphi `$ remains zero. It was also shown that $`𝐉_{tor}=J_\varphi \widehat{\varphi }`$ is zero so long as the poloidal flux surfaces are unperturbed from their initial potential shape. We assume in this second stage that (i) the current is constant in which case $`𝐉_{pol}\times 𝐁_{pol}=0`$ and $`U_\varphi =0`$ and (ii) the poloidal flux surfaces are unperturbed from their initial potential shape in which case $`𝐉_{tor}=0.`$ Thus, during this second stage the magnetic force reduces to

$$𝐉\times 𝐁=𝐉_{pol}\times 𝐁_{tor}=\frac{1}{2\pi }\left(\psi \times \varphi \right)\times \frac{\mu _0I}{2\pi }\varphi =\frac{\mu _0}{8\pi ^2r^2}I^2.$$
(16)
The above force is non-conservative (i.e., $`\times \left(𝐉\times 𝐁\right)`$ is non-zero) and so cannot be balanced by a pressure gradient since a pressure gradient is conservative (i.e., $`\times P=0`$). Thus, it is not possible for equilibrium to be attained in this stage. The only way for an equilibrium to be achieved is for the poloidal profile of the magnetic field to change, which is what happens in the third stage, to be discussed later.

The fact that $`𝐉_{pol}\times 𝐁_{pol}=0`$ means that $`\left(I\times \varphi \right)\times \left(\psi \times \varphi \right)=0`$ which in turn implies that $`I`$ is parallel to $`\psi `$ and so $`I`$ must be a function of $`\psi ,`$ i.e., $`I=I(\psi ).`$ Thus, the poloidal current flows along the poloidal flux surfaces. This is consistent with the well-known Hamiltonian dynamics concept that, because of conservation of canonical angular momentum, particles in a toroidally symmetric geometry cannot make an excursion exceeding a poloidal Larmor radius from a poloidal flux surface \[e.g., see p.207-208 of Bellan(2000)\]. In other words, Hamiltonian mechanics forbids the existence of steady current in the direction normal to a poloidal flux surface.

As sketched in Fig. 1, the poloidal flux function is bulged near $`z=0`$, corresponding to a weaker magnetic field near $`z=0`$ than at the ends $`z=\pm h.`$ This would be the situation if the source currents for the poloidal field were located external to the flux tube and so the middle of the flux tube would be further from the source currents than the ends. Since $`I=I(\psi )`$, the current channel would also be bulged.

Equation (16) implies that the $`z`$ component of the equation of motion is

$$\rho \frac{dU_z}{dt}=\left(𝐉\times 𝐁\right)_z=\frac{1}{8\pi ^2r^2}\frac{I^2}{z}=\frac{}{z}\left(\frac{B_\varphi ^2}{2\mu _0}\right).$$
(17)
This means that there is a force accelerating plasma from regions where $`B_\varphi ^2`$ is strong to regions where $`B_\varphi ^2`$ is weak. Since $`I=I(\psi )`$ and since $`\psi `$ is bulged in the middle, $`B_\varphi `$ must be stronger near $`z=\pm h`$ where the current channel diameter is small than at $`z=0`$ where the current channel diameter is large. There consequently must be an acceleration of plasma from both ends (i.e., $`z=\pm h`$) towards the middle (i.e., $`z=0`$) as shown in Fig. 2. The convergent axial pumping is similar to the “sweeping magnetic twist mechanism” discussed by Uchida and Shibata (1988), but it should not be considered a wave because it involves actual convection of material and not propagation of energy through a material. We note in passing that there could be a few exceptional particles collisionally bouncing back and forth between the approaching fluid elements. These exceptional particles would be accelerated to very high energy by the Fermi acceleration process, and so one would expect to see a tail of energetic particles develop in the vicinity of $`z=0.`$ The Fermi process would thus predict that the most energetic particles would be located around the top of an arched coronal loop and such is indeed what is observed (Feldman, 2002).

## 4 Third stage (stagnation, heating, and straightening)

The flows from both ends converge at the middle and must come to a halt at $`z=0.`$ Convergence of flows means that $`𝐔`$is negative. This has important implications for the magnetic field as can be seen by considering the induction equation toroidal component,

$$\frac{B_\varphi }{t}=r𝐁_{pol}\left(\frac{U_\varphi }{r}\right)r𝐔_{pol}\left(\frac{B_\varphi }{r}\right)B_\varphi 𝐔_{pol}.$$
(18)
We have shown that (i) $`U_\varphi =0,`$ (ii) $`𝐔_{pol}0`$ at the stagnation layer, and (iii) $`𝐔_{pol}`$ is negative. Thus, in the vicinity of the stagnation layer

$$\frac{B_\varphi }{t}=B_\varphi 𝐔_{pol}$$
(19)
showing that $`B_\varphi `$ must increase at the stagnation layer (increase of magnetic field at regions of local flow convergence has been discussed in a more general context by Polygiannakis and Moussas (1999)).

The continuity equation in the vicinity of the stagnation layer gives $`𝐔_{pol}=\rho ^1\rho /t`$ and so in the vicinity of the stagnation layer

$$\frac{1}{B_\varphi }\frac{B_\varphi }{t}=\frac{1}{\rho }\frac{\rho }{t}.$$
(20)
Thus, in the vicinity of the stagnation layer $`B_\varphi `$ increases in proportion to the accumulation of mass. Since $`I`$ is constant during this stage, if $`B_\varphi `$ increases, the radius of the current channel must decrease so as to maintain $`\mu _0I=2\pi rB_\varphi `$ constant. However, since $`I=I(\psi )`$, if the radius of the current channel decreases, then the radius of $`\psi `$ must also decrease, thereby reducing the bulging. Ultimately, the bulge becomes vanishingly small as more and more mass accumulates at the stagnation point and eventually a plasma-loaded, axially uniform flux tube results.

Since stagnation involves conversion of flow velocity into thermal energy, the plasma will be heated at the stagnation layer. The resulting plasma temperature can be estimated from the details of the stagnation process. The axial magnetic force in Eq.(17) involves the gradient of the toroidal field energy density and so the kinetic energy gained by the plasma will be of the order of the change in toroidal field energy density $`B_\varphi ^2/2\mu _0`$ between the ends ($`z=\pm h`$) and the midpoint $`z=0.`$ Because $`B_\varphi ^2/2\mu _0`$ is larger at the ends, the plasma flow kinetic energy is $`\rho U_z^2/2`$ $`\left[B_\varphi ^2/2\mu _0\right]_{z=\pm h}`$ and at the stagnation layer $`z=0`$ this flow energy is converted into heat so $`\rho U_z^2=\left[B_\varphi ^2/\mu _0\right]_{z=\pm h}`$ $`nm_ev_{Te}^2+nm_iv_{Ti}^2.`$

Since the poloidal field is much stronger than the toroidal field and is approximately $`B_z,`$ this means that the flow stagnation causes the plasma to develop a state where

$$\beta =\frac{2\mu _0(nm_ev_{Te}^2+nm_iv_{Ti}^2)}{B_z^2}=2\frac{B_\varphi ^2}{B_z^2}=2\left(\frac{\mu _0I}{2\pi aB_z}\right)^2=\left(\frac{\mu _0I}{\psi }\right)^2\frac{a^2}{2},$$
(21)
where $`a`$ is the radius of the current channel. However, $`\mu _0I=𝑑𝐬\times 𝐁`$ and $`\psi =𝑑𝐬𝐁`$ where the surface integral is over the cross-sectional area of the flux tube. If we define the ratio of poloidal current to poloidal flux as

$$\alpha =\mu _0I/\psi ,$$
(22)
the $`\beta `$ predicted from flow stagnation is

$$\beta =\alpha ^2a^2/2.$$
(23)
Thus, when a current is made to flow along an initially bulging current-free current channel, the current channel will twist up (helicity injection), plasma will be ingested from both ends, accelerated toward the middle where it accumulates and heats up at a stagnation layer. The equilibrium will become straight (filamentary) and have $`\beta =\alpha ^2a^2/2`$ where $`\alpha =\mu _0I/\psi `$ and $`a`$ is the radius of the current channel. To an outside observer the current channel will look field-aligned since the current is axially uniform and appears to be embedded in an axially uniform axial magnetic field. However, oblivious to the outside observer, within the current channel there is a Bennett pinch-like radial force balance between plasma pressure pushing out and $`J_zB_\varphi `$ magnetic force pushing in.

Although the predicted $`\beta `$ is typically very small, its effect is crucial. To see this, consider that in equilibrium $`𝐉\times 𝐁=P`$ so that $`𝐁P=(2\pi )^1\left(\psi \times \varphi +\mu _0I\varphi \right)P=0`$ and so $`P`$ must be a function of $`\psi ,`$ i.e., $`P=P(\psi ).`$ Defining $`\psi _0`$ as the flux on the flux surface where $`P`$ vanishes, we can write $`P(r,z)=(1\psi (r,z)/\psi _0)P_0`$ where $`P_0`$ is the on-axis pressure (i.e., where $`\psi =0`$). We can also write $`\mu _0I(r,z)=`$ $`\alpha \psi (r,z)`$ and so $`𝐉\times 𝐁=P`$ can be written in Grad-Shafranov form as

$$r\frac{}{r}\left(\frac{1}{r}\frac{\psi }{r}\right)+\frac{^2\psi }{z^2}+\alpha ^2\psi =4\pi ^2\mu _0r^2\frac{P_0}{\psi _0}.$$
(24)
If $`\alpha ^2=4\pi ^2\mu _0a^2P_0/\psi _0^2`$ where $`a`$ is the flux tube radius at $`z=0,`$ then the only solution to Eq.(24) satisfying the specified boundary condition that $`P`$ vanishes when $`\psi =\psi _0`$ is the particular solution $`\psi (r,z)=\psi _0r^2/a^2`$. This means that the flux tube must be axially uniform when $`\alpha ^2=4\pi ^2\mu _0a^2P_0/\psi _0^2.`$ Defining $`B_0=\psi _0/\pi a^2`$ as the axial field at $`z=0,`$ it is seen that this condition for axial uniformity corresponds to $`\alpha ^2=4\pi ^2\mu _0a^2P_0/B_0^2\pi ^2a^4`$ or $`\alpha ^2a^2/2=\beta `$ where $`\beta =2\mu _0P_0/B_0^2.`$ This equilibrium has $`J_\varphi =0,`$ so all confinement is provided by the Bennett pinch force $`J_zB_\varphi .`$ The current is purely in the $`z`$ direction, but the magnetic field is helical.

The situation of small but finite $`\beta `$ is substantially different from the case of zero $`\beta `$ because the system is constrained to be axially uniform if and only if $`\beta =\alpha ^2a^2/2.`$ The arguments presented in the discussion of Eqs.(21-23) show that the MHD dynamical pumping tends to produce precisely the situation where $`\beta =\alpha ^2a^2/2,`$ and so it is predicted that MHD dynamical pumping will always cause configurations to tend towards being axially uniform (i.e., filamentary), hot, and dense, and with $`\beta =\alpha ^2a^2/2.`$

The definition of $`\alpha `$ in Eq.(22) is closely related to that used for force-free fields. However, there is an important difference because while Eq.(22) corresponds to having $`\mu _oJ_z=\alpha B_z`$, Eq.(22) makes no statement about any relationship between $`J_\varphi `$ and $`B_\varphi `$. The finite $`\beta `$ equilibrium discussed in the previous two paragraphs is not force-free and involves the radial force balance $`J_zB_\varphi =P/r`$with $`J_\varphi =0;`$this differs from the force-free radial equation $`J_\varphi B_zJ_zB_\varphi =0`$ with $`\mu _oJ_\varphi =\alpha B_\varphi ,`$ $`\mu _oJ_z=\alpha B_z.`$ It is worth noting that the determination of $`\alpha `$ made from vector magnetographs (e.g., Pevtsov et al., 1997) effectively use the definition $`\alpha =\mu _oJ_z/B_z`$ which is equivalent to Eq.(22); these measurements do not provide information on either $`J_\varphi `$ or $`B_\varphi `$ and so do not provide any information on the value of $`\mu _oJ_\varphi /B_\varphi .`$ Thus there is only one definition for $`\alpha ,`$ but its application is different for force-free situations compared to finite-$`\beta `$ situations: for force-free situations $`\alpha `$ gives the ratio of current to flux for both toroidal and poloidal directions whereas for the finite $`\beta `$ situation, $`\alpha `$ refers only to the the ratio of poloidal current to poloidal field.

The prediction that $`\beta =\alpha ^2a^2/2`$ can be compared with the actual observed values of $`\beta `$ in TRACE flux loops. To calculate the predicted $`\beta ,`$ we use the nominal measured flux loop radius $`a=1.6\times 10^6`$ m from Aschwanden et al.(2000) and the nominal measured active region $`\alpha =2\times 10^8`$ m-1 from Fig. 4 of Pevtsov et al. (1997). These parameters give a nominal $`\beta _{predicted}=\alpha ^2a^2/2=5\times 10^4.`$ The observed value $`\beta _{observed}`$ is calculated using a nominal measured density $`n=10^{15}`$m$`^3`$ and a nominal measured temperature $`10^6`$ K \[Aschwanden et al. (2000)\]. In addition a nominal axial magnetic field $`B_z=1.5\times 10^2`$ T is assumed based on the argument that because the flux tube is axially uniform, its axial field must also be axially uniform and so will have the same value as the nominal $`B_z=1.5\times 10^2`$ at the surface of an active region. These parameters give $`\beta _{observed}=2\mu _0n\kappa T/B_z^2=4\times 10^4`$ which is very close to $`\beta _{predicted}`$. If the model were wrong, one would expect no relationship between the predicted and observed $`\beta `$’s, i.e., one would expect a discrepancy of many orders of magnitude between the predicted and the observed $`\beta `$’s.

This model also has implications regarding the brightening typically observed when the axis of a coronal loop starts to writhe and the loop develops a kink instability (sigmoid). Since kink instability occurs when $`\alpha h1`$ and for a long thin flux tube $`a<<h`$, this model predicts that $`\beta =\alpha ^2a^2/2<<`$ $`\alpha ^2h^2/2`$ will still be small even if $`\alpha `$ is increased to the point where $`\alpha h1`$ and kink instability occurs. However, $`\beta `$ will increase as $`\alpha `$ increases and so this model predicts that the loop should brighten in proportion to the writhing of its axis (i.e., in proportion to $`\alpha `$ as $`\alpha h`$ approaches unity).

Finally, we note that Feldman (2002) has recently used purely observational evidence to argue that electric currents with geometry similar to what is discussed here are the means by which the Sun and similar stars produce their coronal activity.

Acknowledgment: Supported by USDOE Grant DE-FG03-97ER54438.

REFERENCES

Aschwanden, M. J., Nightingale, R. W., and Alexander, D. 2000, ApJ 541, 1059

Bellan, P. M., Spheromaks (Imperial College Press, 2000, London)

Chen, F. F., Introduction to Plasma Physics and Controlled Fusion (Plenum 1984, New York)

Feldman, U. 2002, Physica Scripta 65, 1985

Klimchuk, J. A. 2000, Solar Physics 193, 53

Pevtsov, A. A., Canfield, R. C., and McClymont, A. N. 1997, ApJ 481, 973

Polygiannakis, J. M. and Moussas, X. 1999, Plasma Phys. Control. Fusion 41, 967

Uchida Y. and Shibata, K. 1988, Solar Phys. 116, 291


text:
# Budd-Vannimenus theorem for superconductors

\[

## Abstract

The Budd-Vannimenus theorem is modified to the surface of superconductors in the Meissner state. This identity links the surface value of the electrostatic potential to the free energy in the bulk which allows one to evaluate the observed potential without the explicit solution of the charge profile at the surface.

\]

The Budd-Vannimenus theorem  links the value of the electrostatic potential at the surface of a metal to the energy per electron in the bulk. This simple consequence of the Feynman-Hellmann theorem turned out to be very useful in studies of metal surfaces  as it enables one to circumvent demanding calculations of the charge distribution at the surface. We will use this approach to evaluate the surface potential of superconductors within the Ginzburg-Landau (GL) theory.

Perhaps it is useful to start from a historical review. Within the pre-London approaches based on the ideal charged liquid , the electrostatic potential was assumed to balance the Lorentz and the inertial forces acting on currents. Since the current $`j=env`$ flows along the surface and its amplitude falls off exponentially from the surface into the bulk, the electrostatic Bernoulli–type potential, $`e\phi =\frac{1}{2}mv^2`$, falls off on the scale of the London penetration depth. From the Poisson equation, $`\rho =ϵ_0^2\phi `$, one finds that this corresponds to a charge accumulated in the layer penetrated by the magnetic field. Charge neutrality is maintained by the opposite surface charge. This picture has been confirmed by the London theory .

The electrostatic potential equals the Bernoulli potential only at zero temperature when all electrons are in the superconducting condensate. At finite temperatures, a part $`n_n`$ of electrons remains in the normal state, while the rest $`n_s=nn_n`$ contributes to the supercurrent, $`j=en_sv`$. As a consequence, the electrostatic potential reduces to $`e\phi =\frac{n_s}{n}\frac{1}{2}mv^2`$ as shown by van Vijfeijken and Staas . Jakeman and Pike  recovered this result from the time-dependent GL theory. Moreover, their approach describes the surface charge as a space charge localized on the Thomas-Fermi screening length.

The surface charge and the charge accumulated at the London penetration depth form a surface dipole due to which the electrostatic potential in the bulk differs from the potential out of the superconductor. This voltage difference, which like the resistance of the Hall effect is perpendicular to the current, is not accessible by contact measurements but it was observed by a contact-less capacitive pickup .

The first measurements  were not sufficiently accurate and created expectations that the surface voltage includes traces of the pairing mechanism. Rickayzen  used the thermodynamic approach and showed that the potential indeed has a pairing contribution, $`e\phi =\frac{n_s}{n}\frac{1}{2}mv^24\frac{n_n}{n}\frac{\mathrm{ln}T_c}{\mathrm{ln}n}\frac{1}{2}mv^2`$. The first term results from the Lorentz and inertial forces, while the second one reflects the pairing. The pairing term dominates close to the critical temperature $`T_c`$. If $`\phi `$ would be experimentally accessible, one could deduce from it the density dependence of $`T_c`$. The highly accurate measurement by Morris and Brown , however, chilled all the expectations. The electrostatic potential revealed no pairing contribution.

The disagreement between the theory and experiment remained unexplained for a long time and the question of the charge transfer in superconductors was left aside till the discovery of the high-$`T_c`$ materials. It was predicted that in layered materials the superconducting transition induces a charge transfer from Cu02 planes to charge reservoirs. This transfer caused merely by the pairing mechanism has been experimentally confirmed by the positron annihilation, the $`x`$-ray absorption spectroscopy, and the nuclear magnetic resonance.

Apparently, there are two groups of contradictory experimental results. The pairing contribution is absent in the surface potential but it is observed in the bulk. As it was indicated recently , there is an additional surface dipole which cancels the pairing contribution to surface potential seen in capacitive measurements.

The approach presented in Ref.  employs the Budd-Vannimenus theorem, however, it has the same limitation as Rickayzen’s formula – the magnetic field is described within the London’s theory and it has to be weak not to perturb the density of the condensate. In contrast, the measurement of Morris and Brown  explores the entire range of the magnetic fields from weak fields up to the critical value. They found that the potential equals the magnetic pressure $`\phi =\frac{1}{\rho _{\mathrm{lat}}}\frac{1}{2\mu _0}B^2`$, i.e., it is independent of the temperature and the used material enters merely via the charge density $`\rho _{\mathrm{lat}}`$ of its lattice, for both type-I and II superconductors.

Here we modify the Budd-Vannimenus theorem so that it is consistent with the GL theory. Its validity is not restricted to weak magnetic fields. We will show that the surface potential is given by

$$\phi _0=\frac{f_{\mathrm{el}}}{\rho _{\mathrm{lat}}},$$
(1)
where $`f_{\mathrm{el}}`$ is the surface value of the GL free energy (3). With the help of the modified Budd-Vannimenus theorem one can explain the experimental result of Morris and Brown under very general conditions.

The recent implementation of the GL theory has changed the picture of the surface charge. It is not localized on the scale of the Thomas-Fermi screening length as claimed by Jakeman and Pike, but it varies on scales of the GL coherence length . As pointed out in Ref. , if the electrostatic potential is a local functional of the GL wave function, $`\phi [\psi ]`$, the GL boundary condition, $`\psi =0`$ in the direction normal to the surface, implies the zero electric field at the surface, $`E=\phi =\frac{\phi }{\psi }\psi =0`$. Accordingly, no additional surface charge is needed to maintain the charge neutrality.

Since the surface charge is absent, it seems that the surface dipole on the Thomas-Fermi screening length is absent, too. There has to be yet another surface dipole, however. The pairing correlation is weaker on the surface than in the bulk, what results in forces pulling the Cooper pairs inside. Such forces are always balanced by the electrostatic field. The understanding of this effect will require microscopic studies which are not yet feasible. We can merely speculate that the surface dipole is somehow linked to the gap profile at the surface , i.e., it is localized on the BCS coherence length.

To introduce the surface dipole on an intuitive level we first assume that the system is close to the critical temperature. In this regime, the London penetration depth $`\lambda `$ and the GL coherence length $`\xi `$ are much larger than the BCS coherence length $`\xi _0`$. We can then define an intermediate scale $`L`$ such that $`\xi _0L\xi ,\lambda `$, as sketched in Fig. 1. On the scale $`L`$, the GL wave function, the vector potential, and the electrostatic potential change only negligibly, e.g., $`\phi (x)\phi (x0)\phi (0)`$ for $`\xi _0<x<L`$. The extrapolation of the bulk potential towards the surface, $`\phi (0)`$, has to be distinguished from the true surface potential $`\phi _0`$. The difference $`\phi _0\phi (0)`$ is caused by the surface dipole we aim to evaluate.

Close to the critical temperature we can follow the idea of Budd and Vannimenus . Let us assume a virtual compression of the crystal lattice such that the lattice charge density is removed from the surface layer of an infinitesimal width $`\delta L`$. The perturbation of the lattice charge density in the infinitesimal layer $`0<x<\delta L`$ is $`\delta \rho _{\mathrm{lat}}=\rho _{\mathrm{lat}}`$. The compression leads to an increase of the charge density in the layer $`\delta L<x<L`$, where $`\delta \rho _{\mathrm{lat}}=\rho _{\mathrm{lat}}\delta L/L`$, in accordance with the charge conservation.

Let us remind the basic idea of the Feynman-Hellmann theorem. The lattice charge enters the jellium model of metals as an external parameter. If one changes this external parameter, the situation corresponds to doing work on the system, $`\delta 𝒲=S𝑑x\delta \rho _{\mathrm{lat}}\frac{f}{\rho _{\mathrm{lat}}}=S𝑑x\delta \rho _{\mathrm{lat}}\phi `$, where $`f`$ is the density of the free energy including the electrostatic interaction, and $`S`$ is the sample area. According to the Feynman-Hellmann theorem, the change of the electrostatic potential does not contribute to the work up to the first order in $`\delta \rho _{\mathrm{lat}}`$. Now we can proceed with the algebra. We split the integral into three parts. Since $`\delta L`$ is an infinitesimal displacement, the potential in the layer $`0<x<\delta L`$ can be replaced by the surface value $`\phi _0`$. The surface region $`\delta L<x<\xi _0`$ gives a negligible contribution of the order of $`\xi _0/L`$. In the remaining bulk region $`\xi _0<x<L`$, the electrostatic potential is nearly constant and equals $`\phi (0)`$. The work thus reads $`\delta W=S\delta L\rho _{\mathrm{lat}}(\phi (0)\phi _0)`$.

In the same time, the work increases the free energy $``$ of the system, $`\delta 𝒲=\delta `$. We denote the spatial density of the electronic free energy by $`f_{\mathrm{el}}`$. Lower case letters denote free energies per unit volume, i.e., their densities, while the calligraphic upper case letters denote the corresponding total values. The total free energy changes as $`\delta =\frac{(f_{\mathrm{el}}SL)}{(SL)}S\delta L=\left(f_{\mathrm{el}}+n\frac{f_{\mathrm{el}}}{n}\right)S\delta L`$. The first term results from the reduced volume, $`SLS(L\delta L)`$, and the second one from the corresponding increase of the electron density by $`nn(1+\delta L/L)`$. We thus obtain a modification of the Budd-Vannimenus theorem,

$`\rho _{\mathrm{lat}}\left(\phi _0\phi (0)\right)`$ $`=`$ $`f_{\mathrm{el}}n{\displaystyle \frac{f_{\mathrm{el}}}{n}}.`$ (2)
This relation is our final result. It describes the step of the potential at the surface due to the surface dipole in terms of the free energy.

Formula (2) differs from the original Budd-Vannimenus theorem in three points. First, in the original Budd-Vannimenus theorem, the surface potential is related to the potential $`\phi _{\mathrm{}}`$ deep in the bulk. In (2), the extrapolation of the internal potential towards the surface, $`\phi (0)`$, appears instead. Second, in order to cover systems at finite temperature, we use the free energy instead of the ground-state energy. Third, within the original Budd-Vannimenus approach, the electron charge density and the lattice charge density have to be equal because of the charge neutrality. In our approach, the density of electronic charge differs locally from the lattice charge density, $`en\rho _{\mathrm{lat}}`$, due to the charge transfer on the scales $`\xi `$ and $`\lambda `$.

Let us demonstrate how the relation (2) can be used within the GL theory. To this end we introduce the GL free energy

$$f_{\mathrm{el}}=\frac{1}{2m^{}}\left|\left(i\mathrm{}e^{}A\right)\psi \right|^2+\alpha |\psi |^2+\beta |\psi |^4,$$
(3)
where $`\psi `$ is the GL wave function, $`A`$ is the vector potential and $`m^{}=2m`$ and $`e^{}=2e`$ are the mass and the charge of the Cooper pair. The GL parameters $`\alpha `$ and $`\beta `$ depend on the temperature and the electron density $`n=n_n+2|\psi |^2`$. Finally, we add the electromagnetic energy

$$f=f_{\mathrm{el}}+\phi (\rho _{\mathrm{lat}}+en)\frac{ϵ_0}{2}E^2+\frac{1}{2\mu _0}B^2,$$
(4)
with the magnetic field $`B=\times A`$ and the electric field $`E=\phi `$.

Variations of the free energy with respect to its independent variables $`A,\phi ,\psi ,n_n`$ yield the equations of motion in Lagrange’s form

$$\frac{f}{\nu }+\frac{f}{\nu }=0.$$
(5)
For $`\nu =A`$ the variational condition (5) yields the Ampere-Maxwell equation, for $`\nu =\phi `$ the Poisson equation, for $`\nu =\psi `$ the GL equation, and for $`\nu =n_n`$ the condition of zero dissipation,

$$e\phi =\frac{f_{\mathrm{el}}}{n}.$$
(6)
This condition allows one to evaluate the electrostatic potential in the bulk of the superconductor . Of course, one can add any constant to the electrostatic potential. The actual choice (6) simplifies the relations we deal with below. Another convenient choice would be to subtract the value of the potential in the non-magnetic state, i.e., deep inside the superconductor, so that $`\phi `$ approaches zero in the bulk.

Formula (6) does not cover the surface dipole on the scale $`\xi _0`$, therefore at the surface it provides the extrapolated bulk value $`\phi (0)`$. We can thus use (6) to rearrange the Budd-Vannimenus theorem (2) as

$$\rho _{\mathrm{lat}}\phi _0=f_{\mathrm{el}}+\phi (0)\left(\rho _{\mathrm{lat}}+en\right).$$
(7)
Now all terms on the right hand side are explicit quantities which one obtains within the GL theory extended by the electrostatic interaction .

In customary GL treatments, the electrostatic potential and the corresponding charge transfer are omitted. This approximation works very well for magnetic properties since the relative charge deviation, $`\left(\rho _{\mathrm{lat}}+en\right)/\rho _{\mathrm{lat}}`$, is typically of the order of $`10^8`$ to $`10^4`$ leading to comparably small corrections in the GL equation. With the same accuracy one obtains the electronic free energy $`f_{\mathrm{el}}`$. Therefore it is possible to evaluate the surface potential using (1), which follows from (7) if terms proportional to $`\left(\rho _{\mathrm{lat}}+en\right)/\rho _{\mathrm{lat}}`$ are neglected.

So far we have restricted ourselves to systems close to the critical temperature when the validity conditions of the GL theory are well satisfied. In many cases, the GL theory is used beyond the limits of its nominal applicability, however. In these cases, the GL coherence length $`\xi `$ and/or the London penetration depth $`\lambda `$ are comparable to or even shorter than the BCS coherence length $`\xi _0`$. Although the above derivation does not apply to this case, it is possible to show that the Budd-Vannimenus theorem (7) provides physically consistent predictions in these cases, too.

We evaluate now the voltage which develops across a superconducting slab with the magnetic field parallel to the slab. For the slab geometry, the GL equation has an integral of motion, see Bardeen . This integral can be obtained quite generally by the Legendre transformation of the free energy, $`g=f_\nu (f/\nu )\nu `$. Indeed, if the fields $`\nu `$ obey the equations of motion (5), the gradient $`g=0`$ vanishes, i.e. $`g=\mathrm{const}`$. With the help of $`_\nu (f/\nu )\nu =(\mathrm{}^2/m^{})|\psi |^2ϵ_0E^2+B^2/\mu _0`$, we can express the electronic free energy as

$$f_{\mathrm{el}}=g+\frac{\mathrm{}^2}{m^{}}|\psi |^2\frac{ϵ_0E^2}{2}+\frac{B^2}{2\mu _0}\phi (\rho _{\mathrm{lat}}+en).$$
(8)
At the surface, the GL boundary condition demands that $`\psi =0`$ what implies $`E=0`$. The free energy on the surface thus reads

$$f_{\mathrm{el}}=g+\frac{B^2}{2\mu _0}\phi (0)(\rho _{\mathrm{lat}}+en).$$
(9)
From the Budd-Vannimenus theorem (7) then follows $`\rho _{\mathrm{lat}}\phi _0=g+B^2/(2\mu _0)`$, and since $`g=\mathrm{const}`$, we obtain the potential difference across the surface as the difference of left and right magnetic pressures

$$\rho _{\mathrm{lat}}(\phi _0^L\phi _0^R)=\frac{B_L^2B_R^2}{2\mu _0}.$$
(10)
Exactly this value of the voltage has been observed by Morris and Brown .

As noticed already by Bok and Klein , there is a simple argument for the formula (10). The left hand side represents the electrostatic force (per unitary area) on the lattice, $`F_{\mathrm{elst}}=_L^R𝑑xE\rho _{\mathrm{lat}}=\rho _{\mathrm{lat}}(\phi _0^L\phi _0^R)`$. The right hand side is the Lorentz force $`F_{\mathrm{Lor}}=BJ`$ with the mean magnetic field $`B=\frac{1}{2}(B_L+B_R)`$ and the net current $`J=_L^R𝑑xj`$ given by Ampere’s rule, $`B_LB_R=\mu _0J`$. Since the electrostatic field provides the only mechanism by which the force is passed from the electrons to the lattice, the two forces have to be equal, $`F_{\mathrm{Lor}}=F_{\mathrm{elst}}`$, which yields (10).

Turning the argument around, one can view the derivation of the Lorentz force as an alternative proof of formula (7). For the slab geometry, this is more general than the proof based on the Budd-Vannimenus approach, as the approach via the Lorentz force does not rely on the specific order of the characteristic lengths $`\xi _0\xi ,\lambda `$ needed above to define the intermediate length $`L`$. Moreover, the Legendre transformation provides the integral of motion for any modification of the GL free energy which is of the form $`f_{\mathrm{el}}[\psi ,\psi ]`$.

On the other hand, within the Budd-Vannimenus approach the surface dipole is treated as a property of the superconducting condensate, what encourages us to hope that formula (7) or its approximation (1) can be used to obtain the surface potential also for cases when the magnetic field has a component perpendicular to the surface. In particular, we expect that it will be applicable also to the superconductors in the mixed state, especially to evaluate the electric field generated by vortices penetrating the surface .

In conclusion, the Budd-Vannimenus theorem was modified to the surface of a superconductor. It allows one to evaluate the electrostatic potential of the surface from the free energy and the bulk electrostatic potential near the surface. Within approximation (1) one obtains the surface potential from the free energy without the actual knowledge of the bulk potential.

This work was supported by MŠMT program Kontakt ME601 and GAČR 202000643, GAAV A1010312 grants. The European ESF program VORTEX is also acknowledged.


text:
# Constraint on intergalactic dust from thermal history of intergalactic medium

## 1 Introduction

There are metals in the intergalactic medium (IGM) even in low density regions such as Lyman $`\alpha `$ forest (e.g., Cowie et al. 1995; Telfer et al. 2002). Since metal and dust are relating each other, it is sure that dust grains also exist in the IGM. However, the amount of the intergalactic (IG) dust is still quite uncertain although it seems to be not so abundant. We try to put a new constraint on the amount of the IG dust by using the IGM thermal history suggested by the recent observations of Lyman $`\alpha `$ forest (e.g., Schaye et al. 2000).

After the early attempts to estimate the IG extinction (Eigenson, 1949; Humason, Mayall, & Sandage, 1956), Crane & Hoffman (1973) have obtained an upper limit on the amount of the IG dust in terms of the density parameter as $`\mathrm{\Omega }_{\mathrm{dust}}^{\mathrm{IGM}}(z=0)10^4`$, by comparing the observed spectral energy distribution of distant (but $`z0.5`$) giant elliptical galaxies with an average spectrum of similar nearby galaxies (see also Nickerson & Partridge 1971; Takase 1972).

Observations of the redshift evolution of the QSOs spectral slope provide us with another approach to the IG dust (Wright, 1981; Cheng, Gaskell, & Koratkar, 1991). Although Cheng et al. (1991) found no evidence for an appreciable IG reddening up to $`z2`$–3, this may not be so strict constraint because it is difficult to measure the UV slope of QSOs enough precisely owing to the broad and complex iron band superposing on the continua. Recently, we obtain further constraint on the IG dust from observations of high-z supernovae (SNe) (Riess et al., 1998; Perlmutter et al., 1999). An upper limit of the colour excess by the IG reddening is $`E_{BV}_{z0.5}E_{BV}_{z0.05}0.03`$ mag.

By the way, Aguirre (1999) suggests that the IG dust can be “gray” by a selection rule in the transfer of grains from the host galaxies to the IGM. This occurs when small ($`0.1\mathrm{\mu m}`$) grains are destroyed selectively by the thermal sputtering in the hot gas halo of the host galaxies before they reach the IGM. Then, the extinction property becomes “gray” since the relatively large dust survives. Importantly, such “gray” IG dust can avoid the detection by the IG reddening survey like above.

Even if the IG dust is really “gray”, its evidence should be imprinted in the cosmic microwave background (CMB) and infrared background because the dust emits thermal radiation in the wave-band from the far-infrared to submillimetre (submm) (Rowan-Robinson, Negroponte, & Silk, 1979; Wright, 1981). According to the current status, COBE data provides us with a rough upper limit on the IG dust (Loeb & Haiman, 1997; Ferrara et al., 1999; Aguirre & Haiman, 2000). On the other hand, the submm background radiation may give a more strict constraint on the IG dust. The emission from the IG dust expected from the modern cosmic star formation history with the grain transfer mechanism of Aguirre (1999) contributes to a substantial fraction ($``$ 75%) of the measured background radiation in the submm (Aguirre & Haiman, 2000). It is interesting to develop the approach of Aguirre & Haiman (2000) to match a strong constraint by SCUBA whose result of the number count accounts $``$ 90% of the submm background light (Calzetti 2001 and references therein).

In any current status, the amount of the IG dust is still highly uncertain. Hence, it is worth trying to find a new constraint on the amount of the IG dust. In this Letter, it is demonstrated that we can obtain such a new constraint of the IG dust by comparing the observational IGM temperature with a theoretical one if we take into account the dust photoelectric heating in the theoretical model. This is possible because the dust heating is an efficient mechanism in the IGM (Nath, Sethi, & Shchekinov, 1999). Interestingly, the effect of the dust photoelectric heating depends on not only the amount but also the size of dust. Therefore, we can put a limit on the amount of the IG dust as a function of the typical size of the IG grains.

## 2 Photoelectric Heating by Grains

To assess photoelectric effect, we must specify the charge on grain, $`Z_\mathrm{d}`$ (in the electron charge unit), which is given by (Spitzer, 1941):

$$\underset{i}{}R_i+R_{\mathrm{pe}}=0,$$
(1)
where we have assumed a charge equilibrium because a typical charging time-scale is very short ($`10[a/0.1\mathrm{\mu m}]^1`$ yr under the UV background radiation dominated by QSO light given in §3). In the equation, $`R_i`$ is the collisional charging rate by $`i`$-th charged particle, and $`R_{\mathrm{pe}}`$ is the photoelectric charging rate. We assume spherical grains for simplicity.

The collisional charging rate is

$$R_i=\pi a^2Z_is_in_i_{v_{\mathrm{min}}}^{\mathrm{}}\sigma _i(v_i,Z_\mathrm{d},Z_i)v_if(v_i)𝑑v_i,$$
(2)
where $`Z_i`$ is the charge in the electron charge unit, $`s_i`$ is the sticking coefficient, $`n_i`$ is the number density, $`a`$ is the grain radius, $`v_i`$ is the velocity, $`v_{\mathrm{min}}`$ is the minimum velocity required to collide with a grain, $`\sigma _i`$ is the dimensionless collisional cross section depending on both charges and $`v_i`$, and $`f(v_i)`$ is the velocity distribution function. We simply assume $`s_i`$ is always unity. In the above equation, we neglected the effect of the secondary emission.

We can neglect the effect of the “image potential” (Draine & Sutin, 1987) on the collisional cross section because the obtained charges are enough large. If we assume the Maxwellian velocity distribution for the particle, the integral in the r.h.s. of equation (2) is reduced to $`(8k_\mathrm{B}T/\pi m_i)^{1/2}g(x)`$ and $`g(x)=1x`$ for $`Z_\mathrm{d}Z_i0`$ or $`g(x)=\mathrm{exp}(x)`$ for $`Z_\mathrm{d}Z_i>0`$, where $`k_\mathrm{B}`$ is the Boltzmann’s constant, $`T`$ is the gas temperature, $`m_i`$ is the charged particle’s mass, and $`x=e^2Z_\mathrm{d}Z_i/ak_\mathrm{B}T`$.

The photoelectric charging rate is given by

$$R_{\mathrm{pe}}=\pi a^2_{\nu _{\mathrm{min}}}^{\nu _{\mathrm{max}}}Q(a,\nu )Y(a,\nu ,Z_\mathrm{d})\frac{4\pi J_\nu }{h\nu }𝑑\nu ,$$
(3)
where $`Q`$ is the absorption coefficient of grains, $`Y`$ is the photoelectric yield, $`J_\nu `$ is the intensity (averaged for solid angle) of the incident radiation at a frequency $`\nu `$, $`h`$ is the Plank constant, and $`\nu _{\mathrm{max}}`$ is the maximum frequency of the incident radiation. The minimum frequency, $`\nu _{\mathrm{min}}`$, is the threshold photon frequency required for an electron to escape from a grain. That is, $`h\nu _{\mathrm{min}}`$ is equal to the ionization potential, $`IP`$, of a grain: $`IP=W+(Z_\mathrm{d}+1/2)e^2/a`$, where $`W`$ is the work function. Here we neglect quantum effects on $`IP`$ because it is very small for grains with $`a10`$ Å (Weingartner & Draine 2001a, hereafter WD01).

According to WD01, we adopt $`W=4.4`$ eV for graphites and $`W=8.0`$ eV for silicates. For $`Q`$, we adopt the values of “graphite” and “smoothed UV astronomical silicate” calculated by Draine & Lee (1984); Laor & Draine (1993); Weingartner & Draine (2001b). The value of $`Y`$ is estimated, based on the way constructed by WD01, which includes approximately the effect of the energy distribution of photoelectrons and the geometrical enhancement for small grains, and reproduces recent results of the laboratory experiments even for very small grains.

The both equilibrium charges obtained for graphites and silicates are similar when their radii are equal. This is consistent with the result of Nath et al. (1999). In the following discussion, we consider mainly silicate case.

The photoelectric heating rate by a grain with radius $`a`$ is expressed by

$$\gamma (a)=\pi a^2_{\nu _{\mathrm{min}}}^{\nu _{\mathrm{max}}}E_{\mathrm{pe}}(a,\nu ,Z_\mathrm{d})QY\frac{4\pi J_\nu }{h\nu }𝑑\nu ,$$
(4)
and the cooling rate of the electron capture (recombination) by the grain is $`\lambda (a)=(3/2)k_\mathrm{B}T|R_\mathrm{e}|`$, where we have assumed the Maxwellian distribution for the gaseous electrons. We also have defined $`E_{\mathrm{pe}}`$, the mean kinetic energy of photoelectrons that are emitted from a grain with radius $`a`$ and charge $`Z_\mathrm{d}`$ when a photon with energy $`h\nu `$ is absorbed. This is determined from the energy distribution function of the photoelectrons, $`f(E)`$, which is assumed to be a parabolic function introduced by WD01, 111Our $`f(E)`$ is equal to $`f_E(E)=f_E^0(E)/y_2`$ in WD01. and is consistent with a typical energy distribution of photoelectrons in the laboratory experiments. Hence, we estimate

$$E_{\mathrm{pe}}=_0^{E_{\mathrm{max}}}Ef(E)𝑑E=\frac{E_{\mathrm{max}}(E_{\mathrm{max}}2E_{\mathrm{min}})}{2(E_{\mathrm{max}}3E_{\mathrm{min}})},$$
(5)
where $`E_{\mathrm{max}}=h\nu IP`$ is the maximum energy of the photoelectrons and $`E_{\mathrm{min}}=e^2(Z_\mathrm{d}+1)/a`$ is the minimum energy of the photoelectrons appearing on the grain surface (however falling back into the grain). For the parabolic function adopted here, a typical kinetic energy of the photoelectrons is about half of $`E_{\mathrm{max}}`$.

To estimate the total photoelectric heating rate, we specify a characteristic grain size in the IGM instead of the size distribution of the IG grains. This is because the size distribution in the IGM is quite unknown. Here we consider three cases for the typical grain size; 0.001 µm, 0.01 µm, and 0.1 µm. The smallest case of 0.001 µm is based on the size of the primary grains produced in the ejecta of SNe II (Todini & Ferrara, 2001). It is interesting to obtain an information about the size of the IG dust as well as to put limits on the amount. From the analysis of the case of very small IG dust, we can assess the selection rule of the IG dust suggested by Aguirre (1999).

The grain number density is given as $`n(a)=\rho _\mathrm{d}/m(a)`$, where $`m(a)=(4\pi /3)\varrho a^3`$ is the mass of a grain with the radius $`a`$, $`\varrho `$ $`(=2\mathrm{g}\mathrm{cm}^3)`$ is the grain material density, and $`\rho _d`$ is the volume grain mass density, which is given by $`\rho _\mathrm{d}=\mu m_\mathrm{H}n_\mathrm{b}𝒟`$, where $`\mu `$ is the mean atomic weight, $`m_\mathrm{H}`$ is the proton mass, $`n_\mathrm{b}`$ is the baryon number density, and $`𝒟`$ is the dust-to-gas mass ratio. We introduce a relative dust-to-gas ratio, $`\zeta `$, defined as $`𝒟=\zeta 𝒟_{\mathrm{MW}}`$, where $`𝒟_{\mathrm{MW}}=6\times 10^3`$ is the Galactic dust-to-gas ratio (Spitzer, 1978), for convenience below. Therefore, the total photoelectric heating rate and electron capture cooling rate per unit volume are $`\mathrm{\Gamma }_{\mathrm{pe}}=\gamma (a)n(a)`$ and $`\mathrm{\Lambda }_{\mathrm{pe}}=\lambda (a)n(a)`$, respectively.

The total heating rate is approximately proportional to $`a^{0.5}`$ under the situation considered in the next section. This is because the heating rate per a grain roughly has a dependence of $`a^{2.5}`$, i.e., the grain cross section plus an additional $`a`$ dependence caused by $`E_{\mathrm{pe}}QY𝑑\nu `$, whereas the grain number density is proportional to $`a^3`$. Therefore, for a fixed dust mass, the photoelectric effect becomes more important as the grain size is smaller. The cooling rate by the electron capture is always about an order of magnitude smaller than the photoelectric heating rate.

## 3 Dust and Thermal History in IGM

Suppose an ideal fluid element with the mean density of the IGM. The time evolution of its temperature is governed by (e.g., Hui & Gnedin 1997)

$$\frac{dT}{dt}=2HT\frac{T}{X}\frac{dX}{dt}+\frac{2(\mathrm{\Gamma }\mathrm{\Lambda })}{3k_\mathrm{B}Xn_\mathrm{b}},$$
(6)
where $`H`$ is the Hubble constant, $`n_\mathrm{b}`$ is the cosmic mean number density of baryon, $`\mathrm{\Gamma }`$ and $`\mathrm{\Lambda }`$ are the total heating and cooling rates per unit volume, respectively, $`X`$ is the number ratio of the total gaseous particles to the baryon particles, i.e., $`Xn_i/n_\mathrm{b}`$, where $`n_i`$ is the number density of the $`i`$-th gaseous species and we consider H i, H ii, He i, He ii, He iii, and electron. Then, we solve equation (6) coupled with the non-equilibrium rate equations for these gaseous particles. For the rate coefficients, we use those compiled by Cen (1992). The adopted initial conditions are $`T=T_{\mathrm{CMB}}(1+z)`$, H i/H(total)=0.99, He i/He(total)=0.99, He ii/He(total)=0.0099 at redshift $`z=10`$, where $`T_{\mathrm{CMB}}=2.7`$ K is the temperature of the CMB radiation at $`z=0`$. The result is insensitive to the choice of the initial condition if calculation is started at an epoch before H i reionization. We also adopt $`H_0=70`$ km s-1 Mpc-1, $`\mathrm{\Omega }_\mathrm{M}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_\mathrm{b}=0.04`$, and He mass abundance $`Y=0.24`$.

To add to the photoelectric heating and the electron capture cooling by grains, we consider the following heating/cooling mechanisms: For cooling, we use the recombination cooling, collisional ionization/excitation cooling, bremsstrahlung cooling, and Compton cooling compiled by Cen (1992). For heating rates, to take account of the radiative transfer effect (Abel & Haehnelt, 1999), we multiply usual photoionization heating rates by a correction factor of $`(1+C_{\mathrm{RT}}f_i)`$, where $`f_i`$ is the fractional abundance of H i, He i, and He ii relative to total number of H and He. This factor mimics the situation that atoms are ionized by higher energy photons when the optical depth is large (i.e., $`f_i`$ is nearly unity) relative to the optically thin limit. The coefficient $`C_{\mathrm{RT}}`$ is a parameter determined by solving the cosmological radiative transfer. The metal line cooling is negligible because we consider only the case that metallicity is less than 1/100 solar value and temperature is less than a few $`10^4`$ K (see e.g., Sutherland & Dopita 1993).

We need to give the UV background radiation and the cosmic reionization history to obtain the temperature evolution. We adopt, for UV background, $`J_\nu =J_{21}\nu ^\alpha `$, where $`J_{21}`$ is the mean intensity at the Lyman limit of H i in unit of $`10^{21}`$ erg s-1 cm-2 sr-1 Hz-1 and its evolution is given by equation (3) in Kitayama et al. (2001). We assume the spectral index to be $`\alpha =1`$ as a QSO dominated background radiation. For the reionization history, we assume, for simplicity, a sudden reionization of H i (and He i) at $`z=6.0`$ (Becker et al. 2001)222Although Kogut et al. (2003) find a higher-$`z`$i reionization, our constraint of the IG dust dose not change because it is determined mainly from the IGM temperature after the He ii reionization.; there is no UV photon before this redshift. To mimic the He ii reionization at $`z3.4`$ (Theuns et al., 2002), we set the maximum photon energy of the UV background to be 54.4 eV (He ii Lyman limit) for $`z>3.4`$ and to be 1.24 keV333This is determined by the maximum frequency in the data of grain absorption coefficient. for $`z3.4`$ (i.e., a sudden He ii reionization). Other choice of $`\alpha `$ will be discussed later.

We compare our theoretical thermal history at the mean density of the IGM with the observational one obtained by Schaye et al. (2000). They observe the Ly$`\alpha `$ forest clouds with the column density of $`10^{1315}`$ cm-2 (i.e., slightly over density regions), and then they convert the temperature at this density range estimated from the minimum $`b`$ (Doppler) parameters into that at the mean density of the IGM by using the observed equation of state of the IGM. Thus, we can compare both thermal histories directly.

In Figure 1, we show the IGM thermal history obtained for no dust cases and dusty cases of $`\zeta =1/100`$ (silicate grains). For no dust cases (three dotted curves), the calculated thermal histories are consistent with the data points of Schaye et al. (2000) if $`C_{\mathrm{RT}}3.0`$–7.0. For dusty cases (solid \[$`a=`$ 0.1 µm\], dash-dotted \[$`a=`$ 0.01 µm\], and dashed curves \[$`a=`$ 0.001 µm\]), we observe, especially after the He ii reionization, temperature enhancement by the dust photoelectric heating. In Table 1, we summarize the temperature enhancement factors at some redshifts for some sets of dust size and dust-to-gas ratio. All calculations do not include the evolution of the dust-to-gas ratio, for simplicity.

We also present the time evolution of heating rates per a unit volume for the case of $`a=0.01`$ µm and $`\zeta =1/100`$ in Figure 2. The solid, dashed, dotted, and dash-dotted curves are the dust photoelectric, H i, He i, and He ii photoionization heating rates, respectively. These heating rates show sudden change at the He ii reionization epoch ($`z=3.4`$). This is because the maximum photon energy increases from 54.4 eV to 1.24 keV suddenly at that time. We note here that the heating rate by grains is saturated even when the maximum photon energy increases more than about 300 eV in the case of the spectral index of the incident radiation of $`\alpha =1`$. Thus, our results are robust for the choice of the maximum energy after the He ii reionization as long as we adopt the energy larger than 300 eV. In addition, we show the cases of $`a=0.1`$ and 0.001 µm (and $`\zeta =1/100`$) by thin curves in Figure 2. Even if we adopt other sets of $`a`$ and $`\zeta `$, the other curves of the photoionizations do not change so significantly.

Let us discuss each case in more detail. For the cases of $`a=0.01`$–0.1 µm, before the He ii reionization, the photoelectric heating by grains is minor heating mechanism relative to the H i photoionization heating. Thus, the effect of the photoelectric heating on the IGM temperature is very small, i.e., temperature enhancement is less than a few %. After the He ii reionization, the photoelectric heating dominates other heating mechanisms (but by a factor of less than 3) when $`\zeta =1/100`$. However, the temperature enhancement factor is still less than 1.5 (solid and dash-dotted curve in Figure 1). Thus, the thermal histories of these cases are consistent with that obtained by Schaye et al. (2000) when $`C_{\mathrm{RT}}=3.0`$ is adopted.

For the case of very small ($`0.001`$ µm) grains, the photoelectric heating is always dominant when $`\zeta 1/100`$. Especially, this heating exceeds the He ii photoionization heating by a factor of 10 after the He ii reionization. Then the temperature enhancement factor becomes about 2 at $`z=2`$ (dashed curve in Figure 1). Thus, the thermal history in this case may be inconsistent with the data from Schaye et al. (2000). This means that if the typical grain size in the IGM is $`0.001`$ µm, the dust-to-gas mass ratio in the IGM should be much less than 1/100 Galactic value.

## 4 Discussion

In Table 2, we summarize the obtained upper limits of the amount of the IG dust at $`z2`$. Let us compare these limits with the IG dust amount estimated from a possible cosmic star formation history (e.g., Madau, Pozzetti, & Dickinson 1998). Once we assume a cosmic star formation history, the cosmic metal amount can be estimated. Roughly, we find $`\mathrm{\Omega }_{\mathrm{metal}}`$ a few $`\times 10^5`$ at $`z3`$ (e.g., Aguirre 1999). If we postulate $``$ 0.5 as the metal depletion to dust grains, which number is suitable for the Milky Way, we obtain $`\mathrm{\Omega }_{\mathrm{dust}}10^5`$ as the total cosmic dust amount at the redshift. We show in Table 2 the upper limits in terms of the density parameter converted from the dust-to-gas ratios by $`\mathrm{\Omega }_{\mathrm{dust}}^{\mathrm{IGM}}(z)=\zeta 𝒟_{\mathrm{MW}}\mathrm{\Omega }_\mathrm{b}(z)`$. For the larger grain models, we obtain an upper limit, $`\mathrm{\Omega }_{\mathrm{dust}}^{\mathrm{IGM}}7\times 10^6`$ which is comparable to $`\mathrm{\Omega }_{\mathrm{dust}}`$. This means that our upper limit may not deny the possibility that most of the dust grains produced in galaxies before $`z3`$ escape to the IGM as an extreme case. For very small (0.001 µm) grain case, we obtain $`\mathrm{\Omega }_{\mathrm{dust}}^{\mathrm{IGM}}7\times 10^7`$. Thus, a small fraction ($``$ 10%) of such very small grains can escape from the host galaxies to the IGM because of $`\mathrm{\Omega }_{\mathrm{dust}}^{\mathrm{IGM}}/\mathrm{\Omega }_{\mathrm{dust}}0.1`$. The latter conclusion is consistent with the suggestion by Aguirre (1999), who proposes that the IGM is very difficult to be polluted by the very small grains.

The obtained values of the heating rate by the dust photoelectric effect are consistent with those in Nath et al. (1999). It indicates that our result is robust against uncertainties of the adopted efficiency factor $`Q`$, photoelectric yield $`Y`$, and energy distribution function of the photoelectrons, etc., because these are somewhat different between our model and Nath et al. (1999). If we consider only graphite grains, the photoelectric heating rate is found to be about 70% of that of only silicate grains, i.e., silicate grains can be more efficient heating source. The temperature in the graphite case is lower than that in the silicate case, while the decrement is $``$ 10%. These results are also consistent with Nath et al. (1999).

We shall comment some uncertainties of our analysis. First, the effect of the spectral index of $`\alpha `$ is examined. As shown in Zheng et al. (1997), for example, a softer spectrum than $`\alpha =1`$ is still compatible with observations. If $`\alpha =1.5`$, the heating rates by dust at $`z=3`$ decrease by a factor of 0.25 ($`a=`$ 0.1 $`\mu `$m)–0.5 ($`a=`$ 0.001 $`\mu `$m) relatively to our case of $`\alpha =1`$. If $`\alpha =2`$ is used, the heating rates at $`z=3`$ decrease by a factor of 0.1 (0.1 $`\mu `$m)–0.25 (0.001 $`\mu `$m). In these calculations, $`\zeta =1/100`$ and $`C_{\mathrm{RT}}=3.0`$ are assumed. As the dust heating becomes less efficient, a more abundant IG dust is allowed. Then the upper limits for the IG dust amount increase by a factor of 2 (0.001 $`\mu `$m)–4 (0.1 $`\mu `$m) for $`\alpha =1.5`$ and by a factor of 4 (0.001 $`\mu `$m)–10 (0.1 $`\mu `$m) for $`\alpha =2`$ relative to the $`\alpha =1`$ case, because the dust heating rate is proportional to the dust-to-gas ratio linearly. To summarize, a softer spectral index constrains more loosely the amount of the IG dust by a factor.

It might be possible to take very large $`\alpha `$. For example, we consider a background radiation field dominated by only galaxies (e.g., $`\alpha =5`$). In this case, the dust heating rates at $`z=3`$ become 1/100 (0.1 $`\mu `$m)–1/10 (0.001 $`\mu `$m) of those in $`\alpha =1`$. However, the thermal histories with $`\alpha =5`$ never reproduce a steep rise of the observed temperatures at the He ii reionization. In addition, the transfer correction should be smaller as the spectrum is softer (Abel & Haehnelt, 1999). Therefore, the background spectrum is unlikely to be as soft as $`\alpha =5`$ after the He ii reionization. We should assume a harder index.

A caveat against our results is in the correction factor of the transfer effect. Indeed, the effects of the photoelectric heating and the transfer correction can be offset. Thus, we should check this point by solving the cosmological radiative transfer in the future work.

The current work may give a hint for the origin of the large $`b`$ region of the observed IGM clouds. There is a possibility that the large $`b`$ regions are localized. If the very small grains are located at the restricted areas, a systematic heating owing to the dust photoelectrons may explain the large-$`b`$ at the localized region. The principal role of the photoelectric heating in the IGM temperature also indicates that the theoretical equation of state in the IGM should be reconstructed by including the photoelectric heating effect. All of these works are being developed by the authors.

## Acknowledgments

We would appreciate valuable comments by the referee, Dr. A. Aguiree. We are grateful to R. Hirata for continuous encouragement, and also thank H. Hirashita and S. Bianchi for very useful discussion. AKI is supported by the Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists.