text:
# On the Methods of Determining the Radio Emission Geometry in Pulsar Magnetospheres

## 1 Introduction

Since the discovery of pulsars (Hewish et al. 1968) the geometry of radio emission from pulsar magnetospheres was interpreted in terms of emission from purely dipolar magnetic fields (eg. Radhakrishnan & Cooke 1969; Cordes 1978; Lyne & Manchester 1988; Blaskiewicz et al. 1991; Rankin 1993; Gil & Kijak 1993). This assumption is justified by relatively high emission altitudes in comparison with $`R_{\mathrm{ns}}`$ inferred for the radio emission ($`0.01R_{\mathrm{lc}}`$, where $`R_{\mathrm{lc}}=cP/2\pi `$ is the light cylinder radius and $`P`$ is the rotation period of a neutron star with radius $`R_{\mathrm{ns}}`$). Most importantly, however, radio emission from the dipolar magnetic field hopefully can be described by a sufficiently small number of parameters and a limited number of observational parameters. The word “hopefully” reflects a second crucial assumption applied for the radio emission beam: that distinguishable features in pulse profiles (conal components, the outer edges of a profile) are associated with the beam structure (eg. of concentric hollow cones of enhanced radio emission) which in the reference frame corotating with the neutron star (CF) is symmetric with respect to the plane containing the dipole magnetic moment $`\stackrel{}{\mu }`$ and the rotation axis. Without this disputable assumption, the number of parameters required to determine the emission geometry increases considerably and one is left with a multi-parameter theory to be compared with data from which only a few parameters can be deduced. Hereafter, the axial symmetry (in CF) of the radio emission beam is assumed; the problem of its justification is presented in Section 4.1.

To further constrain the parameter space, these two assumptions (I – dipolar magnetic field; II – symmetry of radio beam) are often supplemented with two additional simplifications: III – it is assumed that in the CF the bulk of radio waves is emitted in the direction which is tangent to the local magnetic field at an emission point; IV – identifiable features in pulse profiles (eg. maxima of conal components) observed in a narrow frequency band are interpreted as radiation from a very narrow range of altitudes. The assumption No. III requires relatively large Lorentz factors $`\gamma `$ of radio emitting electrons. Hereafter we assume that the broadening of the pulse features due to the radiation from the low energy electrons is negligible, although it may be significant in reality. Our assumption is fully justified in the case of the most popular radiation mechanism – the curvature radiation, because the electron Lorentz factors must exceed $`100`$ for the curvature spectrum to extend up to the observation frequency ($`300`$ MHz). Curvature radiation from low energy electrons ($`\gamma 10`$) cannot broaden the profile because it does not extend to $`10^2`$ MHz. The model-dependent estimates of $`\gamma `$ cover a large range between $`50`$ and $`10^4`$ (eg. Ruderman & Sutherland 1975; Machabeli et al. 2001; Melrose 2004). For the rather conservative estimate of $`\gamma 100`$ the corresponding angular size of the single electron radiation beam is $`\gamma ^10.5^{}`$, which is an order of magnitude smaller than the typical width of the radio beam corresponding to the mean pulse profile.

With the assumptions I–IV, the radio emission geometry becomes completely determined by four parameters: $`\alpha `$, $`\beta `$, $`r`$, and $`\rho `$. Their meaning is the following: $`\alpha `$ is an inclination of the magnetic dipole with respect to the rotation axis (ie. the angle between the magnetic moment $`\stackrel{}{\mu }`$ and the angular velocity $`\stackrel{}{\mathrm{\Omega }}=2\pi P^1\widehat{z}`$). $`\beta `$ is an “impact angle”, ie. the minimum angle between an observer’s line of sight and $`\stackrel{}{\mu }`$. This parameter is often replaced by $`\zeta =\alpha +\beta `$ – the angle between the observer’s line of sight and the rotational axis. The next parameter, $`r`$, is the radial distance of the radio emission region measured from the center of the neutron star. Hereafter, we will often use its normalized value $`r^{}=r/R_{\mathrm{lc}}`$. The last parameter, $`\rho `$, is the half opening angle of the radio emission cone/beam, ie. it is the angle between the direction of radio emission in CF and $`\stackrel{}{\mu }`$. These four parameters determine completely the emission geometry in the sense that any additional information about the emission region (eg. coordinates of emission point in the dipole frame, with z-axis along $`\stackrel{}{\mu }`$) can be easily derived from textbook equations for the dipole geometry and spherical trigonometry (eg. Radhakrishnan & Cooke 1969; Cordes 1978; Gil et al. 1984).

One desired additional piece of information is an answer to the question “which magnetic field lines does the observed emission come from?”. Following standard conventions, we identify magnetic field lines by $`\theta _{\mathrm{fp}}`$ – the colatitude of foot points of the magnetic field lines at the neutron star surface. With the radio emission altitude $`r`$ and the beam radius $`\rho `$ determined, the value of $`\theta _{\mathrm{fp}}`$ can be easily calculated. The parameter $`\theta _{\mathrm{fp}}`$ is often expressed in terms of a fraction of the polar cap angle: $`\theta _{\mathrm{fp}}^{}=\theta _{\mathrm{fp}}/\theta _{\mathrm{pc}}`$, where $`\theta _{\mathrm{pc}}=\mathrm{arcsin}((R_{\mathrm{ns}}/R_{\mathrm{lc}})^{1/2})`$.

These four quantities are to be deduced from radio data. In most cases, however, a model-dependent analysis of the radio data at a given frequency provides us with only two quantities: $`\beta `$, and $`W(f)`$, where $`W(f)`$ is the apparent width of the radio pulse profile, usually defined by some simple criterion (eg. measured at some fraction $`f`$ of maximum intensity; a wide variety of $`f`$ is employed: $`f=0.0005`$ (Kijak & Gil 2003), $`f=0.1`$ (Blaskiewicz et al. 1991), $`f=0.5`$ (Rankin 1993), $`f=1`$ (Gangadhara & Gupta 2001)). The value of the inclination angle $`\alpha `$ could in principle be determined along with $`\beta `$ from the “rotating vector” model of polarization position angle swings (Radhakrishnan & Cooke 1969). In practice, however, a fit of the model to the observed position angle curve is much less sensitive to $`\alpha `$ than to $`\beta `$ (Rankin 1993). Only in exceptional cases can both of these parameters be derived (Lyne & Manchester 1988; Blaskiewicz et al. 1991; von Hoensbroech & Xilouris 1997). This led Rankin (1990) to a formulation of a method of determining $`\alpha `$ which assumes that the apparent width of the core component does not depend on the impact angle $`\beta `$ and that in all pulsars the core component originates from the same, low altitude. The observed width of the core component becomes then a universal function of the pulsar rotation period $`P`$ and the dipole inclination $`\alpha `$, allowing, thus, to determine the latter parameter.

We are provided, therefore, with just three parameters ($`\alpha `$, $`\beta `$, $`W`$) instead of four. This limitation prompted the development of two kinds of methods to determine the geometry of the radio emitting region: 1) A purely geometric method which assumes that the lowest detectable emission at the leading and at the trailing edge of a radio pulse originates at the last open magnetic field lines, with $`\theta _{\mathrm{fp}}=\theta _{\mathrm{pc}}`$. With the fourth parameter assumed a priori, the method makes it possible to determine the radial position of radio emission $`r`$ from the observed pulse width $`W_0`$ (eg. Cordes 1978; Gil & Kijak 1993; Kijak & Gil 1997; 1998; 2003; Kijak 2001; hereafter we will use the index ‘0’ to refer to the pulse with $`W(f)`$ at the lowest intensity level, ie. practically at $`f=0.00050.1`$). The derived altitudes are a few tens of $`R_{\mathrm{ns}}`$ at observation frequency $`\nu 1`$ GHz. 2) Methods of the second kind are able to derive the fourth observational parameter by a measurement of phase shifts of some profile features with respect to some fiducial points. The methods are relativistic in the sense that the phase shifts are caused by combined effects of aberration and propagation time delays due to the finite speed of light $`c`$ (for brevity, the latter effect will hereafter be called retardation). The second methods are superior to the geometric method in that they do not assume a priori the value of the fourth parameter. However, they must rely on additional assumptions about the radio pulse profile. Gangadhara & Gupta (2001, hereafter GG2001) measure the relativistic phase shift of conal components with respect to the core component, which is assumed to originate from much lower altitude than the cones which surround it. Blaskiewicz et al. (1991, hereafter BCW91) measure the shift of pulse edges with respect to the center (or the “inflection point”) of the position angle curve which is assumed to originate from the same altitude as the emission at the outer wings of the profile. On average, the method of Blaskiewicz et al. (1991) (hereafter BCW method) gives emission radii in rough agreement with that of the geometric method (Gil & Kijak 1993; Kijak & Gil 1997; 1998; 2003). However, the method of GG2001 predicts notably larger values of $`r`$ (by about factor of 2) than the geometric method. Unlike the geometric method (Cordes 1978; Gil & Kijak 1993), both relativistic methods (of GG2001 as well as BCW91) *in principle* make it possible to identify the radio emitting field lines. GG2001 and Gupta & Gangadhara (2003, hereafter GG2003) find $`\theta _{\mathrm{fp}}^{}`$ in the range $`0.220.74`$ for radio emission observed at maxima of conal components. Unfortunately, the method of GG2001 can be applied only for pulsars having unambigously identifiable pairs of conal components (and which possess a core component). So far, this has required application of a window thresholding technique (hereafter WT technique, GG2001) and limited application of the method only to a handful of the brightest objects. On the other hand, the BCW method suffers from a difficulty in finding the center of the position angle curve, and, as we show below, it yields values of $`\theta _{\mathrm{fp}}^{}`$ exceeding 1. For convenience, hereafter we refer to the relativistic methods of GG2001 and BCW91 with the term “relativistic phase shift methods” (RPS methods).

In Section 2 we revise the method of GG2001, which results in a new formula for radio emission altitudes and furnishes the method with new interesting features. As noted in GG2003, the original method of GG2001 yielded radio emission altitudes larger than those derived with the geometric method. Our revision removes part of this discrepancy. In Section 3 we propose to apply the method to the outer edges of averaged pulse profiles. This will hopefully provide a test of the main assumption of the geometric method about the value of $`\theta _{\mathrm{fp}}^{}=1`$ for the beam edge. As an example, we try to perform such a test using the method of BCW91. In Section 4 we discuss in detail the assumptions of the relativistic methods. Our conclusions are in Section 5.

## 2 The modified relativistic phase shift method

