text:
# Understanding the radio emission geometry of multi-component radio pulsars from retardation and aberration effects

## 1 Introduction

Radio emission from a pulsar is believed to originate in the open field line region of the polar cap of the neutron star. The size, shape and location of regions of radio emission in the average profiles of pulsars is thus expected to reflect the arrangement of emission regions in the pulsar magnetosphere. Pulsar average profiles exhibit a great diversity in shape, and their classification based on the number of emission components is a useful starting point to study the emission characteristics of pulsars. Rankin (1993, 1990; and references therein) has carried out such a detailed classification, and has come up with the conclusion that there are two kinds of emission components – core and conal – in pulsar profiles, which result from two distinct types of emission mechanisms. Further, Rankin proposes that the conal components arise from two nested hollow cones of emission, which along with a central core emission region, make up the complete pulsar emission beam. The actual profile observed for a given pulsar depends on the cut that the observer’s line-of-sight makes through this emission beam. From the above work, Rankin also concludes that core radiation originates from very close to the neutron star surface whereas the conal radiation comes from regions higher up in the magnetosphere. The outer cone is postulated to originate higher up in the magnetosphere than the inner cone, but along the same set of field lines. Recent work by Gil et al (2002) which compares frequency dependence of pulsar radiation features for average profiles with those for single pulses, comes up with evidence in support of such a conal beam model. A somewhat different model of the pulsar emission beam is described by Lyne & Manchester (1988) who propose that the emission within the beam is patchy, i.e. the distribution of component locations within the beam is random, rather than organized in one or more hollow cones.

One of the problems in resolving this conflict is related to a correct determination of the total number of emission components for a pulsar and the significance of their arrangement within the pulse window. In a recent paper (Gangadhara & Gupta, 2001; hereafter paper I), we have addressed these issues with the help of a novel method of studying the emission geometry of pulsars. Using the sensitive “window-thresholding” technique (W-T technique hereafter) on single pulse data for a given pulsar, we first identify all possible conal emission components in the profile that can be detected given the sensitivity of the data, and determine accurate locations in pulse longitude for these. From this we determine the total number and arrangement of emission cones in the observed profile. Then, from detection of retardation and aberration effects in the locations of the conal components, we estimate the emission height for each cone. Coupling this with a dipolar field geometry and the available information about the relative orientations of the spin and rotation axes of the pulsar and the line of sight, we are able to identify the field lines from which each cone originates, thereby giving a complete solution for the emission geometry for each cone. Using this technique on single pulse data for PSR B0329+54 at two different radio frequencies, we showed in paper I that there are as many as 9 identifiable components in this pulsar, arranged in the form of 4 concentric rings around the central core component. The emission heights for the cones were found to increase systematically from inner to outer cones for a given frequency, ranging from $``$ 160 to $``$ 1150 km. Also, for a given cone, the height was found to reduce with increasing frequency, providing direct proof of the commonly accepted “radius to frequency mapping” paradigm.

In this paper, we extend our study of emission geometry to another six pulsars. Section 2 describes the observations and data analysis steps. In section 3, we present our results, where we are able to solve for the emission geometry for 5 of the pulsars studied, and have tentative results for the sixth. In section 4 we discuss the implications of our findings for the current understanding of pulsar emission physics.

## 2 Observations and Data analysis

For our study of emission geometry, we selected a set of pulsars with clearly defined properties that make it easy to apply the techniques used for PSR B0329+54 in paper I. These included the following criteria: (i) the presence of multiple emission components (or indications thereof) in the existing profiles; (ii) the presence of a reasonably well identified core component at our frequency of observation, which was the 325 MHz band of the GMRT; (iii) strong enough mean flux such that reasonable signal to noise ratio (SNR) could be expected for single pulse observations with the incoherent array mode of the GMRT. The first two constraints above confine the selected pulsars to the category of “triple” or “multiple” pulsars as defined by Rankin (1993), or equivalently, to the category of “cones with cores” pulsars, as defined by Lyne & Manchester (1988). Of the possible candidates, we avoided pulsars with very wide profiles (e.g. PSR B0826$``$34, PSR B1541+09) as it is not easy to unambiguously identify discrete emission components in these profiles due to the extended emission over a large fraction of the pulse period. The final list of pulsars selected for our study is given in Table 1, along with relevant details. Here, $`\alpha `$ and $`\beta `$ refer to the standard definitions of the angle between the rotation and magnetic axes and the impact angle of the line of sight with respect to the magnetic axis, respectively. The values for these angles are taken from Rankin (1993). These are not very different from the corresponding values given by Lyne & Manchester (1988) – the typical difference for $`\alpha `$ is $``$ a few degrees and that for $`\beta `$ is less than $``$ 0.5 to 1.0 degree.

The observations were carried out at the GMRT, in the 325 MHz band, using the incoherent array mode (see Gupta et al. (2000) for more details about the pulsar modes of operation of the GMRT). The data were obtained by incoherent addition of the dual polarization signals from 12 antennas. The bandwidth used was 16 MHz, divided into 256 spectral channels by the digital back-ends. The raw data were integrated to a time resolution of 0.516 milliseconds before being recorded for off-line analysis, where the data were dedispersed and gated to obtain the single pulse sequence for each pulsar. During this analysis, care was taken to check and flag the data for radio frequency interference signals. For a majority of the pulsars, the centre frequency of the observations was 318 MHz; for one case it was 328 MHz, and 333 MHz for another (see Table 1).

All of the 318 MHz observations were carried out on 10-11 Dec 2000. For most cases, the total number of pulses were obtained from a single observing run. For some pulsars (e.g. PSR B2045$``$16), data from different observing sessions (on different days) was combined (after proper alignment) to obtain the final single pulse sequence with enough pulses for subsequent analysis. The total number of pulses for each pulsar is given in column 7 of Table 1.

The average profile obtained for each pulsar is shown in Figure 1. The time resolution here is the basic sample time interval of 0.516 millisec, except for the case of PSR B2111+46, where the data has been integrated by a factor of 4, to a sample interval of 2.064 millisec. The SNR obtained for the peak of the average profiles is quite good – in the range 200 to 600 – resulting in single pulse peak SNRs of the order to 5 to 10 and better, which are needed for using the W-T technique effectively. All the pulsars show multi-component profiles, with at least 3 or more components easily visible. The classification of each pulsar, according to the scheme proposed by Rankin (1993), is given in col. 6 of Table 1, where the triple (T) pulsars have 3 known components and the multiple (M) pulsars have 5 known components. From published results related to (i) frequency evolution of the profile components and (ii) polarization properties of average profiles (e.g. Rankin 1993, Lyne & Manchester 1988 and references therein), it is clear for most of these pulsars that the central component is the core component. Further, in almost all cases, it is apparent from Figure 1 that the conal components are located asymmetrically with respect to the core component, with the tendency for each trailing conal component to be closer to the core than its leading counterpart. This is consistent with the expectations of retardation and aberration effects due to finite emission height of the conal beams (e.g. paper I).

To carry out a detailed modeling of the geometry of the emission cones for each pulsar, we need to identify all possible emission components in the profiles, that can be determined within the SNR limitations of our data. For this we have used the W-T technique described by us in paper I, where we employed this technique to identify as many as 9 emission components in PSR B0329+54. In this technique, we set a window in a desired pulse longitude range and employ an intensity threshold to select the single pulses that go to make an average profile, which we refer to as a “W-T profile”. We consider all those pulses which have emissions above the threshold within the window. As a result of this averaging of selected pulses, emission components within the window improve in SNR compared to other parts of the profile and are thus more easily detected in the W-T profile. As described in Paper I, we use a few different checks to ensure that the detected components are genuine. Usually, the window size and location are varied by a few bins while checking that the component shows up at the same location (this is not always possible to do for components that are located in close proximity to known, strong components). In addition, the selection threshold is varied over a reasonable range to check the persistence of the component, and to minimize picking up contributions from random noise. These checks are robust enough to avoid detection of spurious components by chance accumulations in the selected window. This is supported by the fact that the W-T technique does not detect any new components when applied to off-pulse regions of the profile.

Using the W-T technique, we have been able to detect new emission components in most of the pulsars studied here, as well as get accurate estimates for the locations of the existing components. The total number of new components detected for each pulsar is given in column 8 of Table 1. Some typical examples of detection of components (new as well as known ones) are illustrated in the sample of W-T profiles shown in Figure 2. Here, Figs. 2(a) & 2(e) show examples of clear and easy detections (usually of already known components); Figs. 2(b),2(d) & 2(g) illustrate detection of new, weak components and Figs. 2(c), 2(f) & 2(i) are cases of detections of new components located in close proximity to known, strong components. In some cases, we were able to increase the confidence of our detections by integrating the data to a larger time constant (2 to 4 times the original sampling rate) and repeating the W-T analysis, albeit with a loss in resolution of the location of the component peak. This was particularly useful in cases where the original pulse profile is relatively broad and well sampled (e.g. PSR B2111+46) and the SNR is marginal for the weaker emission components.

In some cases, we detect likely candidates for new emission components but can not confirm them with enough confidence from our data set. This is due to limitations of either (i) inadequate SNR, (ii) not enough pulses for the weak components to show up reliably often enough or (iii) very close proximity of the suspected component to a known component, especially when the latter is a strong component. In such cases, we label the component detection as “tentative”. We mention this specifically for each pulsar when we discuss the individual results.

To determine accurate locations in pulse longitude for each component of a pulsar, we take the corresponding W-T profile and fit a Gaussian curve to the selected component. The noise contribution to each data point is estimated from the off-pulse noise r.m.s. of the W-T profile. The peak of the best-fit Gaussian is taken as the best estimate of the location of the component and the 1-sigma error estimate on this peak is used to derive the error in the location of the component. For the sample detections shown in Figure 2, we also show the best fit Gaussian to the component detected within the W-T window. For most cases, we get quite good fits, with values for the reduced chi-square $``$ 1.0. The error estimates are generally found to be proportional to the SNR of the fitted component, which is as expected.

