text:
# Comment on “Does the transverse electric zero mode contribute to the Casimir effect for a metal?

## Abstract

Recently J. S. Høye, I. Brevik, J. B. Aarseth, and K. A. Milton \[Phys. Rev. E 67, 056116 (2003)\] proposed that if the Lifshitz formula is combined with the Drude model, the transverse electric zero mode does not contribute to the result for real metals and there arises a linear temperature correction to the Casimir force at small temperatures. The authors claim that in spite of the fact that the Casimir entropy in their approach is negative, the Nernst heat theorem is satisfied. In the present Comment we show that the authors’ conclusion regarding the Nernst heat theorem is in error. We demonstrate also the resolution of this thermodynamic puzzle based on the use of the surface impedance instead of the Drude dielectric function. The results of numerical computations obtained by the authors are compared with those from use of the surface impedance approach which are thermodynamically consistent.

PACS number(s): 05.30.-d, 11.10.Wx, 73.61.At, 77.22.Ch

A recent calculation shows that the transverse electric zero mode does not contribute to the Casimir free energy . This results in an observable linear temperature correction to the Casimir force between parallel plates made of real metal . According to Ref. , the Casimir free energy increases with increasing temperature in some temperature interval, and the Casimir entropy becomes negative in this interval. From the view-point of the authors of Ref. , this, however, does not lead to violation of the Nernst heat theorem. Below we show that on the contrary the transverse electric zero mode does contribute to the Casimir free energy which results in the absence of a linear in temperature correction to the Casimir force at small temperatures. In reality, the Casimir free energy is universally a decreasing function of temperature, and the Casimir entropy is always positive. We demonstrate also that if, as is done in Ref. , the Drude model and the Bloch-Grüneisen formula are used to calculate the Casimir force, the Nernst heat theorem is nesessarily violated.

The correct description of the thermal Casimir force between real metals has assumed great importance due to the recent precision experiments , followed by the prospective applications of the Casimir effect in nanotechnology  and also its use as a test for predictions of fundamental physical theories . Unexpectadly, it was found that the calculations of the thermal Casimir force on the basis of the Lifshitz formula  supplemented by the Drude model ran into serious difficulties . The key question of the controversy is whether the transverse electric zero mode contributes to the Casimir effect in the case of real metals. Boström and Sernelius  have used the Drude model to compute the contributions of both transverse electric and transverse magnetic modes and found that the transverse electric zero mode does not contribute. Later on it was shown  that in this approach the Casimir entropy is negative in some temperature interval and the third law of thermodynamics (the Nernst heat theorem) is violated. On the contrary, the authors of Ref.  claim that the entropy is zero at zero temperature, i.e. the third law of thermodynamics is not violated.

To prove this statement they present in Sec. IV of Ref.  two sets of arguments, analytical and numerical. Both of them start from Eq. (3.4) representing the Lifshitz result for the Casimir free energy per unit area in the configuration of two semispaces separated by a distance $`a`$:

$$\beta F=\frac{1}{2\pi }\stackrel{}{\underset{m=0}{\overset{\mathrm{}}{}}}_{\zeta _m}^{\mathrm{}}\left[\mathrm{ln}\left(1\lambda _m^{TM}\right)+\mathrm{ln}\left(1\lambda _m^{TE}\right)\right]q𝑑q,$$
(1)
where

$`\lambda _m^{TM}=A_me^{2qa},\lambda _m^{TE}=B_me^{2qa},`$ (2)
$`A_m=\left[{\displaystyle \frac{\epsilon _mq\sqrt{(\epsilon _m1)\zeta _m^2+q^2}}{\epsilon _mq+\sqrt{(\epsilon _m1)\zeta _m^2+q^2}}}\right]^2,B_m=\left[{\displaystyle \frac{q\sqrt{(\epsilon _m1)\zeta _m^2+q^2}}{q+\sqrt{(\epsilon _m1)\zeta _m^2+q^2}}}\right]^2.`$ (3)
Here the dielectric permittivity $`\epsilon _m\epsilon (i\zeta _m)`$ is computed on the imaginary Matsubara frequencies, $`\zeta _m=2\pi m/\beta `$, $`\beta 1/T`$, and units with $`k_B=\mathrm{}=c=1`$ are used.

According to the Drude model

$$\epsilon (i\zeta )=1+\frac{\omega _p^2}{\zeta (\zeta +\nu )},$$
(4)
where $`\omega _p`$ is the plasma frequency, $`\nu =\nu (T)`$ is the relaxation parameter. Substituting Eq. (4) into Eq. (2) one obtains $`B_0(0)=0`$, i.e. the transverse electric zero mode does not contribute in the case of the Drude metal in accordance with Refs. .

The authors advance three analytical arguments in support of the validity of the Nernst heat theorem in Ref. , but no direct proof is given. As the starting point, two simple systems containing three harmonic oscillators each are considered. It is shown that for some “subsystem” the entropy may be negative. Needless to say, systems of this type are not realistic models of full-scale Quantum Electrodynamical interaction between the plates made of real metal as they do not reflect the main features of this interaction. Next, as a “real metal” the case where $`\epsilon `$ is independent of $`\zeta `$ is considered. Up to order of magnitude estimation for the entropy is obtained consistent with the Nernst heat theorem. This estimation, however, is not relevant to the Drude metal because the dielectric permittivity (4) depends on frequency. The final analytical argument of Ref.  is for the frequency-dependent $`\epsilon `$. In this case, however, no result is obtained. The authors speculate that “For a real metal obeying the Drude dispersion relation (3) (with $`\nu 0`$),… we expect a $`T^3`$ (or $`T^4`$) correction to the free energy at sufficiently low temperature.” If this were the case, the Nernst heat theorem would be satisfied. We show below that this is, however, not the case.

To find the asymptotic behavior of the free energy, given by Eqs. (1)-(4), on temperature in the region of low temperatures, $`TT_{eff}=1/(2a)`$, one needs the explicit dependence of the relaxation parameter $`\nu `$ on temperature. This is given by the Bloch-Grüneisen law through the linear connection between the relaxation parameter and the static resistivity. As a result, $`\nu T^\alpha `$ where at low temperatures $`\alpha 5`$ whereas at high temperatures $`\alpha =1`$  (the tabulated data for different metals can be found in Ref. ). Note that $`\nu (T)0`$ except of one point $`T=0`$. Performing the perturbation expansion in Eqs. (1)-(4) in powers of three small parameters, $`T/T_{eff}`$, $`\nu /\omega _P`$ and $`\lambda _p/(2\pi a)`$ (we consider separation distances $`a\lambda _p`$ where $`\lambda _p`$ is the plasma wavelength), the following result for the free energy is obtained

$`F={\displaystyle \frac{\pi ^2}{720a^3}}\left(12{\displaystyle \frac{\lambda _p}{\pi a}}+{\displaystyle \frac{18}{5}}{\displaystyle \frac{\lambda _p^2}{\pi ^2a^2}}\right){\displaystyle \frac{\zeta (3)}{16\pi a^3}}\left[\left(1+{\displaystyle \frac{\lambda _p}{\pi a}}\right)\left({\displaystyle \frac{T}{T_{eff}}}\right)^3{\displaystyle \frac{\pi ^3}{45\zeta (3)}}\left(1+2{\displaystyle \frac{\lambda _p}{\pi a}}\right)\left({\displaystyle \frac{T}{T_{eff}}}\right)^4\right]`$ (5)
$`+{\displaystyle \frac{\zeta (3)T}{16\pi a^2}}\left(12{\displaystyle \frac{\lambda _p}{\pi a}}+3{\displaystyle \frac{\lambda _p^2}{\pi ^2a^2}}\right)+{\displaystyle \frac{\nu }{\omega _p}}{\displaystyle \frac{T}{4\pi a^2}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\left[\stackrel{~}{\zeta }_m{\displaystyle _{\stackrel{~}{\zeta }_m}^{\mathrm{}}}{\displaystyle \frac{dy}{e^y1}}+{\displaystyle \frac{1}{\stackrel{~}{\zeta }_m}}{\displaystyle _{\stackrel{~}{\zeta }_m}^{\mathrm{}}}{\displaystyle \frac{y^2dy}{e^y1}}\right].`$ (6)
Here the nondimensional Matsubara frequencies are defined as $`\stackrel{~}{\zeta }_m=2a\zeta _m=4\pi ma/\beta `$. All the details of derivation of perturbation Eq. (6) can be found in Ref. .

It is quite clear from Eq. (6) that the expectations of Ref.  are not borne out . In addition to the higher-order terms in temperature, there is a linear term which leads to a nonzero value of entropy at zero temperature

$$S(T=0)=\frac{F}{T}|_{T=0}=\frac{\zeta (3)}{16\pi a^2}\left(12\frac{\lambda _p}{\pi a}+3\frac{\lambda _p^2}{\pi ^2a^2}\right)<0.$$
(7)
Thus, we confirm the conclusion of Refs.  that in the framework of the Lifshitz formula supplemented by the Drude model with a temperature-dependent relaxation parameter (see Fig. 6 of Ref.) the Nernst heat theorem is in fact violated.

One may mention, as the authors of Ref.  do in the end of Appendix D, that their dependence of $`\nu `$ on $`T`$ neglects the effect of impurities, which give rise to some nonzero residual resistivity at zero temperature . This, however, cannot remedy the situation concerning the inconsistency with the Nernst heat theorem. In fact, the resistivity ratio of a sample can be defined as the ratio of its resistivity at room temperature to its residual resistivity. For pure samples the resistivity ratio may be as high as $`10^6`$ . As an example let us consider $`Au`$ with $`\nu (T=300\text{K})=5.32\times 10^{13}`$rad/s . In this case for the residual value of the relaxation parameter one obtains $`\nu _{res}=5.32\times 10^7`$rad/s. The asymptotic expression (6) is applicable under the condition $`\nu \zeta _1=2\pi T`$ (this condition is automatically fulfiled for the above used dependence $`\nu (T)T^\alpha `$ with $`\alpha 1`$). Thus, with allowance made for impurities, the asymptotic (6) is applicable at temperatures $`T\nu _{res}/(2\pi )=6.5\times 10^5`$K. What this means is that the entropy at temperature $`T=5\times 10^4`$K has a nonzero negative value given by Eq. (7). Physically this is equivalent to the violation of the Nernst heat theorem. We would like to point out also that the usual theory of the electron-phonon interaction, describing the electrons interacting with the elementary excitations of a perfect lattice with no impurities, must satisfy and does satisfy all the requirements of thermodynamics. That is why, the attempt to remedy the violation of the Nernst heat theorem at the expense of impurities is meaningless.

We would now like to address the numerical arguments of Ref.  in support of the validity of the Nernst heat theorem in the Lifshitz theory combined with the Drude model. All the numerical computations in Ref.  are performed for $`Au`$ with $`\omega _p=9.0`$eV, $`\nu (T=300\text{K})=35`$meV. Once again, no direct computations by Eqs. (1)-(4) with $`\nu (T)`$ given by Fig. 6 of Ref.  are done. Instead, when calculating force-temperature relation, the room temperature data for $`\epsilon (i\zeta )`$ or the above value of $`\nu (T=300\text{K})`$ were used ignoring the correct temperature dependence \[e.g. our Eq. (4) combined with Fig. 6 of Ref.  presenting $`T`$-dependence of $`\nu `$\]. All computed force-temperature and force-distance curves demonstrate nonmonotonous behavior (see Figs. 2-4 of Ref. ). This is a counterintuitive effect as the authors of Ref.  themselves recognize. In fact, with an increase of temperature the population of modes and, consequently, force modulus should increase. The nonmonotonous character of the force curves is an artifact. It is accounted for by the absence of the transverse electric zero mode contribution for metals in the formalism under discussion (see below).