The RPS method of Gangadhara & Gupta (2001) applies to pulsars with both core and cone emission, ie. for M and T pulsars in the classification scheme of Rankin (1983). Questionability of the conal pattern of radio emission beam will be discussed in Section 4.1. The model assumes that in the reference frame corotating with the neutron star the cones are symmetric around the narrow core beam. Obviously, in the case of the dipolar magnetosphere the second assumption can be most naturally accounted for by a model of the radio beam consisting of a narrow core beam centered at the dipole axis and a few nested hollow cones of enhanced radio emission surrounding the core. Magnetic field lines corresponding to a given cone have the same $`\theta _{\mathrm{fp}}`$ and the same half-opening angle $`\rho `$ at a fixed altitude. In fact, the second assumption can be made less stringent: it is sufficient for the conal beams to be symmetric with respect to the plane containing the core component and the rotation axis. Thus, the “conal” beams may have eg. elliptical crossection with the core beam located in the plane containing the ellipse center and the rotation axis (ie. not necessary at the ellipse center). The possibility of the elliptical beam is discussed eg. in Mitra & Deshpande (1999).

The method starts with identifying the core component in the radio pulse profile. Then one attempts to identify the same number of conal components on both the leading and the trailing side of the core component. The innermost pair of two conal components located on both sides of the core is identified as emission from the same, innermost hollow cone of the radio beam. The next pair of conal components which bracket/flank the innermost pair is interpreted as emission from the next-to-innermost hollow cone of the radio beam, and so on. For each conal pair the observed phase of the leading component ($`\varphi _\mathrm{l}`$) and the phase of the trailing component ($`\varphi _\mathrm{t}`$) is measured with respect to the core component which defines the phase $`\varphi =0`$. As emphasized by GG2003 in all cases with clear identification of cones, the conal components were shifted toward earlier phases with respect to the core component, ie. $`|\varphi _\mathrm{l}|>\varphi _\mathrm{t}`$. This effect was noticed already by Gil (1985), who proposed the retardation effect to explain the asymmetry. More correctly, GG2001 interpreted the asymmetry as a forward shift of conal components due to combined effects of aberration of conal emission and retardation of core emission. The (negative) relativistic phase shift $`\mathrm{\Delta }\varphi _{\mathrm{obs}}`$ between the midway point of conal pairs and the core component is equal to:

$$\mathrm{\Delta }\varphi _{\mathrm{obs}}=\frac{\varphi _\mathrm{t}|\varphi _\mathrm{l}|}{2}$$
(1)
and the separation between the maxima of conal components is given by

$$W=\varphi _\mathrm{t}+|\varphi _\mathrm{l}|$$
(2)
The observed positions of the two conal components ($`\varphi _\mathrm{l}`$ and $`\varphi _\mathrm{t}`$) measured with respect to the core component are the *two* most desired observational parameters which replace the single $`W`$ parameter. Beyond question, recognizing this point is a great insight of Gangadhara & Gupta (2001).

The next step in the method of GG2001 was to express $`\mathrm{\Delta }\varphi _{\mathrm{obs}}`$ in terms of radial distance of conal radio emission $`r_{\mathrm{em}}`$. In the absence of the aberration and retardation effects (hereafter AR effects), the center (ie. the maximum) of the core component as well as the center of the conal pair would be observed at the same fiducial phase $`\varphi _\mathrm{f}`$ which formally corresponds to emission from $`r=0`$. In general, the radial distance $`r_{\mathrm{em}}`$ of the conal emission is different from the radial distance $`r_{\mathrm{cr}}`$ of the core emission. Therefore, the AR effects shift the center of the conal pair by $`\mathrm{\Delta }\varphi (r_{\mathrm{em}})`$ with respect to $`\varphi _\mathrm{f}`$, whereas the center of the core component is shifted with respect to $`\varphi _\mathrm{f}`$ by a different amount of $`\mathrm{\Delta }\varphi (r_{\mathrm{cr}})`$. The observed shift is the difference between these: $`\mathrm{\Delta }\varphi _{\mathrm{obs}}=\mathrm{\Delta }\varphi (r_{\mathrm{em}})\mathrm{\Delta }\varphi (r_{\mathrm{cr}})`$.

The value of $`\mathrm{\Delta }\varphi `$ is a sum of a shift due to the aberration and a shift due to the retardation: $`\mathrm{\Delta }\varphi =\mathrm{\Delta }\varphi _{\mathrm{ab}}(r)+\mathrm{\Delta }\varphi _{\mathrm{ret}}(r)`$, where the radial distance of radio emission $`r`$ may refer either to $`r_{\mathrm{em}}`$ or $`r_{\mathrm{cr}}`$. The retardation shift is equal to

$$\mathrm{\Delta }\varphi _{\mathrm{ret}}\frac{r}{R_{\mathrm{lc}}},$$
(3)
regardless of $`\alpha `$ and $`\beta `$, at least as long as the lowest order ($`r^{}`$) relativistic effects are considered.

To first order in $`r/R_{\mathrm{lc}}`$ (cf. Appendix A; GG2001), the aberration changes the direction of photon emission by an angle:

$$\eta _{\mathrm{ab}}\frac{v_{\mathrm{rot}}}{c}\frac{r}{R_{\mathrm{lc}}}\mathrm{sin}\zeta $$
(4)
where $`\stackrel{}{v}_{\mathrm{rot}}=\stackrel{}{\mathrm{\Omega }}\times \stackrel{}{r}`$ is the local corotation velocity at the emission point. Following the already published work (Cordes 1978; Phillips 1992), GG2001 assumed that the resulting phase shift $`\mathrm{\Delta }\varphi _{\mathrm{ab}}`$ is equal to $`\eta _{\mathrm{ab}}`$. This is, however, not true in general, because the aberrated emission direction, when rotated by $`360^{}`$ around $`\stackrel{}{\mathrm{\Omega }}`$, does not delineate a great circle on a sphere centered at the star. For this reason, the aberrational shift should be written in the form:

$$\mathrm{\Delta }\varphi _{\mathrm{ab}}\frac{\eta _{\mathrm{ab}}}{\mathrm{sin}\zeta }\frac{r}{R_{\mathrm{lc}}}.$$
(5)
Eq. (5) reflects the fact that any decrease in the aberration angle $`\eta _{\mathrm{ab}}`$ caused by the smaller corotation velocity near the rotation axis is cancelled out by the small circle effect (eq. A8). A rigorous derivation of this approximate result can be found in Appendix A. Thus, the total relativistic shift with respect to $`\varphi _\mathrm{f}`$ is equal to $`\mathrm{\Delta }\varphi 2r^{}`$, whereas the observed shift of the conal pair with respect to the core is equal to $`\mathrm{\Delta }\varphi _{\mathrm{obs}}=2(r_{\mathrm{em}}r_{\mathrm{cr}})R_{\mathrm{lc}}^1`$. The resulting formula for the altitude of the conal emission reads

$$h_{\mathrm{em}}=r_{\mathrm{em}}r_{\mathrm{cr}}\frac{\mathrm{\Delta }\varphi _{\mathrm{obs}}}{2}R_{\mathrm{lc}}.$$
(6)
In the case $`r_{\mathrm{cr}}r_{\mathrm{em}}`$, considered by GG2001:

$$r_{\mathrm{em}}\frac{\mathrm{\Delta }\varphi _{\mathrm{obs}}}{2}R_{\mathrm{lc}}.$$
(7)
This formula gives emission radii smaller by a factor $`a(1+\mathrm{sin}\zeta )/2`$ than those obtained by GG2001 (cf. their eq. 9). As shown in Appendix A, the approximation (7) holds with accuracy of $`10`$% for $`r^{}0.01`$. For objects with moderate or large dipole inclinations ($`\alpha 40^{}90^{}`$) this is a minor correction, however, for nearly aligned rotators the radii become almost two times smaller. The method of GG2001 applies to broad pulse profiles with well resolved conal components. This criterion favours objects with small $`\alpha `$, and, therefore, makes the modification important. GG2003 find particularly large emission radii for two objects, one of which (B2111$`+`$46) has very small inclination angle $`\alpha =9^{}`$. For this object, viewed at the angle $`\zeta =11.4^{}`$, our modification results in emission radii which are smaller by a factor of $`0.6`$ than those derived by GG2003. For the second problematic pulsar (B2045$``$16), viewed at a considerably larger angle $`\zeta =37.1^{}`$, our modification does not reduce the radii significantly. GG2003 report some difficulties with unambiguous identification of conal pairs for this pulsar.
Using the values of $`\varphi _\mathrm{l}`$ and $`\varphi _\mathrm{t}`$ published by GG2001 and GG2003 we recalculated the radio emission altitudes for objects studied therein. The resulting values are given in Table 1 (columns 5 and 6). The errors of $`h_{\mathrm{em}}`$ ($`22`$%) are much larger than the errors of $`\varphi _\mathrm{l}`$ and $`\varphi _\mathrm{t}`$ ($`3`$% as measured by Gangadhara & Gupta), because $`h_{\mathrm{em}}\mathrm{\Delta }\varphi _{\mathrm{obs}}`$ and $`\mathrm{\Delta }\varphi _{\mathrm{obs}}`$ is a very small difference between $`\varphi _\mathrm{t}`$ and $`|\varphi _\mathrm{l}|`$ (eq. 1).

We emphasize that our modification is not barely a rescaling of the method of GG2001, but it furnishes their method with completely new features. Unlike the original equation (9) of GG2001, our formula for the emission radius (eq. 7) *does not* depend on $`\alpha `$ and $`\beta `$, which is a valuable feature given the problems with determining $`\alpha `$. A method of BCW91 is also independent of $`\alpha `$ and $`\beta `$, however, it requires high-quality polarization measurements, and can only be applied to pulsars with well ordered position angle swings. Our modification of the method of GG2001 provides a method for determining $`r`$, which not only can do without the knowledge of $`\alpha `$ and $`\beta `$, but at the same time it does not require the polarization measurements. It is based solely on information contained in M, and T-type profiles. The applicability of different methods to different sets of available data is summarized in Table 2. Since the presence of a core in pulse profiles is often accompanied by a disordered shape of the position angle curve (eg. Rankin 1993), the two methods (ie. the delay-radius relation of BCW91 and our eq. 7) can complement each other. On the other hand, this is a disadvantage, because it makes it difficult to compare the methods.

To identify magnetic field lines from which the conal emission originates, it is necessary to determine the half opening angle $`\rho `$ of the emission cones. This is accomplished with the help of the cosine formula of spherical trigonometry, which connects the values of $`\rho `$, $`\alpha `$, and $`\zeta `$ with the observed separation of conal components $`W(f=1)`$:

$$\mathrm{cos}\rho =\mathrm{cos}\left(\frac{W}{2}\right)\mathrm{sin}\alpha \mathrm{sin}\zeta +\mathrm{cos}\alpha \mathrm{cos}\zeta .$$
(8)
This time the knowledge of $`\alpha `$ and $`\zeta `$ is necessary. Since the formula for $`\rho `$ does not depend on $`r`$ the values of $`\rho `$ calculated by GG2001 and GG2003 are correct.

With $`\rho `$ and $`r`$ determined from eqs. (8) and (7), one can calculate the colatitude of footprints of active magnetic field lines at the star surface with the dipolar formula:

$$\theta _{\mathrm{fp}}=\mathrm{arcsin}\left[\left(\frac{R_{\mathrm{ns}}}{r}\left(\frac{2}{3}\frac{1}{6}\left(x+\left(x^2+8x\right)^{1/2}\right)\right)\right)^{1/2}\right],$$
(9)
where $`x=\mathrm{cos}^2\rho `$ (see Appendix A). In the small angle approximation ($`\rho 1`$) the formula reduces to the well known approximate form:

$$\theta _{\mathrm{fp}}\frac{2}{3}\rho \left(\frac{R_{\mathrm{ns}}}{r}\right)^{1/2}.$$
(10)
Obviously, since the revised method predicts lower emission radii than its original version but the same opening angles $`\rho `$, it must yield larger footprint colatitudes $`\theta _{\mathrm{fp}}`$. Since, $`\theta _{\mathrm{fp}}r^{1/2}`$ the values of $`\theta _{\mathrm{fp}}`$ increase only by a factor of $`b=(2/(1+\mathrm{sin}\zeta ))^{1/2}`$ between $`1`$ and $`1.4`$. Consequently, the modified values of $`\theta _{\mathrm{fp}}^{}=\theta _{\mathrm{fp}}/\theta _{\mathrm{pc}}`$ for objects studied by GG2001 and GG2003 are given in the last column of Table 1. They range from 0.28 to 0.88. The values of $`\theta _{\mathrm{fp}}^{}`$ for the outermost cones cover the range between 0.35 and 0.88. The errors of $`\theta _{\mathrm{fp}}^{}`$ are based only on the errors of $`\varphi _\mathrm{l}`$ and $`\varphi _\mathrm{t}`$ determined by GG2001 and do not include the large uncertainties of $`\alpha `$, and $`\beta `$.

Since our revision rescales the values of $`r`$ and of $`\theta _{\mathrm{fp}}`$ by factors which are the same for a given pulsar, the main findings of GG2001 and GG2003 (ie. higher altitudes for broader cones and lower altitudes for higher radio frequencies) remain valid.

## 3 Comparison of the RPS methods with the geometric method

In addition to the conal components, there is one more feature in radio pulse profiles which for a long time has been believed to be symmetric with respect to the $`(\stackrel{}{\mathrm{\Omega }},\stackrel{}{\mu })`$ plane – the outer edge of the radio beam. The assumption about the symmetry has always been present in the traditional, geometric method of determining radio emission altitudes (eg. Kijak & Gil 1997). Therefore, the method of GG2001 may also be applied to the edges of radio beam, not only to maxima of the conal components. Hereafter, the RPS method of GG2001 with this assumption will be referred to as the “outer edge relativistic phase shift method” (OERPS method). The method requires a measurement of $`\varphi _\mathrm{l}`$ and $`\varphi _\mathrm{t}`$ for the outer edges of an *averaged* pulse profile. Then calculations of the original method of GG2001 can be performed (equations 19) to determine both the altitude and the locations $`\theta _{\mathrm{fp}}^{}`$ of magnetic field lines for the radio emission at the pulse edge. Thus, the OERPS method may be used to test the geometric method which a priori assumes $`\theta _{\mathrm{fp}}^{}(\mathrm{edge})=1`$.

The OERPS method does not require an identification of conal components within the pulse profile. This method will work without the analysis of single pulses (with the WT technique) and it can be applied directly to averaged pulse profiles with large signal to noise ratio. A significant limitation of the OERPS method is that it requires very accurate measurements of the observed positons $`\varphi _\mathrm{l}`$ and $`\varphi _\mathrm{t}`$ of the leading and the trailing edge of the radio pulse profile. These are much less defined than the positions of conal maxima. Therefore, the method may give less accurate results than the original RPS method of GG2001, unless very precise determination of $`\varphi _\mathrm{l}`$ and $`\varphi _\mathrm{t}`$ is achieved. GG2001 and GG2003 found that shifts of outermost cones’ maxima with respect to the core may in some cases be as large as $`3^{}`$, and they were able to measure some shifts with accuracy of $`0.05^{}`$ ie. a few percent. A much larger error (typically $`30`$ %) is inherent in the OERPS method and in the method of BCW91 which must deal with locating the “lowest intensity” points.

Because we have no access to high-quality radio data to apply the OERPS method, instead, below we try to use the BCW method to identify the active magnetic field lines. Like the OERPS method, the method of BCW91 provides the value of $`r`$ in a way which is independent of $`\alpha `$, $`\beta `$, and, most importantly, of $`\theta _{\mathrm{fp}}`$. The lowest intensity width $`W_0`$ (usually measured by BCW91 at the level of $`f=0.02`$ and used to determine the center of a profile), along with $`\alpha `$, $`\beta `$, and $`r`$, provides all necessary information to calculate $`\theta _{\mathrm{fp}}`$ with the help of equations (8) – (9), ie. in the same way as in GG2001.

However, since different components of a pulse profile may originate from different altitudes, the values of $`r`$ derived with some methods may represent an average over several components. To calculate $`\theta _{\mathrm{fp}}`$ for the pulse edge, it is necessary to know the radial distance $`r_{\mathrm{edge}}`$, which refers to the outer edges of the pulse profile. Establishing the relation between the pulse-averaged emission radius $`r`$, and the “edge radius” $`r_{\mathrm{edge}}`$ is an important step in determining $`\theta _{\mathrm{fp}}`$ for the outer edge of the pulse profile. Normally, calculating $`\theta _{\mathrm{fp}}`$ for the pulse edge would require using the equations (8) – (9) with $`r=r_{\mathrm{edge}}`$. An easier way is to use the equation (10) which implies that the ratio of $`\theta _{\mathrm{fp}}`$ corresponding to $`r_{\mathrm{edge}}`$ determined with some method, to $`\theta _{\mathrm{fp},\mathrm{geo}}`$ of the geometric method is:

$$\frac{\theta _{\mathrm{fp}}}{\theta _{\mathrm{fp},\mathrm{geo}}}\theta _{\mathrm{fp}}^{}\left(\frac{r_{\mathrm{geo}}}{r_{\mathrm{edge}}}\right)^{1/2},$$
(11)
because $`\theta _{\mathrm{fp},\mathrm{geo}}=\theta _{\mathrm{pc}}`$ by assumption, which we want to verify. In the above formula $`\theta _{\mathrm{fp},\mathrm{geo}}`$ and $`r_{\mathrm{geo}}`$ refer to the geometric method, whereas $`\theta _{\mathrm{fp}}`$ and $`r_{\mathrm{edge}}`$ refer to another method.<sup>1</sup><sup>1</sup>1For brevity, in eq. (11) we skip the index ‘edge’ at $`\theta _{\mathrm{fp}}`$ and $`\theta _{\mathrm{fp}}^{}`$, which hereafter should be understood as referring to the outer edge of the pulse profile. The above equation holds only when $`r_{\mathrm{edge}}`$ (and $`\theta _{\mathrm{fp}}`$) refer to the outer edge of a pulse profile, just as the values of $`r_{\mathrm{geo}}`$ do. Emission radii provided by the BCW method (hereafter denoted by $`r_{\mathrm{delay}}`$) are calculated under the assumption that all components in a pulse profile originate from the same altitude, ie. $`r_{\mathrm{edge}}=r_{\mathrm{delay}}`$ by assumption. Thus, inserting $`r_{\mathrm{edge}}=r_{\mathrm{delay}}`$ into eq. (11) one can calculate $`\theta _{\mathrm{fp}}^{}`$ predicted by the BCW method for the pulse edge.

We use the values of $`r_{\mathrm{delay}}`$ and $`r_{\mathrm{geo}}`$ calculated for the same data set by BCW91 and presented in their table 3. The values of $`\theta _{\mathrm{fp}}^{}`$ derived with the BCW method for seven clean-cut cases with reasonably small statistical errors, and most probably free from systematic errors, are presented in column 3 of Table 2. Column 4 presents a range of $`\theta _{\mathrm{fp}}^{}`$ allowed by 1 $`\sigma `$ errors derived for $`r_{\mathrm{delay}}`$ and $`r_{\mathrm{geo}}`$ by BCW91 (table 3 therein). The last column shows the level of consistency of the two methods (in $`\sigma `$), ie. the level at which the derived values are consistent with $`\theta _{\mathrm{fp}}^{}=1`$.

According to Table 2, for four cases (B0301$`+`$19 and B0525$`+`$21, both at 0.43 and 1.4 GHz) the BCW method yields $`\theta _{\mathrm{fp}}^{}>1`$, ie. radio emission from the closed magnetic field line region. The inconsistency with $`\theta _{\mathrm{fp}}^{}=1`$ is at $`2\sigma `$ level. The result $`\theta _{\mathrm{fp}}^{}>1`$ is equivalent to $`r_{\mathrm{edge}}<r_{\mathrm{geo}}`$ (cf. eq. 11), and, therefore may result from underestimating $`r_{\mathrm{edge}}`$ and/or from overestimating $`r_{\mathrm{geo}}`$. As noted by BCW91, contrary to the assumption $`r_{\mathrm{edge}}=r_{\mathrm{delay}}`$, the values of $`r_{\mathrm{edge}}`$ may actually be larger than $`r_{\mathrm{delay}}`$, because $`r_{\mathrm{delay}}`$ represents an average over the pulse profile, and inner parts of the profile may originate from lower altitudes than the outer edges of the pulse (eg. Rankin 1993; GG2001). Below we show, however, that even if the central parts of the pulse profiles originated from the star surface, the values of $`r_{\mathrm{edge}}`$ cannot exceed $`2r_{\mathrm{delay}}`$. In the BCW method, the values of $`r_{\mathrm{delay}}`$ are derived from the relative position of the outer edges of the profile and the center of the position angle (PA) curve. Let us assume that the center of the PA curve is determined mostly by the central parts of profiles which may originate from lower altitudes than the beam edge. Let us define a fiducial phase zero as a moment at which an observer detects a light signal emitted from the neutron star center when the dipole axis was located in the plane containing the rotation axis and the observer’s position. The star’s rotation shifts the profile center (determined by the midpoint between the *outer edges* of the pulse) toward earlier phases by $`2r_{\mathrm{edge}}/R_{\mathrm{lc}}`$ with respect to the fiducial phase, where $`r_{\mathrm{edge}}`$ is the radial distance for emission at the edge of the pulse profile. At the same time, the rotation delays the center of the position angle swing by $`2r_{\mathrm{in}}/R_{\mathrm{lc}}`$ with respect to the fiducial phase, where $`r_{\mathrm{in}}`$ refers to the radiation from inner parts of the profile which determine the position of the center of the PA swing, according to our working hypothesis. For the case considered by BCW91 (ie. $`r_{\mathrm{edge}}=r_{\mathrm{in}}`$), the total shift is $`4r_{\mathrm{edge}}/R_{\mathrm{lc}}`$. For the other limiting case $`r_{\mathrm{edge}}r_{\mathrm{in}}`$ the total shift is two times smaller: $`2r_{\mathrm{edge}}/R_{\mathrm{lc}}`$, since the delay of the center of the PA curve with respect to the fiducial phase is negligible. In this extreme case the value of $`r_{\mathrm{edge}}`$ derived from the shift $`\mathrm{\Delta }\varphi _{\mathrm{PA}}`$ between the center of the PA curve and the profile center is given by:

$$r_{\mathrm{edge}}=\frac{\mathrm{\Delta }\varphi _{\mathrm{PA}}}{2}R_{\mathrm{lc}}=2r_{\mathrm{delay}}(\mathrm{for}r_{\mathrm{edge}}r_{\mathrm{in}}).$$
(12)
Therefore, in general the emission radii for the beam edge may be at most two times larger than $`r_{\mathrm{delay}}`$ derived by BCW91 for the case of $`r_{\mathrm{edge}}=r_{\mathrm{in}}`$:

$$r_{\mathrm{edge}}2r_{\mathrm{delay}},$$
(13)
Using the upper limit of $`r_{\mathrm{edge}}=2r_{\mathrm{delay}}`$ in eq. (11), would decrease the values of $`\theta _{\mathrm{fp}}^{}`$ in Table 2 by a factor of $`2^{1/2}`$. In spite of the large error range for $`\theta _{\mathrm{fp}}^{}`$ (column 4 in Table 2), this would make only *two* of the four problematic values consistent with 1 within the level of $`1\sigma `$. The other two values of $`\theta _{\mathrm{fp}}^{}`$ (for B0301$`+`$19 at 1.42 GHz and B0525$`+`$21 at 0.43 GHz) would still remain inconsistent with 1 at $`1\sigma `$ level. Thus, lower emission altitudes for the central parts of profiles are not able to remove the disagreement completely. Indeed, different locations of emission regions for the conal peaks and for the bridge between them have been proposed for the two pulsars based on different spectra and fluctuation properties (Backer 1973; Rankin 1983; BCW91). For B0525$`+`$21, BCW91 suggest the contribution of a quadrupole component of stellar magnetic field as a source of the error in $`r_{\mathrm{delay}}`$. The underestimate of $`r_{\mathrm{delay}}`$ may also be caused by the rotationally-induced magnetic field line sweep-back, since the effect produces a shift of opposite sign to the shift $`\mathrm{\Delta }\varphi _{\mathrm{PA}}`$ found by BCW91. The significance of this effect will be investigated in a forthcoming paper.

The result $`\theta _{\mathrm{fp}}>1`$ (or $`r_{\mathrm{edge}}<r_{\mathrm{geo}}`$) could also arise due to an overestimate of $`r_{\mathrm{geo}}`$. This may be caused by refraction effects, which would broaden the radio beam with altitude (Lyubarskii & Petrova 1998). However, this effect could be compensated (or even dominated) by the unavoidable underestimate of $`r_{\mathrm{geo}}`$ which is inherent to the method of its derivation: the geometric method associates the maximum value of $`\theta _{\mathrm{fp}}^{}=1`$ with the observed pulse width, whereas the true value of $`\theta _{\mathrm{fp}}^{}`$ may well be smaller due to an inactivity of the outer parts of the polar cap and/or due to a limited sensitivity threshold. Therefore, in the case of the rectilinear propagation of radio waves, the values of $`r_{\mathrm{geo}}`$ should be considered as lower limits for $`r_{\mathrm{edge}}`$. In consequence, it is worth noting that BCW91 underestimated their geometric radii $`r_{\mathrm{geo}}`$, because these were calculated for the 10% intensity level ($`f=0.1`$) instead of the lowest intensity level, which BCW91 usually identified with 2% of maximum intensity. Therefore, the values of $`\theta _{\mathrm{fp}}^{}`$ in Table 2 are underestimated, and the actual disagreement of the BCW method with the condition $`\theta _{\mathrm{fp}}^{}1`$ is even larger than given in the last column of Table 2.

For B1914$`+`$13 an upper limit for $`\theta _{\mathrm{fp}}^{}`$ equal to $`0.8`$ is given in Table 2 (based on $`r_{\mathrm{geo}}`$ and $`r_{\mathrm{delay}}`$ from table 3 of BCW91). This is the only case (of single core type – S<sub>t</sub>), for which $`\theta _{\mathrm{fp}}^{}<1`$. Obviously, this could equally well imply the inactivity of the outer part of the polar cap as our inability to detect the outer wings of radio profiles. In either case, the assumption that $`\theta _{\mathrm{fp}}^{}=1`$ at the *observed* profile edge would be invalidated. We emphasize, that like any value of $`\theta _{\mathrm{fp}}^{}`$ in Table 2, the upper limit of $`\theta _{\mathrm{fp}}^{}=0.8`$ may be underestimated because of the relatively high level of intensity (10%) for which the values of $`r_{\mathrm{geo}}`$ were calculated in BCW91.

Thus, by selecting a sample of pulsars studied by BCW91 with the smallest statistical errors, and least probable to be affected by systematic errors, we find that the BCW method implies the radio emission from the region of closed magnetic field lines, (ie. with $`\theta _{\mathrm{fp}}^{}>1`$), at the level $`>2\sigma `$ in some cases. The way out would be to postulate either significant magnetic field deviations from the static dipole shape, or refraction effects (Lyubarskii & Petrova 1998). The distortions of the magnetic field may be caused by the presence of the quadrupolar stellar magnetic field (BCW91), or due to the rotationally induced sweep-back (Shitov 1983).

## 4 Discussion of assumptions of the relativistic phase shift methods

### 4.1 Asymmetry of the radio beam

To provide reliable estimates of $`r`$ the RPS method of GG2001 requires the symmetry of cones of enhanced radio emission with respect to the plane containing the core component and the rotation axis. To calculate $`\theta _{\mathrm{fp}}`$, the core/cones system must be centered at the dipole axis. The OERPS method requires the same symmetry for the outer edge of the radio beam, and the BCW method requires the symmetry of the beam edge with respect to the $`(\stackrel{}{\mathrm{\Omega }},\stackrel{}{\mu })`$ plane, regardless of the position of the core component. All these assumptions are related to the old problem: is the beam shape conal or patchy? The patchy beam would invalidate all these methods. Fortunately, however, statistical analyses of distribution of components within the radio pulse provide arguments for the conal shape (eg. Mitra & Deshpande 1999; Kijak & Gil 2002; Gil et al. 2002). It seems that arguments for patchy structure of the beam (Lyne & Manchester 1988, hereafter LM88) are not equally strong. In addition to the effects described by Gil et al. (1993), the aberration and the retardation may enhance the *apparent* patchiness of the beam. This is because central components of pulse profiles seem to originate from lower altitudes than the outer components (eg. Rankin 1993; GG2001). The aberration and retardation effects shift the relative positions of different components in phase by $`\mathrm{\Delta }\varphi 2\mathrm{\Delta }r^{}`$, where $`\mathrm{\Delta }r^{}`$ is the difference of emission radii between different components in units of $`R_{\mathrm{lc}}`$. Estimating $`\mathrm{\Delta }r^{}0.01`$ one obtains $`\mathrm{\Delta }\varphi 1^{}`$ which amounts to $`13`$% of pulse width observed for pulsars with large dipole inclinations ($`\alpha 90^{}`$). Given that observed separation between components is much smaller than the pulse width (say $`2040`$%), these shifts, different for different objects because of diversity of $`r`$, along with other factors (eg. inexactly determined impact angles $`\beta `$) could easily produce the apparently random distribution of components within the pulse window, ie. in figures like the fig. 12 in LM88.

Other arguments for the conal beam shape are provided directly by the WT technique of GG2001. For 6 (out of 7) objects studied by GG2001 and GG2003, the same number of conal components was found on both sides of the core component, and only B2045$``$16 is a *possible* exception. In all cases, strong (ie. easily identifiable) components on the leading side corresponded to strong components on the trailing side, whereas new weak components, detected with the WT technique on one side of the core were always associated with weak components on the other side. It is natural to interpret this coincidence in terms of the conal, not patchy beam shape. Another interesting implication of this relation is the possibility of applying the method of GG2001 to a larger number of T and M pulsars, without the neccessity of identification of all cones with the WT technique.

Thus, the conal beam shape seems to be consistent with observations. On theoretical grounds, the conal shape is supported by a model of sparks rotating around the magnetic pole due to the $`\stackrel{}{E}\times \stackrel{}{B}`$ drift (Ruderman & Sutherland 1975; Gil & Sendyk 2000). The model, as well as the conal beam hypothesis, is supported by observations of subpulse drifts (Gil & Krawczyk 1997; Deshpande & Rankin 1999; Vivekanand & Joshi 1999).

### 4.2 The reference altitude

All three RPS methods rely on a measurement of a shift between three identifiable points within a pulse profile and make assumptions about relative altitudes of these points. The method of GG2001 provides the *altitude* of the conal radio emission measured with respect to the altitude of the core emission (eq. 6). To identify the active magnetic field lines (ie. to calculate $`\theta _{\mathrm{fp}}^{}`$) it is necessary either to know $`r_{\mathrm{cr}}`$ or to rely on the assumption $`r_{\mathrm{cr}}r_{\mathrm{em}}`$ (as we did in Table 1) which effectively means that the derived values of $`\theta _{\mathrm{fp}}^{}`$ are upper limits.