From our determination of the total number of the emission components, we resolve the emission geometry for each pulsar as consisting of a central core component and one or more cones of emission around it. From our best estimates of the locations of these components, we then solve for the emission heights of the cones, as well as for the transverse location of the emitting field lines on the polar cap, using the technique described in paper I. The polar cap location of the field lines is quantified by the ratio of the distance to the “foot” of the field line (from the magnetic pole) to the distance to the “foot” of the last open field line in the magnetosphere – the parameter $`s_\mathrm{L}^\mathrm{i}`$ of paper I (in this description, the foot of the field line is the point where the field line intersects the neutron star surface).

## 3 Results

The results from our analysis of the 6 pulsars are summarized in Tables 2 to 7. The table for each pulsar gives, for every emission cone that we detect (including the known ones), the following information: $`\varphi _\mathrm{l}^\mathrm{i}`$ & $`\varphi _\mathrm{t}^\mathrm{i}`$, the location (with respect to location of the core component) of the leading and trailing emission components that constitute the cone (columns 2 and 3); $`\eta ^\mathrm{i}`$, the inferred retardation plus aberration angle (column 4); $`\mathrm{\Gamma }^\mathrm{i}`$, the inferred half angle of the conal beam (column 5); $`r_{\mathrm{em}}^\mathrm{i}`$, the computed emission height for the cone (column 6); and $`s_\mathrm{L}^\mathrm{i}`$, the polar cap location of the associated field lines (column 7). The emission heights of the cones are given in kilometers, as well as fraction of the light cylinder radius ($`r_{\mathrm{LC}}`$) in column 6. The errors for the derived quantities (columns 4 to 7) are determined from appropriate propagation of the errors on the estimates of the emission component locations (including that for the core component), and do not include any effects of possible errors in the values of $`\alpha `$ and $`\beta `$.

We first discuss the results for individual pulsars and present a general summary of the results at the end of this section.

### 3.1 PSR B0450$``$18

This pulsar shows a 3 component profile for frequencies below about 1 GHz and is interpreted as a typical triple (core and one conal pair) pulsar (e.g. Rankin 1993). However, this classification is not completely unambiguous. Though the frequency evolution of the components (e.g. Gould, 1994) appears to indicate that the central component is a core, the polarization characteristics do not show clear, unambiguous evidence of core radiation for this component. At our frequency, the profile shows a fairly wide emission region around the central peak, especially on the trailing side (Fig. 1). The W-T analysis readily localises the leading and trailing components in the profile. However, in the central 10 to 15 degrees of pulse longitude we get unclear and conflicting results. The W-T analysis shows evidence for multiple emission peaks in this region (e.g. at $``$4.1$`^{}`$and 3$`^{}`$), but these are generally not found to be stable against variations of window widths and thresholds. Hence, it is difficult to unambiguously identify distinct emission components in the central region, including a clear identification of the reference core component. It is possible that what we are seeing in this region is a somewhat tangential cut through the edge of an inner emission cone, and that the core, if present at all, is a weak component in this region. This would be consistent with the fact that this pulsar has a relatively large value of the impact parameter $`\beta `$ – it is the largest amongst all the pulsars studied here (Table 1). The entry in column 8 of Table 1 for this pulsar reflects our uncertainty in the identification of new emission components.

Using a tentative identification of the core as the peak produced by setting a wide window in the central region, we get the locations of the main leading and trailing emission components with respect to the core as $``$8.79$`^{}`$and +6.76$`^{}`$, respectively. In addition, there is some evidence from our analysis for the presence of weak components in the leading and trailing wings of the profile, but these can not be confirmed with the quality of our data. Solving for a single emission cone for this pulsar, the emission height and polar cap location are obtained as 310 km (1.2% of $`r_{\mathrm{LC}}`$) and 0.56 respectively (see Table 2). We note that the results for this pulsar should be treated with caution, as the emission geometry is probably more complicated than what we present.

### 3.2 PSR B1237+25

The typical average profile for this well studied pulsar clearly shows 5 emission components with the central one known to be the core component. The centers of both cones are clearly offset to earlier longitudes than the core, as would be expected due to retardation and aberration. Our analysis reveals the presence of two more emission components, located between the core and the inner conal ring. Of these two, the component on the leading side (at $``$2.12$`^{}`$) is seen relatively easily in our analysis (Fig. 2b), whereas the one on the trailing side (at +1.56$`^{}`$) is somewhat harder to detect. This is because it is located in the relatively narrower region between the core and the trailing component of the inner cone. We therefore interpret this pulsar as having 3 emission cones around the core component, with each cone clearly showing retardation and aberration effects. Solving for the locations of the emission cones at 318 MHz gives heights of 180, 460 and 600 km (0.3%, 0.7% and 0.9% of $`r_{\mathrm{LC}}`$ respectively) and transverse locations on the polar cap of 0.33, 0.41 and 0.59 (see Table 3).

### 3.3 PSR B1821+05

The average profile for this pulsar shows a core dominated triple component emission beam at 318 MHz (Fig. 1). The conal components, located at $``$13.75$`^{}`$and +10.94$`^{}`$, clearly exhibit retardation-aberration effects. Even though the SNR is not very good, there are indications of extra emission components between the core and the known conal components. With some difficulty, our analysis picks up a component at $``$6.07$`^{}`$from the core and one at +8.2$`^{}`$on the trailing side (Figs. 2c & 2d). In addition, we have marginal detections of two more conal components, at $``$10.25$`^{}`$and +4.61$`^{}`$. These last two are detected with fairly low thresholds in the W-T technique and are, at best, tentative detections, requiring confirmation with higher sensitivity observations. Analysis of the data after further integration by a factor of 2 were used to confirm the presence of these components.

Our final result for this pulsar is a detection of 3 cones for which the emission heights at 318 MHz are obtained as 290, 410 and 570 km (0.8%, 1.2% and 1.6% of $`r_{\mathrm{LC}}`$ respectively), and polar cap locations are obtained as 0.43, 0.57 and 0.63 respectively (see Table 4).

### 3.4 PSR B1857$``$26

This is a well known multi-component (M-type) pulsar, where the typical profile at meter wavelengths readily shows 4 components and has clear indications of a fifth component located close to the leading edge of the central, core component (see Fig. 2e for our detection of this component). Our analysis easily gives the locations for these 5 known components. In addition, we see evidence for two more outer components, located near the leading and trailing edges of the profile, at $``$23.73$`^{}`$and +19.95$`^{}`$with respect to the core (see Fig. 2f for one of these). We label these two as tentative detections, as the threshold and number of contributing pulses are rather small (though the analysis of data integrated by a factor of 2 supports these detections). Further, there are hints of emission components in between the known conal components, but these can not be confirmed with the quality of our data.

We thus have 5 confirmed components (along with 2 tentative ones) for this pulsar, leading to an interpretation of core and 2 emission cones (and a tentative third outer cone). The emission heights (at 318 MHz) for these 3 cones are found to be 220, 480 and 690 km (0.8, 1.7 and 2.4 % of $`r_{\mathrm{LC}}`$) and polar cap locations of 0.61, 0.69 and 0.74 (Table 5).

### 3.5 PSR B2045$``$16

This well known pulsar shows a 3 component profile where the central core component is clearly asymmetrically located with respect to the centre of the conal ring. This is an indication of significant aberration and retardation effects for the emission cone. On using the W-T technique, we are able to detect 3 more emission components, located in the “saddle” regions between the core and the known conal components, at longitudes of $``$7.2$`^{}`$, $``$5.5$`^{}`$and +3.06$`^{}`$with respect to the core (see Fig. 2g for an example of one of these). This gives a total of 5 conal components, leaving us with the somewhat difficult situation of a probable undetected conal component, which needs to be located in the narrow valley between the core and the already known trailing conal component. There are some indications of an emission component at +4.0$`^{}`$, but we can not determine this unambiguously from our analysis, as it is too close to the known component at +4.59$`^{}`$. Lacking this, we have the problem of deciding which of the 3 new conal components is the odd one out. If we take the components at $``$5.5$`^{}`$and +3.06$`^{}`$as forming the second conal ring, the solutions for the two rings are: emission heights at 1240 and 2230 km (1.3 and 2.4 % of $`r_{\mathrm{LC}}`$) and polar cap locations of 0.28 and 0.32 – this is the solution listed in Table 6. Alternatively, taking the components at $``$7.2$`^{}`$and +3.06$`^{}`$as forming the second conal ring changes the emission height estimate for this cone from 1240 km to 2090 km and the polar cap location shifts from 0.28 to 0.22. We believe the first option is more likely to be the real situation.

### 3.6 PSR B2111+46

This pulsar shows a very wide but clearly triple profile with a dominant core component at meter wavelengths (Fig. 1). The location of the conal components is clearly asymmetric with respect to the core, showing retardation-aberration effects. For this pulsar two data sets, one taken at 318 MHz in December 2000 and another taken at 333 MHz in February 2002 were available. Since the data from February 2002 had better SNR and more number of pulses, all the results presented here are from the analysis of this data set. Analysis of the December 2000 data produced results which are similar to those reported here, and well within the errors of the measurements.

Since the profile is very wide, the data were averaged to 4 times the original time resolution of 0.516 millisec, to improve the SNR for the W-T analysis. Besides confirming the location of the main components, our W-T analysis picks up 2 additional conal components, located at $``$18.71$`^{}`$and +14.91$`^{}`$with respect to the core (see Figs. 2h & 2i). In addition, there are hints of another pair of components located at the leading and trailing edges of the profile (at approximately $``$48.3$`^{}`$and +38.7$`^{}`$with respect to the core), but we are unable to confirm them clearly. Hence we confine ourselves to 2 cones of emission around the central core component, for this pulsar. The estimates for the emission heights at 333 MHz are 1360 km and 2080 km (2.8 and 4.3 % of $`r_{\mathrm{LC}}`$), while the values for the polar cap locations are 0.22 and 0.31 (see Table 7).