In an effort to confirm the validity of the Nernst heat theorem in their formalism, the authors of Ref.  compute numerically the force-temperature dependence for dielectrics with constant $`\epsilon =10^2\mathrm{\hspace{0.17em}\hspace{0.17em}10}^3,\mathrm{\hspace{0.17em}10}^4`$, and $`\epsilon =\mathrm{}`$ (Fig. 5 of Ref. ). They consider the horizontal slope near the point $`T=0`$, which is present in all curves, except for $`\epsilon =\mathrm{}`$ and for $`Au`$ as a desirable property. (“If the force had a linear dependence on $`T`$ for small $`T`$ so would the free energy $`F`$, in contradiction with the requirement that the entropy $`S=F/T`$ has to go to zero as $`T0`$.”) The absence of a desirable slope for $`Au`$ is considered as a lack of resolution on the scale of the Figure. It is noted also in  that not only dispersive but also nondispersive curves are nonmonotonous.

All these conclusions depend on the authors’ physically incorrect assumption that there exist dielectrics with arbitrary large constant $`\epsilon `$ and that metals may be considered as a limiting case of such kind dielectrics when $`\epsilon \mathrm{}`$ on the imaginary frequency axis. In actual truth, $`\epsilon `$ can be assumed to be frequency- and temperature-independent only in the case of so called non-polar dielectrics whose atoms or molecules do not have their own dipole moments. The electric susceptibility of non-polar dielectrics arises due to the electronic polarization of atoms and molecules. The values of $`\epsilon `$ for non-polar dielectrics are of order of one only . Large values of $`\epsilon `$ can exist only for polar dielectrics where the partial orientation of permanent dipole moments occurs. But for polar dielectrics $`\epsilon `$ depends strongly on the frequency and temperature. Specifically, their $`\epsilon `$ quickly decreases with the increase of frequency. As a result, at optical and infrared frequencies, which are characteristic for the Casimir effect, the values of $`\epsilon `$ are determined by the electronic polarizability  and cannot exceed several units. Thus, the nonmonotonous force-temperature curves presented in Fig. 5 of Ref.  are not relevant to the resolution of the above puzzle with the Nernst heat theorem.

The reference to the experiment by Bressi et al  where, supposedly, “the observed Casimir forces were lower than those predicted by the traditional (SDM) theory for conducting plates, in cases where the distances were low, $`a0.5\mu `$m” is a misunderstanding. According to Ref. , this reduction effect is apparent from Fig. 4 of  and could experimentally confirm the characteristic temperature variations predicted by the authors. In Fig. 4 of , however, all data are presented only for $`a>0.5\mu `$m and there is no evidence that these data support expectations of .

In view of the above, we arrive at the conclusion that the arguments presented in Ref.  are not correct. The use of the Drude model in combination with the Lifshitz formula actually leads to the violation of the Nernst heat theorem and the other nonphysical features such as a linear temperature correction to the Casimir force at small separations, negative values of entropy and nonmonotonous force-temperature and force-distance relations.

Recently the resolution of these complicated problems was obtained  using an alternative approach to the description of real metals based on the concept of the surface impedance. This approach offers a fundamental understanding of the reason why the Drude model is not compatible with the theory of the thermal Casimir force between real metals. This is due to the fact that the Drude dielectric function (4) is obtained for the frequency region of the normal skin effect, where the electromagnetic oscillations penetrate through the skin layer, and lead to a real current of the conduction electrons. The interaction of the conduction electrons with the elementary excitations of the crystal lattice (phonons) leads to the occurrence of electric resistance and heating of a metal. By contrast, the thermal photons in thermal equilibrium with a metal plates or, much less, the virtual photons giving rise to the Casimir effect cannot, under any circumstances, lead to the initiation of a real current and heating of a metal. Of course, this is strictly prohibited by thermodynamics. Hence the standard concept of a fluctuating electromagnetic field penetrating inside a metal described by the Drude dielectric function (4) fails to describe virtual and thermal photons. As a consequence, the Lifshitz formula in combination with the Drude model leads to the above contradictions with thermodynamics. Note that the dielectric permittivity depending only on frequency is also inapplicable in the frequency domain of the anomalous skin effect (see Ref.  for details).

In contrast to the dielectric permittivity, the surface impedance $`Z(\omega )`$ is defined at all frequencies. The impedance boundary conditions are

$$𝑬_t=Z(\omega )\left[𝑩_t\times 𝒏\right],$$
(8)
where $`𝑬_t`$, $`𝑩_t`$ are the tangential components of electric and magnetic fields, $`𝒏`$ is the internal normal vector to the surface. They take into account the reflection properties of real metal with no consideration of the electromagnetic fluctuations inside of it.

By the use of the surface impedance instead of the Drude model (4), the Lifshitz formula (1) is preserved, but the coefficients $`A_m`$, $`B_m`$ from Eq. (2) should be replaced by

$$A_m^{imp}=\left(\frac{qZ_m\zeta _m}{q+Z_m\zeta _m}\right)^2,B_m^{imp}=\left(\frac{\zeta _mZ_mq}{\zeta _m+Z_mq}\right)^2,$$
(9)
where the impedance is computed on the Matsubara frequencies. Substituting in Eq. (9) the impedance function of the normal skin effect or the anomalous skin effect, one finds $`B_0^{imp}(0)=1`$ . In the region of the infrared optics it follows  $`B_0^{imp}(0)=(\omega _pq)^2/(\omega _p+q)^2`$.

Let us now present several computational results obtained by Eqs. (1), (9) in the framework of the impedance approach in comparison with the results of Ref. . In Fig. 1, the magnitude of the Casimir surface force density $`^T=F/T`$ for gold is computed, in the temperature interval 1 K$`T1200`$K at a separation distance $`a=1\mu `$m. Solid line is calculated in the framework of the impedance approach (separation of $`1\mu `$m corresponds to the domain of the infrared optics), and dashed line is obtained by the approach of Ref.  (i.e. via the Lifshitz formula supplemented by the Drude model with a temperature dependent relaxation parameter). It is clearly seen that the dashed line is nonmonotonous demonstrating the existence of a wide temperature region where force modulus decreases with an increase of temperature (like in Figs. 2, 3 of Ref. ). At the same time, solid line demonstrates the monotonous increase of the modulus of the Casimir force with temperature which is consistent with our expectations on the basis of thermodynamics.

In Fig. 2, the magnitude of the Casimir force density between dielectrics with $`\epsilon =`$const is computed at different temperatures at a separation of $`a=1\mu `$m. Both lines were obtained by the usual Lifshitz formula (1), (2) with $`\epsilon _m=\epsilon =`$const. Solid line is for mica with $`\epsilon =7`$; dashed line faithfully copies the line of Fig. 5 of Ref.  with $`\epsilon =100`$. Solid line demonstrates the monotonous increase of the Casimir force with increasing temperature as is expected from thermodymanics. Dashed line is nonmonotonous because non-polar dielectrics with so high $`\epsilon `$ do not exist in nature due to bounds on the possible values of electronic polarization. From this figure it is clear that the Lifshitz formula for dielectrics is consistent with thermodynamics when the correct input data for $`\epsilon `$ are substituted.

In Fig. 3, the Casimir entropy for gold is plotted as a function of temperature at a separation distance between plates of $`a=1\mu `$m. Solid line is computed in the framework of the impedance approach. Dashed line is obtained by the approach of Ref. , i.e. by the usual Lifshitz formula and Drude model (1)-(4) with a relaxation parameter depending on temperature according to Fig. 6 of Ref. . Evidently, the solid line satisfies all conditions, i.e. positive values of entropy at nonzero temperatures, and the validity of the Nernst heat theorem. By contrast, the dashed line presents negative values of entropy and the violation of the Nernst heat theorem (for two semispaces separated by a gap there is no possibility to introduce some “composite system” whose “subsystem” with a negative entropy these semispaces would be; for this reason the Casimir entropy must be positive which is in fact the case in the impedance approach). The analytical proof of the validity of the Nernst heat theorem in the impedance approach can be found in Ref. .

Now we are in a position to answer the questions raised in the introductory note to this Comment. The answer to the question on whether the transverse electric zero mode contributes to the Casimir effect is yes. As a consequence, there is no linear temperature correction to the Casimir force between metallic plates at small temperatures. The Casimir free energy is always decreasing and the force magnitude increases with increase of temperature. As to the Casimir entropy, it is always positive and obeys Nernst’s theorem as it must be in accordance with thermodynamics. Regarding the doubts not only on the applicability of the Drude model as such, but even more, doubt on the applicability of the fundamental Lifshitz formula, they are carried too far. Although, as discussed above, the Drude model is in fact not appropriate to describe the thermal Casimir effect in the case of real metals, the Lifshitz formula (1) with coefficients (9) in terms of the surface impedance is perfectly well suited to calculate all the quantities of physical interest.


format: text

text:
# Proper motion, age and initial spin period of PSR J0538+2817 in S147

## 1 Introduction

The determination of the birth properties of pulsars is of crucial importance in understanding both the physics of core-collapse supernovae as well as the population and evolution of a radio pulsar. The initial spin period of pulsars, $`P_0`$, is particularly difficult to measure as it requires the knowledge of the pulsar’s age, $`\tau `$, and its spin-down behaviour. It is usually assumed that the observed evolution of the spin-frequency $`\nu =1/P`$ can be described by a power-law $`\dot{\nu }\nu ^n`$, where $`n`$ is the so-called braking index. The initial spin period can then be calculated from

$$P_0=P\left[1\frac{n1}{2}\frac{\tau }{\tau _c}\right]^{\frac{1}{n1}}$$
(1)
where $`\tau _c=P/2\dot{P}`$ is the characteristic age of the pulsar. The characteristic age is a good estimator for the true age of the pulsar, $`\tau _c\tau `$, under the assumptions that $`P_0P`$ and that the spin-down is due to magnetic braking for which $`n=3`$.

Since the use of the characteristic age, rather than the true age, can lead to considerable errors, it is desirable to have an independent age measurement. The age of a supernova remnant (SNR) that originated in the same explosion as the pulsar, can serve as such an estimator but the only pulsar for which the age of an associated SNR is clearly known, is the Crab pulsar. From the observation of the explosion in A.D.1054 and a measured braking index of $`n=2.51(1)`$, the initial spin period is computed to be $`P_0=`$19 ms (Lyne et al. 1993). It may also be possible to measure a proper motion of a pulsar associated with a SNR. If this transverse motion is directed away from the center of the SNR, then this is strong evidence that the pulsar is genuinely associated with the remnant. It also allows one to determine the age of both the pulsar and the SNR by comparing the present offset from the center with its speed. As SNRs typically fade away after $`100,000`$ yr, pulsars genuinely associated with SNRs are necessarily young. Such pulsars, however, often show rotational instabilities in the form of glitches and/or timing noise (e.g. Lyne et al. 1995) which makes the measurement of proper motion via timing observations usually a difficult task. Interferometric measurements offer a solution but so far it has only been possible for PSR B1951+32 in CTB 80 for which an initial spin period of $`P_0=27(6)`$ ms was derived (Migliazzo et al. 2002).

