{'fraction_non_alphanumeric': 0.03400635793453046, 'fraction_numerical': 0.02866225620757797, 'mean_word_length': 4.493298093260336, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Explaining the subpulse drift velocity of pulsar magnetosphere within the space-charge limited flow model

We try to explain the subpulse drift phenomena adopting the space-charge limited flow (SCLF) model and comparing the plasma drift velocity in the inner region of pulsar magnetospheres with the observed velocity of drifting subpulses. We apply the approach described in a recent paper of van Leeuwen&Timokhin (2012), where it was shown that the standard estimation of the subpulse drift velocity through the total value of the scalar potential drop in the inner gap gives inaccurate results, while the exact expression relating the drift velocity to the gradient of the scalar potential should be used instead. After considering a selected sample of sources taken from the catalog of Weltevrede, Edwards&Stappers (2006) with coherently drifting subpulses and reasonably known observing geometry, we show that their subpulse drift velocities would correspond to the drift of the plasma located very close or above the pair formation front. Moreover, a detailed analysis of PSR B0826-34 and PSR B0818-41 reveals that the variation of the subpulse separation with the pulse longitude can be successfully explained by the dependence of the plasma drift velocity on the angular coordinates.


INTRODUCTION
The research on pulsar magnetospheres started 45 years ago from the pioneering paper of Goldreich & Julian (1969), where it was shown that the magnetic field, together with the fast rotation of the pulsar, generates strong electric fields tending to pull out charged particles from the surface of the neutron star and to form the plasma magnetosphere. However, the very next arising question -how many particles will actually leave the surface of the star under the action of this force -is still a subject of scientific debate. In fact, the magnetic field which generates the pulling out electric field also leads to the substantial increase in the cohesive energy of the surface charged particles (positive ions in larger degree than electrons), making the outer layer of the star very dense and strongly bound. Medin & Lai (2007) showed that for each magnetic field intensity there exists a critical surface temperature, above which the particles are able to freely escape from the surface of the star. The magnetic field of pulsars is typically inferred from observations under the assumption that it is strictly dipolar, thus providing B d = 2 × 10 12 PṖ × 10 15 G (here P is the period andṖ is the period derivative of the pulsar). Using this formula, it turns out that the majority of pulsars satisfy the condition for the free particle outflow. However, in many works, starting from Ruderman & Sutherland (1975) it was suggested that the magnetic field close to the surface of the pulsar should have a multipole structure, with the surface magnetic field several orders of magnitude larger than the estimated B d . This idea is supported by X-ray observations (Zhang, Sanwal & Pavlov 2005;Kargaltsev, Pavlov & Garmire 2006;Pavlov et al. 2009), and several studies have been performed searching for a mechanism of generation and maintenance of such small-scale strong magnetic fields [see Geppert, Gil & Melikidze (2013) and references therein].
The amount of charged particles extracted from the surface of the star by the rotationally induced electric field is a key aspect of any pulsar magnetosphere model. The model of Ruderman & Sutherland (1975), for example, assumes that no particle leaves the pulsar surface and there is a vacuum gap formed above the star with a huge difference in the scalar potential between the bottom and the top, i.e., ∼ 10 12 V. According to this idea, the gap will be periodically discharged and rebuilt, making it intrinsically "non-stationary". On the contrary, in the "stationary" space-charge limited flow (SCLF) model of Arons & Scharlemann (1979) there is a free flow of charged particles from the pulsar surface. The partially screened gap model of Gil, Melikidze & Geppert (2003) combines some features of both the previous models and is based on the assumption that the vacuum gap discharge will lead to the back flow of charged particles bombarding and heating the surface of the star, causing the thermal ejection of ions from the surface and partially screening the original accelerating electric field.
The choice among different models of pulsar magnetosphere should be made by accurate comparison of their predictions with the results of observations. A very interesting phenomena serving as a diagnostic tool is the subpulse drift. Pulsar radio emission comes to us in the form of pulses, which may look very different among each other and typically consist of individual subpulses. Although the average pulse profile is very stable and represents a unique fingerprint of each pulsar, a range of successive pulses plotted on top of each other in a so-called pulse stack, quite often shows an organized phase shift of the subpulses forming drift bands. This phenomenon was reported for the first time in Drake & Craft (1968), while the first systematic analysis of "drifting" pulsars dates back to Backus (1981) and Rankin (1986). So far, the largest statistical study of the phenomenon has been presented in the works of Weltevrede, Edwards & Stappers (2006) and Weltevrede, Stappers & Edwards (2007), who considered a sample of 187 pulsars, 55% of which show drifting subpulses. Usually the subpulse drift bands are characterized by the horizontal separation between them in subpulse longitude, P2, and the vertical separation in pulse periods, P3. The subpulse behavior of the individual pulsars may be rather complex and demonstrate smooth or abrupt change of drift direction, phase steps, longitude and frequency dependence of the separation P2, or even presence of subpulses drifting in the opposite directions at the same time. For the graphical representation of the periods P2 and P3 we refer the reader to Fig. 1 of Weltevrede, Edwards & Stappers (2006), where one can find several examples of the different subpulse behavior of individual sources.
The vacuum gap model, and especially the partially screened gap model, has been widely used to explain the subpulse drift phenomena (Gil & Sendyk 2000;Melikidze, Gil & Pataraya 2000;Gil, Melikidze & Geppert 2003;Gupta et al. 2004;Bhattacharyya et al. 2007;Gil et al. 2008), while for a long time the SCLF model has been regarded unable to account for it. However, recent analytical (van Leeuwen & Timokhin 2012) and numerical (Timokhin 2010b;Timokhin & Arons 2013) progresses have shown that the door can be left open even for the SCLF model. The main goal of this paper is to explain the subpulse drift velocity in the framework of the relativistic SCLF model, without addressing the question of the generation mechanism of the plasma features responsible for the appearance of the subpulses, while trying to compare the results with the available observational data.
The plan of the paper is the following. In Sect. 2 we give a brief review of the models used so far to explain the subpulse phenomena and motivate our choice to concentrate on the SCLF model. In Sect. 3 we present the basic equations to explain the subpulse drift velocity in the framework of the SCLF model. In Sect. 4 we consider a set of pulsars from the catalog of Weltevrede, Edwards & Stappers (2006) with coherently drifting subpulses, and try to deduce in which regions of the pulsar magnetosphere the SCLF model would predict the plasma with the observed velocities. In Sect. 5 we focus on two specific sources, PSR B0826-34 and PSR B0818-41, trying to account for their phenomenology. Finally, Sect. 6 is devoted to the summary of the results obtained and to the conclusions.

A BRIEF SURVEY OF EXISTING MODELS
The first theoretical explanation of the subpulses was provided by Ruderman & Sutherland (1975), who associated the subpulses to the spark discharges of the vacuum gap above the pulsar surface. In its original form, the model applied to the pulsars with anti-parallel angular velocity Ω and magnetic moment µ, and assumed that the charged particles (positive ions) are tightly bound to the surface of the star and cannot be pulled out by the rotationally induced electric field. This requirement leads to the formation of a vacuum gap in the region where the magnetic field lines are open and with a potential drop between the top and the bottom of the order of 10 12 V for typical pulsar parameters. Due to the presence of strong curved pulsar magnetic fields, the gap will be unstable and periodically discharged by the photon induced pair creation process. The discharges will build up plasma columns, which are subject to the E × B drift in the electromagnetic field of the magnetosphere. Ruderman & Sutherland (1975) showed that, unless the potential drop of the gap is completely screened, the plasma columns will not exactly corotate with the star but always lag behind the rotation of the star, and this is responsible for the visual drift of the subpulses along the pulse longitude. The sparks are assumed to form rings and the so-called tertiary periodicity P4 is the time needed for the spark carousel to make one full rotation around the magnetic axis. Although in this model the subpulses cannot outrun the rotation of the star, due to the effect of aliasing (Gil, Melikidze & Geppert 2003;Gupta et al. 2004), the apparent velocity of the subpulses may be both positive (from earlier to later longitudes) and negative (from later to earlier longitudes). Ruderman & Sutherland (1975) estimated the subpulse drift velocity to be proportional to the full potential drop across the gap, resulting in excessively large values of the drift velocity compared with the observed ones. Later on, Gil & Sendyk (2000), Melikidze, Gil & Pataraya (2000), Gil, Melikidze & Geppert (2003) generalized this model to account for arbitrary inclination angles χ between Ω and µ, and modified it to allow for the partial outflow of ions and electrons from the surface of the star, forming partially screened gap instead of the pure vacuum. It was argued that the favorable conditions for the spark discharge persist even if the original vacuum gap is screened up to 95% or more, making the velocity of drifting subpulses consistent with the observed values. However, the partially screened gap model requires surface values of the magnetic field of the order of 10 14 G, much larger than those deduced when the magnetic field is dipolar, 1 i.e. ∼ 10 12 G. The partially screened gap model has been used in a number of works to describe the subpulses of specific pulsars as well as their X-ray emission (Gupta et al. 2004;Bhattacharyya et al. 2007;Gil et al. 2008) and it has typically revealed a strong predictive power.
In Clemens & Rosen (2004); Rosen & Clemens (2008); Rosen & Demorest (2011), the drifting subpulses are instead explained by non-radial oscillations of the surface of the star. This model gives a very natural explanation to the subpulse phase shift, relating it to the intersection of the observer's line of sight with the nodal line. The empirical model of Wright (2003) relates the formation of drifting subpulses to the interaction between electron and positron beams traveling up and down between the inner and the outer gaps of the pulsar magnetosphere. Kazbegi, Machabeli & Melikidze (1991) and later Gogoberidze et al. (2005) proposed a model where the subpulses are formed due to the modulation of the emission region by large-scale "drift waves", generated by oppositely directed curvature drifts of electrons and positrons. Finally, Fung, Khechinashvili & Kuijpers (2006) proposed a model to explain drifting subpulses that is based on the diocotron instability in the pair plasma on the open field lines. Detailed description of all these models may be found in the review of Kuijpers (2009) and references therein.
Recently van Leeuwen & Timokhin (2012) have shown that the order of magnitude estimation of the subpulse drift velocity used in Ruderman & Sutherland (1975), and subsequent works, can be replaced by a more precise expression, relating the velocity to the radial derivative of the potential instead of the absolute value of the potential drop. This simply comes from the fact that the subpulse drift velocity is determined by the E × B drift of the distinct plasma features in the magnetosphere, while the electric field responsible for this drift is given by the gradient of the scalar potential. The final expression for the subpulse drift velocity in degrees per period is (van Leeuwen & Timokhin 2012) whereṼ ≡ V /∆Vvac is the scalar potential normalized to the potential drop ∆Vvac = ΩBr 2 pc /2c between the rotation axis and the boundary of the (small) polar cap, rpc is the radius of the polar cap, B is the value of the magnetic field (which van Leeuwen & Timokhin (2012) assumed to be constant across the polar cap), ξ ≡ θ/Θ is the polar angle normalized to the colatitude of the polar cap boundary Θ, c is the speed of light. As emphasized in van Leeuwen & Timokhin (2012), this velocity, being the true plasma drift velocity in the open magnetic field line region, should introduce periodic modulation to the observed radio emission from pulsars, independently on the particular emission generation mechanism. We stress that this periodic modulation (which is exactly the tertiary periodicity, mentioned before) can be found even in the spectra of those pulsars which don't reveal regularly drifting subpulses Gil, Melikidze & Zhang 2007). The same is seen for many pulsars listed in Weltevrede, Edwards & Stappers (2006) and Weltevrede, Stappers & Edwards (2007).
The main reason for which the SCLF model has been so far regarded unable to explain the subpulse drift phenomena is that it has not a prescribed mechanism for the formation of spark-like features, which are invoked to explain the subpulses in the vacuum and in the partially screened gap models. However, recent progress in numerical simulations of pulsar magnetosphere have shed new light on this subject. One of the main tendencies in the numerical studies of pulsar magnetosphere is to consider it as a global object, with the different regions closely interlinked and interdependent. For example, self-consistent simulations of pair cascades in the polar cap region (Timokhin 2010b;Timokhin & Arons 2013) show that the cascade behavior is mostly determined by the global magnetospheric current density and that periods of plasma generation are interleaved with quiet periods, both for the vacuum gap and for the SCLF regimes. Based on that, van Leeuwen & Timokhin (2012) proposed the idea that the distinct emitting features in the inner magnetosphere, responsible for the appearance of the subpulses, may be caused by global current filaments, similar to those observed in auroras.
These arguments motivated us to adopt the approach of van Leeuwen & Timokhin (2012) for studying the possible subpulse behaviour within the SCLF model. In the rest of the paper we do not consider the problem of the subpulse generation. Rather, we concentrate on the velocity given by Eq. (1), calculated from the analytical expressions for the scalar potential in the relativistic SCLF model, which are available in the literature (see Muslimov & Tsygan (1992); Harding & Muslimov (1998, 2002). The main issue that we try to address is whether it is possible, in the framework of the SCLF model, to explain the observed subpulse drift velocities, and if so, to infer in which part of the magnetosphere they should be produced. In the second part of the paper we will instead apply our arguments to two specific sources.

SUBPULSE DRIFT VELOCITY IN THE FRAMEWORK OF THE SCLF MODEL
Scharlemann, Arons & Fawley (1978) and Arons & Scharlemann (1979) were the first to show analytic solutions for the scalar potential in the vicinity of the pulsar polar cap and in the framework of the SCLF model. In their analysis, the accelerating electric field parallel to the magnetic field of the pulsar is due to the curvature of magnetic field lines and to the inertia of particles. Later, Muslimov & Tsygan (1992) have shown that, due to the effect of dragging of inertial frames in general relativity, it is possible to obtain accelerating electric fields which are two orders of magnitude larger than those normally expected. This approach has further been elaborated in Harding & Muslimov (1998, 2002. For convenience, in this subsection we present the main results found by Muslimov & Tsygan (1992) as well as the expressions for the subpulse drift velocity that we obtained, i.e., using Eq. (1). In general relativity the dipole-like magnetic field in the exterior spacetime close to the surface of a slowly rotating neutron star described by the metric ds 2 = −N 2 c 2 dt 2 + N −2 dr 2 +r 2 dθ 2 + r 2 sin 2 θdφ 2 − 2ωr 2 sin 2 θdtdφ (2) is given by the expressionŝ Here the spherical coordinates (r, θ, φ) are used with the polar axis oriented along the magnetic moment of the pulsar,r = r/R, R is the radius of the neutron star, B0 = 2µ/R 3 is the value of the magnetic field at the pole, N = (1 − 2GM/rc 2 ) 1/2 is the lapse function of the metric, G is the gravitational constant, M is the mass of the star, ω is the frequency of dragging of inertial frames, ε = 2GM/Rc 2 is the compactness parameter, while the function f (r) is given by The polar angle of the last open magnetic field line Θ is equal to where is the polar angle of the last open magnetic field line at the surface of the star. The scalar potential Φ in the polar cap region of the inner pulsar magnetosphere is obtained from the solution of the equation where ρGJ is the Goldreich-Julian charge density and ρ is the actual charge density in the open field line region of the magnetosphere. The boundary conditions of this Poisson equation are (i) Φ = 0 at the surface of the star and along the last open magnetic field lines, and (ii) E = 0 at large distances from the star. Very close to the surface of the star (r − 1 << 1) the solution reads while at the distances Θ0 <<r − 1 << c/ΩR one gets Here Φ0 = ΩB0R 2 /c, κ ≡ εβ, while β = I/I0 is the stellar moment of inertia in units of I0 = M R 2 . The parameter ξ = θ/Θ changes from 0 to 1 inside the polar cap region, Jm is the Bessel function of order m, ki andki are the positive zeroes of the Bessel functions J0 and J1, arranged in ascending order. Moreover, These results allow one to find the plasma drift velocity in the polar cap region of the magnetosphere with the electric field E = −∇Φ and the magnetic field (3)-(4). One can easily show that the largest contribution to the azimuthal drift in the corotating frame of the star is due to the term −(cE θ Br/B 2 )φ, which, after proper transformations (see subsection 2.2 of van Leeuwen & Timokhin (2012) for the details), leads to the subpulse drift velocity in degrees per period as given by (1). The final expressions for the drift velocity, obtained from (8) and (9) and For inclination angles χ < 90 • (except for the almost orthogonal pulsars) the scalar potential in the polar cap region is positive, has a maximum close to the magnetic axis and goes to zero at the last open magnetic field lines, so that the value of ωD is negative almost everywhere. From the point of view of observations it means that the SCLF model predicts negative drift (from larger to smaller longitudes) in case of the outer line-of-sight geometry and positive drift in case of the inner line-of-sight geometry. In the rest of our work we claim that the velocities (12) and (13) represent the true drift velocities of whatever features are responsible for the subpulses in a specific portion of the magnetosphere [see also van Leeuwen & Timokhin (2012)]. Moreover, the expressions (12) and (13) predict the longitude dependent (not constant) apparent drift velocity of the subpulses along any observer's line of sight, unless the inclination angle of the pulsar is exactly zero and the line of sight is exactly concentric with the magnetic field axis. In section 5 we will use this fact to explain the longitude dependence of the subpulse separation in case of two individual pulsars.  Weltevrede, Edwards & Stappers (2006) for which the inclination angle χ is known from previous studies. a Value taken from Rankin (1993). b Value taken from Lyne & Manchester (1988).

COMPARISON WITH THE OBSERVED VELOCITIES OF THE DRIFTING SUBPULSES
Equations (12) and (13) give the subpulse drift velocity within the SCLF model. In order to check whether these expressions predict numbers in agreement with observations, we have used the data from the catalog of Weltevrede, Edwards & Stappers (2006) and Weltevrede, Stappers & Edwards (2007). These authors present the results of the observations of 187 pulsars in the northern hemisphere at wavelengths of 21 cm and 92 cm. Pulsars revealing the drifting subpulses phenomenon are divided into three classes, depending on the character of their Two-Dimensional Fluctuation Spectrum (2DFS). Coherent drifters (marked as Coh) have narrow pronounced feature in their 2DFS spectra, meaning that P3 has a stable value through the observations. Diffuse drifters of the classes Dif and Dif * have a broader feature in 2DFS spectra, which for Dif pulsars is clearly separated from the alias borders of the spectra (P/P3 = 0 and P/P3 = 0.5), while for Dif * pulsars is not (see Weltevrede, Edwards & Stappers (2006) for more details and examples). For our purposes we considered the coherent drifters from Tab. 2 of Weltevrede, Edwards & Stappers (2006) (corresponding to the observations at 21cm), selecting only those with known inclination angle χ. The values of the inclination angles were taken from Rankin (1993) and, if absent there, from Lyne & Manchester (1988). The resulting sample is reported in Table 1. When two values of the periods P2 and P3 were given in Weltevrede, Edwards & Stappers (2006), we chose the first one.
We assumed all pulsars to have the typical numbers for compactness ε = 0.4, κ = 0.15, and stellar radius R = 10 6 cm. A clear picture of drifting subpulses is observed when the line of sight of the observer grazes the emission cone, so that it is reasonable to take ξ = 0.9. One may also notice that the second term in (12) and (13), containing the dependence on φ, is smaller than the first term (it depends on a higher degree of the small angle Θ0) and plays a role mostly for the pulsars with large inclination angles. Hence, for the purposes of this subsection we fixed 2 φ = π. Under these assumptions, the drift velocities (12) and (13) for each individual pulsar depend only on the radial coordinater. So, by solving numerically the equation for each pulsar of Table 1, we can find the altituder0 − 1 of the plasma features that are responsible of the subpulses. When solving the equation (14) in the low altitude approximation we took the first 30 terms of the expansion (8), which reduces the error to less than one percent.
In Table 1 we have pulsars with both positive and negative subpulse drift velocities, and no preferred direction of the drift (sign of P2) was found in Weltevrede, Edwards & Stappers (2006) and Weltevrede, Stappers & Edwards (2007). For comparison, we report in the table also the values of the impact angle β, taken from Rankin (1993) and Lyne & Manchester (1988), which is the angle of the closest approach between the magnetic axis and the line of sight. As we already mentioned, the SCLF model predicts negative drift for the outer line of sight (positive β) and positive drift for the inner line of sight (negative β). However, from Table 1 we don't see a correlation between the signs of β and P2. According to the vacuum/partially screened gap model, the discrepancy between the predicted and the observed direction of the drift is usually explained in terms of aliasing (Gupta et al. 2004), but an alternative explanation is possible. Expressions (12) and (13), depending on the inclination angle χ, give negative drift velocities for pulsars with inclination angle in a range from zero to almost 90 • and positive velocities for the inclination angles from ∼ 90 • to π (pulsars with inclination angles close to 90 • need to be considered separately). This suggests that pulsars with positive subpulse drift velocities may have Ω · B < 0, while those with negative velocities have Ω · B > 0 (provided that the line of sight geometry is outer in all cases).
One may assume that in reality any of the three proposed explanations for the visible direction of the drift (aliasing, line of sight, inclination angle) may be applicable depending on the individual pulsar properties. There are pulsars showing different sign of the drift velocity in different modes, or in the same mode, like B0826-34 or J0815+09, which obviously cannot be explained by the change of the inclination angle. However, keeping this in mind as a possible scenario, we have performed the following analysis.
In case when the subpulse drift velocity of the pulsar is positive, we change the inclination angle χ reported in Table 1 into π − χ in order to solve the equation (14), getting a meaningful value for the altitude (r > 1). The pulsar B2045-16 in Table 1 breaks this assumption, having opposite sense of the subpulse drift for the observations at 21 cm and 92 cm. At least six more pulsars among those reported by Weltevrede, Edwards & Stappers (2006) and Weltevrede, Stappers & Edwards (2007) show a similar phenomenology. However, all of them are diffuse drifters (class Dif * ) at 92 cm (including B2045-16 itself), suggesting that aliasing is very likely to occur for them. Changing the inclination angle from χ to π − χ in our calculations for the majority of pulsars (apart from the almost orthogonal ones, where the term with sin χ in expressions (12), (13) starts to play a role) is essentially analogous to considering the inner line of view instead of the outer one. If we simply put negative sign in front of all drift velocities from Table 1, we get the same values ofr0, solving Eq. (14). The results of our analysis are schematically represented in Fig. 1. The red points represent the values of the altitude above the surface of the star in units of stellar radii (r0 − 1) obtained from Eq. (14), with the error bars calculated from the errors of P2 and P3 given in Table 2 of Weltevrede, Edwards & Stappers (2006) (at 21 cm). The green points represent the same quantity obtained from the observations at 92 cm (Weltevrede, Stappers & Edwards (2007)). The pulsars are arranged in ascending order of the period. Black stars show the angular radii of the polar caps at the surface of the star, Θ0, for comparison. Pulsars which are absent in the 92 cm catalog are marked with "no data". Arrows indicate those values of the altitude too large to fit in the plot.
For the sake of comparison, we have also computed the altitudes as given by the heights of the pair formation front (PFF) following Hibschman & Arons (2001); Harding & Muslimov (2001, 2002. The expressions for the PFF location in the framework of the SCLF model are given in the Appendix. We recall that the main processes determining the height of the PFF in pulsar magnetosphere are curvature radiation and inverse Compton scattering (resonant as well as non-resonant). In Fig. 1 the shaded blue regions indicate the PFF altitude from Hibschman & Arons (2001) for the curvature radius of the surface magnetic field lines starting from the values ∼ R (lower boundary of the shaded zones) up to the values of the standard dipolar magnetic field (upper boundary of the shaded zones). The blue points show the PFF altitudes from Harding & Muslimov (2001, 2002. The temperature of the star is taken to be T = 2 × 10 6 K. Even though our approach fails for some of the considered pulsars, we find it very promising that for most of them the altitude corresponding to the observed value of the subpulse drift velocity lies within the range or slightly above the PFF altitudes. It is a well-known fact that radio emission in pulsar magnetosphere originates at heights of the order of tens of stellar radii (Kramer et al. 1997;Kijak 2001;Kijak & Gil 2003). However, it was assumed in Ruderman & Sutherland (1975) and subsequent works on drifting subpulses that the geometrical pattern and drift velocity of the subpulses are determined by the distribution of sparks within the gap, i.e., very close to the surface of the star. In this respect we think that the heights that we have obtained, so close to the PFF values, are not occasional and may reflect the fact that, whatever is the mechanism of subpulse generation in the SCLF model, they are likely to form in the vicinity of the pair formation front.
Our analysis was not meant to consider systematically all the sources in the catalog of Weltevrede, Edwards & Stappers (2006), but rather it was aimed at demonstrating the potential ability of the model in accounting for the observations. While we are confident that the main conclusion of this subsection may be applicable also to other pulsars, a systematic study of all the sources may become necessary.

DISCUSSION OF SPECIFIC SOURCES
The expressions for the subpulse drift velocity (12) and (13) derived in the framework of the SCLF model naturally contain a dependence on ξ and φ, and they predict different velocities for different regions of the polar cap, in contrast to the estimations of the vacuum gap model. In this subsection we attempt to exploit these additional degrees of freedom to explain the variability of the subpulse velocities along the pulse longitude in case of two specific pulsars, i.e., PSR B0826-34 and PSR B0818-41. Although not included in the catalog of Weltevrede, Edwards & Stappers (2006), both of them have been repeatedly investigated at several observing frequencies, and, since they have wide profiles allowing to track several subpulse drift bands at a time, they can be regarded as ideal test cases.

Basic parameters
The pulsar B0826-34, with spin P = 1.8489 s andṖ = 1.0 × 10 −15 , has an unusually wide profile, extending through the whole pulse period. The pulsar emits in its strong mode for 30% of the time. For the rest of the time, the pulsar stays in the weak mode, with an average intensity of emission which is ∼ 2% of the emission of the strong mode (Esamdin et al. 2005;Serylak 2011;Esamdin et al. 2012). Because of its weakness, the very existence of the weak mode was confirmed only very recently and for a long time it was thought to be a null pulsar (Durdin et al. 1979;Biggs et al. 1985;Bhattacharyya, Gupta & Gil 2008).
The average pulse profile of B0826-34 consists of the main pulse (MP) and the interpulse (IP), separated by regions of weaker emission. The intensity of the MP is much larger than the intensity of the IP at the frequencies 318 MHz and 606 MHz, while at frequencies larger than ∼ 1GHz the IP starts to dominate. The MP itself has a double peaked structure with a separation between the peaks decreasing at higher frequencies, following the common trend described by the radius-to-frequency mapping model of Kijak & Gil (2003). A detailed description of the average profile evolution with frequency may be found in Gupta et al. (2004), Bhattacharyya, Gupta & Gil (2008).
Additional relevant physical parameters are those related to the observing geometry of PSR B0826-34, i.e., the values of the inclination angle χ and of the impact angle β. Usually the values of these angles are determined by fitting the polarization profile of the pulsar. According to the "rotating vector model" of Radhakrishnan & Cooke (1969), the polarization angle of the pulsar radio emission is equal to where l is the pulse longitude, related to the azimuthal coordinate φ by means of standard theorems of spherical geometry (Gupta et al. 2004) 3 as However, in many cases this method does not give a unique value for the inclination and for the impact angles, but rather a wide range of possible combinations (Miller & Hamilton 1993). For example, early estimations of Biggs et al. (1985) for PSR B0826-34 based on the polarization measurements suggested a large range for χ and β with the best fit of χ = 53 • ±2 • and β+χ = 75 • ±3 • , a fact which does not agree with the large width of the profile. In Gupta et al. (2004) these angles were estimated from the polarization information 4 together with the frequency evolution of the profile and were found to lie in the range 5 1.5 • χ 5.0 • and 0.6 • β 2.0 • .

Subpulse phenomenology
The drifting subpulses are seen almost along the whole range of longitudes, showing from 5-6 up to 9 visible tracks at a time. The character of the observed subpulse drift is irregular, the apparent velocity reveals an oscillatory behavior, changing sign with a periodicity of the order of tens to several hundred periods of the pulsar. The values of the apparent subpulse drift velocities were measured in a number of papers for different observing frequencies (Gupta 3 In the corresponding formula of Gupta et al. (2004) the azimuthal coordinate with respect to the magnetic axis is denoted by σ, while φ is used to denote the pulse longitude, which may introduce a confusion when making comparison with our results. However, we preferred to keep the notation φ for the azimuthal angular coordinate as in Muslimov & Tsygan (1992). 4 We recall that the measurements of linear polarization are sometimes contaminated by the orthogonal polarization mode switching. Single pulse polarization observations may be necessary for the reconstruction of the mode-corrected polarization angle swing (Gil & Lyne 1995).  Table 2. In all cases the declared velocity range is not symmetric with respect to zero, having a positive average velocity (see Figs. 5 and 6 of Esamdin et al. (2005) for an example of a positive average drift). This phenomenology suggests the interpretation that the average drift velocity reflects some "undisturbed" characteristic of the pulsar, while the oscillatory change of the velocity is determined by some sort of perturbation (see the discussion at the end of this subsection for possible explanations). Another interesting property of PSR B0826-34 is the longitude dependence of the period P2. Esamdin et al. (2005) reports the values for P2 between 26.8 • and 28 • in the MP region (average 27.5 • ) and between 19 • and 23.5 • in the IP region (average 22.2 • ). In order to explain this behaviour the authors proposed a model of spark carousel consisting of two rings of 13 sparks each, with the separation between the sparks in the outer ring (responsible for the MP) 27.5/22.2 ≈ 1.2 larger than in the inner ring (responsible for the IP). In Gupta et al. (2004) P2 was found to vary between 21.5 • and 27 • with the mean value of 24.9 • . As a possible explanation of the observed phenomena the authors proposed a scheme where the ring of sparks is centered not around the dipolar axis of the pulsar magnetic field, but around the "local magnetic pole", shifted with respect to the global dipole one. This agrees well with the models suggesting that the pulsar magnetic field has a multipole structure near the surface of the neutron star (Gil & Sendyk 2000), which arises due to the dynamo mechanism in the newborn stars (Urpin & Gil 2004) or, more probably, due to the Hall drift (Geppert, Rheinhardt & Gil 2003;Geppert, Gil & Melikidze 2013).

Analysis of the subpulse drift
We start our analysis of the subpulse drift by estimating the altitude above the surface of the star, which would correspond to the average subpulse drift velocities of PSR Sec. 4]. Since the average drift velocities are positive at all the observing frequencies, while our model predicts negative values of the velocity for small angles χ and outer line of sight, we are lead to infer that PSR B0826-34 is anti-parallel, with χ close to 180 • . Assuming for the moment χ = 180 • and ξ = 0.5, we have computed from Eq. (14) the altitudes corresponding to the values of the average drift velocities and we have reported them in the last column of Table 2. The results depend very weakly on the chosen ξ and χ, provided the latter is close to 180 • . For comparison, when the temperature of the polar cap is T = 3 × 10 6 K, the altitude of the PFF h HA lies in the range 0.0035 − 0.036 for the different values of fρ and h HM = 0.062 (for the definitions of h HA and h HM see the Appendix). On the contrary, when the temperature is T = 4 × 10 6 K, the corresponding values for the h HA are 0.0026 − 0.029, with the same h HM . Hence, the values ofr0 − 1 that we have obtained lie within the corresponding range of values deduced through the PFF approach, and they are typically closer to its lower boundary.
A closer look at the resultingr0 − 1 shows that they satisfy the frequency scaling suggested by the radius-to-frequency mapping, i.e., rem ∼ ν −0.26±0.09 (Kijak & Gil 2003), where rem is the height of the radio emission generation and ν is the observed frequency. To demonstrate it, we have fitted the four data points of Table 2 by the function Aν −0.26 , getting A = 0.026. In Fig. 2 we plot the four considered data points with the error bars calculated under the assumption that the accuracy of the measurements of the apparent subpulse drift velocity is of the order of ∼ 0.05 • /P together with the obtained fit (solid blue line). The two dashed lines correspond to the fits A1ν −0.26+0.09 , with A1 = 0.014, and A2ν −0.26−0.09 , with A2 = 0.046. As we already mentioned, it is commonly accepted that the generation of radio emission takes place in the upper region of the magnetosphere at altitudes ≈ (10 − 60) R (Kramer et al. 1997;Kijak 2001;Kijak & Gil 2003). Our results tend to confirm the existence of a radius-tofrequency scaling, although the lack of more robust observational data does not allow us to draw stronger conclusions. Fig. 3 is devoted to the illustration of the geometry of PSR B0826-34. In the upper panel we show a sketch of the polar cap, whose boundary is represented with a black circle. At the altituder 0 − 1 = 0.0044, the polar cap has an angular radius Θ ∼ 0.5 • [cf., Eq. (7)]. It should be stressed that different observations of the angular size of the polar cap do not provide consistent conclusions. According to Lyne & Manchester (1988), for instance, Θ ∼ 13 • P −1/3 /2, which, for PSR B0826-34, gives Θ ≈ 5.3 • , i.e., an order of magnitude larger compared to what we have found. On the contrary, the angular size of the polar cap deduced from Xray observations is much smaller and close to the values that we have obtained through Eq. (7).
Another relevant quantity is the trajectory of the observer's line of sight, which is given by (Manchester & Taylor 1977) Lyne & Manchester (1988) 6 . The lower panel of Fig. 3 reports the polarization profile computed through Eq. (15). It should be compared with Fig. 6d of Biggs et al. (1985) and with the upper panel of Fig. 1 of Gupta et al. (2004) 7 .
An additional quantity that our model can potentially reproduce is the ratio between the velocities (12) at φ = 0 and φ = π along the line of sight, which, for the considered geometry, is equal to 1.22. This value is very close to the observed ratio of the subpulse separation P2 in the MP to the subpulse separation P2 in the IP, namely 27.5/22.2 ≈ 1.239 (Esamdin et al. 2005). This suggests us that P2 is proportional to the velocity ω at any given point of the polar cap. Intuitively, if one imagines the carousel of sparks (or any other feature responsible for the subpulse phenomena) moving around the polar cap with the longitude dependent velocity, it seems plausible that in the regions with lower velocity the subpulse tracks will tend to look closer, while in the regions with larger velocity they will look farther away from each other. Thus, in our model we associate the region around φ = 0 with the MP, while the region around φ = π with the IP. Concerning the value of P3, 6 This makes our picture for the observed geometry a bit different from Fig. 10 of Esamdin et al. (2005), where the angle β is apparently larger than χ. 7 Note that the plots of the polarization angle and of apparent drift velocity in this subsection are shifted in longitude in order to make them easier comparable with the corresponding plots in the literature, associating zero of the pulse longitude with the bridge region before the IP. In our analysis the zero of the azimuthal angle and of the pulse longitude is instead associated with the peak of the MP. The shift is taken to be 105 • to match the distance between the second zero of the position angle curve and the end of the pulse from the top panel of Fig.1 of Gupta et al. (2004).  which, we recall, also enters the expression for the apparent drift velocity, we could not find it in the literature for PSR B0826-34, possibly because the subpulse tracks are irregular. However, from the observations of other pulsars we know that the values of P3 are essentially the same at different observing frequencies for a given pulsar (Weltevrede, Edwards & Stappers 2006;Weltevrede, Stappers & Edwards 2007;Bhattacharyya, Gupta & Gil 2009), suggesting that this subpulse parameter is more "stable". We therefore argue that the longitude dependence of P2 is explained in terms of the longitude dependence of the drift velocity of the features that generate the subpulse, while P3 could represent a specific characteristic property of the individual pulsar.
As a further step, we have tried to account for the longitude dependence of P2 reported in Gupta et al. (2004) and our results are reported in Fig. 4. The solid blue curve shows the drift velocity (12) projected on the trajectory of the line of sight across the polar cap as a function of the pulse longitude l. The value ofr taken for this plot corresponds to the average drift velocity 0.55 • /P , as from Gupta et al. (2004). The grey and black sets of data points have been obtained in the following way. We first started from the ob-  Figure 4. Comparison of the variation of the plasma drift velocity along the pulse longitude (solid blue curve) with the corresponding variation of the period P 2 of PSR B0826-34, measured in Gupta et al. (2004). In the absence of the measured P 3 for this pulsar, the values of P 2 from the upper panel of the Fig.8 of Gupta et al. (2004) are divided on 44 (grey data points) and shifted on 50 • (black points), to bring them in better visual correspondence (these parameters are not a result of fit, but estimated by eye). served values of P2 given in the top panel of Fig. 8 of Gupta et al. (2004) along the pulse profile and shifted them for 105 • to bring in correspondence with our longitude scale. Then, since the value of P3 is not available for PSR B0826-34, we divided the data values (including the error bars) by a number (44 in this specific case) which brings the data points and the blue curve in better visual overlap, and we plotted them in grey color. The black points are the same as the grey ones, but shifted by 50 • in longitude. Although the value of this shift is also chosen "ad hoc" and not the result of a fitting procedure, there is a positive uncertainty in our choice of the origin of the longitude (see the footnote 7), as well as a possible longitude shift between the 318MHz data of Gupta et al. (2004) and the 1520MHz data of Esamdin et al. (2005), that we have chosen as a reference. In spite of the observational uncertainties, the correspondence that we have found between the data points and the analytical curve is very promising and needs further investigation and comparison with more data.
Finally, as already mentioned before, the measured subpulse drift velocity of PSR B0826-34 reveals an irregular behavior on the timescales of tens to hundreds of period, for which a firm explanation is still lacking. Gupta et al. (2004) explained these variations within the partially screened gap model by invoking small fluctuations of the polar cap temperature around the mean value 8 (Gil, Melikidze & Geppert 2003). On the contrary, van Leeuwen & Timokhin (2012) argued that the potential drop in the polar cap region may be determined not only by the local physical conditions, but by the global structure of the magnetosphere (Li, Spitkovsky & Tchekhovskoy 2012;Kalapotharakos et al. 2012;Timokhin 2010a,b). As a result, the long term changes of the observed drift rate may be related to the evolution of the magnetospheric current density distribution, for example, due to switches between metastable magnetosphere configurations (Timokhin 2010a).
The alternative explanation that we propose is that the observed variations are related to stellar oscillations. In our preceding research (Morozova, Ahmedov & Zanotti 2010;Zanotti, Morozova & Ahmedov 2012) we have studied the influence of the non-radial stellar oscillations on the scalar potential of the polar cap region of the magnetosphere. The oscillation velocity at the stellar surface modulates the linear velocity of the pulsar rotation, introducing a new term in the charge density, in the scalar potential and in the accelerating electric field above the surface of the star. We also shown that oscillations may increase the electromagnetic energy losses of the pulsar, causing its migration above the death-line in the P −Ṗ diagram. Taking into account that PSR B0826-34 is located relatively close to the death-line (τ = 3 × 10 7 yr, B d = 1.4 × 10 12 G) and that it is very intermittent, staying in the ON state for only 30% of the time, it is likely that this pulsar is visible only when it oscillates, which may also determine the character of variation of the observed subpulse velocity. 9 5.2 PSR B0818-41

Basic parameters
The pulsar B0818-41, with the main parameters P = 0.545, τ = 4.57 × 10 8 yr, B d = 1.03 × 10 11 G, was discovered during the second pulsar survey (Manchester et al. 1978;Hobbs et al. 2004). The width of the average pulse profile of PSR B0818-41 is close to 180 • with a pronounced subpulse drift along a wide range of longitudes. The typical subpulse drift pattern (see Fig. 2 of Bhattacharyya et al. (2007) at the frequency of 325 MHz) consists of an inner region with slower apparent drift velocity, surrounded by an outer region with larger intensity of subpulses and steeper drift. Multifrequency observations of Gil (2009) at 157, 244, 325, 610 and1060 MHz show that at lower frequencies the subpulses become weaker and may be seen only in the outer regions at 244 MHz and only in the trailing outer region at 157 MHz. The observing geometry, i.e., the values of the inclination angle and of the impact angle, is not uniquely determined for this pulsar and the polarization profile admits several interpretations. However, based on the average polarization behaviour, Qiao et al. (1995) concluded that the inclination angle of PSR B0818-41 should be small. Unique nulling properties of PSR B0818-41 are studied in Bhattacharyya, Gupta & Gil (2010).

Subpulse phenomenology
The value of P3 = 18.3 ± 1.6P was found in Bhattacharyya et al. (2007), observing at the frequency 325 MHz, by means of the fluctuation spectrum analysis, and the same value was confirmed later in Bhattacharyya, Gupta & Gil (2009) for all other observing frequencies. In the inner part of the subpulse drift region, where several (typically 3 to 4) subpulse tracks are observed within one pulse, the value of P2 was found from the second peak of the autocorrelation function to be 17.5 • ± 1.3 • (Bhattacharyya et al. 2007). In the outer regions of the profile the value of P2 is larger (already from the visual inspection of the subpulse tracks) and not easily measurable, because typically no more than one subpulse per pulse is seen 9 The observed periodicity of hundreds of seconds may be reached by core g-modes of the neutron star (McDermott, van Horn & Hansen 1988). in these regions. Estimations for the different observing frequencies and different parts of the profile can be found in Bhattacharyya, Gupta & Gil (2009), while the average value of P2 may be taken around 28 • (Bhattacharyya et al. 2007).
As indicated in Gupta et al. (2004), the measured values of the periods P2 and P3 are not necessarily equal to the true intrinsic ones. The value of P3 may be affected by aliasing, which starts to play a role when P3 < 2P . The value of P2 is affected by the finite time required for the line of sight to traverse the polar cap as well as by the difference between the longitude l along the pulse (in which we measure P2) and the azimuthal coordinate φ around the magnetic axis. However, if we assume that there is no aliasing, the correction to the measured value of P2 due to the finite traverse time is given by a factor of 1/[1 + P2/(360 • P3)] [derived from the equation (5) of Gupta et al. (2004)], which in our case is ∼ 0.997. As long as we are not concerned with the structure of the carousel as a whole and we don't consider the possibility of aliasing in our calculations, we assume everywhere that the measured values of P2 and P3 are equal to the intrinsic ones.

Analysis of the subpulse drift
As in the previous subsection, we start the analysis from the determination ofr, corresponding to the measured subpulse drift velocity of PSR B0818-41. Choosing for this estimation χ = 0 and ξ = 0.5 and using the other known parameters of the pulsar, we get r0 = 1.011 for ωD = −0.956 • /P (this value ofr0 depends very weakly on the chosen ξ and χ, provided the inclination angle is small). For comparison, assuming that the temperature of the polar cap is T = 3 × 10 6 K, the altitude of the PFF h HA lies in the range 0.013 − 0.106 for the different values of fρ and h HM = 0.182, while for the temperature T = 4 × 10 6 K the corresponding values for the h HA are 0.0097 − 0.0796 with the same h HM . So, the estimated value for the altitude lies close to the lower boundary obtained through the PFF approach. Fig. 5 is devoted to the illustration of the geometry of PSR B0818-41, with the inclination angle between the magnetic and the rotational axes χ = 0.34 • and the impact angle β = 0.51 • . The considered geometry corresponds to an outer line of sight, which, within the SCLF model, naturally results in negative value of the subpulse drift velocity. The chosen inclination and impact angles are much smaller than those suggested before in the literature 10 . These small values are required by the fact that the size of the polar cap is very small Θ = 0.94 • at the considered altitude. However, one may notice that the G-2 geometry of Bhattacharyya, Gupta & Gil (2009), which gives the best fit to the polarization angle profile of PSR B0818-41, effectively corresponds to an outer line of sight with χ = 180 • − 175.4 • = 4.6 • and β = 6.9 • , so that sin β/ sin χ = 1.5. Based on this, we chose the inclination and impact angles of our geometry to satisfy β/χ = 1.5 in order to match the polarization profile of the pulsar. The polarization angle, calculated through the "Rotating Vector Model", is plotted in the lower panel of the Fig. 5 as a function of the longitude and reproduces well the observational data [cf., Fig. 6 of Bhattacharyya, Gupta & Gil (2009)].
We can use the geometry of Fig. 5 to model the observed behavior of the subpulse tracks. Starting from the trailing end of the  profile (because the drift velocity is negative) we evolve the azimuthal coordinate φ with time using the projected velocity (18). At each time step we calculate the corresponding pulse longitude using relation (16) and add a new subpulse track every P3 = 18.3P . The resulting pattern of the subpulse motion for 200 pulses as a function of pulse longitude is shown in Fig. 6 through blue solid lines, which may be directly compared to Fig. 2 of Bhattacharyya et al. (2007). The solid vertical lines indicate the region, where the subpulse drift is actually observed. We estimated its width as 120 • , corresponding to the pulse longitudes 140−260 • of Fig. 2 of Bhattacharyya et al. (2007). Note that the center of the profile there is slightly shifted with respect to 180 • . When compared to the observational data, we find that the SCLF model is able to reproduce the curved subpulse tracks reasonably well.
The pattern of the subpulse tracks obtained with our model is symmetric by construction. However, according to the observations of Bhattacharyya, Gupta & Gil (2009), the drift bands of the trailing outer region of the pulse appear to be steeper than the drift bands of the leading outer region. Although a clean explanation for this effect is still missing, we argue that some degree of asymmetry can be due to the effects of retardation, aberration and refraction of the signal in the outer magnetosphere (Gangadhara & Gupta 2001;Gupta & Gangadhara 2003;Petrova 2000;Weltevrede et al. 2003).
We emphasize that the original spark model of Gil & Sendyk (2000), on the basis of very general arguments, predicts that the pulsar polar cap should be densely filled by equidistant equalsize sparks. On the contrary, the carousel pattern proposed later in Gupta & Gil (2009) has larger and wider separated sparks in the outer ring in order to explain the observed subpulse behavior (compare Fig. 1 of Gil & Sendyk (2000) with Fig. 10 of Esamdin et al. (2005)). We find it interesting that, within the SCLF model, it is possible to explain the observed curved subpulse tracks by means of the velocity variations only, without breaking the assumption that the features responsible for the subpulses are equal in size and equidistant. This supports the argument, already presented in the previous subsection, according to which the variations of P2 with the pulse longitude, observed for many pulsars, may be completely explained by the variability of the subpulse drift velocity across the polar cap, while the value of P3, which seems to be independent of the pulse longitude and even of the observing frequency, should reflect an intrinsic characteristic property of the individual pulsar. One interesting observation made in Bhattacharyya, Gupta & Gil (2009) is that the leading and trailing outer regions of the pulse profile maintain a unique phase relationship, with the maximum of the energy in the trailing component being shifted in time by ∼ 9P with respect to the maximum of the energy in the leading component. Based on this observation, the authors propose an elegant solution to the aliasing problem, arguing that the considered shift may not be explained without aliasing and suggesting the model of 20 sparks ring with first-order alias and a true drift velocity of 19.05 • /P as the most plausible description of the system. As an alternative, we propose that the position of the picks of the pulse profile is not strictly determined by the position of the sparks in the outer ring, but modulated by the outer regions of the magnetosphere. One may assume, as it is customary, that the major radio emission mechanism of the pulsar is due the formation of the secondary plasma from the energetic photons emitted by the primary particles. Near the axis this process is negligible due to the large curvature radius of the magnetic field lines, while on the edges of the polar cap region the acceleration potential itself drops to zero. This produces the "hollow cone" distribution of the secondary plasma in the magnetosphere above the PFF (see Petrova (2000), Weltevrede et al. (2003), where this model is used for the study of magnetospheric refraction, also Fung, Khechinashvili & Kuijpers (2006)). We therefore believe that such a distribution modulates the emission profile and changes the position of the maxima with respect to the position of the outer ring of sparks. However, this issue is beyond the scope of the present paper and we leave it for a future study.

CONCLUSIONS
The phenomena of drifting subpulses are typically explained by resorting to the vacuum or to the partially screened gap models of pulsar magnetospheres. For example, the partially screened gap model allows for the formation of a spark carousel due to discharges of the large potential drop through the inner polar gap and these sparks are thought to be responsible for the appearance of subpulses. However, it has been recently shown by van Leeuwen & Timokhin (2012) that the expression used for the estimation of the subpulse drift velocities, both in the vacuum and in the partially screened gap model, is not accurate enough.
On the other hand, considering the pulsar magnetosphere as a global object, one can propose alternative mechanisms for the formation of distinct emitting features representing subpulses and in this paper we have reconsidered the ability of the Space-Charge Limited Flow model to explain this phenomenology. The SCLF model provides analytical solutions for the scalar potential in the polar cap region of the pulsar magnetosphere in case of free outflow of the charged particles from the surface of the star. Hence, the drift velocity of subpulses along the pulse can be interpreted in terms of the plasma drift velocity, which in turn depends on the gradient of the scalar potential rather than on its absolute value.
After considering a selected sample of sources taken from the catalog of Weltevrede, Edwards & Stappers (2006), we have found the following conclusions: • the SCLF model predicts the subpulse drift velocities compatible to the observed ones at heights above the surface of the star close to the pair formation front; • the angular dependence of the plasma drift velocity in the SCLF model provides a natural explanation for the variation of the subpulse separation P2 along the pulse. In particular it may explain the curved subpulse driftbands of PSR B0818-41 and the difference between the values of P2 in the main pulse and interpulse of PSR B0826-34.
These results suggest that the role of the SCLF model in explaining the drifting subpulse phenomena has been underestimated, calling for additional investigations and systematic comparisons with all available observations.