### 3.7 Summary of results

Our main results can be summarized as follows:

(i) All the 6 pulsars show clear signatures of an asymmetry of the cones with respect to the core radiation, in that the centers of the cones occur at earlier pulse longitudes, with respect to the location of the core. For a given pulsar, the magnitude of the effect increases from inner to outer cones.

(ii) Interpreting the shift of the cone centers as the effect of retardation and aberration on the emitted conal beam leads to reasonable values for the emission heights. These range from about 200 km to 2200 km over all the 6 pulsars. In terms of the light cylinder distance, this range is from 0.3% to 4.3% of $`r_{\mathrm{LC}}`$. For a given pulsar, the heights are found to increase systematically from the inner to the outer cones (though in some cases the sizes of the error bars prevent this from being made a firm conclusion).

(iii) Combining the emission height estimates with a dipole field line geometry results in estimates of the transverse polar cap locations of the foot of the field lines that each cone is radiated from. These are found to vary from about 0.22 to 0.74 of the radial distance to the last open field line from the magnetic axis. Furthermore, this value is found to increase systematically from inner to outer cones for a given pulsar (though again, for some cases, the sizes of the errors reduce the significance of this conclusion).

## 4 Discussion

There are several interesting aspects of pulsar emission geometry that can be addressed by our study. Firstly, it should be clear from this work (and paper I) that it is not sufficient to study only the average profiles to understand the emission geometry. We have shown that a careful analysis of good quality single pulse data can provide valuable extra information about the emission geometry. It would be almost impossible to detect some of the emission components we have detected, using other techniques such as fitting multiple gaussians to the average profile. Here, the window-thresholding method used by us comes in as a handy tool to find components that are weak and/or emit intermittently, or those which are located close to other strong components. Using this, we find that most triple and multiple component pulsars can be shown to have 5 or 7 components. It is then likely that most pulsars, when observed with sufficient quality data, are going to show multiple emission cones around the central core component.

Further, a careful identification of the emission components appears, almost invariably, to support the picture of a core beam plus multiple, concentric conal beams of emission, in that the conal components turn out to be arranged in equal numbers on either side of the core. What is even more interesting is that in all cases, there is an asymmetry in the location of the core with respect to the centers of the cones, and that too always in the same sense of having the cone centers advanced with respect to the core component location. Though we have presented the detailed analysis for only a handful of pulsars, this trend is easily visible in a vast majority of pulsars that have clearly defined core and conal emission components. A check of the data from different sources in the literature (e.g. Gould (1994); Hankins & Rickett (1986); Weisberg et al. (1999); the EPN database) reveals several clear cases: PSRs B0105+68, B1700$``$32, B1737+13, B1826$``$17, B1916+14, B1946+35, B2002+31 & B2003$``$08 (detailed studies of several such pulsars have been initiated and will be reported upon shortly). Thus it would seem that this is a fairly universal trend among pulsars that show clear core and conal emission components, and our interpretation of this as being due to retardation-aberration effects in the magnetosphere, is a reasonable one.

The emission heights that we obtain for the different pulsars have very plausible values. For PSRs B0450$``$18, B1237+35, B1821+05 and B1857$``$26, we obtain emission heights in the range of 200 to 700 km, while for PSRs B2045$``$16 and B2111+46, the values range from 1200 to 2200 km. Typical estimates of emission heights reported in literature are a few 100 km for most pulsars (e.g. Rankin (1993), Kijak & Gil (1998)). Rankin (1993) estimates the emission heights for the cones from the measured half-power widths of the conal component pairs and making the assumption that the cones are radiated from near the last open field lines of the polar cap. This yields 1 GHz height estimates of about 130 and 220 km for the two cones. Our height estimates will naturally be larger than these as we find that the emission is located on field lines that are well inside of the last open field line: a cone of the same angular width radiated from inner field lines has to come from a greater height in the magnetosphere. Specifically, for PSRs B1237+25 and B1857$``$26, Rankin(1993) gives estimates of emission heights for the cones at a range of frequencies. Comparing with the values at the nearby frequency of 270 MHz (Table 3 of Rankin, 1993), we find that our estimates for the same cone at 318 MHz are more than twice as large. If, however, the estimates of Rankin (1993) are corrected for emitting field line with $`s_\mathrm{L}^\mathrm{i}0.6`$ (as postulated in Mitra & Rankin (2002)), then the results are in better agreement with ours. Our emission height estimates are closer to those obtained by Kijak & Gil (1997,1998), though still somewhat larger.

We note that the height estimates for PSRs B2045$``$16 and B2111+46 are somewhat unusually large. Correspondingly, these two pulsars also have rather small values for the polar cap locations of the radiating field lines, as compared to the lower emission height pulsars. The reason for these larger heights is not clear. As a fraction of the light cylinder radius distance, however, all these heights are still quite small, being $``$ 5% of $`r_{\mathrm{LC}}`$ at the largest. This also means that effects of magnetic field sweepback (see paper I) can safely be ignored for all these cases.

Strictly speaking, all our height estimates are heights of the cones relative to that of the core radiation, as the retardation-aberration of the cone centers is measured with the core component as reference. Since the core radiation is thought to originate very close to the neutron star surface (e.g. Rankin, 1993 & 1990), these are good estimates for the actual emission heights of the cones above the neutron star surface (the extra height of 10 km for the radius of the neutron star is much smaller than our typical error bars). However, if this assumption about the height of the core emission beam is not correct, then the height of the core emission has to be added to the height estimates we have reported for the cones – in this sense, our results can be taken as lower limits to the actual emission heights. The other effect of a finite core height will be to cause the inferred polar location of the cones to become somewhat smaller (see Eqn. 15 of paper I), with the fractional change diminishing progressively from inner to outer cones. For example, an emission height of 100 km for the core radiation of PSR B1237+25 would cause the polar cap locations of the cones to reduce to 0.26, 0.37 and 0.55 from their present values of 0.33, 0.41 and 0.59. In any case, it is quite clear from our results that the core emission originates at heights that are substantially ($``$ 100s of km) less than that for the cones – otherwise we would not see any relative retardation-aberration effects between the core and cones.

Another interesting feature of the results is that they are quite self-consistent with radio emission originating within the open field line region of the magnetosphere, in that we do not obtain results for polar cap locations that are greater than 1.0. This can happen in principle, as the emission height estimate is completely independent of the width of the cone, whereas both these parameters are relevant for the calculation of the polar cap location (see Eqns. 9 and 15 of paper I). Polar cap locations greater than 1.0 can result if the emission height estimates are less than the current results by factors of 4 to 16. In a sense then, this argues that our interpretation of component location asymmetries as being due to retardation-aberration effects is a sensible interpretation.

Perhaps the most interesting result from this work is that the different cones for a pulsar originate on different sets of field lines in the magnetosphere. Clearly, all the conal radiation does not originate near the last open field line region. In fact, it appears that a significant fraction of the polar cap region may be active in the generation of the multiple conal beams. Furthermore, our results indicate that this activity is confined along distinct and concentric rings on the polar cap. Generally, it is thought (e.g. Ruderman & Sutherland (1975), Cheng & Ruderman (1980)) that this activity may be produced by spark discharges in vacuum gaps formed immediately above the neutron star surface. For such a model, our results indicate that the sparks are confined to concentric rings on the polar cap. We note that there are recent reports of inference of such patterns of circulating sparks for pulsars with conal-single geometries showing drifting sub-pulses (e.g. Deshpande & Rankin (1999)). Our results tend to support a similar conclusion for multi-component pulsars with more central traverses of the observing line-of-sight. For our pulsars, we can actually constrain the radii of these rings, via the $`s_\mathrm{L}^\mathrm{i}`$ parameter. For a 10 km radius neutron star with a period of P seconds, the radial size of the polar cap is just $`145\mathrm{P}^{0.5}`$ meters. For three concentric conal rings anchored to field lines at $`s_\mathrm{L}^\mathrm{i}`$ = 0.33, 0.41 and 0.59 (as in PSR B1237+25), this translates to ring radii of 41, 51 and 73 meters. Similar numbers ($``$ a few to several tens of meters) are obtained for other pulsars also. It is indeed remarkable that we can infer the presence of features of this size on the polar cap.

Gil & Sendyk (2000) have recently proposed a model for pulsar radiation that invokes concentric rings of sparks to explain the observed properties of pulsar radiation. In their theory, the separation between concentric rings of sparks, very close to the neutron star surface, is of the order of h, the height of the vacuum gap. They estimate the value of h to be of the order of 10 meters for pulsars with period $``$ 1 second, which is quite similar to the typical separation of rings that we infer. For example, for PSR B1237+25, the separation between first and second rings is 10 meters and that between the second and third is 22 meters. For PSR B1821+05, these numbers are 23 and 12 meters respectively; for PSR B1857$``$26, they are 15 and 9 meters. Thus our results make interesting connections with some of the existing theoretical models and support the idea of concentric rings of sparks circulating on the polar cap as being the cause of the radio radiation.

Although the activity that is the ultimate source of the radio radiation that we observe may be occurring very close to the neutron star surface, the final radiation that we see appears to originate at significant heights in the magnetosphere. Further, we have shown that the inner cones are emitted at lower altitudes and outer cones at progressively higher altitudes. In addition, it is generally accepted that the same cone is emitted at higher heights at lower frequencies. Clearly, all this is pointing to the requirement of the appropriate combination of parameters at specific points in the magnetosphere for the plasma disturbance to produce a beam of radiation that can be seen along our line of sight through the magnetosphere. This could be a combination of physical and geometrical parameters that contrive to give the right conditions at specific points in the magnetosphere. The localisation of the emitting regions achieved by the kind of work reported here and in paper I should help in better constraining the parameter space and thus help the theories trying to find the correct explanation for the generation of the radio waves.