In this work we present a timing proper motion measurement for the 143-ms pulsar J0538+2817 which was found within the boundaries of G180.0$``$1.7, also called S147 (Anderson et al. 1996). The SNR S147 has a prominent shell structure with a radius of $`\theta =83(3)`$ arcmin (Sofue et al. 1980). With an estimated age of $`100`$ kyr (e.g. Kundu et al. 1980) it is considered to be one of the oldest well-defined SNRs in the Galaxy, although other authors derived much younger ages (e.g. 20 kyr, Sofue et al. 1980). Therefore, S147 has been studied rather extensively at radio frequencies (e.g. Fürst et al. 1982). Very recently, X-ray observations with CHANDRA have revealed a structure interpreted as a pulsar wind nebula (Romani & Ng 2003) while XMM-Newton observations revealed pulsed X-ray emission from the surface (McGowan et al. 2003). In the following we describe our timing observations and data analysis which leads to a proper motion measurement which clearly associates the pulsar with SNR S147. We hence obtain accurate estimates of the age of the pulsar and derive its initial spin period. The results are finally compared to the recent X-ray observations.

## 2 Observations and Data Analysis

The observations were made with the 100-m Effelsberg radiotelescope from April 1994 at 1410 MHz and with the 76-m Lovell telescope at Jodrell Bank from March 1996 at 606 MHz and 1400 MHz. At both telescopes, two circularly polarized signals were mixed down to intermediate frequencies, detected and incoherently de-dispersed in hardware using filterbanks before sub-integrations of 15 sec (Effelsberg) and 60 sec (Jodrell Bank) were written to disk for off-line processing. All data were time-stamped with clock information provided by local H-maser clocks which were later synchronized to UTC by using signals from the Global Positioning System (GPS) satellites. Details of the observing systems can be found in Anderson et al. (1996) and Hobbs et al. (2003a).

The data from both telescopes were first processed to a common time resolution of 238.3 $`\mu `$s before being subjected to the same template matching procedure which used identical templates to determine the pulse times-of-arrival (TOAs). Since PSR J0538+2817 shows mode-changing, exhibiting two distinct profiles which differ in the relative height of the two prominent components (Anderson et al. 1996), we applied a proven technique that uses two different, carefully created templates for the given modes (see e.g. Stairs et al. 2000). Any possible remaining effects due to the occurrence of mixed-mode profiles are accounted for by adopting relatively large minimum errors of $`\mathrm{\Delta }t180\mu `$s (Effelsberg) and $`\mathrm{\Delta }t130\mu `$s (Jodrell Bank), respectively. The success of this procedure is notable by the fact that no arbitrary clock offsets between the telescopes or profile modes were needed in the further timing analysis.

In a first step of the timing analysis using the DE200 planetary ephemerides, the dispersion measure (DM) was determined by fitting for DM and a simple spin-down model to data obtained with the Lovell telescope at 606 and 1400 MHz over a small period of time. The DM value was then held fixed for the subsequent analysis. Like many other young pulsars, PSR J0538+2817 shows long-term timing noise visible in the timing residuals. Hobbs et al. (2003a) developed a new technique to remove such timing noise, allowing one to reliably separate the signature due to proper motion. Applying this technique, the transverse proper motions have been determined for more than 300 pulsars and excellent agreement is found for those pulsar where interferometric measurements are available (Hobbs et al. 2003a,b). We applied this technique, before fitting to a spin-down model which now included proper motion and a second period derivative. Iterative tests were made to check the robustness of the obtained solution which also studied the effects of omitting and including different parts of the data in the analysis. We do not believe that the non-zero second period derivative is due to magnetospheric braking, therefore we do not attempt to calculate a braking index, $`n`$ (see Hobbs et al. 2003a).

As the pulsar lies close to the ecliptic with an ecliptic latitude, $`\beta `$, of only $`\beta =4^{}.9`$, position and proper motion measurements in the latitudinal direction are necessarily much less accurate than those in ecliptic longitude. In order to minimize covariances between the astrometric parameters, fits were made in ecliptic coordinates using the software package TEMPO111http://www.atnf.csiro.au/research/pulsar/timing/tempo, resulting in post-fit residuals shown in Fig. 1. The final spin and astrometric parameters are presented in Table 1. Quoted uncertainties are derived from twice the formal TEMPO error and standard Monte-Carlo simulations. Procedures for the latter are detailed in Lange et al. (2001) and results obtained for proper motion are shown in Fig. 2.

## 3 Results and Discussion

As can be seen from Table 1 and Fig. 2, we have obtained a proper motion measurement for PSR J0538+2817. Whilst the movement in ecliptic longitude is measured to high significance, the uncertainty for the proper motion in latitude is much larger. With a probability of 78%, the proper motion is positive in latitudinal direction, making the pulsar moving in the right quadrant on the sky to be consistent with a movement away from the center of S147. We measure a position angle of P.A.$`=311_{}^{}{}_{56}{}^{+28}`$. The obtained proper motion is consistent with a comparison of the interferometric position obtained by Anderson et al. (1996) and our present timing position.

Using the measured values, we can compute the pulsar’s motion in the past and compare it to the location of the center of the SNR. Fitting a circular shape to the radio contours of S147 over a frequency range from 430 MHz to 4750 MHz (Kundu et al. 1980, Sofue et al. 1980, Angerhofer & Kundu 1981, Fürst et al. 1982, Fürst & Reich 1986) we determine the SNR center to be at $`\lambda _{SNR}=85^{}.57(1)`$ and $`\beta _{SNR}=4^{}.44(1)`$ which agrees well with the center determined by Sofue et al. (1980) at 4750 MHz only. We mark the central position on the 2.7-GHz map obtained by Fürst & Reich (1986) shown in Fig. 3. We also mark the position of the pulsar 30 kyr ago as computed from the measured proper motion ($`\lambda _{30k}=85^{}.57(3)`$, $`\beta _{30k}=4^{}.5(5)`$). It is clear that the previous position of the pulsar agrees very well with the center of the SNR, strongly suggesting that the pulsar was born in the same explosion that created S147, about 30,000 years ago.

Assuming that the pulsar was born at the center of S147, we use the offset of the pulsar from this position, $`\mathrm{\Delta }\mathrm{\Theta }=2.18(4)\times 10^6`$ mas and the measured proper motion to determine a true age of pulsar and remnant of $`\tau =33_9^{+17}`$ kyr. The uncertainty in this age is dominated by the error in the proper motion measurement in latitudinal direction. In order to derive more accurate estimate, we can use the offset and motion in longitudinal direction only. This results in an age of $`\tau _\lambda =30\pm 4`$ kyr, confirming that the kinematic age of the pulsar is dramatically smaller than the characteristic age of $`\tau _c=618`$ kyr. The only assumption made in deriving this age is that the pulsar was born in the center of the SNR.

Anderson et al. (1996) already discussed the possible association of PSR J0538+2817 with S147 in detail and concluded that an association is plausible. They based their arguments on the proximity of the pulsar to the SNR center and the consistent distance estimates for both pulsar and SNR. The distances estimated for S147 ranging from of 0.8 to 1.6 kpc are well consistent with the dispersion measure distance of the pulsar, i.e. $`d=1.2`$ kpc as derived from the NE2001 model (Cordes & Lazio 2002). At this distance, our proper motion measurement yields a transverse speed of $`v=385_{130}^{+260}`$ km s-1 which is in excellent agreement with mean observed velocity of pulsars (Lyne & Lorimer 1993). Given the position of the pulsar within the SNR boundaries and a normalized angular distance of only $`\delta =\mathrm{\Delta }\mathrm{\Theta }/\mathrm{\Theta }=0.43(2)`$ away from the center, and its location in the Galactic anti-center region where we find only a rather sparse population of known SNRs and pulsars, an association seemed indeed very likely. This is now confirmed by the pulsar’s movement away from the center.

Further independent evidence is available that the characteristic age is much larger than the true age of the pulsar. Firstly, while some authors estimate a blast wave age of S147 in a range from 80 kyr to 200 kyr (e.g. Kundu et al. 1980), Sofue et al. (1980) estimated an age of only 20 kyr. Even though these age estimates for SNRs depend also on the density of the ambient interstellar medium and are known to be highly uncertain (e.g. Fürst & Reich 1986), all values are lower than $`\tau _c`$, and the latter estimate agrees indeed very well with our findings. Given this young age of the SNR, the well defined shell structure appears less surprising. If the explosion occurred in a low-density, hot stellar wind cavity blown up by the progenitor star, the expansion will not be describable by a Sedov-phase but will be free until it reaches the cavity boundaries (E. Fürst, private communication). At a distance of 1.2 kpc, the observed SNR radius corresponds to $`30`$ pc, implying an expansion velocity of 1000 km s-1 in the free expansion phase. Future CO observations may be able to reveal the wind cavity.

Secondly, we have access to polarization information for the radio emission of PSR J0538+2817 obtained by Mitra et al. (2003), who measured a rotation measure of RM=$`7\pm 12`$ rad m-2. The pulsar exhibits an extremely high degree of polarization of 92(2)%. While this is not uncommon, it is usually found in young sources (e.g. Morris et al. 1981), again supporting a young age of the pulsar.

Finally, we can compare our results to recent X-rays observations. Romani & Ng (2003) reported the discovery of a faint nebula surrounding the pulsar. They interpret this nebula as an equatorial torus, supporting the association of pulsar and S147. In their calculation they assumed a pulsar age of 100 to 200 kyr, but we can derive consistent results with our smaller derived age of $`\tau =30`$ kyr with a somewhat smaller ISM density. The orientation of their modelled torus with a symmetry axis located at position angle $`\mathrm{\Psi }=154(6)^{}`$ (measured N through E) and $`\xi =100(6)^{}`$ (into plane of the sky) is consistent with the position angle swing measured for our polarization data. Applying a rotating vector model (Radhakrishnan & Cooke 1969), we find a best fit for a magnetic inclination of $`\alpha 95^{}`$ and an impact angle of $`\sigma 2^{}`$ although the uncertainties are considerable. These findings support Romani & Ng’s interesting conclusion derived from the torus geometry that the velocity direction appears to be nearly aligned with the rotation axis of the pulsar. However, with a pulse width of about $`40^{}`$ and such a geometry, one may expect to observe also emission from the opposite magnetic pole. This is not the case, but the emission beam may not be circular and/or the emission beam could be patchy (Lyne & Manchester 1988). We also note that the bilateral symmetry axis of S147 appears to be very nearly parallel to both the symmetry axis of the torus and the proper motion direction of the pulsar. While the remnant is somewhat asymmetric about the orthogonal axis, the center of the explosion could in principle deviate from the SNR center which is determined by the sharp circular rim. However, it is reasonable to assume that the explosion center is located on the symmetry axis between the determined SNR center and the current pulsar position. In such a case, the age of the pulsar would be even younger, and our conclusions would be essentially unaffected.