There is a very simple geometric argument against the hypothesis that the core components originate from the very vicinity of the star surface: With the only exception of B1237$`+`$25, for all other objects studied by GG2001 and GG2003, and listed in our Table 1, the impact angles $`\beta `$ are *larger* than the half opening angle of the last open field lines at the surface of the neutron star with radius $`R_{\mathrm{ns}}=10`$ km. Had the central components (identified as cores) originated close to the surface, our line of sight should miss them. Therefore, the impact angles determined for these objects from fitting the position angle swing are in clear disagreement with the assumptions of the method of determining $`\alpha `$ proposed by Rankin (1990). Assuming rectilinear propagation of radio waves, this discrepancy can be interpreted only in two ways: 1) either the radial distance for the core emission is much larger than $`R_{\mathrm{ns}}`$, or 2) the impact angles $`\beta `$ are systematically overestimated. In the first case, the dipole inclination angles $`\alpha `$ derived with Rankin’s method would be underestimated, at least as long as the width of the core beam is assumed to reflect the width of the open field line region at a given altitude. For the time being, it seems to be a matter of personal preference which of the possibilities (either 1. or 2.) is the case – Rankin (1990) assumed that the values of $`\beta `$ are in error. Regardless of the choice, however, either $`\alpha `$ (based on the assumption $`r_{\mathrm{cr}}R_{\mathrm{ns}}`$) or $`\beta `$ *must* be in error, since they contradict each other. Therefore, the values of $`\theta _{\mathrm{fp}}^{}`$ given in the last column of our Table 1 should be treated with caution.

PSR B0450$``$18, studied by GG2003 has the largest impact angle ($`\beta 4^{}`$) among all the objects listed in Table 1. The half opening angle of the open field line region at the surface of the pulsar is equal to $`1.5\theta _{\mathrm{pc}}1.68^{}`$. To exhibit the estimated width of $`8^{}`$ (Rankin 1993), the central component in the pulse profile of this object would have to originate from radial distance $`r_{\mathrm{cr}}67`$ km, assuming the uncertain values of $`\alpha =24^{}`$, $`\beta =4^{}`$ and $`\theta _{\mathrm{fp}}^{}=1`$ from Rankin (1990). This value is larger than the error of $`h_{\mathrm{em}}`$ estimated for this object in Table 1 and should be added to $`h_{\mathrm{em}}`$ to obtain $`r_{\mathrm{em}}`$. Had the core component not filled in the open field line region, the value of $`r_{\mathrm{cr}}`$ would be even larger.

The RPS method of BCW91 assumes that radio emission at the edge of a profile originates from *the same* altitudes as the emission which determines the center of the position angle curve, in which case the predicted shift is equal to $`4r_{\mathrm{delay}}/R_{\mathrm{lc}}`$. As noted above this assumption, if not satisfied, would lead to an underestimate of emission altitudes for the outer edge of a profile. As shown in Section 3, however, the values of $`r_{\mathrm{edge}}`$ cannot exceed $`2r_{\mathrm{delay}}`$.

## 5 Conclusions

We have modified the relativistic method of GG2001 of determining the radio emission geometry within pulsar magnetospheres. Our modification results in a method of determining the radio emission altitudes which does not depend on viewing geometry nor does it require polarization measurements. According to this method, the altitude of the radio emission region in units of $`R_{\mathrm{lc}}`$ is equal to a half of the relativistic phase shift of a pair of conal components with respect to the core component. We propose to extend application of the revised method to the outer edge of radio pulse profiles to identify magnetic field lines associated with the edge. This may provide a test of the geometric method of determining $`r`$, based on a measurement of a pulse width at the lowest intensity level.

As noted by GG2003, the radio emission radii derived with their original method for the outermost cones were notoriously larger than the radii derived for the pulse edge with the geometric method by Kijak & Gil (1997; 1998). Our revision has removed part of this discrepancy. However, since the geometric method assumes that the observed radio emission fills in the open field line region entirely ($`\theta _{\mathrm{fp}}^{}=1`$), the emission radii derived with the geometric method should be considered as lower limits for $`r_{\mathrm{edge}}`$ (provided the refraction effects can be neglected). Therefore, other methods of determining $`r_{\mathrm{edge}}`$ for the pulse edge should always yield $`r_{\mathrm{edge}}r_{\mathrm{geo}}`$.

We determined the values of $`\theta _{\mathrm{fp}}^{}`$ for the pulse edge at the 10% level using the relativistic method of Blaskiewicz et al. (1991). This has revealed that in individual cases the method often implies emission from the region of closed magnetic field lines. BCW91 associated this with the fact that the emission radii $`r_{\mathrm{delay}}`$ determined with their method may be smaller than $`r_{\mathrm{edge}}`$ because inner parts of pulse profiles can originate from lower altitudes. We have shown that even accounting for this, the value of $`r_{\mathrm{edge}}`$ cannot exceed $`2r_{\mathrm{delay}}`$, which is still not sufficient to solve the problem. Apparently, other effects, like the rotation-induced sweep-back of the magnetic field lines, or refraction effects may be responsible for the discrepancy.

JD thanks V. S. Beskin for comments on pulsar magnetosphere. We are grateful to J. Gil and R. Gangadhara for their valuable comments on the manuscript. We appreciate the comments of the anonymous referee who helped us to improve the paper significantly. This work was performed while JD held a National Research Council Research Associateship Award at NASA/GSFC. The work was also supported by the grant 2P03D.004.24 (JD and BR).

## Appendix A Derivation of equations (4), (5), and (9).

Let the radio waves in the reference frame corotating with the star be emitted in the direction $`\stackrel{}{k}^{}`$ at angle $`\zeta ^{}`$ with respect to the rotation axis. Because of the aberration, in the inertial observer frame the emission direction will be $`\stackrel{}{k}\stackrel{}{k}^{}`$, and the observer will have to be located at an angle $`\zeta `$ with respect to $`\stackrel{}{\mathrm{\Omega }}`$ to detect the radiation (see Fig. 1). In general, $`\zeta `$ differs from $`\zeta ^{}`$ (eg. Kapoor & Shukre 1998). The unit vectors of the emission directions $`\stackrel{}{k}^{}`$ and $`\stackrel{}{k}`$ are related by the aberration formula:

$$\stackrel{}{k}=\frac{\stackrel{}{k}^{}+[\gamma +(\gamma 1)(\stackrel{}{\beta }_{\mathrm{rot}}\stackrel{}{k}^{})/\beta _{\mathrm{rot}}^2]\stackrel{}{\beta }_{\mathrm{rot}}}{\gamma (1+\stackrel{}{\beta }_{\mathrm{rot}}\stackrel{}{k}^{})}$$
(A1)
where $`\stackrel{}{\beta }_{\mathrm{rot}}`$ is the local corotation velocity in units of the speed of light $`c`$, and $`\gamma =(1\beta _{\mathrm{rot}}^2)^{1/2}`$ (eg. Dyks & Rudak 2003). Neglecting the terms of the order of $`\beta _{\mathrm{rot}}^2`$ and higher, the above formula reduces to:

$$\stackrel{}{k}(1+\stackrel{}{\beta }_{\mathrm{rot}}\stackrel{}{k}^{})\stackrel{}{k}^{}+\stackrel{}{\beta }_{\mathrm{rot}}.$$
(A2)
Eq. (10) implies that in the small angle approximation (ie. for emission from within the open field line region at altitudes $`rR_{\mathrm{lc}}`$) the value of the term $`\stackrel{}{\beta }_{\mathrm{rot}}\stackrel{}{k}^{}`$ is limited as follows

$$\stackrel{}{\beta }_{\mathrm{rot}}\stackrel{}{k}^{}\beta _{\mathrm{rot}}\mathrm{cos}\left[\frac{\pi }{2}\frac{3}{2}(r^{})^{1/2}\right]=\beta _{\mathrm{rot}}\mathrm{sin}[1.5(r^{})^{1/2}](r^{})^{3/2},$$
(A3)
where $`r^{}=r/R_{\mathrm{lc}}1`$ is the radial distance of emission. Neglecting this term, one obtains

$$\stackrel{}{k}\stackrel{}{k}^{}+\stackrel{}{\beta }_{\mathrm{rot}}$$
(A4)
with accuracy of the order of $`(r^{})^{1/2}`$ (ie. with 10% error for radio emission at 1% of $`R_{\mathrm{lc}}`$; hereafter we will use $`r^{}=0.01`$ in all estimations). Therefore, to the first order in $`r^{}`$, the aberration angle $`\eta _{\mathrm{ab}}`$ between $`\stackrel{}{k}`$ and $`\stackrel{}{k}^{}`$ is approximately equal to $`\beta _{\mathrm{rot}}`$. By applying the cosine theorem of spherical trigonometry to the spherical triangle ABC in Fig. 1:

$$\mathrm{cos}\zeta =\mathrm{cos}\zeta ^{}\mathrm{cos}\beta _{\mathrm{rot}}+\mathrm{sin}\zeta ^{}\mathrm{sin}\beta _{\mathrm{rot}}\mathrm{cos}\frac{\pi }{2}$$
(A5)
and neglecting the $`\beta _{\mathrm{rot}}^2`$ terms and higher, one finds that

$$\zeta \zeta ^{},$$
(A6)
ie. the aberration does not change the colatitude of photon emission direction measured with respect to $`\stackrel{}{\mathrm{\Omega }}`$ (to the order of $`\beta _{\mathrm{rot}}`$).

Application of the sine theorem to the ABC triangle:

$$\frac{\mathrm{sin}\mathrm{\Delta }\varphi _{\mathrm{ab}}}{\mathrm{sin}\eta _{\mathrm{ab}}}=\frac{\mathrm{sin}(\pi /2)}{\mathrm{sin}\zeta }$$
(A7)
gives, for $`\mathrm{\Delta }\varphi _{\mathrm{ab}}1`$, and $`\eta _{\mathrm{ab}}1`$, the relation:

$$\mathrm{\Delta }\varphi _{\mathrm{ab}}\frac{\eta _{\mathrm{ab}}}{\mathrm{sin}\zeta },$$
(A8)
used in eq. (5) (the minus sign in eq. (5) has been inserted to account for the fact that larger azimuths correspond to earlier detection phases). Eq. (A4) implies that $`\eta _{\mathrm{ab}}\beta _{\mathrm{rot}}`$, ie.:

$$\eta _{\mathrm{ab}}r^{}\mathrm{sin}\theta =r^{}\mathrm{sin}(\zeta +\delta ),$$
(A9)
where $`\theta `$ is the colatitude of the emission point (measured from the rotation axis) and $`\delta =\theta \zeta `$ so that $`|\delta |3^1|\beta |2^1(r^{})^{1/2}1`$, with $`\beta `$ being the impact angle. Thus, for $`|\delta |1`$ and $`|\delta |\zeta `$, one obtains