## 5 Conclusions

Using the window-thresholding technique to analyse single pulse data from observations in the 325 MHz band of the GMRT, we have conducted a detailed analysis of the emission geometry of 6 radio pulsars that have prominent, multiple-component profiles. We have been able to detect new emission components in almost all of the pulsars. In some cases, we report tentative detections of weak components that need further confirmation studies. It thus looks likely that many multiple component pulsars may have at least 5 to 7 components. The total number and location of the emission components we detect clearly supports a picture of multiple conal beams around a central core beam for these pulsars (and are a pointer to the possibility that most pulsars may have such an emission beam structure).

Further, for all cases of conal emission, we find that the cone centre is offset to earlier longitudes with respect to the location of the core, and this effect increases in magnitude from the innermost to the outermost cones for every pulsar. It appears this effect may be present in most pulsars that have clear core and conal emission geometries. We interpret this effect as being due to retardation and aberration of the conal emission beams in the magnetosphere of the pulsar. From this, we are able to solve for the emission geometry and estimate the emission height for each cone and also the polar cap location of the field lines associated with the cone, following the technique developed in Gangadhara & Gupta (2001). We find emission height estimates in the range of 200 to 700 km for four pulsars – PSRs B0450$``$18, B1237+25, B1821+05 and B1857$``$26 – and 1200 to 2200 km for two others – PSRs B2045$``$16 and B2111+46. In terms of the light cylinder distance, the emission altitudes range from 0.3% to 4.3%. The estimates for the polar cap location of the radiating field lines lines come out to be in the range 0.22 to 0.74 of the distance to the last open field line, from the magnetic axis. Further, for every pulsar with multiple cones, the emission altitude and the polar cap location are found to increase systematically from the innermost to the outermost cone.

Mapping the locations of the field lines associated with each cone, on to the neutron surface, gives concentric rings with radii of the order of a few to several tens of meters and spacings between adjacent rings of about 10 to 20 meters. Our results thus tend to support a model of concentric rings of sparks produced in the vacuum gap region just above the neutron surface (e.g. Gil & Sendyk (2000)).

The fairly precise localization of emission regions that we have reported should provide additional information for constraining the parameter space for theories and models of radio emission mechanisms.

We thank P. Ahmadi and the staff of the GMRT for help with the observations. The GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. We wish to thank the referee for useful suggestions that helped improve the quality of the results, and J. Gil and T. Hankins for comments on the manuscript. In addition, YG would like to acknowledge the support of the Arecibo Observatory of the NAIC for a sabbatical visit, during which the bulk of the work on this manuscript was done.


format: text

text:
# THE FRACTAL DIMENSION OF AN OIL SPRAY

## 1 Introduction

The study of the geometric complexity in the break-up region of the liquid-gas interface in sprays represents a topic of special interest due to the fact that the physical process of atomization111We use here the common, albeit incorrect, sense for the word. According to its etymology, atom is an indivisible entity, and thus, to atomize really means to make indivisible, a meaning which is exactly the opposite to the intention when commonly used. is intimately related to many commercial and technological applications . The use of aerosol visualization analysis includes fields as diverse as the application of pesticides and insecticides to crops, hand-held sprays for painting or drying, fuel injection in combustion engines, liquid fuel spraying in the power generation industry, spray coating of pharmaceutical products, nuclear core cooling, and commercial filling processes, to mention just a few , and the studies are carried out with the efficient application of liquids, fuels or oils, in mind. In our study, the emphasis is set on the details of the relation between some properties of the fluid, like homogeneity or terminal speed, and the distance to the spray nozzle.

The use of shadowgraph images as a visualization tool is not new; their use in complexity studies however, has not been thoroughly exploited. Since our images are obtained with a high-resolution shadowgraph technique, they have proven specially pertinent for the implementation of the detailed analysis presented elsewhere for experiments in transport by fluids .

## 2 Experimental Setup

Among the different optical techniques that are used for the study of particles, liquids or gases in motion, the shadowgraphy stands out as an inexpensive and powerful technique . Its main advantage is to highlight the relative refractive index of bodies, and thus to allow the visualization of objects that would be impossible to attain with the conventional photographic techniques. This is particularly useful when the objects to photograph are transparent instead of opaque, i. e.: air and water, or transparent fluids with different refractive indexes.

In the shadowgraph technique the object is illuminated from its back through an optical system designed specifically for such purpose. In this way, the light that strikes the photographic film is the one that passes through, without interacting with the object; the light that interacts with the object is absorbed or deviated aside due to a difference in its refractive index and does not influence the film. Figure 1 shows a schematic view of the shadowgraph system used for this work.

The atomization process occurs as a result of the interaction between the oil in its liquid state and the surrounding air, and involves several stages through which the oil becomes an aerosol. These stages can be clearly observed in any of the pictures. Initially, as soon as the oil is forced to leave the atomizer nozzle with a velocity that has both axial and radial components, the fluid is turned into thin laminar waves that are gradually converted by the aerodynamic forces into thin ligaments. Downstream, these ligaments are broken up into a cloud of small droplets that continue to move with a terminal average velocity of several meters per second.

Our photographs correspond to the spray of transparent oil produced by a pressure swirl atomizer: the nozzle orifice is $`1`$ mm in diameter and produces a hollow cone spray with an amplitude of $`45^{}`$ (Fig. 2). The atomization pressure, obtained with a common gear pump, is $`4\times 10^5`$ Pa. The fluid is a blend of industrial oils designed to match the characteristics of diesel fuel which produces an impressive and complex droplets formation spray .

The measurement of the terminal speed was carried out with a technique preceding the Particle Image Velocimeter (PIV) based on a double photographic exposure. The break-up region is optically magnified (three times in our case) and the ISO $`400`$ b&w professional film is imprinted by a double flash from an Argon jet stabilized spark gap which lasts for $`300`$ ns and is shot again after $`0.1`$ ms; this creates a double image whose characteristic and identifiable features are separated by a certain distance that, together with the time lapse between flashes and the optical magnification, allow for the determination of the droplets speed. The mean speed thus measured at the break-up region, $`30`$ mm from the nozzle, has a value of $`15.5`$ m/s. Although no magnification exists in the photographs used for this work, the type of film employed allows for a photographic magnification of up to twenty times without any sharpness loss. Finally, in order to analyze their structure, the photographs were digitalized with a resolution of $`150`$ dpi in an 8-bit gray-scale comprising $`254`$ tonalities of gray.

## 3 Measurements

We have used the well-know box-counting method to calculate the fractal dimension of the spray contour: a grid of size $`\epsilon `$ is superimposed on the black and white images and then, the boxes of this grid that intersect any part of the image are counted to provide the number $`N`$; the process is repeated reducing the size of the grid and thus, the discrete function $`N(\epsilon )`$ is generated. The initial size is usually a one pixel grid  and the number of pixels of size $`1/\epsilon `$ is increased until the grid size reaches a certain previously determined value (this value is usually given by the smallest possible pixel size). Finally, one is interested in the limit

$$D=\underset{\epsilon 0}{lim}\frac{\mathrm{ln}N(\epsilon )}{\mathrm{ln}(1/\epsilon )}$$

which defines the fractal dimension $`D.`$ For the analysis of the contour of the oil aerosol, examples shown in Figs. 3 and 4, we varied the size of the grid from $`1`$ to $`1,300`$ pixels (the longitude of the image in pixels) and 25 different images of the spray were taken under the same physical conditions; within the attainable accuracy, the corresponding results produce all the same value.

The box-counting method is designed for monochromatic images  and thus one has to be very careful when analyzing either color or varying intensity images. This is due to the fact that the features from the image whose intensity is below a certain value, called the threshold level, will be lost when transforming the image to a black and white one . In our case, this problem is represented by some parts of the image being less opaque than others due to light extinction; some thin filaments are less opaque than their thicker similes. Even though every precaution was taken to adjust the depth of field to the size of the spray, some small contribution to this intensity problem may arise from a lack of perfect focusing.mTo avoid complications coming from this difficulty, our analysis was carried out by filtering each image at every gray tonality from the $`0^{th}`$ to the $`254^{th}`$ level of the digitalized images, i. e., the first filter allows only for the black parts to be taken into account, and then the intensity is diminished one digital step to incorporate gradually the less opaque features of every image. If this filtering process is not carried out, the simple transformation of the images to black and white ones produces the equivalent of filtering every image at half its maximum intensity value, Fig. 6. However, as shown in the following section, the filtering process just described reveals the rich fractal structure of the oil spray over a much wider range of intensity values or gray tonalities.

## 4 Results

The fractal dimension of the contour of the oil spray when the images are only filtered at half their maximum intensity is $`1.57\pm 0.03`$ with a linear correlation coefficient of $`0.98`$; the result is illustrated in Fig. 6. The whole spectrum of the fractal dimension of the contour however, is shown in Fig. 7, where the diameter of the open circles indicates the uncertainty at each intensity value and corresponds to the standard deviation of the data in all the $`25`$ analyzed photographs. To illustrate the fractal structure and the self-similarity of the contour at different intensities, the graph in Fig. 7 includes the fractal dimension value obtained from the filtering process starting at the level in which the image becomes a single black rectangle up to the level of a single white one. The graph displays the fractal structure for the whole digital intensity range; a small star indicates the already mentioned fractal dimension at half-intensity filtering ($`126th`$ digital level). A linear fit with a correlation coefficient of $`0.997`$ is also shown, the slope and ordinate are: $`(9.969\pm 0.154)\times 10^4`$ and $`1.436\pm 0.002`$, respectively. After the $`249^{th}`$ digital level, the filtering leaves no tractable structure in the image; the gradual decay of the fractal dimension observed in Fig. 7, from the $`249^{th}`$ to the $`254^{th}`$ digital level of gray, is therefore due to a tenuous gray intensity present in the background of our photographs.