Equally intriguing is the comparison of our results with those by McGowan et al. (2003) who detected pulsed X-ray emission from PSR J0538+2817. Their data are well explained by blackbody radiation from a heated polar cap. However, using the dispersion measure distance the derived temperature is significantly higher than predicted by standard cooling theories. This discrepancy can be reduced by using atmospheric fits, but with a pure-H non-magnetized atmosphere McGowan et al.’s result still falls above the expected temperature. In their analysis McGowan et al. used the characteristic age of the pulsar, and the high temperature indeed suggests that the true age of the pulsar is much smaller. Using our derived age of $`\tau =30`$ kyr, the temperature from the atmospheric fit falls well below the standard cooling curve. However, this fit also produces a distance which is much smaller than the dispersion measure estimate. A proper modelling of the observed X-ray spectra needs to include the surface magnetic field of $`B=7.3\times 10^{11}`$ G, which may change the results considerably (e.g. Pavlov et al. 2001, Zavlin & Pavlov 2002). This is important since with such a low temperature one would be forced to consider the presence of exotic cooling processes. Interestingly, a similar conclusion has been reached recently by Slane et al. (2002) for PSR J0205+6449 associated with 3C58 (but see also Yakovlev et al. 2002).

Pavlov et al. (2002) studied CHANDRA data of the neutron star 1E 1207.4$``$5209 which is likely to be associated with SNR PKS 1209$``$51/52. They also find a characteristic age which is much larger than that estimated for the SNR, suggesting a long birth period of the neutron star. Romani & Ng already already pointed out that the initial spin period of PSR J0538+2817 is likely to be large and close to the present value. Using our age estimate we derive an initial spin period of $`P_0=139.6`$ ms ($`n=3`$), which is insensitive to the actual choice of the braking index as the ratio of $`\tau /\tau _c=0.05`$ is very small (e.g. $`P_0=139.6`$ ms for $`n=10`$ and $`P_0=139.8`$ ms for $`n=0.5`$). It has been attempted to estimate the initial spin period for seven radio pulsars (see Migliazzo et al. 2002 and references therein) but apart from the results for the Crab (based on the known age of the SNR and pulsar) and PSR B1951+32 in CTB 80 (based on kinematic age derived from proper motion), all estimates rely on less certain ages estimated for associated SNRs. It appears that the initial period of those pulsars is $`P_060`$ ms. An exception is PSR J1124$``$5916 with $`P_090`$ ms (Camilo et al. 2002). All estimated birth periods are already significantly larger than what is expected from core-collapse theory of massive stars and it appears difficult to explain even spin periods of a few tens milliseconds (see Heger et al. 2003 for a recent review). The estimated initial spin period for PSR J0538+2817 is much larger still, placing strong constraints on the origin of birth kicks imparted on neutron stars as discussed by Romani & Ng (2003).

In summary, we have measured the proper motion of PSR J0538+2817 which clearly associates the pulsar with the SNR S147. From the separation of the pulsar from the SNR center we determine an age of $`\tau =30\pm 4`$ kyr, making the pulsar significantly younger than is indicated by its characteristic age of $`\tau _c=618`$ kyr. This implies a large initial spin period of $`P_0=139`$ ms. The implied possibility of exotic cooling should be revisited after a magnetized atmosphere has been fitted to the recent X-ray data.

## Acknowledgements

We thank Ernst Fürst for useful disucssions and the 2.7-GHz map of S147. We also thank Werner Becker and Slava Zavlin for helpful discussions and are grateful to an anonymous referee for his comments.


text:
# A data Grid testbed environment in Gigabit WAN with HPSS

## I Introduction

For the Atlas Japan collaboration, International Center for Elementary Particle Physics, University of Tokyo(ICEPP) will build a “Tier-1” regional center and High Energy Accelerator Research Organization(KEK) will build a “Tier-2” regional center for the Atlas experiment of the Large Hadron Collider (LHC) project at CERN. The two institutes are connected by the SuperSINET which is an ultrahigh-speed network for Japanese academic researches. On the network a Grid test bed was constructed to study requisite functionality and performance issues for the tiered regional centers.

High Performance Storage System (HPSS) with high density digital tape libraries could be a key component to handle petabytes of data produced by Atlas experiment and to share such data among the regional collaborators. HPSS parallel and concurrency data transfer mechanisms, which support disk, tape and tape libraries, are effective and scalable to support huge data storage. This paper describes about integration of HPSS into a Grid architecture and the performance measurement of HPSS in use over a high-speed WAN.

## II Test bed system

The computer resources for the test bed were installed in ICEPP and KEK site. One Grid server in each site and HPSS servers in KEK were connected with 1-Gbps Ethernet through the SuperSINET. All resources including network were isolated from other users and dedicated for the test. Figure 1 and Table 1 shows our hardware setup.

Three storage system components were employed. A disk storage server shared its host with the Grid server each at KEK and ICEPP. The remaining HPSS software components were used in the KEK Central Computer system. The HPSS data flow is depicted in Fig. 2. The HPSS Servers including core servers, disk movers, and tape movers are tightly coupled by an IBM SP2 cluster network switch.

In the case of original pftp(parallel ftp with Kerberos authentication) server performance measurement, pftpd was run in the core HPSS server. In the case of GSI-enabled HPSS server which will be mentioned in the section IV, pftpd was run in the same processors with the disk mover. The disk movers were directly connected to the test bed LAN through their network interface cards. Two HPSS disk movers were dedicated to the test.

NorduGrid middleware ran on the Grid servers. Other computing elements (CE) acted as a Portable Batch System (PBS) PBS  that was not required to install with the NorduGrid middleware.

The NorduGridnordugrid  is a pioneer Grid project in Scandinavia that added upper layer functionality, which is necessary to HEP computing, on the Globus tool kit. The middleware was simple to understand and offered functionality sufficient for our test bed study.

Table 2 shows the versions of middleware used in the test bed.

## III HPSS performance over high-speed WAN

### III.1 Basic network performance

Before end-to-end measurement, basic Gigabit Ethernet performance between IBM HPSS servers and a host at ICEPP through the WAN as well as a host on the KEK LAN was measured using netperf netperf . It is shown in Figure 3 as a function of the TCP buffer size of the client. Round Trip Time (RTT) averaged was 3 to 4 ms. The network quality of service was quite good and almost free from packet loss ($`<0.1\%`$). In this measurement, maximum TCP window size in HPSS server had 256kB (the buffer size of 256kB optimized to IBM SP2 switching network). The clients at both KEK and ICEPP had 64MB. Due to rather slower clock-speed proceessors on the HPSS servers the maximum raw TCP transfer performance was limited below 1Gbps. As seen in the graph, network access performance through both LAN and WAN became almost equivalent and saturated beyond 0.5MB buffer size.

Figure 4 shows the network performance with the number of simultaneous transfer sessions through the WAN. In the situation where TCP buffer size was 100KB, up to 4 parallel simultaneous stream sessions improved network throughput. Using greater buffer size than 1MB, multiple stream sessions did not improve the aggregate network transfer speed.

### III.2 HPSS client API performance

Figure 5 shows data transfer speed by using the HPSS client API and comparison between access from LAN and over WAN. The transfer was from/to the disk of HPSS disk-mover disk to/from client host memory. The transferred file size was 2GB in all case. Disk access speed in the disk mover was 80MB/s. It shows that even with a larger API buffer size in the client API, WAN access speed was about a half of LAN access both for reading and writing from/to HPSS server.

To increase HPSS WAN performance in future tests, the newer pdata protocol provided in HPSS 4.3 can be employed. This will improve pget performance. To get the same effect on pputs, the pdata-push protocol provided in HPSS 5.1 is required.

The existing mover and pdata protocols are driven by the HPSS mover with the mover requesting each data packet by sending a pdata header to the client mover. The client mover then sends the data. This exchange creates latency on a WAN. The pdata-push protocol allows the client mover to determine the HPSS movers that will be the target of all data packets when the data transfer is set up. This protocol eliminates the pdata header interchange and allows the client to just flush data buffers to the appropriate mover. The result is that the data is streamed to the HPSS mover by TCP at whatever rates it can be delivered by the client side mover and written to the HPSS mover devices.

### III.3 pftp client - pftpd server transfer speed

Figure 6 shows data transfer speed by using HPSS pftp from HPSS disk mover to client /dev/null dummy device. Again as in the previous HPSS client API transfer, even with a pftp buffer size of 64MB, access speed from WAN was about a half of LAN access. In addition, enabling single file transfer with multiple TCP stream by using the pftp ‘pwidth’ option was not effective in our situation. In our server layout, two disk mover hosts each had two RAID disks. Therefore, up to 4 concurrent file transfers could effect higher network utilization and overall throughput, and was so seen in WAN and LAN access case. In the same figure (Fig. 6) data transfer speed was shown from HPSS disk mover to the client disks which had writing performance of 35-45MB/s. Though disks both in server and client hosts had the access speed exceeding 30MB/s and also network transfer speed exceeded 80MB/s, overall transfer speed dropped to 20MB/s. It is because these three resources were not accessed in parallel but in series.

Figure 7 shows elapsed time for accessing data in tape library. Thanks to HPSS functionality and an adequate number of tape movers and tape drives, the data I/O throughput scaled with the number of concurrent file transfers. However, since the library had only two accessors and could load upto two tape cassetes to drives simultaneously, in the case where data in more than three off-drive tapes is required to access, the throughput goes down.

Comparison of writing to HPSS disk mover from client over WAN and LAN is shown in Fig. 8. In the figure, ‘N files $``$ N files’, for example, means that ‘reading’ N files simultaneously at client and ‘writing’ N files to the server. In our setup, HPSS server had 4 independent disks but client had only one. Reading multiple files in parallel from a single disk at client side degrades the aggregate performance by contention of disk heads.

## IV GSI-enabled pftp

GridFTPgridftp  is a standard protocol for building data GRID and supports the featues of Grid Security Infrastructure (GSI), multiple data channels for parallel transfers, partial file transfers, third-party transfer and reusable and authenticated data channels.

The pftp and ftp provided with HPSS software was not required or designed to support data Grid infrastructure. For future releases, HPSS Collaboration Members have introduced data Grid pftp requirements and the HPSS Technical Committee (TC) has convened a Grid Working Group to propose a development plan. As an interim and partial HPSS data Grid interface solution, the HPSS Collaboration is distributing the GSI-enabled pftp developed by Lawrence Berkeley National Laboratory (LBL). The HPSS TC is also working with the GridFTP development project underway at Argonne National Laboratory.

To acquire an HPSS data Grid interface necessary for our test bed, we requested and received a copy of latest version of GSI-enabled pftp. The protocol itself is pftp but it supports GSI-enabled AUTH and ADAT ftp-command.

Table 3 shows commands in each FTP protocol. While GSI-enabled pftp and GridFTP have different command set for parallel transfer, buffer management and Data Channel Authentication (DCA), the base command set is common. Fortunately unique functions of each protocol are optional and the two protocols are able to communicate. Installing and testing the GSI-enabled pftp, we proved that the GSI-enabled pftp daemon from LBL could be successfully accessed from gsincftp and globus-url-copy with no dcau option(standard globus client utilities). From NorduGrid, the server was accessible as well. The below is a sample XRSL(Extended Resource Specification Language) which utilize GSI-enabled pftp server as a storage element(SE) of the NorduGrid.

A sample XRSL

```
&(executable=gsim1)
(arguments=’’-d’’)
(inputfiles=
("Bdata.in"
"gsiftp://dt05s.cc:2811/hpss/manabe/data2"))
(stdout=datafiles.out)
(join=true)
(maxcputime="36000")
(middleware="nordugrid")
(jobname="HPSS access test")
(stdlog="grid_debug")%
(ftpThreads=1)

```