$$\eta _{\mathrm{ab}}r^{}\mathrm{sin}\zeta \left(\mathrm{cos}\delta +\frac{\mathrm{sin}\delta }{\mathrm{tan}\zeta }\right)r^{}\mathrm{sin}\zeta ,$$
(A10)
which proves eq. (4). The neglected term $`\mathrm{sin}\delta /\mathrm{tan}\zeta `$ becomes large for the nearly aligned geometry. In general, ie. allowing for the maximum possible range of impact angles $`\beta `$ between $`\pm 1.5(r^{})^{1/2}`$, the resulting error in $`\eta _{\mathrm{ab}}`$ could be equal to 10% for $`\zeta `$ as large as $`30^{}`$. In the context of the RPS models discussed in this paper, however, the angle $`\zeta `$ corresponding to the 10% error will be much smaller. This is because the line of sight must traverse nearly through the center of the open field line region for the core/cone structure to be discernible, ie. $`\beta \epsilon 1.5(r^{})^{1/2}`$ with the parameter $`\epsilon `$ considerably smaller than 1. Adopting a reasonable value for $`\epsilon `$, eg. $`\epsilon 0.2`$, one finds that the error of $`\eta _{\mathrm{ab}}`$ exceeds 10% when $`\zeta 6^{}`$.

Inserting eq. (A10) into eq. (A8) gives $`\mathrm{\Delta }\varphi _{\mathrm{ab}}r^{}`$ and proves eq. (5) The independence of $`\mathrm{\Delta }\varphi _{\mathrm{ab}}`$ on $`\zeta `$ results directly from the fact that the decrease in the aberrational angle $`\eta _{\mathrm{ab}}`$ for the line of sight approaching the rotation axis (eq. A10) is cancelled out by the “small circle” effect (eq. A8). Thus, we find that *as long as the emission direction in the CF does not deviate much from the meridional plane, the aberrational phase shift does not depend on viewing geometry (ie. neither on $`\alpha `$ nor on $`\zeta `$)*. The simplicity of the delay-radius relation derived in BCW91 has the same origin.

We now turn to the derivation of eq. (9). Let $`\stackrel{}{\mu }`$ be the magnetic moment of the dipole, and let $`\stackrel{}{k}^{}`$ be the emission direction tangent to the local magnetic field in the CF. The equation of the dipolar magnetic field lines is: $`r^1\mathrm{sin}^2\theta _\mathrm{m}=\mathrm{const}`$, where $`\theta _\mathrm{m}`$ is the colatitude of the emission point in the reference frame with $`\widehat{z}`$ axis along $`\stackrel{}{\mu }`$. This equation implies that

$$\theta _{\mathrm{fp}}=\mathrm{arcsin}\left[\left(\frac{R_{\mathrm{ns}}}{r}\mathrm{sin}^2\theta _\mathrm{m}\right)^{1/2}\right],$$
(A11)
where $`R_{\mathrm{ns}}`$ is the radius of the neutron star. The angle $`\rho =\mathrm{}(\stackrel{}{\mu },\stackrel{}{k}^{})`$ is given by

$$\mathrm{cos}\rho =\frac{\stackrel{}{\mu }\stackrel{}{k}^{}}{|\stackrel{}{\mu }|}=\frac{3\mathrm{cos}^2\theta _\mathrm{m}1}{\left(1+3\mathrm{cos}^2\theta _\mathrm{m}\right)^{1/2}}$$
(A12)
By calculating $`\theta _\mathrm{m}`$ from eq. (A12) and inserting into (A11) one obtains eq. (9). The uncertainty of $`R_{\mathrm{ns}}`$, present in eq. (A11), does not affect the ratio $`\theta _{\mathrm{fp}}/\theta _{\mathrm{pc}}`$, as long as $`\theta _{\mathrm{fp}}1`$. The measurement errors of $`\varphi _\mathrm{l}`$ and $`\varphi _\mathrm{t}`$ result in the error of $`\theta _{\mathrm{fp}}`$ mainly through the inverse square root dependence on $`r`$. The error of the cones’ separation $`W`$ (which enters the formula (9) through eq. 8) is negligible. This is because $`\mathrm{\Delta }\varphi _{\mathrm{obs}}`$ (and, thereby, $`r`$) is calculated as a *small difference* between $`\varphi _\mathrm{t}`$ and $`|\varphi _\mathrm{l}|`$ (with $`\varphi _\mathrm{t}|\varphi _\mathrm{l}|`$) whereas $`W`$ is a sum of these (cf. eqs. 1 and 2 and the comments in Sec. 2). The dependence $`\theta _{\mathrm{fp}}r^{1/2}`$ implies that the fractional error of $`\theta _{\mathrm{fp}}`$ ($`13`$%) is approximately two times smaller than the fractional error of $`r`$ (cf. columns 5 and 7 in Table 1).


format: text

text:
# The Intrinsic Time Scale of Transient Neuronal Responses

## Abstract

In a generic neuron model, we present the linear response theory for the firing rate in response to both time dependent input currents and noise amplitudes. In both cases the signal transmission is strongly attenuated for frequencies above the stationary firing rate. For high frequencies both the mean input and the noise transmission function decay as $`\omega ^2`$, independent of model details. Our results indicate that previously suggested mechanisms for near instantaneous transmission of information are not consistent with the spike generation mechanism of real neurons.

In the brain information is processed through a hierarchy of neural layers. Based on the observation that human subjects can successfully distinguish complex visual stimuli in only a few hundred milliseconds Thorpe , it has been argued that the processing time of individual layers must be very small. Recently two different mechanisms have been proposed for the practically instantaneous relaying of information at the single neuron level. They utilize the fact that in vivo cortical neurons exhibit substantial subthreshold fluctuations of their membrane potentials, due to the large number of presynaptic neurons DestexhePare . The first mechanism assumes that information is transferred by a time dependent noise amplitude SigmaModulation1 ; SigmaModulation2 , the second assumes information transfer via a modulation of the mean input in the presence of temporally correlated background noise Gerstner ; Fourcaud . These studies, which are based on a seminal work by Knight Knight , have been conducted on integrate-and-fire models. These models, however, are highly idealized models of cortical neurons. Their main ingredient is a fixed voltage threshold. Each time the membrane potential reaches this threshold, a spike is said to be emitted. Real neurons, however, as well as biophysically realistic conductance based neuron models do not have a fixed voltage threshold and are rather comparable to excitable systems. This naturally raises the question whether the response properties found in simple threshold models are preserved if one considers models with a more realistic spike generating mechanisms. In this letter we address this question by a combination of analytical and numerical techniques considering a generic Type-I neuron model, the $`\theta `$-Neuron GutkinErmentrout . We show that in a realistic regime the response to a transient stimulus, both in the mean input as well as in the noise amplitude is *not* instantaneous as predicted by the integrate-and-fire models. Moreover, we show that the maximum transmission frequency is approximately determined by the stationary firing rate. Thus our results suggest that the $`\theta `$-Neuron reproduces the dynamical behavior of cortical neurons faithfully, while still being tractable analytically.

Neurons have been divided into several classes based on the kind of excitability they exhibit. In our study we focus on Type-I excitability which can be encountered in many neurons withing the cortex HodgkinHuxley . In models of Type-I neurons, repetitive firing typically emerges via a saddle-node bifurcation. The normal form of this bifurcation is equivalent to a phase oscillator, the $`\theta `$-Neuron GutkinErmentrout . It has been shown that it has the same mechanism of excitability and super-threshold behavior as cortical neurons HanselMato . Its dynamics is:

$$\tau \dot{\theta }=(1\mathrm{cos}\theta )+I(t)(1+\mathrm{cos}\theta ),$$

where $`\tau `$ is the time constant and $`I(t)`$ represents the total synaptic input to the neuron. Each time the oscillator crosses the point $`\theta =\pi `$ a spike is said to be “fired”. For constant inputs this model exhibits a periodic firing regime for $`I>I_c=0`$ with firing rate $`\nu =\sqrt{I}/\left(\pi \tau \right)`$, and an excitable regime for $`I<0`$. We analyze the response of an ensemble of such neurons, which is given by the ensemble averaged firing rate $`\nu (t)`$. The input current $`I(t)`$ is decomposed,

$$I(t)=I_0+\sigma \sqrt{\tau }z(t),$$

into a mean current $`I_0`$ and a noise term which reflects the fluctuations induced by the synaptic inputs Tuckwell . The noise is modeled by an Ornstein-Uhlenbeck process with a correlation time $`\tau _c`$:

$$\tau _c\frac{dz}{dt}=z+\eta (t),$$

where $`\eta (t)`$ is Gaussian white noise. In the limit $`\tau _c0`$, $`z(t)`$ becomes white noise. In the following, we have chosen the time constant $`\tau =0.25\text{ms}`$ which results in a spike duration of about $`1\text{ms}`$ as found in real neurons. $`I_0`$ and $`\sigma `$ have been chosen to give realistic $`\theta `$-correlation times of approximately $`10\text{ms}`$ and firing rates in the range $`120\text{Hz}`$.

The state of an ensemble of such neurons is described by a probability density function $`P(\theta ,z,t)`$. Its dynamics is determined by the Fokker-Planck equation:

$$_tP=\widehat{L}P,$$
(1)
with the operator,

$`\widehat{L}`$ $`=`$ $`\tau ^1_\theta \left(\left(1\mathrm{cos}\theta \right)+\left(I_0+\sigma \tau ^{1/2}z\right)\left(1+\mathrm{cos}\theta \right)\right)`$
$`+\tau _c^1_zz+{\displaystyle \frac{1}{2}}\tau _c^2_z^2,`$
periodic boundary conditions in the $`\theta `$-direction and natural boundary conditions in the $`z`$-direction. The firing rate $`\nu (t)`$ is identical to the total probability current through the line $`\theta =\pi `$:

$$\nu (t)=2\tau ^1_{\mathrm{}}^{\mathrm{}}P(\pi ,z,t)𝑑z.$$

In the limit of temporally uncorrelated input, $`\tau _c0`$, the rate can be calculated analytically Risken ; LindnerBulsara  which gives:

$$\nu _{\text{WN}}^1=\frac{4\tau \sqrt{\pi }}{\sigma }_0^{\mathrm{}}𝑑y\mathrm{exp}\left\{\frac{4}{\sigma ^2}\left(\frac{y^6}{3}+I_0y^2\right)\right\}$$

To investigate how the neuron responds to time dependent synaptic inputs it is important to consider both the response to a time dependent mean input current $`I=I_0+ϵe^{i\omega t}+\sqrt{\tau }\sigma z`$ and a time dependent noise amplitude $`I=I_0+\sqrt{\tau }\left(\sigma +ϵe^{i\omega t}\right)z`$. Inserting the modulated input current into Eq. (1) and expanding $`P(\theta ,z,t)=P_0(\theta ,z)+ϵ\stackrel{~}{P}_\omega (\theta ,z)e^{i\omega t}+\mathrm{}`$ gives in linear order in $`ϵ`$ Risken :