## 5 Conclusions

The shadowgraph technique allows to obtain high-definition images of the whole oil spray which are ideal for the analysis of the dynamics of the aerosol and even of its smallest details since the image can be magnified several tenth times without noticeable sharpness losses; a higher density of dots is equivalent to a smaller pixel size which enables the analysis of smaller regions in the image.

For complexity studies of similar systems through fractal dimension analysis , it may be possible to use b&w photographs since the strong linear relation between intensity and threshold assures the existence of an interval (centered at the b&w intensity value) where the fractal dimension will remain linear in relation to the intensity. However, the analysis carried out for the contour intensity clearly indicates that, in this and similar cases, it is not suitable to directly digitalize the images in a purely black and white format since then all the rich structural information contained in the gray tonalities is irremediably lost.

This work has been partially supported by the IIE (México), and DGAPA-UNAM (project IN101100); AACP acknowledges financial support from DGEP-UNAM


text:
# Stable two-dimensional solitary pulses in linearly coupled dissipative Kadomtsev-Petviashvili equations

## 1 Introduction and derivation of the model

Localized structures, such as solitary pulses (SPs), play a dominant role in many conservative and dissipative nonlinear physical systems. As is commonly known, in conservative systems SPs are supported by balance between nonlinearity and dispersion , while in dissipative models, such as the Ginzburg-Landau equations, it must be supplemented by the balance between losses and gain .

An important example of a one-dimensional (1D) model that combines conservative and dissipative effects is a mixed Kuramoto - Sivashinsky (KS) – Korteweg - de Vries (KdV) equation, which was first introduced by Benney  and is therefore also called the Benney equation. This equation finds various applications in plasma physics, hydrodynamics and other fields . SPs are, obviously, important objects in systems of this type ; however, they cannot be completely stable objects in the Benney equation, as the zero solution, which is a background on top of which SPs are to be found, is linearly unstable in this equation due to the presence of the linear gain (however, if the dispersion part of the Benney equation is large enough, the growing perturbation, moving at its group velocity, does not actually overlap with the SP and therefore does not destroy it, see Refs.  and references therein). A stabilized version of the Benney equation was recently proposed in Ref. . It is based on the KS – KdV equation for a real wave field $`u(x,t)`$, which is linearly coupled to an additional linear equation of the diffusion type for an extra real field $`v(x,t)`$, that provides for the stabilization of the zero background:

$`u_t+uu_x+u_{xxx}v_x`$ $`=`$ $`\alpha u_{xx}\gamma u_{xxxx},`$ (1)
$`v_t+cv_xu_x`$ $`=`$ $`\mathrm{\Gamma }v_{xx}.`$ (2)
Here, $`\alpha `$, $`\gamma `$ and $`\mathrm{\Gamma }`$ are coefficients accounting for the gain and loss in the $`u`$–subsystem and loss in the $`v`$–subsystem, respectively, and $`c`$ is a group-velocity mismatch between the fields.

It was shown both analytically (by means of the perturbation theory) and numerically in Ref.  that the system (1), (2) gives rise to a completely stable SP, as well as to stable bound states of SPs, in a broad parametric region. As a matter of fact, Eqs. (1) and (2) furnish the first example of a model of the KS type that gives rise to fully stable pulses, and they can be easily observed in experiment (see a detailed discussion of physical realizations of the model – first of all, in terms of liquid films flowing down under the action of gravity – in the next section).

The liquid films and other systems in which the observation of stable SPs is expected (e.g., double-front flames, see below) are two-dimensional (2D) media, therefore a relevant issue is to introduce a physically meaningful 2D version of equations (1) and (2) and seek for stable 2D pulses in the generalized model. This is the objective of the present work. It will be demonstrated that the 2D model which will be derived here is a generic one for a number of different applications. The results will directly point at a new type of 2D pulses that can be observed experimentally in a straightforward way. Besides that, a possibility of the existence of stable 2D pulses may help to understand the phenomenon of turbulent spots, see, e.g., Refs.  and references therein.

The paper is organized as follows. The 2D model is derived in detail in section 2, starting with a particular physical problem, viz., a downflowing liquid film carrying a surfactant, and then proceeding to the generic character of the model. Parallel to the derived model, we will also consider its counterpart which differs by the form of 2D dissipative terms, in order to demonstrate that basic results are insensitive to particular features of the model (which is relevant to show, even if the model is generic). In section 3, we consider stability of the zero solution in the 2D model, which, as well in the 1D case, is a necessary condition for the full stability of SPs. In section 4, an analytical perturbation theory for SPs is developed by treating the gain and loss constants as small parameters. To this end, first a family of exact 2D soliton solutions of the zero-order system (the one without the gain and loss terms) is obtained, following the pattern of the well-known “lump” solitons of the KP-I equation. Then, using the balance equation for the net field momentum, similar to how it was done in the 1D model , we demonstrate that the combination of the gain and loss terms may select two (or none) stationary pulses out of the continuous soliton family existing in the zero-order system; in the case when two stationary pulses are found, it is very plausible that the one with the larger value of its amplitude is stable. In section 5, we present results of direct numerical simulations of the full 2D model, which completely confirm the analytical predictions, i.e., the existence of stable 2D localized SPs. The paper is concluded by section 6.

## 2 The model

The physical meaning of the model based on Eqs. (1) and (2), and its 2D generalization developed below, can be understood in terms of a particular application to a thin downflowing liquid layer with a surfactant trapped on its surface. As is well established (see a review ), in the 1D case the evolution of the flow velocity field $`v(x,t)`$ in the layer is governed, in appropriately chosen units, by the KS–KdV (Benney) equation, $`u_t+uu_x+u_{xxx}=\alpha u_{xx}\gamma u_{xxxx}`$. In this equation, the gain $`\alpha `$ is induced by gravity, and the loss parameter $`\gamma `$ is proportional to the fluid’s viscosity coefficient, while the left-hand side is generated by the Euler’s equations (for an irrotational flow), exactly the same way as in the classical derivation of the KdV equation. If the surfactant is distributed on the surface of the layer with a density $`c+v(x,t)`$, where $`c`$ and $`v`$ are, respectively, its constant and small variable parts ($`\left|v\right|c`$), the gradient of $`v`$ creates, through the variation of the surface tension, an additional force $`v_x`$ which drives the flow, hence the Euler’s equation adds the term $`v_x`$ to the right-hand side of the KS–KdV equation (the coefficient in front of this term may be scaled to be $`1`$), so that the equation takes precisely the form (1). Further, the evolution of the surfactant density is governed by an obvious advection-diffusion equation: $`v_t+\left[u\left(c+v\right)\right]_x=\mathrm{\Gamma }v_{xx}`$, where $`\mathrm{\Gamma }`$ is the surface diffusion constant. With regard to the condition $`\left|v\right|c`$, the latter equation takes the form (2).

In the 2D case, we consider a quasi-1D (weakly 2D) flow of the film, with the $`y`$–scale much larger than that along the $`x`$ axis. In other words, if the wave is taken as $`\mathrm{exp}\left(iKx+iQy\right)`$, the small wavenumbers are ordered so that

$$QK^2.$$
(3)
Then, according to the classical derivation , the KdV part of Eq. (1) is replaced either by the KP-I equation,

$$\left(u_t+uu_x+u_{xxx}\right)_x=u_{yy},$$
(4)
or by the KP-II equation, which is

$$\left(u_t+uu_x+u_{xxx}\right)_x=u_{yy},$$
(5)
the coordinate $`y`$ being properly rescaled. The choice between Eqs. (4) and (5) is determined by the sign of the 2D correction to the dispersion; in particular, the capillarity gives rise to KP-I. The difference between the KP-I and KP-II equations is that, although both of them have quasi-1D (i.e., $`y`$-independent) soliton solutions that reduce to the usual KdV solitons, only in the KP-II equation this soliton is stable against $`y`$-dependent perturbations. On the other hand, the KP-I equation has stable 2D solitons (the so-called “lumps”), which are weakly (non-exponentially) localized in both $`x`$ and $`y`$, see below; the KP-II equation does not have 2D solitons.

The next step is to accordingly generalize, for the 2D situation, the dissipative term in Eq. (1). Dissipative generalizations of the KP equations were introduced in some works, see, e.g., Refs. . Those generalizations follow the pattern of the arrangement of the KP equations proper: if one starts from a corresponding 1D equation containing dissipative terms \[for instance, Eq. (1)\], which is written as $`\mathrm{something}=0`$, its 2D counterpart is $`\left(\mathrm{something}\right)_x=\pm u_{yy}`$. The accordingly modified Eqs. (1) and (2) then take the form

$`\left(u_t+uu_x+u_{xxx}v_x+\alpha u_{xx}+\gamma u_{xxxx}\right)_x`$ $`=`$ $`\pm u_{yy},`$ (6)
$`\left(v_t+cv_xu_x\mathrm{\Gamma }v_{xx}\right)_x`$ $`=`$ $`\pm v_{yy}.`$ (7)
The second equation in this system can be simplified, as the usual ordering of the partial derivatives, adopted in course of the quasi-1D derivation , implies that $`v_{yy}`$ is a small quantity of a higher order than $`v_{xxx}`$, hence $`v_{yy}`$ may be dropped, so that Eq. (7) remains one-dimensional. As a result, this version of the 2D system takes the form

$`\left(u_t+uu_x+u_{xxx}v_x+\alpha u_{xx}+\gamma u_{xxxx}\right)_x`$ $`=`$ $`\pm u_{yy},`$ (8)
$`v_t+cv_xu_x\mathrm{\Gamma }v_{xx}`$ $`=`$ $`0.`$ (9)
However, while keeping Eq. (9) in the 1D form is quite acceptable, the way the dissipative and gain terms in Eq. (8) were made two-dimensional was formal, being not based on any physical argument. Moreover, it will be shown in the next section that, unlike its 1D counterpart (1) and (2), this 2D model cannot produce any stable solitary pulses, as its zero solution is always unstable.