In the performance measurement with 2GB file being accessed from ’pftp client’, GSI-enabled pftp server and normal kerberos pftp server had equivalent elapsed data transfer time in any situation. Accessing from ‘Grid-FTP client’, GSI-enabled pftp server and normal pftp server, as well, had equivalent transfer time in usual. However, in the case where multiple disk movers were utilized and accessed data and GSI enabled-pftpd server resided in separated disk movers, transfer speed was halved. Figure 9 shows aggregate transfer speed over the number of independent simultaneous file transfer and shows the situation. After investigating the detailed communication between client and server, we found the differnece behaviour of the two servers. In original pftp where pftpd running in HPSS core server, data path is directly established between pftp client and disk mover. On the other hand, GSI-enabled pftp, data flow was from disk mover, via pftpd to client host. When the disk mover and pftpd server do not reside in the same host, two successive network transfer are incurred.

## V summary

ICEPP and KEK configured NorduGrid test bed with HPSS storage server over high speed GbE WAN. Data access performance was measured with several system configurations in comparison between LAN and WAN access. From that, we found that network latency affected data transfer speed with HPSS pftp and client API. The “GSI-enabled pftpd” developed by LBL was successfully adapted to the interface between Grid infrastructure and HPSS.

Our paper is a report on work-in-progress. Final results require that the questions relative to raw TCP performance, server/client protocol traffic, and pftp a protocol be further evaluated; that any necessary modifications or parametric changes be acquired form our HPSS team members; and that measurements be taken again. Further understanding of the scalability and the limitation of multi-disk mover configurations would be gained by measuring HPSS network utilization and performance using higher performance network interfaces adapters, system software and infrastructure, and processor configurations.


text:
# I Introduction

## I Introduction

ALS is the most common motor neuron disease in human adults that is characterized by selective motor neuron death Cleveland01 ; Beckman01 ; Rowland01 ; Siddique91 . In approximately 10% of the cases of ALS, the disease is inherited, thus called Familial ALS (FALS) Cleveland01 . Mutations in the cytoplasmic enzyme SOD1 were identified as the primary cause of approximately 20% of FALS cases Rosen93 ; Gurney94 . More than 90 mutations have been identified so far. The FALS mutations are scattered throughout the primary sequence and three-dimensional structure of SOD1, which exists as a homo-dimer in native state (each monomer is 153 amino acids in length) Gaudette00 . One hypothesis to explain the toxic gain-of-function of the mutants is that the toxicity is derived from intracellular aggregates and/or failure of SOD1 degradation, which is supported by the observation that in both mice and cell culture models, death of mutant neurons is preceded by formation of cytoplasmic aggregates containing mutant SOD1 Johnston00 ; Bruijn98 ; Durham97 . In addition, SOD1 knockout mice do not develop motor neuron disease Reaume96 . Aggregation or misfolding is therefore a characteristic of SOD1 mediated FALS. Toxicity may arise through (i) aberrant chemistry due to misfolded mutants Bromberg99 ; Cleveland99 ; Elliot01 , (ii) saturation of essential cellular machinery such as chaperones and proteasome components Tran99 ; Bruening01 , or (iii) oligomeric forms of the aggregate may themselves be toxic to the cells, as has been suggested for other neurodegenerative diseases Walsh02 .

Aggregation of proteins is mediated by a variety of factors including native state stability Chiti00 ; Canet02 , $`\beta `$-sheet propensity Ding02b , net charge and overall hydrophobicity Chiti02a . In addition, folding and aggregation are kinetically competing processes Chiti02b ; Ding02b ; Dobson99 . Thus, it is crucial to study the folding thermodynamics and kinetics of SOD1 in order to understand the causes of SOD1 misfolding, aggregation and toxicity. A microscopic picture of SOD1 folding and/or aggregation dynamics of SOD1 is not available from either theory or experiment. We reconstruct the SOD1 folding mechanism from molecular dynamics simulations.

Neither the structure of SOD1 aggregates nor the mechanism of aggregation is fully understood. Several experimental studies have characterized the folding thermodynamics of SOD1 and the associated FALS mutants Hallewell91 ; McRee90 ; Parge92 ; Lepock90 ; Stroppolo00 ; Mei92 ; Rodriguez02 . Lepock et al. Lepock90  and Rodriguez et al. Rodriguez02  measured the stability of mutants in differential scanning calorimetry studies of the irreversible unfolding of the distinctly metallated species of the enzyme, and showed that mutants are destabilized relative to the wild type. In addition, a crystal structure of the G37R mutant Hart98  shows higher atomic displacement parameters (B-factors) for the side-chains compared to the wild type, but the backbone conformation is not significantly different, thus indicating greater molecular flexibility in some portions of the structure. It was recently shown that the for some mutants, the apo-state (without metals) of the mutants is markedly destabilized compared to that of the wild-type, while the holo-state (with metals) of both mutants and wild-type is not significantly affected Lindberg02 . It is, therefore, possible that the apo-state of SOD1 — the subject of our studies — is implicated in misfolding and aggregation.

Most Cu, Zn superoxide dismutases, including human SOD1, have been shown to undergo irreversible aggregation on exposure to temperatures higher than their respective melting temperatures. The presence of free cysteines — those that do not form disulfide bonds in the native state — is known to be involved in this phenomenon. It is well-established that improper disulfide bond formation, concomitant with cysteine oxidation, possibly enhanced by metals in SOD1, is a cause of heat-denatured aggregation Perry87 . In the case of bovine SOD1, substitution by site-directed mutagenesis of free cysteine residues was found to greatly increase the reversibility of denaturation without substantially affecting the conformational stability McRee90 ; Lepock90 . Yeast SOD1 Lepock85  and E. coli SOD1 Battistoni98 , neither of which contain free cysteines, show partial and full reversibility respectively upon denaturation. However, SOD1 from P. leiognathi also denatures irreversibly despite having no free cysteines Bourne96 . Clearly, the presence of free cysteines is one of the crucial factors (but not the only factor) responsible for irreversible denaturation. One plausible reason for the effect of free cysteines is that the formation of improper disulfide bonds makes SOD1 more prone to kinetic traps. Therefore, we also study the effect of disulfide bonds on SOD1 folding reversibility.

An atomic resolution simulation of the folding process using traditional molecular mechanics force-fields is difficult by direct computational approaches because of the vast dimensionality of the protein conformational space Karplus94a . Simplified models such as the Gō model go81 ; Abe81  provide a powerful alternative to study folding because of their ability to simulate folding on reasonable time scales and to reproduce the basic thermodynamic and kinetic properties of proteins both on- Shimada01 ; Clementi00 ; Abkevich94  and off-lattices Dokholyan98b ; Dokholyan02b . In the Gō model, the energy of the protein is expressed as a sum of pairwise native contact energies. A native contact exists between amino acids if they are closer to each other than a cut-off distance in the native state. Thus, in the Gō approximation, folding can be construed as a transition from a state with no or a small number of native contacts to a state with all native contacts. Despite their apparent simplicity, native topology-based approaches have yielded results in agreement with experiments, in particular for small proteins ($`<100`$ amino acids) Borreguero02 ; Ding02a  reflecting the underlying simplicity in the folding of these proteins. However, a Gō model of longer proteins may lead to discrepancy in protein folding dynamics between simulation and experiment because sequence-specific interactions may significantly alter the protein folding dynamics. Thus, we develop methodology to incorporate sequence specific information of the relative contribution of various amino-acids into the Gō model of SOD1.

## II Results

We perform DMD simulations using the unscaled Gō model, to determine the thermodynamic properties of a model of SOD1 (as described in Methods). The temperature dependence of average potential energy is shown in Fig. 1a. There is a sigmoidal increase in the potential energy with increasing temperature. At low temperatures, SOD1 is present mostly in its native state, as evidenced by the small ($`2.83`$ Å) root-mean-square deviation, $`RMSD`$, from the native state (Fig. 1b), and the small fluctuations in the radius of gyration, $`R_g`$, around its native state value (14.2 Å) (data not shown). At high temperatures, SOD1 is unfolded with the $`RMSD`$ from the native state being greater than 40 Å, demonstrating the loss of any structural similarity with the native state. Also, $`R_g`$ is approximately three times its value in the native state (data not shown), indicating that the average distance between any two amino acids is approximately three times that in the native state.

In a typical trajectory (Fig. 2a) at the temperature corresponding to the mid-point of the abrupt change in potential energy (Fig. 1a) dependence on temperature (defined as the the folding transition temperature $`T_F`$), we observe three distinct states: folded, unfolded and intermediate. A histogram of potential energies of protein conformations near $`T_F`$ is, therefore, tri-modal (Fig. 2c), instead of the expected bi-modal for a two-state protein, indicating the presence of a folding intermediate. The detection of the intermediate state is in contrast with the experimental observation of the two-state folding dynamics of the SOD1 monomer Stroppolo00 . Thus, the unscaled Gō model for SOD1 does not accurately reproduce the folding dynamics of SOD1. Next, we calculate the pairwise interaction energies in the native state of SOD1 monomer and then employ the scaled Gō model, to uncover SOD1 folding thermodynamics and kinetics.

A few contacts make a large thermodynamic contribution to SOD1 native state stability. We decompose the stability of SOD1 into pairwise contributions (as described in Methods) that are plotted in the lower triangle in Fig. 3 (the energy map). In the energy map, some contacts, corresponding to both short-range and long-range interactions, have significantly stronger CHARMM interaction energies (from three to five times) than average. We find that there are 16 long-range interactions that have large (about five times larger than the average) interaction energies: Lys3-Glu21, Lys3-Glu153, Lys30-Glu100, Ser34-Asp96, Lys36-Asp92, Glu40-Lys91, Glu49-Arg115, Arg69-Glu77, His71-Asp83, Arg79-Asp83, Arg79-Asp101, His80-Asp83, Asp101-Val103, Glu121-Ser142, Asp125-Lys128, Glu133-Lys136. The residues involved in the strong long-range contacts are highlighted in Fig. 4 and are involved in contacts between oppositely charged amino acid pairs on the surface of the protein. Out of the 16 identified interactions, 6 have been found to be directly disrupted in patients with FALS as a result of point mutations (listed in Table I) in the residues making these contacts 222http://www.alsod.org. In most cases, the mutations correspond to the change of a charged amino acid for an oppositely charged or a neutral one, suggesting an important role these key contacts play in SOD1 folding.

Scaling thermodynamically important contacts modulates folding kinetics in the Gō model. We test the importance of the identified amino acids in DMD simulations by scaling their contact strength in the Gō potential. We increase the depth of the interaction well of the identified key contacts to five times their original value in the unscaled Gō model ($`b=5`$). Using this scaled Gō model, we perform DMD simulations to calculate the thermodynamic properties of the SOD1. We find that the potential energy, $`RMSD`$ and $`R_g`$ (data not shown) exhibit the same sigmoidal trend as for the unscaled Gō model, but the folding transition temperature increases to approximately 1.1 from its earlier value of 1.0. Remarkably, a typical trajectory near the folding transition temperature ($`T_F=1.15`$) shows (Fig. 2b) that the intermediate state is no longer populated and the protein only exists in either the folded or the unfolded state (Fig. 2b and d) near $`T_F`$. Thus, by strengthening these key contacts we reproduce the experimentally observed two state folding of SOD1 monomer Mei92 ; Stroppolo00 . As control, we randomly select several combinations of sixteen contacts and scale them by the same procedure. In such control simulations, we do not observe such effect on the folding kinetics (Fig. 2b), thus demonstrating that the effect of strengthening the identified contacts is specific. Importantly, this result uncovers a crucial connection between the folding thermodynamics and kinetics of SOD1: in our model, the residues that contribute the most to the SOD1 thermodynamic stability of the native state are also crucial for keeping the folding kinetics two-state.