$$i\omega \stackrel{~}{P}_\omega (\theta ,z)=\widehat{L}\stackrel{~}{P}_\omega f(\theta ,z),$$
(2)
with $`f(\theta ,z)=\tau ^{\text{-1}}_\theta (1+\mathrm{cos}\theta )P_0(\theta ,z)`$ for modulated input currents and $`f(\theta ,z)=\tau ^{1/2}_\theta \left(1+\mathrm{cos}\theta \right)zP_0(\theta ,z)`$ for modulated noise amplitudes. The formal solution of eq. (2) is:

$$\stackrel{~}{P}_\omega (\theta ,z)=e^{i\omega t}_{\mathrm{}}^te^{(tt^{})\widehat{L}}f(\theta ,z)e^{i\omega t^{}}𝑑t^{}.$$
(3)
This integral can be solved in terms of eigenfunctions $`P_k(\theta ,z)`$ of the operator $`\widehat{L}`$:

$$\lambda _kP_k(\theta ,z)=\widehat{L}P_k(\theta ,z),$$
(4)
with the associated eigenvalues $`\lambda _k`$. Because of the applied boundary conditions the spectrum $`\left\{\lambda _k\right\}`$ is discrete, and since detailed balance is not fulfilled, eigenvalues $`\lambda _k`$ and the corresponding eigenfunctions $`P_k(\theta ,z)`$ are either real or form complex conjugate pairs. Moreover, the eigenfunctions are typically not orthogonal. However, an orthonormal set of basis functions $`\left\{|\varphi _k\right\}`$ can be constructed from them by Gram-Schmidt orthogonalization (in Dirac notation) with expansion coefficients $`a_l^k`$:

$$|\varphi _k(\theta ,z)=\underset{l=1}{\overset{k}{}}a_l^kP_l(\theta ,z).$$

Inserting $`_k|\varphi _k\varphi _k|=\mathrm{𝟏}`$ into Eq. (3) yields,

$`\stackrel{~}{P}_\omega (\theta ,z)`$ $`=`$ $`e^{i\omega t}{\displaystyle \underset{k,lk}{}}b_ka_l^k{\displaystyle _{\mathrm{}}^t}e^{i\omega t^{}+(tt^{})\widehat{L}}P_l(\theta ,z)`$
$`=`$ $`{\displaystyle \underset{k,lk}{}}{\displaystyle \frac{a_l^kb_k}{i\omega \lambda _l}}P_l(\theta ,z),`$
with $`b_k=\varphi _k|f`$. The rate response is then given by:

$`\nu (t)`$ $`=`$ $`2\tau ^1{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left(P_{\text{st}}(\pi ,z)+\stackrel{~}{P}_\omega (\pi ,z)e^{i\omega t}\right)𝑑z`$
$`=`$ $`\nu _0+\nu _1(\omega )e^{i(\omega t+\varphi (\omega ))}.`$
The limit $`\omega \mathrm{}`$ can be treated analytically. Since $`f(\theta ,z)`$ vanishes at $`\theta =\pi `$, the modulus of $`\stackrel{~}{P}`$ at this point has to be proportional to $`\omega ^2`$:

$$\left(i\omega \widehat{L}\right)\stackrel{~}{P}_\omega (\pi ,z)e^{i\omega t}=f(\pi ,z)=0.$$

An expansion in $`\omega ^1`$ reveals that the relative phase of $`\nu _1(\omega )`$ is $`\pi `$ in this limit. We would like to stress that this decay is universal and does not depend on model details. It is only due to the insensitivity to external inputs at the point where a spike is emitted, i.e. $`f(\pi ,z)=0`$.

The eigenvalues and eigenfunctions of the time-independent operator $`\widehat{L}`$ determine $`\nu _1(\omega )`$ completely. We computed them using a matrix representation of Eq. (4) obtained by expanding $`P_k(\theta ,z)`$ into a complete set of orthonormal functions:

$`P_k(\theta ,z)`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}a_{n,m}\left(2^{n+1}\sqrt{\pi /2\tau _c}m!\right)^{1/2}`$
$`e^{in\theta }H_m(\sqrt{2\tau _c}z)e^{z^2\tau _c}.`$
Here $`H_m(z)`$ are the Hermite polynomials AbramowitzStegun . Inserting this into Eq. (4), multiplying from left with $`\left(2^{n^{}+1}\sqrt{\pi /2\tau _c}m^{}!\right)^{1/2}e^{in^{}\theta }H_m^{}(\sqrt{2\tau _c}z)e^{z^2\tau _c}`$ and integrating over the entire domain leads to the following eigenvalue problem:

$$\lambda a_{n,m}=\underset{n^{},m^{}}{}L_{n,m;n^{},m^{}}a_{n^{},m^{}}.$$
(5)
Because the Fokker-Planck operator has only two Fourier components in the $`\theta `$-direction and is polynomial in the $`z`$-direction, $`𝐋`$ is sparse with:

$`L_{n,m;n^{},m^{}}`$ $`=`$ $`i\tau ^1(1+I_0)n\tau _c^1m`$
$`L_{n,m;n^{}\pm 1,m^{}}`$ $`=`$ $`(2\tau )^1i(1I_0)n`$
$`L_{n,m;n^{},m^{}1}`$ $`=`$ $`in\sigma \left(4\tau \tau _c\right)^{1/2}\left(m+1\right)`$
$`L_{n,m;n^{},m^{}+1}`$ $`=`$ $`in\sigma \left(4\tau \tau _c\right)^{1/2}m`$
$`L_{n,m;n^{}\pm 1,m^{}\pm 1}`$ $`=`$ $`in\sigma \left(16\tau \tau _c\right)^{1/2}(m+1)`$
$`L_{n,m;n^{},m^{}2}`$ $`=`$ $`\left(\sqrt{2\tau \tau _c}\right)^1\sqrt{(m+1)(m+2)}`$
We solve (5) numerically using the Arnoldi-method Arnoldi , a high performance iterative algorithm. Figure 1 shows the spectrum and the stationary density for two different correlation times $`\tau _c`$ together with the spectrum in the white noise limit. The eigenvalues are arranged in a series of wedges. The tip of each wedge is located at integer multiples of $`\tau _c^1`$. The different wedges account for the decay of excitations in the $`z`$-direction, whereas the eigenvalues within a wedge account for the decay (real part) and oscillation (imaginary part) of excitations in the $`\theta `$-direction. In the limit $`\tau _c0`$ the spectrum exhibits only one wedge (open circles). For $`\tau _c=10\text{ms}`$ (upper plot), the interaction between eigenvalues leads to deviations from the white noise limit only for strongly damped modes. The deviations become more pronounced only for large values of $`\tau _c=50\text{ms}`$ (lower plot). In both regimes the stationary density strongly deviates from a separable density.

Examples of the response amplitude $`\nu _1(\omega )`$ and the phase $`\varphi (\omega )`$ are depicted for different values of $`\nu _0`$ in Fig. 2 for both mean input and noise stimulation. In the case of mean input modulation, the linear response amplitude exhibits a resonance maximum at approximately the frequency of its stationary firing rate and then decays rapidly to zero. The response phase starts at zero and then drops to $`\pi `$. For a modulation in the noise amplitude the behavior is similar, except for additional resonances at higher frequencies, which are, however, strongly damped.

Whereas in general all eigenfunctions and eigenvalues contribute to the rate response, the relatively simple Lorentz-like shape of the response functions suggests that one frequency effectively dominates the neurons’ behavior. We observed that this cut-off frequency is given by the imaginary part of the first excited eigenvalue of the first wedge for the current response, and by the corresponding eigenvalue of the second wedge for the noise response. This is demonstrated by the green dashed lines in Fig. 2. These are the sum of two Lorentzians with their maxima at the positive and negative frequency given by the imaginary part of the second eigenvalue in the corresponding wedge. The reason for this is apparent from for the decomposition of $`f(\theta ,z)`$ into eigenfunctions of the Fokker-Planck operator $`\widehat{L}`$:

$$f(\theta ,z)=\underset{l}{}\alpha _lP_l(\theta ,z).$$

It has no contribution from the stationary density, since $`f(\theta ,z)𝑑\theta =0`$, and for mean input modulations has a maximum contribution for the eigenfunctions corresponding to second eigenvalue in the first wedge. For noise modulations it has a maximum contribution for the eigenfunction corresponding to the second eigenvalue in the second wedge. The imaginary parts of the second eigenvalues in the first two wedges are, however, almost identical for the parameters used.

The cut-off frequency for subthreshold mean inputs $`I_0<0`$ is depicted in Fig. 3 together with the dependence on the noise correlation time $`\tau _c`$. The cut-off frequency increases for increasing values of $`I_0`$ but always stays below the stationary rate $`\nu _0`$. Increasing the noise correlation time at $`I_0=0`$ shifts the cut-off frequency to slightly larger values compared to $`\nu _0`$.

Thus our analysis demonstrates that in a realistic regime responses much faster than $`\nu _0^1`$ are strongly damped, since the transmission function decays as $`\omega ^2`$. We would like to conclude this letter with a comparison of the dynamic behavior of the $`\theta `$-neuron with a biophysical realistic conductance based neuron WangBuszaki  as well as with the classical leaky integrate-and-fire (LIF) model.

Figure 4 shows the step response of the three models for an identical current correlation, identical initial and final firing rates and similar membrane potential correlation times. The dynamics of the $`\theta `$-neuron and the conductance based model agree well, their response time is about $`10\text{ms}`$. The dashed line shows the linear response result, which predicts a slower response but is of the same order of magnitude. The LIF model on the other hand responds practically instantaneously. This is due to the fact that in the LIF model high input frequencies are not substantially damped, i.e. the response amplitude does not decay for large frequencies. This is impossible in the $`\theta `$-neuron in which the response amplitude always decays as $`\omega ^2`$. We would like to stress that this decay is a universal property due to the insensitivity to input currents at the point where a spike is fired and is independent of model details. Thus, although one observes that the response times of both, the $`\theta `$-neuron and the LIF model decrease with increasing $`\tau _c`$, the mathematical origin and nature of this dependence is very different. Whereas in the $`\theta `$-neuron it is a consequence of the dependence of the eigenvalues on $`\tau _c`$, in integrate-and-fire models they result from the voltage threshold.

In conclusion, we presented the linear response theory for the firing rate of the $`\theta `$-neuron in response to both time dependent input currents and time dependent noise amplitudes. For an effective numerical treatment we derived a sparse matrix representation of the Fokker-Planck operator. Using the eigenvalues and eigenfunctions of this operator, we showed that the transmission amplitude is in both cases strongly damped for frequencies above a cut-off frequency. In a wide range of parameters this cut-off frequency is always below the mean firing rate. We showed that the response behavior agrees well with the dynamics of a conductance based model neuron and is different from the behavior of the LIF model. Our results indicate that the $`\theta `$-neuron, although simple, captures well the dynamical properties of real neurons. They also reveal that previously proposed mechanisms of practically instantaneous transmission of information are incompatible with the spike generating mechanism of real neurons. Our findings suggest that the $`\theta `$-neuron reproduces the dynamical behavior of cortical neurons with Type-I excitability faithfully, while still being tractable analytically.

We acknowledge useful discussion with M. Bethge, P. Latham, A. Morrison, K. Pawelzik, M. Timme and C. v. Vreeswijk. This research was supported in part by the National Science Foundation under Grant No. PHY99-07949.


text:
# Quasi-Langmuir-Blodgett Thin Film Deposition of Carbon Nanotubes

## Abstract

The handling and manipulation of carbon nanotubes continues to be a challenge to those interested in the application potential of these promising materials. To this end, we have developed a method to deposit pure nanotube films over large flat areas on substrates of arbitrary composition. The method bears some resemblance to the Langmuir-Blodgett deposition method used to lay down thin organic layers. We show that this redeposition technique causes no major changes in the films’ microstructure and that they retain the electronic properties of as-deposited film laid down on an alumina membrane.

Carbon nanotubes have emerged as materials of fundamental importance and great application potential due to their exceptional electrical, mechanical, and thermal properties . Various proposals exist for their incorporation into devices  in single tube or thin film architectures. It has recently been found that networks of nanotubes can act as conducting channels in field effect transistors (FETs) . In addition, such films could be used in fault tolerant chemical or biological sensors , thermal heat shunts, as well as in measurements of fundamental nanotube properties in cavity and optical experiments . For such applications the preparation of uniform flat films is of paramount importance.

Although a number of methods do exist it is still far from a trivial task to prepare thin relatively uniform films. The deposition of such films from solution is difficult as nanotubes have very poor solubility in typical solvents without the use of surfactants. Strong intertube attractions, violent hydrophobicity, and low solubility at moderate concentrations all fight against typical wet chemistry techniques in making uniform films. Even if nanotubes can be suspended at a low concentrations under certain conditions , simple air drying of well suspended nanotubes on substrates results in flocculation when the local concentration approaches the solubility limit. Surfactants that make the nanotubes compatible with aqueous dispersions may be inappropriate for applications that require pure nanotubes. It is possible to deposit from dilute suspensions onto filtering membranes, but unless other steps are taken then one is constrained to use the filter as a substrate. This may be inappropriate for applications where one wants to apply gate voltages to the film (in FETs) or use optically transparent substrates. Under some situations, nanotubes can be deposited with spin coating, but for thin films ($`<`$1 $`\mu `$m), it is difficult to get adequate uniformity with such a method. A method has been found to create aligned thin films from evaporation of Triton X suspensions at the air-substrate-suspension triple line, but this technique seems to be limited to collodial dispersions and short nanotubes .

Very thin nanotube films can be grown using new chemical vapor deposition (CVD) based methods . However, this latter process requires innovative catalyst deposition and high temperatures - not compatible with CMOS technology - for the nanotube growth. These are barriers for a nanotube-CMOS integration process, in particular in the case of substrates that are not heat tolerant.

In this work we describe a method of laying down thin uniform films of carbon nanotubes on substrates of arbitrary composition that has a number of advantages over the above described techniques. A dilute nanotube solution is deposited onto an alumina membrane. The volume below the membrane is then back-filled with a fluid (typically deionized water) immiscible with the suspending liquid; the nanotube film can then be floated on an aqueous layer as a ’raft’. Upon drawing a substrate through the free liquid surface a thin uniform coating of nanotubes can be transferred.

Our method bears some resemblance to the Langmuir-Blodgett thin film deposition used to create uniform layers of organic molecules. In the Langmuir-Blodgett technique a hydrophobic group enables a thin film to lay on a free surface of a subphase compound (typically water). Although true Langmuir-Blodgett thin films have been made with nanotubes embedded in a surfactant matrix suspended on top of an aqueous subphase and then pulling the substrate through the surface , such deposition methods are not useful if the final nanotube product need be pure.

This technique has a number of important advantages over simple air drying of a liquid suspension of nanotubes. As detailed above, the liquid phase interaction between nanotubes in solution results in large flocculation effects as a nanotube suspension dries. Air drying gives totally unsuitable results, where nanotubes clump themselves in concentrated 0.1-mm ’piles’. These kind of effects even come into play when spin-coating nanotubes and prevent very thin films from being deposited uniformly.

High pressure carbon monoxide (HiPCO) process single wall nanotubes were obtained from Carbon Nanotechnologies Inc and additionally purified based on a procedure derived from Chiang et al. . Nanotubes, held in an alumina crucible, were allowed to react in moist air (obtained by bubbling ambient air through an immersed frit in room temperature water) at 225 degree Celsius for 18h (using a tubular oven; flow = 0.1 L/min). The remaining solid was suspended in concentrated HCl (37%, ACS grade Aldrich) and immersed in an ultrasound bath for 30 min. The nanotubes were then filtered, rinsed, and washed with water until a neutral pH was reached and then dried in a vacuum oven.

The purified nanotubes were ultrasonically dispersed for 30 minutes in 10:1 ratio of 1,2-dimethylbenzene (ortho-xylene) and 1,2-dichlorobenzene. The nanotube concentration in solution was $``$ 3 mg/L. It is a matter of some debate whether a true solution of nanotubes versus a suspension can be obtained. As we don’t make a distinction, we estimate the saturation concentration of nanotubes in our solvent to be approximately $``$15 mg/L as defined in Ref.  which is significantly greater than our working concentration.

A specified amount of solution (for the below images - an amount sufficient to give 8.3$`\mu `$g/cm2 which gives an $`1\mu `$m thick film) is deposited onto an alumina membrane (Whatman 0.02 $`\mu `$m pore size). A vacuum can then be engaged which allows one to remove the liquid smoothly and uniformly on a time scale short enough that flocculation and large-scale structures do not have time to form. We henceforth refer to the film formed on the alumina membrane that employs it as a substrate as ’as-deposited’. As-deposited films were also measured for comparison purpose. After the liquid has been removed, pumping continues and the film is allowed to dry for some additional time. Exact numbers are difficult to give; at our vacuum pumping rate, we dried the films for $``$ 1 minute after the excess liquid had been removed. This leaves a nanotube network with no excess fluid but which is still somewhat damp with the suspension fluid.

After deposition on the alumina membrane the chamber underneath the film is backfilled with deionized water. Our apparatus is constructed such that this water can be allowed to wash up and over the deposited nanotube film. We have found that the wetting of the alumina membrane and intended deposition substrate is aided by the addition of $``$2% isopropanol to the deionized water, although under some circumstances this caused the nanotube layer to ’shrivel’. As the fluid washes over the filter, parts of the film, still wet with the suspending fluid, will lift off the substrate and float on top of the water in the form of a single or possibly multiple large ’rafts’. The intended final deposition substrate (glass cover slides in the below images) can be slipped through the free surface and under the floating rafts, whereby the vacuum can be reengaged and the water removed, redepositing the nanotube film directly on top of the new substrate. If the as-deposited film did not immediately lift up, subsequent water washes typically dislodged it. Using the above method, we have succeeded in laying down large scale relatively uniform films with areas of up to 3.3cm2 as seen in Fig. 1.

The nanotube suspending solution must be chosen carefully. We chose to use a hybrid mixture to optimize the necessary features of high nanotube solubility, a specific gravity of less than 1, and immiscibility in water. Two of the best solvents for nanotubes 1,2-dichlorobenzene (solubility up to 95 mg/L) and chloroform ($``$ 30 mg/L)  have relatively high specific gravities (1.3 and 1.48 respectively). Such solvents are inappropriate as it is essential to our technique that the solvent have a density less than that of water, otherwise the ’raft’ will not float. A number of other candidate solvents (for instance tetrahydrofuran and dimethylformamide) are miscible in water and hence inappropriate. For this reason we chose to use a 10:1 ratio of 1,2-dimethylbenzene (ortho-xylene) and 1,2-dichlorobenzene. The solution combines the features of reasonably high nanotubes solubility (estimated to be $``$ 15mg/L), low specific gravity (0.92g/ml), and immiscibility in water.

One may have concerns regarding adverse effects the films may suffer during their immersion and redeposition. For instance the extreme hydrophobicity of the nanotubes might cause them to clump into thicker ropes upon exposure to water despite being reasonably well dispersed prior to alumina membrane deposition. That this is not a concern can be seen in Figs. 2a and 2b. Here scanning electron microscope (SEM) micrographs of films, both ones as-deposited and with the water immersion and redeposition steps, show that the films appear to consist of well separated $``$6 nanotube wide ’ropes’. There is essentially no difference in their microstructure indicating that the film morphology does not suffer any gross effects due to water immersion.

Temperature dependent DC resistivity measurements support the inference that there are no subtle structural changes which may influence the film’s transport properties. Data were taken on a number of nanotube films, both as-deposited and also ones immersed in water and then redeposited on glass slides as described above. The measurements were done $`via`$ the standard 4-probe technique with currents of 1 $`\mu `$A under 1 atm. of helium gas. Electrical contact was made via Epotek silver epoxy. In Fig. 3, is shown the resistivities of two representative films. The displayed data have been corrected for geometrical factors and this is given in ohms per square unit. Both the overall magnitude of the resistivity and its temperature dependence are unaffected by the redeposition process. The very small differences observed are within the typical variability of as-deposited films.

In conclusion we have described a method to deposit thin films of carbon nanotubes on substrates of arbitrary composition. We have shown that the structural and electronic properties of the films are essentially unaffected by this redeposition method, opening the way for the incorporation of such films into nanoscale electronic devices. While CVD growth of carbon nanotubes will most probably remain the method of choice for nanotube-silicon integration, the method described above may have promise in fabricating thin films on surfaces that cannot withstand the CVD growth environment. As such, Quasi-Langmuir-Blodgett thin film deposition may become an effective way for integration of networks into transparent or plastic substrates, which would essential for flexible electronics and display applications.

The authors would like to thank N.S. Armitage and R. Crane for various insightful conversations and S. Kwan for help in the SWNT purification.