In order to derive a physically relevant form of the dissipative part of the $`u`$–equation, one should resort to the standard procedure that derives a *generic* set of dissipative terms in the quasi-1D situation (in the context of convective flows) in the Newell-Whitehead-Segel (NWS) equation . This equation gives a simple prescription, which, as well as the KP equations, is based on the ordering (3) of the longitudinal and transverse wavenumbers: the longitudinal dissipative term $`\gamma u_{xxxx}`$ must be supplemented by its transverse counterpart $`\gamma u_{yy}`$ \[the scaling of the transverse coordinate $`y`$ in the NWS equation is precisely the same which cast the KP equation in the standard form (4) or (5); the identity of the scalings is not accidental, being a consequence of the fact that both KP and NWS equations are generic ones in the quasi-1D geometry, the former one in the class of dispersive equations, and the latter among dissipative equations\]. Thus, the proper form of the 2D system is

$`\left(u_t+uu_x+u_{xxx}v_x+\alpha u_{xx}+\gamma u_{xxxx}\gamma u_{yy}\right)_x`$ $`=`$ $`\pm u_{yy},`$ (10)
$`v_t+cv_xu_x\mathrm{\Gamma }v_{xx}`$ $`=`$ $`0.`$ (11)
In fact, exactly the same combination of dissipative terms as in Eq. (10) has been earlier derived in asymptotic equations governing nonlinear waves on thin downflowing liquid films  and in two-fluid flows .

Note that Eqs. (10) and (11) does not contain any additional free parameter in comparison with the 1D system (1), (2). As a matter of fact, this is another consequence of the fact that both the dispersive and dissipative parts of the system were derived for the generic quasi-1D case.

A special form of a 2D quasi-isotropic (rather than quasi-1D) generalization of the Benney equation was also derived, which, for instance, describes Rossby waves in a rotating atmosphere  (see also Refs. ):

$$u_t+uu_x+\mathrm{\Delta }u_x+\alpha u_{xx}+\gamma \mathrm{\Delta }^2u=0,$$
(12)
where $`\mathrm{\Delta }_x^2+_y^2`$, hence $`\mathrm{\Delta }^2`$ in Eq. (12) is a fourth-order isotropic dissipative operator. Following the pattern of Eq. (12), a quasi-isotropic generalization of Eqs. (1) and (2) may be introduced, replacing the term $`u_{xxxx}`$ in Eq. (8) by $`\mathrm{\Delta }^2`$, which leads to a system

$`\left(u_t+uu_x+u_{xxx}v_x+\alpha u_{xx}+\gamma \mathrm{\Delta }^2u\right)_x`$ $`=`$ $`\pm u_{yy},`$ (13)
$`v_t+cv_xu_x\mathrm{\Gamma }v_{xx}`$ $`=`$ $`0.`$ (14)
Note that all the 2D systems introduced above conserve two masses,

$$M=_{\mathrm{}}^+\mathrm{}_{\mathrm{}}^+\mathrm{}u(x,y)𝑑x𝑑y,N=_{\mathrm{}}^+\mathrm{}_{\mathrm{}}^+\mathrm{}v(x,y)𝑑x𝑑y$$
(15)
(which indeed have the meaning of masses in the application to liquid-film flows, in the cases of both the single layer with surfactant and two layers with the lower one dominated by viscosity).

Thus, Eqs. (10) and (11) can be derived in a consistent way as a system describing the downflow of a liquid viscous film carrying a surfactant. Actually, the derivation outlined above clearly suggests that this model is a generic one for weakly 2D systems combining dispersion, gain, and viscosity. In particular, a derivation involving more technicalities and following the lines of Ref.  shows that the same model applies to a downflow of a two-layered liquid film in the case when the substrate layer is dominated by viscosity. A physically different example may be a double-front flame propagating in a combustible gaseous mixture, in the well-studied case when the combustion involves two consecutive reactions (see Refs.  and references therein). In the case when a single-flame front is unstable, it is well known that its evolution is governed by the KS equation proper . It is known too that a situation with one front stable and one unstable is possible, which is described by a linearly coupled system, consisting of a KS equation and the one tantamount to Eq. (13). Dispersion, which is missing in the KS equation proper, can be induced by a background shear flow tangent to the flame , but detailed derivation of the full model for this case is beyond the scope of this work.

As concerns the model (13), (14), we consider it as a semi-phenomenological one, which may apply to cases which are “more isotropic” than those obeying the condition (3). We will study this model parallel to the physical one (10), (11) in order to see if qualitative results are sensitive to details of a given model. Accordingly, the systems (10), (11) and (13), (14) will be below referred to as physical and phenomenological ones, respectively. In the analysis presented in the following sections, we focus on the case of the KP-I type, i.e., with the upper sign in Eqs. (8), (10) and (13), as only in this case one may expect the existence of nontrivial 2D pulses, while the models of the KP-II type may only extend the SP found in Ref.  into a quasi-1D ($`y`$–independent) pulse in 2D.

## 3 The stability of the zero solution

As previously mentioned, completely stable SPs can only exist in a system whose trivial solution, $`u=v=0`$, is stable, therefore our first objective is to analyze this necessary stability condition. We substitute into the corresponding linearized equations a 2D perturbation in the form $`u\mathrm{exp}\left(ikx+iqy+\lambda t\right)`$, $`v\mathrm{exp}\left(ikx+iqy+\lambda t\right)`$, where $`k`$ and $`q`$ are arbitrary real wave numbers of the perturbation, and $`\lambda `$ is the corresponding instability growth rate \[note that, as all the equations that we are going to consider are written in the scaled form, $`k`$ and $`q`$ are not assumed to be specially small, unlike $`K`$ and $`Q`$ in Eq. (3)\]. This leads to a linearized dispersion equation

$$\left[k\left(\lambda ik^3\alpha k^2+\gamma k^4\right)iq^2\right]\left(\lambda +ick+\mathrm{\Gamma }k^2\right)+k^3=0$$
(16)
for the formal model (8), (9), or

$$\left[k\left(\lambda ik^3\alpha k^2+\gamma (k^2+q^2)^2\right)iq^2\right]\left(\lambda +ick+\mathrm{\Gamma }k^2\right)+k^3=0$$
(17)
for the phenomenological model (13), (14), or

$$\left[k\left(\lambda ik^3\alpha k^2+\gamma (k^4+q^2)\right)iq^2\right]\left(\lambda +ick+\mathrm{\Gamma }k^2\right)+k^3=0$$
(18)
for the physical model (10), (11). The stability condition states that both solutions of the quadratic equations (16), (17), or (18) must satisfy the inequality

$$\mathrm{Re}\left[\lambda (k,q)\right]0$$
(19)
at all the real values of $`k`$ and $`q`$.

As it was already mentioned, the zero solution in the formal model (8), (9) is always unstable, which can be shown as follows: in the case when $`k`$ is small, while $`q`$ is $`O(1)`$, then two roots of Eq. (16) can be expanded as

$$\lambda _1(k)=ick\mathrm{\Gamma }k^2+\mathrm{},\lambda _2=iq^2/k+\alpha k^2+\mathrm{}.$$
(20)
Obviously, the second root in Eq. (20) yields instability.

The zero solution may be stable in the physical and phenomenological models. Although direct check of the condition (19) for Eqs. (17) and (18) at all real values of $`k`$ and $`q`$ with the four free parameters is a formidable algebraic problem, it is possible to link the stability condition for the physical system to that which was studied in detail for the 1D system (1), (2) in Ref. . An algebraic transformation shows that, if the condition $`\mathrm{Re}\left[\lambda (k)\right]0`$ holds at all real values of $`k`$ in Eq. (8) of Ref. , then the inequality $`\mathrm{Re}\left[\lambda (k,q)\right]0`$ is true at all real values of $`k`$ and $`q`$ in Eq. (18), or, in other words, the stability of the zero solution in the 1D case guarantees its stability in the 2D case for the physical model. However, rather than following formal algebra, it is easier to understand this result from Fig. 1, which shows a 3D plot of the instability growth rate $`\mathrm{Re}\lambda `$ vs. $`k`$ and $`q`$ for a set of typical values of the parameters. Obviously, when $`k`$ is large, the stability is secured by the higher-order dissipative term in Eq. (10). The most dangerous case is that when $`k`$ and $`q`$ are relatively small. It can be seen from Fig. 1 that the growth rate, considered as a function of $`q`$, attains its maximum value at $`q=0`$ if $`k`$ is very small. As $`k`$ increases, two additional local maxima of the growth rate appear at nonzero $`q`$, and they may be greater than the local maximum at $`q=0`$. However, since the local maximum at $`q=0`$ becomes more negative as $`k`$ increases, the two side maxima remain negative too. Thus, the zero solution to Eqs. (10) and (11) is stable in the same parameter region in which it was found to be stable in the 1D system (1), (2) in Ref. .

In the phenomenological model, the zero solution is also stable in a certain parametric region. However, no simple relation of the stability condition to that in the 1D system can be found in this model.

## 4 The perturbation theory for two-dimensional solitary pulses

Both the physical and phenomenological models reduce to a zero-order system by setting $`\alpha =\gamma =\mathrm{\Gamma }=0`$, while keeping an arbitrary value of $`c`$. This zero-order system is conservative, consisting of the KP-I equation coupled to an extra linear one,

$$\left(u_t+uu_x+u_{xxx}v_x\right)_x=u_{yy},v_t+cv_x=u_x.$$
(21)
Looking for a solution to Eqs. (21) in the form of a soliton traveling at a constant velocity $`s`$ in the $`x`$-direction, so that