Improper disulfide bonds cause kinetic trapping of SOD1. Human SOD1 is known to aggregate irreversibly on heat denaturation Rodriguez02  in vitro. One of the reasons for the irreversibility is improper disulfide bond formation between free cysteines. There are four cysteine residues in the SOD1 monomer, at positions 6, 57, 111 and 146 of which residues 57 and 146 form a disulfide linkage in the native state. We study the effect of non-native disulfide bond formation by starting from a fully unfolded conformation, generated at T=2.3 ($`T_F`$), and cooling the system to low temperature T=0.4 ($`T_F`$) slowly, while allowing the possibility of interactions between any of the four cysteines to form improper disulfide bond(s). We model disulfide bonds by an attractive potential ($`b=2`$, $`\gamma =1`$ in Eq. (2)). As a control, we repeat the annealing procedure without these non-native disulfide bonds ($`b=1`$, $`\gamma =1`$ in Eq. (2)). We observe that out of 20 trajectories with improper disulfide interactions, 4 result in SOD1 structures that are trapped in non-native conformations while all of the 20 control trajectories fold to the correct native structure on the same time scale, implying that occurrence of improper disulfide linkage increases the number of kinetic traps for the SOD1 chain en route to its native state. The above observation is statistically significant (the probability of observing this result by chance: $`p`$-value $``$0.01). The increased number of kinetic traps makes intermediate states more populated and interaction of these intermediate states may make SOD1 prone to aggregation. We conclude that free cysteine residues play an important part in the kinetics of folding, and possibly, aggregation.

## III Discussion

Experimental evidence suggests that the unfolding of dimeric superoxide dismutase occurs in a three state equilibrium Mei92 ; Stroppolo00 .

$$D2MU$$
(1)
where D, M and U are SOD1 dimer, monomer and unfolded state respectively. In order to understand the effect of mutations on SOD1 folding, it is important to study the effect of mutation on both processes in Eq. (1). Here, we study the second process ($`MU`$) and report thermodynamically and kinetically important residues for SOD1 monomer folding.

Using our hybrid method, which incorporates molecular mechanics energetics into the Gō model, we identify the key residues involved in the folding of the protein. Surprisingly, the key residues are found on the surface of the protein and not in its hydrophobic core (Fig. 4). By comparing the identified residues with a comprehensive list of mutations associated with FALS 333http://www.alsod.org, we find that, 6 out of the 16 contacts identified in our simulations are found to be implicated in the disease. In Table I, we list those FALS mutations, that disrupt contacts identified from our simulations. We find that these mutations replace a charged residue in the wild-type protein to a neutral or oppositely charged residue in the mutant, thereby disrupting the electrostatic interaction. The presence of mutations disrupting these contacts further underscores the importance of these key contacts and indicates that the disruption of even one of these contacts may potentially result in protein misfolding.

There have been no experimental studies characterizing SOD1 monomer folding kinetics. However, other members of the immunoglobulin-like (Ig-like) superfold family Murzin95 , to which SOD1 belongs, are experimentally well-characterized. Clarke and colleagues have established that the folding of pathways of Ig-like proteins share some common features Clarke99 . They observe that for five structurally similar Ig-like proteins, the folding rates are correlated with thermal stabilities, suggesting a similar folding pathway for all members of the family. Thus, for Ig-like proteins, the kinetic importance of residues in folding ($`\varphi `$-values Fersht95 ) is largely determined by the topology of the fold. We, therefore, expect that the residues important for SOD1 folding are in structurally equivalent positions to other members of the Ig-like super-fold family, whose folding kinetics is experimentally well-studied Clarke99 . We choose Tnf3, the fibronectin type III domain of human tenascin, for which, based on $`\varphi `$-values, eight residues in the hydrophobic core were identified as part of the folding nucleus Hamill00 . We identify residues in SOD1 that are structurally equivalent to the nucleus residues in Tnf3 by performing a structural alignment of the two proteins using DALI Holm95  and list these residues in Table II. We discover that two of the residues corresponding to this putative nucleus correspond to the strongly interacting key residues identified by our simulations, while three others are nearest neighbors of other residues we identify. This is significant because (i) mutations in the nucleus residues, being in proximity to the identified key residues, may disrupt the network of interactions of key residues and affect folding kinetics and (ii) in the Gō model entropic contributions are included as effective energetic interactions. Therefore, the identified key residues may be kinetically crucial because of the entropic contribution of the corresponding region of the structure (which is reflected by high $`\varphi `$-values of the neighboring nucleus residues). Therefore, kinetic evidence further supports our hypothesis that residues important for SOD1 native thermodynamics are important for folding kinetics.

Richardson and Richardson have pointed out Richardson02  that the interactions between edges of $`\beta `$-strands are crucial to regulate protein aggregation propensity. Naturally occurring $`\beta `$-sheet proteins prevent aggregation by protecting the edges of $`\beta `$ strands through a multitude of fold-specific mechanisms. In particular, for Ig-like $`\beta `$-sandwiches, strategically placed side-chain charges at the edges and loop-crossing are dominant mechanisms which shield the edges of the sheets. A recently solved crystal structure of mutant apo-SOD1 supramolecular assembly Elam03  identified three regions of the SOD1 structure where non-native gain of interaction was observed leading to amyloid like arrangement of SOD1 in crystals: the cleft between strands S5 and S6 (residues 83-100), the zinc loop (residues 65-78) and the electrostatic loop (residues 125-142). It was noted that disorder in the zinc and electrostatic loops promoted non-native interactions with the S5-S6 cleft, which is the interface between the two $`\beta `$-sheets of SOD1 and hence a likely site for edge-to-edge aggregation Richardson02 . We find that most of the key contacts for proper folding are made by residues in the three regions identified above. In particular, we see that contacts on both the edge-strands S5 and S6 play a crucial role in folding kinetics. The contacts made by residues Asp92, Asp96 and Lys91 which form the turn in the S5-S6 hairpin are crucial to seal the interface between the two $`\beta `$-sheets and likely prevent gain of non-native interactions, as shown in Fig. 5. Thus, our simulations constitute another piece in the growing body of evidence which implicates edge protection as the dominant mechanism by which natural proteins avoid aggregation. We postulate that the disruption of key interactions is a plausible scenario for mutant mediated aggregation of SOD1. The re-organization of the S5-S6 edge may be the dominant feature which distinguishes the aggregation-prone misfolded state of SOD1, from the native state. The disruption of charge-charge interactions of the side-chains of the identified key residues, which stabilize the $`\beta `$-barrel architecture of SOD1, might expose the edges of $`\beta `$-strand main-chain to interact with other structures, thus leading to aggregation. Wang and Hecht Wang02  have shown that strategically designed mutations in $`\beta `$-strand overhangs disfavor aggregation of normally aggregating proteins, further supporting our hypothesis.

We show that the presence of oxidisable cysteines enhances the kinetic trapping during folding. In the presence of metals, cysteine oxidation is enhanced suggesting a plausible cause for the misfolding of SOD1. However, it has recently been suggested Lindberg02  that misfolding of SOD1 likely takes place in the apo-enzyme state, when the metals have not yet been delivered to SOD1. In the reducing environment of the cell, cysteine-mediated misfolding/aggregation of the apoenzyme is, therefore, not a likely mechanism of SOD misfolding/aggregation. However, our observations do suggest a generic mechanism for SOD1 misfolding especially for mutations that do not affect the native state conformational stability Rodriguez02 . In the scenario similar to the one we have presented for free cysteines, other mutations may affect the folding kinetics of SOD1 by introducing non-native interactions along the folding pathway, leading to kinetic trapping without substantially affecting the native state stability.

## IV Methods

The unscaled Gō model. We use the Discrete Molecular Dynamics (DMD) algorithm Dokholyan98b ; Zhou97 ; Smith97  to study the folding thermodynamics of SOD1. DMD has recently been used to study the folding kinetics Borreguero02 ; Ding02a ; Zhou99a ; Zhou99b  and aggregation of proteins Ding02b ; Smith01 . We model the SOD1 chain by the beads-on-a-string model developed by Ding et al. Ding02a  with beads corresponding to all $`C_\alpha `$ and $`C_\beta `$ atoms and constraints between neighboring beads to mimic real protein flexibility. We study the folding of one monomer derived from the crystal structure of SOD1 (Protein data bank access code: 1SPD). We use the Gō potential go81  to model the interaction energy, $`V_{ij}`$, between $`C_\beta `$ atoms ($`C_\alpha `$ for Gly) $`i`$ and $`j`$ of SOD1:

$`V_{ij}=\{\begin{array}{cc}+\mathrm{},\hfill & \text{ }|r_ir_j|r_0\text{ }\hfill \\ \gamma u_{ij},\hfill & \text{ }r_0<|r_ir_j|r_1\text{ }\hfill \\ 0,\hfill & \text{ }|r_ir_j|>r_1\text{ }\hfill \end{array},`$ (5)
where $`|r_ir_j|`$ is the distance between atoms $`i`$ and $`j`$. The parameters $`r_0`$ and $`r_1`$ are the hard-core diameter (set to 3.25 Å) and the cut-off distance (7.5 Å) between $`C_\beta `$ atoms ($`C_\alpha `$ for glycine) in the native state respectively. We assign attractive interaction for native contacts ($`\gamma =1`$), defined to exist between pairs of residues whose $`C_\beta `$ atoms ($`C_\alpha `$ for glycine) are closer than 7.5 Å in the native state, and repulsive interaction ($`\gamma =1`$) for non-native contacts, defined to exist between $`C_\beta `$ atoms ($`C_\alpha `$ for glycine) which are farther than 7.5 Å in the native state.

The depth of the attractive square-well $`u_{ij}`$ in Eq. (5) determines the strength of the attractive interaction between two amino acids. First, we consider the case where all contacts are equally strong, i.e. the depth of each attractive well is $`u`$. We call this case the unscaled Gō model and use it to perform DMD simulations of SOD1. We find that our model undergoes a collapse transition to its native state through a meta-stable intermediate (see Results) in contrast with the experimental observation of two-state folding dynamics of SOD1 monomer Mei92 ; Stroppolo00 . A possible reason for the existence of meta-stable intermediates is that in the unscaled Gō model, all native contacts have equal strength. While such an approximation holds for smaller proteins (such as Src SH3 Borreguero02 ) with relatively low number of contacts ($`162`$), SOD1 monomer has 495 native contacts and the assumption of equal contact strength overestimates the population of intermediate species. Thus, in order to include sequence specificity of native interactions, we increase the strength of few specific native contacts according to an all-atom model of SOD1 described below by making the depth of the attractive well of some chosen native contacts $`bu_{ij}`$, where $`b`$ is the scale factor. We call this the scaled Gō model. By scaling Gō potential according to an all atom model of the protein, we effectively include the sequence information.

The scaled Gō model. To determine the relative strength of each individual native contact in SOD1, we use an effective energy function to decompose the total energy of the protein as a sum of pairwise terms Paci02 :

$$E=\underset{i}{}\underset{j>i}{}E_{ij},$$
(6)
where the summation is taken is over all pairs of residues of SOD1 and $`E_{ij}`$ is the energy of interaction between residues $`i`$ and $`j`$. To determine the terms $`E_{ij}`$, we perform simulations with the molecular mechanics force field CHARMM Brooks83  to construct an all atom model of the SOD1 monomer in the native state, starting from its crystal structure (PDB accession code:1SPD). Recently, higher resolution crystal structures of apo-SOD1 have appeared in the literature Strange03 , but our results are not qualitatively affected by starting with the higher resolution structure. We add hydrogen atoms and remove any steric clashes occurring in the structure by energy minimization. The interactions which contribute the most to the relative contact interaction strength are the non-bonded interactions and the solvent-induced interaction screening Paci02 :

$$E_{ij}=E_{ij}^{nonbonded}+G_{ij}.$$
(7)
We compute $`E_{ij}^{nonbonded}`$ as the sum of the Lennard-Jones and electrostatic interaction energies. The solvation free energy, $`G_{ij}`$, measures the screening of the pairwise interaction due to polarization of the solvent.

We calculate $`G_{ij}`$ using an implicit description of the solvent Sharp94 . We approximate the solvent as a continuous dielectric medium with a dielectric constant ($`ϵ=80`$) of water and the protein as a continuous medium with dielectric constant ($`ϵ=1`$) of vacuum. We compute the solvation energy by solving the Poisson-Boltzmann equation using the finite difference method Nicholls91  available in CHARMM. Following Dominy et al. Dominy02 , in order to decompose the solvation energy into pairwise contributions, we eliminate all charges in the protein except those corresponding to residue $`i`$ and calculate its “self” solvation energy $`U_i`$. Similarly, we compute $`U_j`$ and $`U_{ij}`$ when charges corresponding to both $`i`$ and $`j`$ are present. Then, the solvent mediated screening of the interaction is:

$$G_{ij}=U_{ij}[U_i+U_j].$$
(8)
So, we calculate the effective pairwise energy contribution of a contact in the folded state.

The $`E_{ij}`$ values that we thus obtain correspond to the folded state of the protein. In order to calculate the pairwise contribution of residues to the stability, we subtract the value of the contact energy in the unfolded state from the $`E_{ij}`$ values obtained above. The unfolded state is an ensemble of conformations. It has been noted that there is considerable native like structure in the unfolded state Zagrovic02 ; Shortle01 ; Elcock99 , meaning that the average environment of a contact is likely the same as in the folded state, but this contact is formed less frequently. Therefore, we model the unfolded state by simulating the unscaled Gō model near $`T_F`$ and estimating the frequency of native contact formation. The total contribution of a given contact to the stability can be written as:

$$\mathrm{\Delta }E_{ij}=E_{ij}(1p_{ij})$$
(9)
where $`E_{ij}`$ is the contact energy in the folded state and $`p_{ij}`$ is the probability of forming the contact in the unfolded state. The energy and contact maps are plotted in Fig. 3.

## V Conclusions

We present the folding thermodynamic and kinetic analysis of SOD1 using a hybrid molecular dynamics method. We show that our simple model qualitatively reproduces the folding of SOD1 observed in experiments. We find that interactions at the edges of SOD1 $`\beta `$ sheets and the zinc, electrostatic and cross-over loops are the key interactions that modulate the folding of SOD1 suggesting that mutants destabilizing these interactions may make the enzyme more prone to aggregation. In particular, we identify the S5-S6 cleft (residues 83-100) as a likely site for reorganization to yield an aggregation prone misfolded state, an observation consistent with recent experimental evidence Elam03  and the negative design paradigm Richardson02 . We also provide a plausible explanation for the phenomenon of increased reversibility of folding upon mutation of free cysteine residues. We find good agreement between the positions identified by us to be crucial for folding kinetics and those identified by experiments Hamill00  on other members of the super-fold family (Ig-like proteins) to which SOD1 belongs. Further, we find that mutations at six of the residue contacts we identify are implicated in FALS.

## VI Acknowledgements

We thank S. V. Buldyrev for his help with DMD simulations, B. N. Dominy for help with CHARMM calculations and J. M. Borreguero, L. J. Hayward, B. Kuhlman and Z. Xu for helpful discussions. We acknowledge the support of the UNC-CH Research Council Grant. SDK acknowledges the support of UNC Molecular and Cellular Biophysics Training Program and Scholars for Tomorrow Fellowship of UNC Graduate School.


text:
# Selective Transport and Mobility Edges in Quasi-1D Systems with a Stratified Correlated Disorder

## Abstract

We present analytical results on transport properties of many-mode waveguides with randomly stratified disorder having long-range correlations. To describe such systems, the theory of 1D transport recently developed for a correlated disorder is generalized. The propagation of waves through such waveguides may reveal a quite unexpected phenomena of a complete transparency for a subset of propagating modes. We found that with a proper choice of long-range correlations one can arrange a perfect transparency of waveguides inside a given frequency window of incoming waves. Thus, mobility edges are shown to be possible in quasi-1D geometry with correlated disorder. The results may be important for experimental realizations of a selective transport in application to both waveguides and electron/optic nanodevices.

During last few years there is a burst of interest to the problem of localization-delocalization transition in systems with correlated disorder (see, e.g.,  and references therein). This fact is due to the possibility to observe an anomalous transport in 1D models with random potentials. In particular, it was shown  that specific long-range correlations give rise to a complete transparency of electron waves for given energy intervals. Experimental realization of such potentials for single-mode waveguides with delta-like scatters  has confirmed theoretical predictions.

The subject of wave propagation through disordered many-mode waveguides is important both from the theoretical viewpoint and for experimental applications such as optic fibers, remote sensing, radio wave propagation, shallow water waves, etc (see, for example, ). It has also a direct relation to the problem of electronic transport in mesoscopic conducting channels. So far, main results in this field have been obtained for random potentials of white-noise type. However, there are many physical situations where this assumption is not correct. Therefore, the understanding of generic properties of transport in the models with correlated disorder is important for the modern theory of electron/wave propagation. It should be also stressed that existing experimental technics allow for the construction of systems with sophisticated scattering potentials resulting in anomalous transport properties .

In this paper we study transport properties of quasi-1D waveguides with a stratified bulk disorder that has long-range correlations. Our main interest is in exploring the possibility of constructing such random potentials that result in frequency windows of a perfect propagation of waves.

In what follows we consider a plane waveguide (or conducting wire) of width $`d`$, stretched along the $`x`$-axis. The $`z`$-axis is directed perpendicularly to the waveguide so that its lower edge is $`z=0`$ and the upper one is $`z=d`$. The waveguide is assumed to have a stratified (along the $`x`$-axis) region of the length $`L`$ ($`|x|<L/2`$). As a physically plausible model for such a disorder we consider random potential $`V(x)`$ with the zero average, $`V(x)=0`$. The angular brackets here and below stand for the statistical average over realizations of the random function $`V(x)`$. Its correlator $`V(x)V(x^{})=𝒲(|xx^{}|)`$ is assumed to have a maximal value at $`x=x^{}`$ and to decrease with an increase of $`|xx^{}|`$. For statistical treatment of the problem it is also important to assume that the characteristic scale of the decrease of $`𝒲(|xx^{}|)`$ is much less than the length $`L`$ of the scattering region.

We shall characterize transport properties of the waveguide by its transmittance $`T(L)`$ that within the linear response theory is defined by the Kubo’s formula ,

$$T(L)=\frac{4}{L^2}𝑑\stackrel{}{r}𝑑\stackrel{}{r^{}}\frac{𝒢(\stackrel{}{r},\stackrel{}{r^{}})}{x}\frac{𝒢^{}(\stackrel{}{r},\stackrel{}{r^{}})}{x^{}}.$$
(1)
Here the integration with respect to $`\stackrel{}{r}=(x,z)`$ runs over the disordered region. Since the potential $`V(x)`$ does not depend on $`z`$, the retarded Green function has the form,

$$𝒢(\stackrel{}{r},\stackrel{}{r^{}})=\frac{2}{d}\underset{n=1}{\overset{N_d}{}}\mathrm{sin}\left(\frac{\pi nz}{d}\right)\mathrm{sin}\left(\frac{\pi nz^{}}{d}\right)G_n(x,x^{}).$$
(2)
Here $`N_d=[kd/\pi ]`$ is the total number of propagating waveguide modes (conducting channels) determined by the integer part $`[\mathrm{}]`$ of the ratio $`kd/\pi `$. The wave number $`k`$ is equal to $`\omega /c`$ for a classical scalar wave of frequency $`\omega `$, and to the Fermi wave number for electrons.

By substituting Eq. (2) into Eq. (1), we come to the Landauer’s formula ,

$$T(L)=\underset{n=1}{\overset{N_d}{}}T_n(L),$$
(3)
with $`T_n(L)`$ being obtained from Eq. (1) in which the integration is performed over $`x`$, $`x^{}`$ only, and $`𝒢`$ is replaced by $`G_n`$. The mode Green function $`G_n(x,x^{})`$ is entirely determined by the 1D equation with the random potential $`V(x)`$ and lengthwise wave number $`k_n=\sqrt{k^2(\pi n/d)^2}`$,

$$\left[\frac{d^2}{dx^2}+k_n^2V(x)\right]G_n(x,x^{})=\delta (xx^{}).$$
(4)
From Eq. (3) one can see that the total transmittance $`T(L)`$ of a quasi-1D stratified structure is expressed as a sum of partial transmittances $`T_n`$. Every transmittance $`T_n`$ describes the transparency of the corresponding $`n`$th propagating mode. In such a way we reduce the transport problem for the quasi-1D disordered system to the 1D random model (4) of wave scattering in every $`n`$th channel. Eq. (4) can be solved by one of the methods developed in the transport theory of one-dimensional disordered systems, such as the perturbative diagrammatic method of Berezinskii , invariant imbedding method , or the two-scale approach .

The main result is that the average partial transmittance $`T_n`$ as well as all its moments, are described by the universal function that depends on one parameter $`\mathrm{\Lambda }_nL/L_{loc}(k_n)`$ only (one-parameter scaling). The quantity $`L_{loc}(k_n)/4`$ is, in fact, the backscattering length for $`n`$-th propagating mode. The inverse value $`L_{loc}^1(k_n)`$ is equal to the corresponding Lyapunov exponent that can be found in the transfer matrix approach, so that the quantity $`L_{loc}(k_n)`$ is the localization length  associated with a given 1D channel.

Since general expression for the mode transmittance has quite complicated form (see, e.g., Ref. ), here we refer to limit cases only. Specifically, for large localization length, $`\mathrm{\Lambda }_n=L/L_{loc}(k_n)1`$, the transmittance $`T_n`$ exhibits the ballistic behavior and the corresponding $`n`$th normal mode is practically transparent,

$$T_n14\mathrm{\Lambda }_n.$$
(5)
On the contrary, the transmittance is exponentially small when the localization length is much less than the length of the waveguide, $`\mathrm{\Lambda }_n1`$,

$$T_n\frac{\pi ^{5/2}}{16}\mathrm{\Lambda }_n^{3/2}\mathrm{exp}(\mathrm{\Lambda }_n).$$
(6)
This implies strong wave localization in a given $`n`$th channel. Since in this case the transmittance $`T_n`$ is not a self-averaged quantity, it is more reasonable to deal with its average logarithm, $`\mathrm{ln}T_n=4\mathrm{\Lambda }_n`$.

Thus, in order to describe transport properties of a quasi-1D waveguide with 1D bulk disorder, the value of localization lengths $`L_{loc}(k_n)`$ should be known. From the solution of the equation (4) one can obtain ,