$$u(x,y,t)=u(\xi ,y),v(x,y,t)=v(\xi ,y),\mathrm{with}\xi xst,$$
(22)
we immediately conclude that, as well as in the 1D case, for such a solution we have

$$v(\xi ,y)=\left(cs\right)^1u(\xi ,y).$$
(23)
With regard to the relation (23), an *exact* solution to Eqs. (21) giving a 2D weakly localized soliton (“lump”) is

$$u(\xi ,y)=\frac{24\left(1+css^2\right)}{cs}w(\zeta ,z),$$
(24)
where

$$z\frac{1+css^2}{cs}y,\zeta \sqrt{\frac{1+css^2}{cs}}\xi ,w(\zeta ,z)\frac{z^2\zeta ^2+3}{\left[z^2+\zeta ^2+3\right]^2}.$$
(25)
This solution exists in the case

$$\frac{1+css^2}{cs}>0.$$
(26)
Note that, in contrast to the lump soliton solutions of the KP-I equation proper, which may only have positive velocities, the solitons (24) may also move in the negative direction, provided that the velocity satisfies the condition (26). In the particular case $`c=0`$, this condition means that either $`s>1`$, or $`1<s<0`$.

We have checked by direct simulations of Eqs. (21) that all the lump-type SP solutions are stable within the framework of the unperturbed equations (21). We have also simulated catch-up and head-on collisions between two lump pulses with different velocities. The results clearly show that the collisions are *inelastic*, hence, on the contrary to the KP-I equation, the conservative system (21) is not an exactly integrable one. An example of the head-on collision between two lump-type SPs with velocities $`s=2.0`$ and $`s=0.8`$ is shown in Fig. 2. The inelastic character of the interaction is obvious.

In fact, the collision shown in Fig. 2 is not quite generic, as it is generated by an initial configuration in which the mass $`N`$ of the $`v`$–component, see Eq. (15), is almost zero ($`N`$ vanishes exactly if the initial velocity of the second soliton is$`0.756`$ instead of $`0.8`$, which is the case in Fig. 2). In this case, a resonance probably happens, resulting in generation of an additional pair of solitons, which is suggested by Fig. 2 \[especially, by the panel 2(d)\]. Far from this resonant case, a typical collision (see an example in Fig. 3) is also strongly inelastic, but in this case the result may be regarded, in the first approximation, as fusion of two solitons into one.

To check the accuracy of the results produced by direct simulations of the 2D equations, including those displayed in Figs. 2 and 3 and in the next section, the simulations were repeated increasing the number of points. For instance, Figs. 2 and 3 were obtained in a 2D domain of $`256\times 256`$ points; repeating the same simulations with $`512\times 512`$ points, we had obtained the same pictures, without any visible difference.

Next, we restore the terms with the coefficients $`\alpha `$, $`\gamma `$ and $`\mathrm{\Gamma }`$ as small perturbations, with the aim to predict, at the first order of the perturbation theory, which solutions will be selected from the family of the lump solitons (24) found in the zero-order system (21). As it is suggested by the analysis of the 1D system , it is convenient to choose the relative velocity, $`\delta cs`$, as a parameter of the family.

As well as in the 1D case, we select stationary pulses, using the balance condition for the net field momentum,

$$P=\frac{1}{2}_{\mathrm{}}^+\mathrm{}𝑑x_{\mathrm{}}^+\mathrm{}𝑑y\left[u^2(x,y)+v^2(x,y)\right],$$
(27)
which is a dynamical invariant of the conservative (zero-order) system \[the masses (15) cannot be used for this purpose, as they remain dynamical invariants in the full dissipative models\]. The dissipative and gain terms give rise to the following exact balance equations for the momentum:

$$\frac{dP}{dt}=_{\mathrm{}}^+\mathrm{}𝑑x_{\mathrm{}}^+\mathrm{}𝑑y\left[\alpha u_x^2\gamma (u_{xx}^2+u_y^2)\mathrm{\Gamma }v_x^2\right]$$
(28)
in the physical system (10), (11), and

$$\frac{dP}{dt}=_{\mathrm{}}^+\mathrm{}𝑑x_{\mathrm{}}^+\mathrm{}𝑑y\left[\alpha u_x^2\gamma (u_{xx}+u_{yy})^2\mathrm{\Gamma }v_x^2\right]$$
(29)
in the phenomenological system (13), (14).

Substituting the unperturbed exact solution (24) into the equilibrium condition, $`dP/dt=0`$, which follows from Eqs. (28) or (29), we cast the equilibrium condition in the form

$$C_1\left(\alpha \frac{\mathrm{\Gamma }}{\delta ^2}\right)\left(c\delta +\frac{1}{\delta }\right)C_2\left(c\delta +\frac{1}{\delta }\right)^2=0,$$
(30)
where, for both the physical and phenomenological systems,

$$C_1_{\mathrm{}}^+\mathrm{}𝑑z_{\mathrm{}}^+\mathrm{}𝑑\zeta \left(\frac{w}{\zeta }\right)^2,$$
(31)
and the constant $`C_2`$ is defined differently for the two models:

$$C_2^{(\mathrm{phys})}=_{\mathrm{}}^+\mathrm{}𝑑z_{\mathrm{}}^+\mathrm{}𝑑\zeta \left[\left(\frac{^2w}{\zeta ^2}\right)^2+\left(\frac{w}{z}\right)^2\right],$$
(32)
$$C_2^{(\mathrm{phen})}=_{\mathrm{}}^+\mathrm{}𝑑z_{\mathrm{}}^+\mathrm{}𝑑\zeta \left(\frac{^2w}{\zeta ^2}+\frac{^2w}{z^2}\right)^2.$$
(33)
Numerical computation of the integrals yields $`C_1=0.34708`$, $`C_2^{(\mathrm{phys})}=0.75984`$, and $`C_2^{(\mathrm{phen})}=0.85770`$. Then, the equilibrium equation (30) yields the following cubic equation for $`\delta `$ (where $`\stackrel{~}{\alpha }\alpha /\gamma `$, $`\stackrel{~}{\mathrm{\Gamma }}\mathrm{\Gamma }/\gamma `$):

$$\delta ^3+\left(0.46\stackrel{~}{\alpha }c\right)\delta ^2\delta 0.46\stackrel{~}{\mathrm{\Gamma }}=0,$$
(34)
$$\delta ^3+\left(0.41\stackrel{~}{\alpha }c\right)\delta ^2\delta 0.41\stackrel{~}{\mathrm{\Gamma }}=0;$$
(35)
for the physical and phenomenological models, respectively. Note that these equations have the same general form as the one which selects equilibrium values of the parameter $`\delta `$ in the 1D system investigated in Ref. .

Physical roots of Eq. (34) or (35) are defined as those which not only are real, but also satisfy the condition (26). Physical roots select particular lump-type SPs, from the family of lump solitons of the zero-order system, that remain steady pulses in the presence of the perturbations. As well as in the 1D model, the existence of at least two physical roots is necessary for one pulse to be stable (the other one, with a smaller value of the amplitude, is then automatically unstable) . Numerical solution of Eqs. (34) and (35) verifies that both equations have exactly two physical roots in a broad parametric region, while three physical roots can never be found. Interestingly, the velocities corresponding to these two physical roots are both positive, i.e., a stable lump SP cannot move in the negative $`x`$–direction.

Regions in the parameter plane ($`\alpha ,\mathrm{\Gamma }`$) in which exactly two physical roots have been found are shown, for $`\gamma =0.05`$, in Figs. 4 and 5 for the physical model, with $`c=0`$ and $`c=1`$, respectively, and in Fig. 6 for the phenomenological model with $`c=1`$. These regions are bounded by two solid lines in Figs. 4 through 6. Beneath the lower solid line, Eq. (34) or (35) has three real roots, but no more than one of them is physical, and above the upper solid line, two physical roots bifurcate into a pair of complex ones, leaving no physical roots.

The dashed line is the border above which the zero solution is stable in the 1D model (1), (2), hence it is also stable in the 2D physical model, according to the results reported in the previous section. As is seen in Fig. 4, in the parameter plane ($`\alpha ,\mathrm{\Gamma }`$) of the physical system with $`c=0`$, there is no overlapping between regions where the zero solution is stable, and where the perturbation theory selects two physical solutions for the 2D pulses, but a narrow overlapping region is found in the same model for $`c=1`$, see Fig. 5. It is expected that stable lump-type SPs exist inside this narrow region, which is confirmed by direct simulations, as reported in the next section. Despite the existence of an overlapping region in Fig. 6, we have to check the condition (26) for Eq. (17), i.e., the stability of the zero solution in the phenomenological model, point by point, since the stability of the zero solution in the 1D model does not guarantee the stability of the zero solution in the phenomenological system, see the previous section.

The same analysis can be performed for values of the loss parameter $`\gamma `$ other than $`0.05`$, which was fixed in the above consideration. The results show that the variation of $`\gamma `$ produces little change in terms of the expected SP stability region. On the other hand, the group-velocity mismatch $`c`$ affects the stability region significantly. We had numerically found that there is no stability region in both physical and phenomenological model for $`c0`$. Then, by inspecting the region of negative $`c`$, it was found that there is a critical value $`c_{\mathrm{cr}}`$ such that stable pulses become possible for $`c<c_{\mathrm{cr}}`$. In the physical model, $`c_{\mathrm{cr}}0.8`$.

To conclude the analytical consideration, it is necessary to stress that the stability conditions obtained here are only necessary ones. Obtaining a full set of sufficient stability conditions might be possible within the framework of a rigorous spectral analysis of the full system linearized around the pulse solution, which was done (in a numerical form) for the 1D Benney (KS – KdV) equation in Refs. . Extending this analysis to the two-component 2D case is an extremely difficult problem, which is beyond the scope of the present work. Nevertheless, numerical results reported in the next section strongly suggest that the necessary stability conditions, that were obtained above by means of the simple analytical methods, are, most probably, sufficient for the stability of the 2D pulses.

## 5 Direct simulations of the two-dimensional solitary pulses

To check predictions produced by the above analysis for both the physical and phenomenological models, against direct numerical simulations, we have upgraded the implicit pseudo-spectral scheme, which was developed previously for the 1D model (1), (2) (see Appendix in Ref. ), to a scheme capable of dealing with both 2D models. It it relevant to mention that results obtained with a rather coarse grid, $`\mathrm{\Delta }x=\mathrm{\Delta }y=0.5`$, and relatively large time steps, $`\mathrm{\Delta }t=0.1`$, remained virtually unchanged if reproduced with an essentially finer grid and smaller time step. The initial conditions were taken as the lump-soliton solutions of the zero-order system, but with arbitrarily modified values of the amplitude, in order to check whether strongly perturbed pulses relax to stable ones, i.e., whether the stable pulses are attractors.

Collecting data produced by the systematic simulations, it has been found that, for the physical model, stable lump-type SPs do exist and are stable *everywhere* inside the stability region in the ($`\alpha ,\mathrm{\Gamma }`$) parameter plane which was predicted by the analytical perturbation theory, see Fig. 5. Moreover, all the stable pulses were found to be strong attractors indeed. A typical lump-type SP with $`\alpha =0.2`$, $`\gamma =0.05`$, $`c=1.0`$, and $`\mathrm{\Gamma }=0.55`$ is displayed in Fig. 7. The initial-state pulses selected as mentioned above definitely relax to this single stationary lump-type SP, provided that the initial amplitude $`A_0`$ of their $`u`$-component exceeds $`1.5`$. For instance, starting with the initial amplitudes $`A_0=1.96`$ and $`A_0=13.33`$ \[the corresponding values of the unperturbed-solitons’ velocity are $`s_0=0.8`$ and $`s_0=2.0`$, see Eq. (24)\], a lump-type SP develops with the values of the amplitude $`A_{\mathrm{num}}=7.11`$ and $`A_{\mathrm{num}}=7.16`$, and velocities $`s_{\mathrm{num}}=1.32`$ and $`s_{\mathrm{num}}=1.33`$, respectively, by the moment of time $`t=400`$. Meanwhile, the analytical prediction for the amplitude and velocity of the presumably stable pulse, given by Eq. (34) for the same values of the parameters, is $`A_{\mathrm{anal}}7.38`$ and $`s_{\mathrm{anal}}1.35`$. Thus, the numerical results match the theoretical prediction well.

On the other hand, if the initial amplitude is too small, e.g., $`A_0=0.89`$ ($`s_0=0.7`$), the pulse decays to zero, which is natural too, as the stable zero solution has its own attraction basin. Note that for the second (smaller) steady-state pulse, which is expected to play the role of the separatrix between the attraction basins of the stable pulse and zero solution, the perturbation theory predicts, in the same case, the amplitude $`\stackrel{~}{A}_{\mathrm{anal}}1.26`$, so it seems quite natural that the initial pulses with $`A_0=1.96`$ and $`A_0=0.89`$ relax, respectively, to the stable pulse and to zero.

For the phenomenological model, no definite stability region for the pulses has been found above; instead, the stability of the zero solution at each point inside the region bounded by the solid lines in Fig. 6 should be checked separately. For instance, the zero solution is found to be stable at the parameter point $`(\alpha =0.2,\mathrm{\Gamma }=0.55)`$ in Fig. 6. At this point, the two physical roots of Eq. (35) are $`2.073`$ and $`1.786`$. The amplitudes of the corresponding solitons predicted by the perturbation theory are $`\stackrel{~}{A}_{\mathrm{anal}}4.72`$ and $`\stackrel{~}{A}_{\mathrm{anal}}1.81`$, respectively; the one with the larger amplitude may be stable according to the general principle formulated above. Simulations demonstrate that, indeed, there is a stable 2D pulse with a shape very close to the predicted one, all the initial states in the form of the conservative-model lump, whose amplitude exceeds $`2.0`$, relaxing to this stable SP. Figure 8 shows the shapes of the $`u`$– and $`v`$–components of the SP at $`t=400`$, generated by the initial configuration taken as the soliton of the zero-order system (21) with the amplitude $`13.33`$. The amplitude of the thus obtained $`u`$-components is about $`5.14`$, the discrepancy with the analytical prediction being less than $`9\%`$.

It was also found that all the initial lump pulses whose amplitude is less than $`1.6`$ decay to zero. This fact complies well with the expectation that the second (smaller-amplitude) steady-state pulse, predicted by the perturbation theory, whose amplitude is $`1.81`$, ought to be the separatrix between the attraction basins of the stable pulse and zero solution in the phenomenological system.

The simulations also show that, even if the zero background is unstable, lump-type SPs may be quite stable when the integration domain (supplemented by periodic boundary conditions in both $`x`$ and $`y`$) is not very large. An explanation given in Ref.  to a similar “overstability” effect found in the 1D model, which is based on suppression of nascent perturbations by the pulse periodically traveling around the domain, applies to the 2D case as well.

Numerical simulations of interactions between two or more stable lump-type SPs were also carried out. Unlike the 1D system (1), (2), in which stable bound states of two and three pulses were easily found  (and unlike a physically relevant 2D model of a different type, which follows the pattern of the Zakharov-Kuznetsov equation  and will be considered elsewhere), bound states of pulses have *not* been found in both the physical and phenomenological systems based on the KP-I equation. Actually, this fact may be benign for experimental observation of the 2D pulses, as there will be no competition with multi-humped structures that might render the picture much more complex.

## 6 Conclusion

In this work, we have extended the 1D stabilized Kuramoto - Sivashinsky system to the 2D case. The 2D model is, quite naturally, based on the Kadomtsev-Petviashvili (KP) equation supplemented by loss and gain terms. A model with the dissipative terms of the Newell-Whitehead-Segel type was derived in a consistent way for the weakly 2D case. The derivation was outlined for a particular physical system, a downflowing liquid film carrying a surfactant that diffuses along its surface. It was argued that, in fact, the derived model is generic, therefore it also applies, for instance, to the description of interfacial waves in a two-layered flowing film, and double-front flames in a gaseous mixture. Another model, which may apply to a more isotropic situation, was put forward as a semi-phenomenological one. Both models consist of the generalized KP-I equation with gain and loss terms, linearly coupled to an extra equation of the advection-diffusion type. The additional linear equation stabilizes the system’s zero solution, thus paving the way to the existence of completely stable 2D localized solitary pulses, which are objects of an obvious physical interest. Treating the losses and gain as small perturbations, and employing the balance equation for the net field momentum, we have found that the condition of the equilibrium between the losses and gain may select two steady-state 2D solitons from their continuous family, which was found, in an exact analytical form, in the absence of the loss and gain (the exact solitons are similar to the “lump” solutions of the KP-I equation). When the zero solution is stable and, simultaneously, two lump-type pulses are picked up by the balance equation for the momentum, the pulse with the larger value of the amplitude is expected to be stable in the infinitely long system, while the other pulse must be unstable, playing the role of a separatrix between attraction domains of the stable pulse and zero solution. These predictions have been completely confirmed by direct simulations for both the physical and phenomenological systems. Another noteworthy finding is that, unlike their 1D counterpart, both KP-based 2D systems do not generate stable bound states of pulses.

Figure Captions

Fig. 1. The instability growth rate $`\mathrm{Re}\lambda `$ for the zero solution in the physical model \[Eqs. (10) and (11)\] vs. the longitudinal and transverse perturbation wave numbers $`k`$ and $`q`$. The values of the parameters are $`\alpha =0.2`$, $`\gamma =0.05`$, $`c=1.0`$, and $`\mathrm{\Gamma }=0.55`$.

Fig. 2. An example of the inelastic collision between two stable lump-type solitons with different velocities, simulated within the framework of the zero-order conservative system (21) with $`c=0`$. The $`u`$\- and $`v`$-fields prior to the collision (at $`t=0`$) are shown in panels (a) and (b), and after the collision (at $`t=40`$) – in panels (c) and (d). The initial velocities of the two solitons are $`s_1=2.0`$ and $`s_2=0.8`$.

Fig. 3. The inelastic collision in the conservative system (21) with $`c=0`$ in the case when the initial velocities of the two solitons are $`s_1=2.0`$ and $`s_2=0.5`$. The panels have the same meaning as in Fig. 2, with a difference that the final state is displayed for $`t=45`$.

Fig. 4. The stability region of the zero solution, and the region inside which two physical roots exist for Eq. (34) with $`\gamma =0.05`$ and $`c=0`$. The zero solution is stable above the dashed line, and inside the region bounded by two solid lines, Eq. (34) produces two physical solutions.

Fig. 5. The expected stability region for the solitary pulses in the parametric plane ($`\alpha `$,$`\mathrm{\Gamma }`$) for $`\gamma =0.05`$ and $`c=1.0`$. The solid and dashed lines have the same meaning as in Fig. 4.

Fig. 6. The stability region of the zero solution of the 1D model (1), (2) and the region inside which two physical roots exist for Eq. (35) when $`\gamma =0.05`$ and $`c=1`$. The solid and dashed lines have the same meaning as in Fig. 4.

Fig. 7. A stable lump-type solitary pulse in the physical model, found as a result of long numerically simulated evolution in the case $`\alpha =0.2`$, $`\gamma =0.05`$, $`c=1.0`$, and $`\mathrm{\Gamma }=0.55`$. The panels (a) and (b) show established shapes of the $`u`$ and $`v`$ fields, respectively.

Fig. 8. A stable lump-type solitary pulse in the phenomenological model, found the same way as the pulse shown in Fig. 7 for the same values of the parameters. The panels (a) and (b) have the same meaning as in Fig. 7.