$$L_{loc}^1(k_n)=\frac{W(2k_n)}{16k_n^2},$$
(7)
where $`W(k_x)`$ is the Fourier transform of the binary correlator $`𝒲(|xx^{}|)`$ for the scattering potential,

$$𝒲(|xx^{}|)=_{\mathrm{}}^{\mathrm{}}\frac{dk_x}{2\pi }\mathrm{exp}[ik_x(xx^{})]W(k_x).$$
(8)
The main feature of the expression (7) for $`L_{loc}(k_n)`$ is its dependence on the mode index $`n`$. One can see that the larger $`n`$ is, the smaller is $`L_{loc}(k_n)`$ and, consequently, the stronger is the coherent scattering within this mode. This dependence is quite strong due to the squared wave number $`k_n`$ in the denominator of Eq. (7). Evidently, with an increase of the index $`n`$ the value of $`k_n`$ decreases. An additional dependence appears because of the power spectrum $`W(2k_n)`$. Since the correlator $`𝒲(|xx^{}|)`$ is a decreasing function of $`|xx^{}|`$, the numerator $`W(2k_n)`$ increases with $`n`$ (it is a constant in the case of the delta-correlated potential only). Therefore, both the numerator and denominator contribute in the same way for the dependence of $`L_{loc}(k_n)`$ on $`n`$. As a result, we arrive at the concept of hierarchy of mode localization lengths,

$$L_{loc}(k_{N_d})<L_{loc}(k_{N_d1})<\mathrm{}<L_{loc}(k_2)<L_{loc}(k_1).$$

Thus, a remarkable phenomenon arises. On the one hand, the concept of one-parameter scaling holds for any of $`N_d`$ channels whose partial transport is characterized solely by the value of $`\mathrm{\Lambda }_n`$. On the other hand, this concept turns out to be broken for the total waveguide transport. Indeed, due to the revealed hierarchy of $`L_{loc}(k_n)`$, the total transmittance (3) depends on the whole set of scaling parameters $`\mathrm{\Lambda }_n`$. This fact is in contrast with quasi-1D bulk-disordered models, for which all transport properties are described by one parameter only. In this sense, our model is similar to quasi-1D waveguides with rough surfaces, for which the similar hierarchy of attenuation lengths was recently found .

The interplay between the hierarchy of $`L_{loc}(k_n)`$ and the one-parameter scaling for $`T_n`$ leads to that, in general, quasi-1D stratified structures exhibit three transport regimes.

(I) If the largest of the localization lengths, $`L_{loc}(k_1)`$, is much less than the waveguide length $`L`$,

$$L_{loc}(k_1)L,$$
(9)
all the propagating modes are localized and the waveguide is non-transparent.

(II) On the contrary, when the smallest localization length $`L_{loc}(k_{N_d})`$ is much larger than $`L`$,

$$L_{loc}(k_{N_d})L,$$
(10)
all the propagating modes are open ($`T_n1`$) and the waveguide has almost perfect transparency. The total ballistic transmittance in this case is equal to the total number of the propagating modes, $`T(L)=N_d`$.

(III) The intermediate situation arises when $`L_{loc}(k_1)`$ of the first mode is larger, while $`L_{loc}(k_{N_d})`$ of the ”last” $`N_d`$-th mode is smaller than the waveguide length $`L`$,

$$L_{loc}(k_{N_d})LL_{loc}(k_1).$$
(11)
In this case a very interesting phenomenon of the coexistence of ballistic and localized transport occurs. Namely, while lowest modes are in the ballistic regime, highest modes are strongly localized.

As a demonstration, let us consider waveguides with a large number of conducting channels, $`N_d=[kd/\pi ]kd/\pi 1`$, and with potentials having the widely used Gaussian correlator,

$`𝒲(|x|)`$ $`=`$ $`W_0k_0\pi ^{1/2}\mathrm{exp}\left(k_0^2x^2\right),`$ (12)
$`W(k_x)`$ $`=`$ $`W_0\mathrm{exp}\left(k_x^2/4k_0^2\right).`$ (13)
It is convenient to introduce two parameters

$$\alpha \frac{L}{L_{loc}^w(k_1)}=\frac{W_0L}{16k^2},\delta \frac{L_{loc}^w(k_{N_d})}{L_{loc}^w(k_1)}=\frac{2\{kd/\pi \}}{(kd/\pi )}1,$$
(14)
where $`L_{loc}^w(k_1)`$ and $`L_{loc}^w(k_{N_d})`$ refer, respectively, to the largest and the smallest localization lengths in the limit case of the white-noise potential ($`k_0\mathrm{}`$), and $`\{kd/\pi \}`$ is the fractional part of the parameter $`kd/\pi `$.

By applying the above discussed approach, one can find that for Gaussian correlations (13) all the propagating modes are localized when $`\alpha \mathrm{exp}(k^2/k_0^2)`$. This inequality is stronger than that valid for the white-noise case, $`\alpha 1`$. The intermediate situation occurs when

$$\mathrm{exp}(k^2/k_0^2)\alpha \delta \mathrm{exp}(\delta k^2/k_0^2).$$
(15)
One can see that the longer range $`k_0^1`$ of the correlated disorder is, the simpler are the conditions (15) of the coexistence of ballistic and localized transport.

Finally, the waveguide is almost perfectly transparent in the case when $`\alpha \delta \mathrm{exp}(\delta k^2/k_0^2)`$. In particular, at any given value of $`\alpha 1`$ (fixed values of the waveguide length $`L`$, wave number $`k`$ and the disorder strength $`W_0`$) the ballistic transport of all the conducting channels can be realized by a proper choice of the stratified correlated disorder. From Eq. (15) one can also understand under what conditions the coexistence of ballistic and localized modes can be achieved for the standard white-noise case.

The fundamentally different situation arises when the stratified medium has specific long-range correlations. To show this, we would like to note that the localization length $`L_{loc}(k_n)`$ of any $`n`$th conducting channel are entirely determined by the dependence $`W(k_x)`$, see Eq. (7). Therefore, if $`W(2k_n)`$ abruptly vanishes within some interval of wave number $`k_n`$, then $`L_{loc}(k_n)`$ diverges and the corresponding propagating mode appears to be fully transparent for any length of the waveguide. Thus, it is naturally to ask how to construct such random profiles of stratified medium that result in a complete transparency of waveguides within any desired frequency region. The answer to this question can be found from a generalization of the methods developed in Refs. . Specifically, having a desirable dependence for the randomness power spectrum $`W(k_x)`$, one should obtain the function $`\beta (x)`$ whose Fourier transform is $`W^{1/2}(k_x)`$. Then, the random potential $`V(x)`$ can be constructed as a convolution of white noise $`Z(x)`$ with the function $`\beta (x)`$,

$$V(x)=_{\mathrm{}}^{\mathrm{}}𝑑x^{}Z(xx^{})\beta (x^{}).$$
(16)
We emphasize that the transition between localized and ballistic wave/electron transport can be arranged in an abrupt way at any given point inside the allowed interval for $`k_n`$. In order to achieve such a transition, it is convenient to take $`W(k_x)`$ with a discontinuity at the desired point. In other words, the random potential $`V(x)`$ should be of a specific form with long-range correlations along the waveguide.

To show how to arrange a sharp transition, let us take a random potential $`V(x)`$ with the following power spectrum $`W(k_x)`$,

$$W(k_x)=W_0\mathrm{\Theta }(|k_x|2k_0),$$
(17)
where $`\mathrm{\Theta }(x)`$ is the unit-step function and the characteristic wave number $`k_0>0`$ is a correlation parameter to be specified. According to Eq. (16), the random potentials having such correlator can be constructed as

$$V(x)=W_0^{1/2}\left[Z(x)_{\mathrm{}}^{\mathrm{}}𝑑x^{}Z(xx^{})\frac{\mathrm{sin}(2k_0x^{})}{\pi x^{}}\right].$$

Now one can see that in the case under consideration all low modes with wave numbers $`k_n>k_0`$ have finite localization lengths while for high modes with $`k_n<k_0`$ the localization lengths diverge,

$$L_{loc}^1(k_n>k_0)=W_0/16k_n^2,L_{loc}^1(k_n<k_0)=0.$$
(18)
Remarkably, in contrast to the potentials with Gaussian correlations (see above), all propagating modes with

$$n>N_{loc}=[(kd/\pi )(1k_0^2/k^2)^{1/2}]\mathrm{\Theta }(kk_0),$$
(19)
exhibit the ballistic behavior. Since their mode transmittance $`T_n=1`$, such modes form a coset of completely transparent channels. As for other propagating modes with low indices $`nN_{loc}`$, they remain to be localized for large enough length size $`L`$ when the condition (9) holds, i.e. for $`W_0L/16k_1^21`$.

The expression (19) determines the total number $`N_{loc}`$ of localized modes and total number $`N_{tr}=N_dN_{loc}`$ of completely transparent modes. Since localized modes do not contribute to the total transmittance (3), the latter is equal to the number $`N_{tr}`$ of completely transparent modes and do not depend on the waveguide length $`L`$,

$$T=[kd/\pi ][(kd/\pi )(1k_0^2/k^2)^{1/2}]\mathrm{\Theta }(kk_0).$$
(20)
We remind that square brackets stand for the integer part of the inner expression.

For $`k_0k`$ the number of localized modes $`N_{loc}[(kd/\pi )(1k_0^2/2k^2)]`$ is of the order of $`N_d`$. Consequently, the number of transparent modes $`N_{tr}`$ is small, or there are no such modes at all. Otherwise, if $`k_0k`$, the integer $`N_{loc}[\sqrt{2}(kd/\pi )(1k_0/k)^{1/2}]`$ turns out to be much less than the total number of waveguide modes $`N_d`$, and the number of transparent modes $`N_{tr}`$ is large. When $`k_0>k_1`$, the number $`N_{loc}`$ vanishes and all modes become fully transparent. In this case the correlated disorder results in a perfect transmission of waves. One can see that the value $`k_0`$ is, in essence, the total mobility edge that separates the region of complete transparency from that where lower modes are localized.

In conclusion, we have studied quasi-1D waveguides with a stratified disorder that has long-range correlations. We have shown that with a proper choice of correlations one can arrange a complete transparency of transverse channels in a given frequency window of incoming waves. It is worthwhile to emphasize a non-monotonic step-wise dependence of the total transmittance (20) on the wave number $`k`$. Specifically, within the region $`k<k_0`$ for fixed values of $`k_0`$ (when $`W_0Ld^2/16\pi ^2(k_0d/\pi )^21`$), all propagating modes are transparent and the transmittance exhibits step-wise increase with an increase of $`k`$. Each step up arises for an integer value of $`kd/\pi `$ when new conducting channel emerges. This effect is similar to that known to occur in quasi-1D ballistic (non-disordered) structures (see, e.g., ).

On the other hand, for $`k>k_0`$ the transmittance reveals a step down decrease due to successive localization of low modes. In contrast with the steps up, every step down occurs at the local mobility edge of the corresponding channel, when the second term in Eq. (20) takes the integer value. In general, these values do not coincide with the integer values of $`kd/\pi `$, thus resulting in a new kind of step-wise dependence for the transmittance. Finally, for very large values $`k^2W_0L/16`$, the transmittance starts to increase again, due to a successive delocalization of the modes. Our results can be used for fabrications of electron or optic nanodevices with a selective transport.

This research was supported by CONACYT grant 34668-E, and by BUAP grant II-104G02 (México).