text:
# The Parkes Multibeam Pulsar Survey – II. Discovery and Timing of 120 Pulsars

## 1 INTRODUCTION

The Parkes multibeam survey is an on-going survey of a $`10\mathrm{°}`$-wide strip along the Galactic plane ($`|b|<5\mathrm{°}`$ and $`l=260\mathrm{°}`$ to $`l=50\mathrm{°}`$). The survey aims to detect a large sample of pulsars for population studies (more than 600 have been discovered so far) as well as find individual objects of interest including binary pulsars, young pulsars, EGRET source associations and high magnetic field pulsars. The survey uses a 13-beam receiver on the 64-m Parkes radio telescope, receiving two polarisations per beam over a 288-MHz bandwidth centred on 1374 MHz. The survey’s receiver system, data acquisition and control procedures, offline data analysis and calculations of sensitivity are published in Manchester et al. (2001) (hereafter Paper I) along with the discovery and timing parameters for 100 pulsars with declinations $`<35\mathrm{°}`$. Timing observations of the majority of new pulsars with declinations $`>35\mathrm{°}`$ are made using the Jodrell Bank 76-m Lovell radio telescope. Positional and rotational parameters for 120 high-declination pulsars are given in this paper along with corresponding derived parameters and averaged pulse profiles111The data release policy and results for the survey may be found at http://www.atnf.csiro.au/research/pulsar/pmsurv.. The parameters for a further 150 pulsars are given in Bell et al. (in preparation), hereafter Paper III.

In order to facilitate follow-up observations, it is desirable to determine a better position estimate than that available from the survey data. For this reason, a technique known as ‘gridding’ is carried out on newly discovered pulsars before long-term timing begins. This technique, introduced in Paper I, is described in detail here. We also describe the Jodrell Bank observing system, summarise some interesting properties of the sample of 370 pulsars published here and in Papers I and III, and highlight two pulsars: PSR J1734$``$3333, a young pulsar with a high surface magnetic field and PSR J1830$``$1135, a pulsar with a 6-s period.

## 2 TIMING OBSERVATIONS AND ANALYSIS

### 2.1 Gridding

The centre of the discovery beam defines the initial position given to a pulsar candidate. If the pulsar is offset from the centre of the $`14\stackrel{}{.}4`$ beam then it will be detected with a lower signal-to-noise ratio (S/N) than for an equivalent observation centred on the pulsar. Furthermore, if the pulsar is to be observed with the narrower beam of the Lovell telescope then the S/N may be further reduced, hence the most accurate position estimate available should be used until the precise position is determined from timing. The use of a more accurate position maximises the S/N and therefore reduces the amount of telescope time required to obtain a pulse arrival time. Improved knowledge of the position also reduces the density of observations required through the year to obtain a coherent solution without period ambiguities in the timing residuals.

In order to improve the position estimate, a ‘gridding procedure’ is carried out using the Parkes telescope. The original discovery position is re-observed (with a shorter integration time of a few minutes — the exact time depending on the S/N in the original discovery) along with four other integrations of similar durations at positions $`9^{}`$ to the North, South, East and West. The peak of a Gaussian beam shape (with a half-power width of $`14\stackrel{}{.}4`$) fitted to the S/Ns in these consecutive observations determines the pulsar position with an uncertainty substantially smaller than a single beamwidth. Fig. 1a contains a scatter plot of the offset between the initial candidate position and the precise position determined after one year of timing observations for all of the pulsars with gridding observations (many of the earliest candidates were not ‘gridded’). Many pulsars lie outside the half-power beamwidth: the root-mean-square (rms) positional uncertainty is $`7\stackrel{}{.}2`$. However, the rms positional error of $`1\stackrel{}{.}6`$ for the differences between the gridded and precise positions (Fig. 1b) is well within the beamwidth, and the typical area of uncertainty has been reduced by more than an order of magnitude.

### 2.2 Timing at Jodrell Bank

For most of the timing observations, we use the Lovell telescope at Jodrell Bank at a frequency near 1400 MHz. We employ a dual-channel cryogenic receiver system sensitive to the two hands of circular polarisation. Each hand of polarisation is down-converted, fed through a multichannel filterbank and digitized. Until August 1999, the filterbank consisted of $`2\times 32\times 3`$ MHz channels centred on 1376 MHz. Subsequent observations use a filterbank with $`2\times 64\times 1`$ MHz channels centred on 1396 MHz.

In the observing procedure, data are dedispersed and folded on-line according to the pulsar’s dispersion measure (DM) and topocentric period. The folded pulse profiles for each polarisation are stored for subsequent analysis. The duration of each observation is between 6 and 36 minutes (depending upon the flux of the pulsar) and is divided into a number of sub-integrations, each lasting between one and three minutes. During offline data processing, the polarisations are combined to produce the total intensity and the sub-integrations are added to produce a final averaged profile for the observation. Pulse times-of-arrival (TOAs) are determined by cross-correlating the profile with a template of high S/N. These TOAs are corrected to the solar system barycentre using the Jet Propulsion Laboratory DE200 solar-system ephemeris \[Standish 1982\]. The positional and rotational parameters for the pulsar are then obtained by model-fitting the TOAs using tempo222See http://pulsar.princeton.edu/tempo..

The search code produces an initial estimate of a pulsar’s DM, and observations at Jodrell Bank use a hardware dedisperser giving an averaged pulse profile at a particular DM. In order to make a more precise measurement, we use either the discovery or the gridding observations that used the 288-MHz bandwidth available at Parkes. The band is divided into four sub-bands and the DM is obtained by fitting to the TOAs of the pulses determined from each sub-band.

The flux density scale at Jodrell Bank is calculated by switching a noise diode on and off between sub-integrations. The noise diode is calibrated using standard flux calibrators such as 3C295. The pulsar flux density for a single observation is then determined as the mean flux density in the profile averaged over the pulse period.

Included in this paper are ten pulsars with low flux densities and with declinations $`>35\mathrm{°}`$ that were timed at the Parkes telescope because of the higher sensitivity afforded by the use of the larger bandwidth. The observing and data-reduction techniques used at Parkes are described in Paper I. A few pulsars were observed with both instruments, but for consistency, only Jodrell Bank data were used to produce the pulse profiles. Similarly, nine pulsars were timed on a monthly basis with the 305-m Arecibo telescope, using a $`2\times 128\times 0.0625`$ MHz filterbank centred on 1400 MHz. These pulsars were observed for 5-minute or shorter integrations, with the polarisations summed and the data folded on-line using the known DM and the predicted topocentric pulse period. TOAs were determined using a standard template as for the Jodrell Bank observations.

## 3 DISCOVERY AND TIMING OF 120 PULSARS

The two tables giving positional, discovery and flux parameters (Table 1) and rotational parameters (Table 2) take the same form as the equivalent tables in Paper I. In summary, Table 1 lists the positional properties of each pulsar: the pulsar’s name, J2000 right ascension and declination from the timing solution and the corresponding Galactic coordinates. This table also contains discovery information: the beam in which the pulsar was detected, the radial angular distance of the pulsar from the beam centre and the S/N of the pulse profile in the discovery observation. The median flux density averaged over all the observations included in the timing solution and the pulse widths at 50 per cent and 10 per cent of the peak of the mean pulse profile are also included (the 10 per cent width is only measured for pulsars with high S/N).

Table 2 contains the solar-system barycentric pulse period $`P`$, period derivative $`\dot{P}`$, the epoch to which the period refers to, the number of TOAs used in the timing solution, the MJD range covered by the timing observations, the final rms value for the timing residuals and the pulsar’s DM. A non-deterministic period second derivative has also been fitted to a few pulsars, indicated by asterisks in Table 2, which show significant amounts of timing noise.

Data sets for each pulsar have been folded at twice and three times the tabulated period in order to confirm that this represents the fundamental period for the pulsar. The TOAs range from October 1997 to March 2002 with the shortest data span being 1.0 years and the longest 4.2 years.

This paper includes data for four binary pulsars which have previously been published: PSRs J1740$``$3052 \[Stairs et al. 2001\], J1810$``$2005, J1904$`+`$0412 \[Camilo et al. 2001\] and J1811$``$1736 \[Lyne et al. 2000\]. For completeness, their binary parameters are summarised in Table 3. PSR J1837$``$0604, a young energetic pulsar, has been discussed in D’Amico et al. (2001) and between MJDs 51305 and 51894, a glitch occurred in PSR 1806$``$2125 with a magnitude $``$2.5 times greater than any previously observed event \[Hobbs et al. 2002\].

For each pulsar, the profiles used in the determination of the TOAs for the timing analysis were averaged to form mean pulse profiles at 20 cm (Fig. 2). For each profile, the pulsar’s name, pulse period and DM are given along with the effective time resolution of the profile, and pulsars timed at Parkes or Arecibo are indicated. There are several small features in the baselines of some pulsar profiles which may be other pulse components or may be radio frequency interference; higher S/N observations are required to resolve this issue.

## 4 DISCUSSION

Table 4 lists derived parameters for the 120 pulsars. The first column contains the pulsar’s name, followed by $`\mathrm{log}_{10}`$ of the characteristic age, $`\tau _c=P/(2\dot{P})`$ in years, the surface dipole magnetic field strength, $`B_s=3.2\times 10^{19}(P\dot{P})^{1/2}`$ in gauss and the rate of loss of rotational energy, $`\dot{E}=4\pi ^2I\dot{P}P^3`$ in erg s-1, where a neutron star with moment of inertia $`I=10^{45}`$ g cm2 is assumed. The next two columns give the pulsar distance, $`d`$, computed from the DM assuming the Taylor & Cordes (1993) model for the Galactic distribution of free electrons, together with the corresponding Galactic $`z`$-height. The uncertainty in distance from this model is generally far greater than the precision with which the distance is quoted. The final column in Table 4 gives the radio luminosity $`L_{1400}S_{1400}d^2`$.

A $`P`$$`\dot{P}`$ diagram is shown in Fig. 3. The pulsars published here and in Papers I and III are overlaid on previously known pulsars. Anomalous X-ray pulsars and soft $`\gamma `$-ray repeaters are included along with lines of constant surface magnetic field strength, a pulsar emission ‘death line’ and a theoretical boundary separating radio loud and radio quiet pulsars.

### 4.1 The sample of 370 pulsars

A full population analysis will follow the completion of the survey. Here, we briefly highlight some of the salient features of the sample of pulsars presented in this paper and in Papers I and III. To compare the Parkes multibeam pulsars with previously known pulsars, Paper I displayed histograms of period, DM and flux density. We continue this comparison with histograms of pulse widths, ages, magnetic fields and rate of loss of rotational energy for the 370 pulsars (Fig. 4). There is a clear similarity between the measured widths of the profiles of these pulsars and of previously known pulsars (Fig. 4a). The characteristic ages for our sample are lower than previously measured characteristic ages (Fig. 4b). The observed age difference likely arises because we search along the plane of the Galaxy, where pulsars are born and where older pulsars have drifted from. This is also indicated by the results of the Swinburne intermediate latitude survey which uses the same observing system and has discovered many intermediate age to old pulsars \[Edwards et al. 2001\]. Lorimer et al. (1995) suggest that young pulsars have flatter radio spectra than older pulsars. Our high observing frequency may, therefore, be an additional reason for the discovery of so many young pulsars. However, this difference in spectral indices is not confirmed by Malofeev (1996) or Maron et al. (2000).

The possibility of a relationship between characteristic age and luminosity was also studied; the results are shown in Fig. 5a which is consistent with a constant radio luminosity for ages $`40`$ kyr–20 Myr, although there seems to be a significant reduction in luminosity for even older ages. The surface magnetic field histogram (Fig. 4c) displays two features: the lack of pulsars with low magnetic fields (reflecting the absence of millisecond and moderately recycled pulsars in our sample) and the relatively large number of high magnetic field pulsars, due to an excess of long period pulsars. A graph of surface magnetic field versus luminosity (Fig. 5b) is suggestive of the possibility that high magnetic field pulsars have a higher luminosity than their lower-field counterparts. A plot of 1400-MHz luminosity versus $`\dot{E}`$ for our sample (Fig. 5c) shows no indication of any trend, however, these figures do not take any geometrical beaming effects into account.

The histogram of the rate of loss of rotational energy (Fig. 4d) contains a component representing highly energetic pulsars ($`\dot{E}>10^{36}`$ erg s-1). Three (PSRs J1809$``$1917, J1828$``$1101 and J1837$``$0604) are both young ($`10^4<\tau _c<10^5`$ yr) and energetic. PSR J1913+1011 is slightly older ($`\tau _c=170`$ kyr), but has similar parameters to PSR B1951+32 in the supernova remnant CTB 80. These ‘Vela-like’ pulsars are discussed in Paper III.

The average pulse profiles (Fig. 2) are diverse: the majority contain a single component, $``$ 20% contain multiple components and three (PSRs J1806$``$1920, J1828$``$1101 and J1913$`+`$0832) clearly have interpulses. Eight (PSRs J1734$``$3333, J1735$``$3258, J1809$``$2004, J1811$``$1736, J1818$``$1541, J1828$``$1101, J1833$``$0559 and J1837$``$0604) have features suggestive of scattering tails. Although in general pulsars showing evidence of scattering also have high DM, we find no simple relationship between DM and the width of the scattering tail (see Paper I).

### 4.2 PSR J1734$``$3333

Camilo et al. (2000) discuss PSR J1814$``$1744, whose inferred surface magnetic field strength, $`B_s=5.5\times 10^{13}`$ G, is the largest of any known radio pulsar. PSR J1734$``$3333 has a similar field strength of $`5.2\times 10^{13}`$ G (and has a very small age, $`\tau _c=8`$ kyr). Both these values were calculated using the standard expression333Note that occasionally in the literature the magnetic field is calculated with an additional factor of two. In that case the field represents the maximum surface dipole field, as opposed to the value at the equator. $`B_s=3.2\times 10^{19}(P\dot{P})^{1/2}`$ G. Taking these fields at face value, these pulsars thus have surface magnetic field strengths that are larger than the ‘quantum critical field’, $`B_c=4.4\times 10^{13}`$ G, at which the cyclotron energy equals the electron rest-mass energy. Baring & Harding (2001) propose that photon splitting may inhibit pair creation at these high magnetic fields and therefore hypothesize the existence of a boundary separating radio-loud and radio-quiet pulsars. This boundary is indicated on Fig. 3 assuming that the emission occurs from the surface of the neutron star and that the photons propagate parallel to the local field. PSRs J1734$``$3333 and J1814$``$1744 lie in the predicted radio quiescence region. However, as detailed in Baring & Harding (2001) the boundary moves up the $`P`$$`\dot{P}`$ diagram for emission away from the neutron star’s surface. For an emission radius of 1.5 times the neutron star’s radius, the boundary lies above all the known radio pulsars. We note that given the trend toward pulse narrowing with increasing spin period \[Gould 1994\], the fact that no radio counterparts of AXPs, which also lie above this boundary, have been detected, may be due to their radio beams not intersecting our line of sight \[Gaensler et al. 2001\].

The success of this survey at finding new pulsars with fields close to and above the critical field will help improve our understanding of the delineation between radio-loud and radio-quiet pulsars.

### 4.3 PSR J1830$``$1135

After the discovery of PSR J2144$``$3933, a radio pulsar with an 8.5-s period \[Young, Manchester & Johnston 1999\], there has been much interest in detecting long-period pulsars. PSR J1830$``$1135, with a period of 6-s, is the second slowest known radio pulsar. It lies on the $`P`$$`\dot{P}`$ diagram at the extreme right-hand edge of the ‘normal’ pulsars, but due to a relatively high surface magnetic field of $`10^{13}`$ Gauss, lies well above the ‘death line’ (Fig. 3). Although the magnetic field is below the critical field, PSR J1830$``$1135 (and PSR J1814$``$1744) lies roughly in the neighbourhood of some AXPs (shown as open crosses in Fig. 3). Pivovaroff, Kaspi & Camilo (2000) infer, using PSR J1814$``$1744, that the X-ray emission from AXPs depends upon more than the implied surface magnetic field strength. The reasons for this are, however, not well understood.

Among the 370 published Parkes multibeam pulsars, only PSR J1830$``$1135 has a period greater than 6 seconds. Although the search software limits candidates to periods below 5 seconds (see Paper I), strong long-period pulsars with narrow pulse profiles, such as PSR J1830$``$1135, may be detected from their second or third harmonics. The probability of detecting a pulsar is also dependent upon the width of the pulsar beam, which decreases with increasing period (e.g. Rankin 1983 and Gould 1994). Although determination of the physical beam size requires polarisation information, the small angular pulse width of $`4\mathrm{°}`$ for PSR J1830$``$1135 is compatible with current understanding.

### 4.4 Conclusion

The sample of 370 pulsars reported here and in Papers I and III represents approximately half of the total number of new pulsars that we expect to discover on completion of this survey. Details of the remaining pulsars will be published in due course as the data from the timing observations become available. Once complete, the survey will provide a large and relatively unbiased dataset which may be used to study many aspects of pulsar physics.

## Acknowledgements

We gratefully acknowledge the technical assistance provided by George Loone, Tim Ikin, Mike Kesteven, Mark Leach and the staff at the Parkes and Jodrell Bank Observatories toward the development of the Parkes multibeam pulsar system. We thank R. Bhat, D. Lorimer, M. McLaughlin and P. Freire for assistance with observations. DJM and GH acknowledge the receipt of post-graduate studentships from the UK Particle Physics and Astronomy Research Council. IHS received support from NSERC and Jansky postdoctoral fellowships. FC is supported by NASA grants NAG 5-9095 and NAG 5-9950. VMK is an Alfred P. Sloan Research Fellow and was supported in part by a US National Science Foundation (NSF) Career Award (AST-9875897) and by a Natural Sciences and Engineering Research Council of Canada grant (RGPIN 228738-00). The Parkes radio telescope is part of the Australia Telescope which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. The Arecibo Observatory is part of the National Astronomy and Ionosphere Center, which is operated by Cornell University under a cooperative agreement with the US NSF.


format: text

text:
# Pulsed radiation from neutron star winds

## 1 Introduction

Recently, the non-axisymmetric “striped pulsar wind” investigated by Coroniti (1990) and Michel (1994) has been reexamined (Kirk & Lyubarsky, 2001). A striped wind is a structure similar to a Parker spiral, produced by the outward radial advection of magnetic field lines with foot points anchored on the surface of the rotating neutron star. If the field lines originate from the poles of an obliquely rotating dipole, there exists a region around the equatorial plane where the polarity of the field at a fixed radius reverses at the rotation period, producing stripes of alternating magnetic field direction. Magnetic energy is dissipated into particle energy if reconnection occurs at the stripe boundaries. Coroniti (1990) proposed that this kind of dissipation proceeds at a rate sufficient to maintain the thickness of the current sheet approximately equal to the gyro-radius of the heated particles. Michel (1994) used essentially the same criterion, phrased in terms of the velocity of the current carriers. These early papers concluded that observations of the Crab Nebula, which indicate conversion of Poynting flux to particle-borne energy flux within $`10^9r_\mathrm{L}`$, (where $`r_\mathrm{L}`$ is the radius of the light cylinder) could be explained by this mechanism. However, Lyubarsky & Kirk (2001) point out that reconnection is accompanied by acceleration of the wind. Both processes occur on the same timescale in the Coroniti/Michel approach, which is also the timescale on which the plasma expands. As a result, time dilation in the accelerated flow significantly reduces the dissipation rate. In the case of the Crab, conversion of Poynting flux to particle-borne energy flux is not achieved until $`10^{12}r_\mathrm{L}`$ — too slow to account for the observations \[the $`\sigma `$-problem; see Melatos (1998)\].

However, our present understanding of the physics of the dissipation process does not exclude the possibility that it proceeds much more rapidly than suggested by Coroniti and Michel. Recent investigations of reconnection in a relativistic pair plasma indicate that the timescale is related to the crossing time of a fast magnetosonic wave (Zenitani & Hoshino, 2001). In the case of a transsonic wind, the process could then be completed within a few pulsar periods, giving the plasma insufficient time to accelerate to high Lorentz factor. If rapid reconnection in fact occurs, it will be accompanied by radiation losses (Melatos, 1998) that may give rise to an observable signature. In this letter we adopt this hypothesis. We assume that dissipation is triggered at a surface $`r=r_0`$ well outside the light cylinder in the supersonic MHD wind of a pulsar and takes place rapidly. Using two simple models of the electron distribution we show (i) that for the winds of young pulsars (e.g., the Crab) the predicted emission is pulsed, with, in general, a pulse and interpulse that are not symmetrically spaced in phase and (ii) that the predicted synchrotron luminosity in the pulses agrees to order of magnitude with that observed in the high energy emission of the Crab pulsar.

Our calculations suggest that the high energy pulses emitted by young pulsars originate as synchrotron emission in the pulsar wind, far outside the light cylinder. This is in marked contrast with conventional models of high-energy pulse production at an inner or outer gap \[for a review see Harding (2001)\] and also with magnetospheric synchrotron based models (Machabeli et al, 2000; Crusius-Wätzel et al, 2001). The current sheet just beyond the region of closed field lines in the magnetosphere has been proposed as the source of high energy radiation from the Crab pulsar (Lyubarsky, 1996). But, in common with the others, this model relies for pulse production on the corotation with the star of a beam of radiation emitted from within the light cylinder.

## 2 Geometry of the emission region

Any radial, relativistic flow containing a periodic modulation of the emissivity is likely to appear pulsed to an observer because of two effects. Firstly, the strong Doppler boosting of approaching parts of the flow means that only a small cone of the flow is visible, propagating within an angle of roughly $`1/\mathrm{\Gamma }`$ with respect to the line of sight, where $`\mathrm{\Gamma }=(1\beta ^2)^{1/2}`$ is the Lorentz factor and $`c\beta `$ the radial speed. The spread in arrival times at the observer of photons emitted at the same radius but from different parts within this cone is roughly $`r/(c\mathrm{\Gamma }^2)`$. Secondly, if emission is restricted to a range $`\mathrm{\Delta }r`$ in radius, the fact that the photons move only slightly faster than the flow when seen in the laboratory frame means that the spread in arrival times is roughly $`\mathrm{\Delta }r/(c\mathrm{\Gamma }^2)`$. Thus, for a wind modulated at the period $`P=2\pi r_\mathrm{L}/c`$ of a rotating star, the approximate condition for the observation of pulsed emission is

$`\mathrm{\Delta }r,r`$ $``$ $`\mathrm{\Gamma }^2r_\mathrm{L}`$ (1)
If this is fulfilled, the observed pulse shapes appear similar to the pattern of modulation of the flow. In particular, for a flow containing the corrugated spiral current sheet of an oblique rotator (Fig. 1), there should in general be two components (pulse and interpulse) which are symmetrically placed if the pulsar is viewed from the equatorial plane, but are asymmetrical when viewed from higher latitudes. Denoting by $`\alpha `$ the angle between the rotation axis and the magnetic axis, and by $`\zeta `$ the angle between the rotation axis and the viewing direction, these pulses appear only for $`|90^\mathrm{o}\alpha |<\zeta <180^\mathrm{o}|90^\mathrm{o}\alpha |`$.

To improve on the estimate (1) and calculate the expected pulse shapes we adopt the following model of the volume emissivity $`ϵ`$ of the wind:

* The emission region is located on the current sheet in the supersonic MHD wind of a rotating neutron star with an oblique, split-monopole magnetic field.
* The dependence on frequency $`\nu `$ of the volume emissivity is a power law: $`ϵ\nu ^a`$.
* The emission switches on abruptly as plasma in the sheet crosses the surface $`r=r_0`$. Thereafter the volume emissivity decays as a power law in radius: $`ϵr^{2q}`$.

Outside the light cylinder, assumption (i) should be a good approximation, independent of the actual configuration of the magnetic field close to the pulsar surface. The assumption that emission occurs only at the exact position of the current sheet is an idealisation — in reality, the entire flow may be modulated, making the observed pulse more complex. In this case, our computations provide the Green’s function with which the pattern is to be folded to find the observed pulse shape. If the wind carries enough particles to ensure MHD behaviour at least initially \[for an alternative approach see Michel & Li (1999)\], a simple solution is available in the ultra-relativistic limit (Bogovalov, 1999). In it, the velocity is radial and constant and the magnetic field is almost toroidal. The thin current sheet in such a wind is described by the equation

$`r`$ $`=`$ $`r_\mathrm{s}(\theta ,\varphi ,t)`$ (2)
$`=`$ $`\beta r_\mathrm{L}\left[\pm \mathrm{arccos}\left(\mathrm{cot}\alpha \mathrm{cot}\theta \right)+ct/r_\mathrm{L}\varphi +2n\pi \right]`$
(Kirk & Lyubarsky, 2001), where ($`r`$,$`\theta `$,$`\varphi `$) are polar coordinates, $`t`$ is the observer-frame time and $`n`$ is an integer.

Assumption (ii) ensures the pulse shape is independent of frequency. In reality, the electron distribution may exhibit breaks and cut-offs, which translates into different pulse shapes at different frequencies. Our assumption, however, ensures that all parts of the flow can contribute to the emission at a given frequency, making it harder to obtain sharply defined pulses.

The dependence of the observed flux $`F`$ on $`t`$, $`\alpha `$ and the direction of observation $`\widehat{n}`$ is:

$`F`$ $``$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}dt^{}{\displaystyle _{r_0}^{\mathrm{}}}dr{\displaystyle _{90^\mathrm{o}\alpha }^{90^\mathrm{o}+\alpha }}d\theta {\displaystyle _0^{2\pi }}d\varphi `$ (3)
$`r^qD^{2+a}\delta \left[rr_\mathrm{s}(\theta ,\varphi ,t^{})\right]\delta \left(t^{}t+\widehat{n}r/c\right)`$
where $`D=1/\left(1\beta \widehat{r}\widehat{n}\right)`$ is the conventionally defined Doppler factor, apart from a (constant) factor of $`\mathrm{\Gamma }`$. Choosing $`a=0`$ and $`q=3`$, which is appropriate for the synchrotron radiation of a fixed number of particles that undergo predominantly adiabatic losses in the expanding plasma, we present sample pulses for various Lorentz factors in Fig. 2. For the higher Lorentz factors this figure also displays the effect of a faster fall-off of the emissivity towards higher radii ($`q=6`$), mimicking the effect of an accelerating wind. The asymmetry in time of the pulses is produced by our assumption of a sudden switching on of the emission and is reduced by a more rapid switching off.

The parameters of this figure have been chosen to correspond to the Crab pulsar. The radio emission of this object is complex. The shape of the high energy pulses is repeated in the radio, suggesting a common emission site, which in our model would be the wind. But there are additional radio components (Moffatt & Hankins, 1996), one of which displays ‘core emission’ properties (Rankin, 1990). This suggests emission from close to the polar cap, with the magnetic axis passing within a few degrees of the line of sight. The core emission is in the precursor of the main pulse, which constrains the relative phase of the polar cap emission and wind emission in this interpretation. However, a complete solution linking the inner magnetosphere and the wind structure would be required to check this point. Nevertheless, core emission indicates that the inclination angle $`\alpha `$ of the magnetic axis to the rotation axis roughly equals the angle $`\zeta `$ between the rotation axis and the line of sight. The latter is determined by X-ray and optical observations of the ‘torus’ (Aschenbach & Brinkmann, 1975; Hester et al, 1995) to be $`60^\mathrm{o}`$. In Fig. 2 the emission has been assumed to switch on at $`r_0=30r_\mathrm{L}`$ and pulses appear for $`\mathrm{\Gamma }>5`$, in rough agreement with the estimate (1). At higher Lorentz factors, the pulses sharpen. From Eq. (2) it can be seen that the emission splits naturally into two contributions, corresponding to different signs of the term $`\mathrm{arccos}(\mathrm{cot}\alpha \mathrm{cot}\theta )`$. In the first two plots of Fig. 2 these contributions are plotted separately as dotted lines. Although the detailed pulse shapes will depend on the structure of the reconnecting sheets in the Crab wind, as well as the details of the switching on mechanism and the amount of acceleration of the wind, the separation of the two peaks is not sensitive to these unknowns, being determined by the location of the sheets. This separation in phase equals $`\mathrm{arccos}(\mathrm{cot}\alpha \mathrm{cot}\zeta )/180^\mathrm{o}`$, corresponding to $`0.39`$ for the Crab pulsar, in good agreement with observation \[e.g., Kanbach (1998); Kuiper et al (2001)\].

For the Vela pulsar, conflicting interpretations of the polarisation sweep and the X-ray morphology give $`\alpha 71^\mathrm{o}`$, $`\zeta 65^\mathrm{o}`$ (Radhakrishnan & Deshpande, 2001) or $`\alpha \zeta 53^\mathrm{o}`$ (Helfand et al, 2001). These imply separations of the gamma-ray (sub-)pulses of $`0.45`$ and $`0.31`$ of a period respectively, so that the observed spacing of $`0.4`$ marginally favours the former interpretation.

## 3 Luminosity estimates

Electrons heated in the current sheets will enter the surrounding magnetic field and emit synchrotron radiation. An estimate of the power radiated can be found by assuming a monoenergetic electron population of Lorentz factor $`\gamma `$ ($`1`$) and requiring pressure equilibrium between the hot plasma and the surrounding magnetic field. In the delta-function approximation for the synchrotron emissivity, the luminosity at frequency $`\nu `$ into a solid angle $`\mathrm{\Omega }`$ is

$`\mathrm{d}L/\mathrm{d}\nu \mathrm{d}\mathrm{\Omega }`$ $``$ $`{\displaystyle _{r_0}^{\mathrm{}}}drr^2\mathrm{\Delta }\mathrm{\Gamma }n_\mathrm{h}^{}P_{\mathrm{s}.\mathrm{p}.}\delta \left(\nu \nu _0\right)`$ (4)
(Skjæraasen & Kirk, 2001), where $`\mathrm{\Delta }`$ ($`1`$) is the fraction of a wavelength of the wind pattern occupied by hot electrons of proper density $`n_\mathrm{h}^{}`$, $`P_{\mathrm{s}.\mathrm{p}.}=(4/3)c\sigma _\mathrm{T}\gamma ^2(B^2/8\pi )`$ is the synchrotron power emitted by a single electron, and $`\nu _0`$ is the characteristic frequency of emission, measured in the observer frame ($`\sigma _\mathrm{T}`$ is the Thomson cross section). Pressure equilibrium implies $`n_\mathrm{h}^{}\gamma mc^2/3=B^2/8\pi `$, where $`B^{}`$ is the magnetic field in the comoving frame and $`m`$ the electron mass. Well outside the light cylinder, but before energy release begins, i.e., at $`r_\mathrm{L}<r<r_0`$, the density and magnetic field scale as $`r^2`$ and $`r^1`$ respectively, and can be characterised by the values $`n_\mathrm{L}`$ and $`B_\mathrm{L}`$ obtained by extrapolating back to $`r=r_\mathrm{L}`$. Provided the flow does not accelerate significantly whilst emitting synchrotron radiation, these same scalings, which imply $`\gamma `$ independent of $`r`$, can be used in the integration in Eq. (4). This is convenient for an estimate of the luminosity, but is sensitive to the assumption that the wind does not accelerate upon reconnection.

If we identify the ratio of Poynting flux to kinetic energy (at the light cylinder) $`\sigma _\mathrm{L}=B_\mathrm{L}^2/(4\pi \mathrm{\Gamma }n_\mathrm{L}mc^2)`$, define an average fraction of hot particles $`\overline{\mathrm{\Delta }}=r_0_{r_0}^{\mathrm{}}dr\mathrm{\Delta }/r^21`$, and use the standard relation $`B_\mathrm{L}=2.9\times 10^8\dot{P}^{1/2}P^{5/2}`$G, the luminosity becomes

$`{\displaystyle \frac{\mathrm{d}\widehat{L}}{\mathrm{d}\mathrm{\Omega }}}`$ $`=`$ $`3.9\times 10^7\left({\displaystyle \frac{\sigma \overline{\mathrm{\Delta }}r_\mathrm{L}\dot{P}}{\mathrm{\Gamma }^3r_0P^4}}\right)(P\mathrm{in}\mathrm{secs})`$ (5)
where $`\mathrm{d}\widehat{L}/\mathrm{d}\mathrm{\Omega }`$ is the synchrotron luminosity per steradian divided by the wind luminosity per steradian. Similarly, the characteristic frequency is

$`h\nu _0(r)`$ $`=`$ $`2.6\times 10^6\left(r_\mathrm{L}/r\right)\sigma ^2\dot{P}^{1/2}P^{5/2}\mathrm{MeV}.`$ (6)
The spectrum of this model with monoenergetic electrons extends to a maximum frequency given by Eq. (6) with $`r=r_0`$ and is flat: $`\mathrm{d}L/\mathrm{d}\nu \mathrm{d}\mathrm{\Omega }\nu ^0`$. However, this aspect is sensitive to expansion and acceleration of the plasma after reconnection. If, in a more realistic scenario, the electron density is a power law $`N(\gamma )\gamma ^p`$, with $`p>2`$, then most of the energy resides in the particles with lowest $`\gamma `$ and the above estimate (6) refers to the photons emitted by these particles. The emission at frequencies $`\nu <\nu _0(r_0)`$ is dominated by the same low energy particles radiating at larger radius in a weaker magnetic field, and is determined by the expansion and acceleration of the hot plasma. In the absence of acceleration, the flat spectrum derived above prevails. However, the intrinsic spectrum, which is a power law $`\mathrm{d}L/\mathrm{d}\nu \mathrm{d}\mathrm{\Omega }\nu ^a`$, with $`a>0.5`$, is revealed at frequencies $`\nu >\nu _0(r_0)`$, where particles radiating at $`r_0`$ dominate. In the case of the Crab pulsar, the optical to X-ray emission is consistent with a flat spectrum (Shearer & Golden, 2001), which steepens to $`a>1`$ at a photon energy of approximately $`1`$MeV (Kanbach, 1998; Kuiper et al, 2001). Equation (6) then implies $`\sigma =1.1\times 10^4\sqrt{r_0/r_\mathrm{L}}`$. Observations suggest $`\mathrm{d}\widehat{L}/\mathrm{d}\mathrm{\Omega }10^3`$, which, from Eq. (5), is consistent with a mildly supersonic flow, $`\mathrm{\Gamma }250`$, that starts to radiate at $`r_0=30r_\mathrm{L}`$. This satisfies the lower limit of $`5r_\mathrm{L}`$ (Bogovalov & Aharonian, 2000) on the radius at which Poynting flux can be converted to particle-borne energy flux.

## 4 Discussion

Current models of the high energy emission from rotation powered pulsars fall into two groups: “polar cap” and “outer gap” models \[see Harding (2001)\]. Each locates the emission region in the corotating magnetosphere of the star. The predicted pulse shapes depend sensitively on the uncertain geometry of the magnetic field in this region, usually assumed dipolar. In contrast, in our model, pulses are emitted outside the corotation region. If the pulsar drives a supersonic, MHD wind, the field geometry there approaches a simple asymptotic solution (Bogovalov, 1999), determining the basic properties of the pulses.

Pacini (1971) and Pacini & Salvati (1983, 1987) suggested that particles close to the light cylinder emit synchrotron radiation and are responsible for the optical pulses of the Crab and other pulsars, a theory which appears to be in reasonable agreement with observation (Shearer & Golden, 2001). Our model is similar in that the radiation mechanism is synchrotron emission and the scaling of Eq. (5) is close to that originally given by Pacini (1971). It goes further, however, by specifying for the emission region a precise geometrical structure and location.

It is known that the dilution of the plasma in the wind must lead to non-ideal MHD behaviour (Usov, 1975; Michel, 1982) and so can be responsible for triggering the emission. The position at which this occurs depends on the initial concentration of charges in the current sheet and is taken as a free parameter. The major uncertainty in the model is the speed with which dissipation proceeds. Our computations implicitly assume rapid dissipation (over a scale small compared to the radius). In this case, it is possible to make a rough estimate of the synchrotron luminosity which is in agreement with observations. It is conceivable that at least part of the radio emission could be produced in the wind region, as is suggested by the similarity of the pulse profiles at all frequencies seen in the Crab pulsar. However, until a candidate coherent mechanism can be identified, this connection remains very speculative. Nevertheless, the wind scenario is a viable alternative to current theories of gamma-ray, X-ray and optical pulses from rotation-driven pulsars and is capable of making testable predictions about the pulse profile.

###### Acknowledgements.

This work is a collaboration of the TMR Network “Astroplasmaphysics” of the European Commission, contract FMRX-CT-98-0168. Y.G. is supported by a Marie Curie Fellowship from the European Community, contract no. HPMFCT-2000-00671 under the IHP programme.


text:
# Surpassing the standard quantum limit for optical imaging using non-classical multimode light

## Abstract

Using continuous wave superposition of spatial modes, we demonstrate experimentally displacement measurement of a light beam below the standard quantum limit. Multimode squeezed light is obtained by mixing a vacuum squeezed beam and a coherent beam that are spatially orthogonal. Although the resultant beam is not squeezed, it is shown to have strong internal spatial correlations. We show that the position of such a light beam can be measured using a split detector with an increased precision compared to a classical beam. This method can be used to improve the sensitivity of small displacement measurements.

It has long been known that optical measurements are ultimately limited in their sensitivity by quantum noise, or shot noise, of the light. For more than a decade the usage of nonclassical light has provided ways of improving the sensitivity beyond this standard quantum limit Bachor . For example, squeezed light has been used to improve interferometric Interf  and absorption Absor  measurements. However, these improvements can only be applied to signals that correspond to the time modulation of light, as they rely on the temporal quality of the light. On the other hand, many applications require spatial measurement of light. While improvements for spatial applications based on non-classical light have been proposed theoretically Seng ; Fabre , no experimental demonstration has yet been shown to work with continuous wave light. The challenge is to create strong spatial correlations within a laser beam, rather than the temporal correlation typically found in non-classical light sources kolobov . While some experiments involving sub-Poissonian VCSELS operating in a transverse multimode regime exhibited a non-random spatial distribution of the quantum noise VCSELS , no spatial correlation was observed within the produced beam. Here we present the first successful experimental demonstration of a spatially ordered light source and a measurement of the spatial modulation of a laser beam position to below the standard quantum limit in the continuous wave regime.

This experimental work builds on theoretical work done on non-classical multimode states of light Lugiato2 . Such states display strong spatial correlations, and their productions have been the subject of extensive studies in recent years kolobov3 . In particular, the process of parametric down-conversion in a nonlinear optical medium has been extensively studied, as it produces “twin photons” which are quantum correlated both temporally and spatially. Such strong spatial quantum correlations in the plane perpendicular to the direction of propagation are produced in spontaneous down-conversion Teich  and in multimode transverse optical parametric oscillators Lugiato . Nevertheless, to our knowledge, there has been no experimental demonstration of quantum correlations with a multimode transverse light in the continuous wave regime.

Precision optical imaging using CCD cameras or photodetector arrays is required in many areas of science, ranging from astronomy to biology. Ultimately, the performance of optical imaging technology is limited by quantum mechanical effects. Of particular importance, as far as applications are concerned, is the measurement of image displacements, for example, the position of a laser beam. Techniques that rely on determining the position of a laser spot include atomic force microscopy AFM , measurement of very small absorption coefficients via the mirage effect Boccara  and observation of the motion of single molecules biologie . These measurements are usually performed as shown in Fig. 1.

The beam is incident on a split detector that delivers two currents proportional to the light intensity integrated over the two halves ($`x<0`$ and $`x>0`$) of the image plane. If the beam is initially centered on the detector, the mean value of the photocurrent difference is directly proportional to the relative displacement $`d`$ of the beam with respect to the detector. With a classical, shot noise limited laser source, the smallest displacement that can be measured (with a signal-to-noise ratio of one) is shown to be Fabre

$$d_{\mathrm{SQL}}=\frac{\sqrt{N}}{2I(0)}.$$
(1)
Here $`N`$ is the total number of photons recorded by the two detectors during the measurement time, and $`I(0)`$ is the local density of photons (photons per unit transverse length) at the position of the boundary between the two detectors. For a TEM00 Gaussian beam with radius $`w_0`$, the minimum measurable displacement is found to be

$$d_{\mathrm{SQL}}=\sqrt{\frac{\pi }{8}}\frac{w_0}{\sqrt{N}}.$$
(2)
For maximum focusing of the Gaussian beam, $`w_0=\lambda `$, and we obtain $`d_{\mathrm{SQL}}\lambda /\sqrt{N}`$, which is the absolute minimum displacement of a physical system that can be measured with classical beams AFM . Equation (2) shows that a more powerful laser, or a longer measurement time, gives increased measurement precision. However, in many applications these alternatives are simply not practical. In the case of atomic force microscopy, for example, excessive laser power ultimately leads to radiation pressure noise FM . For biological applications, large laser power may damage the samples under investigation and an increased integration time leads to loss of bandwidth. This is the motivation for looking for alternative methods of increasing measurement precision.

The limit of equation (2) can be surpassed only using multimode non-classical light. Let us consider a beam of light with an electric field distribution given by $`E(x)`$. We can build an orthonormal basis of the transverse plane $`\{u_i\}`$ such that $`u_0=E(x)/E(x)`$ is the first vector; $`u_1`$ is a “flipped” mode, given by $`u_0(x)`$ for $`x<0`$ and $`u_0(x)`$ for $`x>0`$ (see Fig. 2); and the other modes are choosen in order to form a basis. In this basis, the mean field of our light lies only in the first mode $`u_0`$ but, a priori, all of the modes contribute to the quantum noise. In order to determine the relevant modes of our measurement, we consider the interference quantities between two modes on each half of the split detector:

$`I_{x<0}(u_i,u_j)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^0}u_i^{}(x)u_j(x)𝑑x`$
$`I_{x>0}(u_i,u_j)`$ $`=`$ $`{\displaystyle _0^+\mathrm{}}u_i^{}(x)u_j(x)𝑑x.`$ (3)
Then the interference quantities relevant for a total measurement (sum of the two photodetectors) and a differential measurement (difference of the two photodetectors) can be written:

$`I_{\mathrm{sum}}(u_iu_j)`$ $`=`$ $`I_{x<0}(u_i,u_j)+I_{x>0}(u_i,u_j)`$
$`I_{\mathrm{diff}}(u_iu_j)`$ $`=`$ $`I_{x<0}(u_i,u_j)I_{x>0}(u_i,u_j)`$ (4)
One can then show that for any transverse mode $`u_i`$,

$`I_{\mathrm{sum}}(u_iu_1)`$ $`=`$ $`I_{\mathrm{diff}}(u_iu_0).`$ (5)
Since all $`u_i`$, for $`i2`$, are orthonormal to $`u_1`$ (i.e. $`I_{\mathrm{sum}}(u_iu_1)=0`$), equation (5) demonstrates that these modes have a zero overlap integral with $`u_0`$ in a differential measurement. It can then be shown that only $`u_1`$, which has a non-zero overlap integral with $`u_0`$, has to be considered along with $`u_0`$ in the noise calculation Fabre ; Delaubert .

We note that the modes $`u_0`$ and $`u_1`$ have perfect interference visibility as shown by their complete overlap integral for the differential measurement, ie. $`I_{\mathrm{diff}}(u_0u_1)=1`$. In this regard, the measurement is analogous to a perfect homodyne measurement with a beam splitter. The two modes are equivalent to the two input beams of a beamsplitter and the two halves of the multimode beam are equivalent to the two outputs. Therefore, similar to a homodyne measurement, the noise on the differential measurement is completely cancelled when the flipped mode is occupied by a perfect squeezed vacuum, with the squeezed quadrature in phase with the coherent field of the $`u_0(x)`$ mode. Conversely, the same result is also obtained when the mode profiles of the squeezed and the coherent fields are interchanged. In order to avoid the effect of losses, we have chosen a squeezed vacuum mode $`u_0(x)`$. We would like to stress that this simplified explanation can be applied only because we have conveniently identified the two relevant transverse modes of the measurement. However, contrary to a homodyne measurement, the entire measurement is performed using a single beam. Furthermore, a more general analysis is not limited to only two-mode beams.

The experimental setup is shown in Fig. 3. A stable Nd:YAG 700 mW laser provides a cw single mode beam at 1064 nm. A part of this beam is sent to a locked MgO:LiNbO3 frequency doubling cavity. The 532 nm output of the frequency doubler is used to pump a degenerate optical parametric amplifier (OPA) that produces a stable 10 $`\mu `$W squeezed beam in the TEM00 mode at 1064 nm. The noise reduction of the OPA output is measured to be 3.5 dB. Details of this squeezing system may be found in Buchler . The flipped mode, $`u_1(x)`$, is produced by sending the remaining part of the initial 1064 nm laser beam through a specially designed phase plate. This phase plate consists of two birefringent half-wave plates, one rotated by 90 with respect to the other, forming the two halves $`x<0`$ and $`x>0`$ of the transverse plane. These elements introduce a phase shift of 180 between the field amplitudes of the two halves. The squeezed output from the OPA is required to be superimposed onto the flipped mode with minimal loss. This is achieved by using a beam splitter that reflects 92% of the squeezed state and transmits 8% of the coherent state. The reflected output is then sent to a quadrant InGaAs detector (EPITAXX 505Q) with quantum efficiency greater than 90%. Only two of the four quadrants, of dimensions 500$`\mu `$m $`\times `$ 500$`\mu `$m each and with a dead zone between the pixels of $`25\mu `$m, are used in this experiment. A lens of focal length 30mm is used to image the phase plate on the detector plane and to counteract the diffraction of the flipped mode, which undergoes an abrupt phase change and therefore contains high spatial frequency components.

Fig. 4 shows the different noise levels monitored as a function of time when the relative phase between the coherent state and the squeezed state is chosen for maximum noise reduction. Due to the high stability of the various servo-loops in the experimental setup, the actively locked operation of the setup can be kept for hours. The noise measured on the sum of the two halves (Fig. 4a), i.e. on the total beam, coincides with the shot noise level for the conditions of this experiment, as expected from the coherent beam, which is not affected by the presence of a squeezed vacuum in an orthogonal mode. The noise measured on each individual half (Fig. 4c) is reduced by 1.08 $`\pm `$ 0.06 dB below the quantum noise limit. The fact that the intensity noise on each half of the beam is below the quantum noise limit, whereas the whole beam is at shot noise, shows the strong non-classical characteristic of this multimode beam. This is corroborated by the experimental data of Fig. 4b, which gives the noise on the intensity difference between the two halves at 2.34 $`\pm `$ 0.05 dB below the quantum limit. The results suggest that the beam is made of two strongly quantum correlated parts, indicating that a significant amount of spatial correlation has been created among the photons. With the measured noise reduction in the squeezed vacuum and a perfect setup (i.e. a perfect phase plate and a perfect mode-matching between the two transverse modes), one would expect 2.5 $`\pm `$ 0.2 dB of noise reduction on the difference between the two pixels. This demonstration is, to our knowledge, the first experiment in which spatial quantum effects have ever been observed in a bright beam of light.

This spatial noise correlation can now be used to improve the precision of displacement measurements in the image plane. For practical reasons, we have choosen to induce the displacement only in the coherent mode, before the mixing on the beamsplitter. However this displacement is of the order of the nanometer, which is several order of magnitude smaller that the relevant precision for the mode-matching of the two transverse modes, and the theoretical prediction for the measurements is the same as if the displacement was done on the total beam. In order to produce a small controllable beam displacement in the frequency range of the previous measurements, we use two electro-optic modulators (EOMs) driven at 4.5 MHz. Fig. 3 shows that EOM2 is slightly tilted with respect to the propagation of the light beam. When a voltage is applied across EOM2, a change in refractive index is induced and the transmitted beam experiences a parallel transverse displacement measured at about 3nm/V. We introduce a modulation at 4.5 MHz as signal for our displacement measurement which can be easily distinguished from the low frequency beam displacements induced by mechanical or acoustic vibrations. Apart from the parallel displacement, EOM2 will also introduce an unwanted phase modulation on the transmitted beam which is detrimental to our measurement. EOM1 of Fig. 3 is therefore used to compensate for this introduced phase modulation. When correct gains are chosen for both modulators, the transmitted beam will not have any phase or amplitude modulation and is only left with pure transverse displacement modulation. Figure 5 shows the differential signal monitored by a spectrum analyzer when the light beam undergoes a displacement modulation with an amplitude of 2.9Å. With a resolution bandwidth of 100 kHz, our setup recorded a modulation peak in the Fourier spectrum. Fig. 5a shows the trace when vacuum instead of the squeezed vacuum is used in mode $`u_0(x)`$. Thus this noise floor gives the standard quantum limit in such a displacement measurement. The signal-to-noise ratio (SNR) of this measurement is 0.68. When the two-mode non-classical beam is utilized in the measurement (Fig. 5b), we obtain a SNR of 1.20. This gives an improvement of the displacement measurement sensitivity by a factor of 1.7. The result is in agreement with the theoretical value calculated with the noise reduction reported in the previous paragraph. Similar measurements have been performed with a 10 kHz resolution bandwidth (and therefore a longer measurement time) and the results show the same improvement of the SNR.

Our results demonstrate that multimode non-classical states of light can be utilized to improve the optical measurement of small displacements. The noise floor of displacement measurements can actually be reduced to below the standard quantum limit. Of particular relevance are the potential usage of multimode squeezed light in atomic force microscopy and biological microscopy. Though our experimental demonstration is restricted to one-dimensional displacement measurements, it can be extended to two-dimensional displacement measurements with more complex forms of multimodal non-classical light.

###### Acknowledgements.

We would like to thank L. Lugiato and M. Kolobov for many enlightening discussions and CSIRO, Sydney, for the manufacturing of the special phase plate. This work is funded by the European project IST-2000-26019 ”Quantum images”, the Centre National de la Recherche Scientifique and the Australian Research Council.


text:
# Measurement of Γ(𝐾_𝑆→𝜋⁺𝜋⁻(𝛾))/Γ(𝐾_𝑆→𝜋⁰𝜋⁰)

## Abstract

We have measured the ratio $`R_S^\pi =\mathrm{\Gamma }(K_S\pi ^+\pi ^{}(\gamma ))/\mathrm{\Gamma }(K_S\pi ^0\pi ^0)`$ with the KLOE detector at the DAΦNE $`e^+e^{}`$ collider. This measurement is fully inclusive with respect to the $`\pi ^+\pi ^{}\gamma `$ final state. The sample of over $`10^6`$ two-pion decays of tagged $`K_S`$ mesons allows a statistical error as low as $`0.1\%`$ to be obtained. The accuracy is limited by systematic uncertainties, which are estimated primarily from data. We find $`R_S^\pi =2.236\pm 0.003_{\mathrm{stat}}\pm 0.015_{\mathrm{syst}}.`$

, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , The most recent measurement of $`R_S^\pi =\mathrm{\Gamma }(K_S\pi ^+\pi ^{})/\mathrm{\Gamma }(K_S\pi ^0\pi ^0)`$ dates back to 1976. The authors quote an accuracy of 4.3% . This and other results of similar accuracy determine the current world-average value of $`R_S^\pi `$, $`2.197\pm 0.026,`$ the fractional error on which is 1.2%. Moreover, the averaging of the existing results is somewhat questionable, since the various experiments have not clearly described their procedures for handling radiative $`\pi ^+\pi ^{}`$$`\gamma `$ decays. It is currently of interest to improve the accuracy on $`R_S^\pi `$, with particular attention to an accurate determination of what part of the radiative spectrum is included in the quoted result . The validity of the $`\mathrm{\Delta }I=1/2`$ rule, empirically established in $`K`$, $`\mathrm{\Lambda }`$, $`\mathrm{\Sigma }`$, and $`\mathrm{\Xi }`$ decays, has never been satisfactorily understood. This lack of understanding is closely connected with the difficulty in relating $`\mathrm{}(ϵ^{}/ϵ)`$ to the CKM matrix elements . The value of $`R_S^\pi `$ also enters into the theoretical calculation for the ratio of branching ratios for the decays $`K_L\mu ^+\mu ^{}`$ and $`K_L\pi ^+\pi ^{}`$ .

We describe the first measurement of $`R_S^\pi `$ to achieve an accuracy below 1%, which moreover fully includes $`\pi ^+\pi ^{}`$$`\gamma `$ radiative decays. This should allow the extraction of the isospin $`I=0`$ and $`I=2`$ decay amplitudes.

The data were collected with the KLOE  detector at DAΦNE, the Frascati $`\varphi `$-factory . DAΦNE is an $`e^+e^{}`$ collider which operates at a center of mass energy of $``$1020 MeV, the mass of the $`\varphi `$-meson. In the following, we use a coordinate system with the $`z`$ axis along the nominal beam line, the $`x`$ axis in the horizontal plane toward the machine center, the $`y`$ axis vertical, and the origin at the center of the collision region. Equal-energy positron and electron beams collide at an angle of 25 mrad, producing $`\varphi `$-mesons nearly at rest: $`p(\varphi )_{y,z}=0`$ and $`p(\varphi )_x=12.5\mathrm{MeV}/c.`$ $`\varphi `$-mesons decay $``$34% of the time into nearly collinear $`K^0`$$`\overline{K}^0`$ pairs. Because $`J^{PC}(\varphi )=1^{},`$ these pairs are in an antisymmetric state, so that the final state is always $`K_S`$$`K_L`$. Two photon intermediate states can lead to $`K_S`$$`K_S`$ and $`K_L`$$`K_L`$ final states, but this contamination is negligible . All of the above implies that the detection of a $`K_L`$ guarantees the presence of a $`K_S`$ of given momentum and direction. This fact can be exploited to identify $`K_S`$-mesons independent of their decay mode. We refer to this technique as $`K_S`$ tagging.

The KLOE detector consists of a large cylindrical drift chamber, DC, surrounded by a lead-scintillating fiber sampling calorimeter, EMC. A superconducting coil surrounding the calorimeter provides a 0.52 T magnetic field.

The drift chamber , 4 m in diameter and 3.3 m long, has 12 582 all-stereo tungsten sense wires and 37 746 aluminum field wires. The chamber shell is made of carbon fiber-epoxy composite and the gas used is a 90% helium, 10% isobutane mixture. These features maximize transparency to photons and reduce $`K_LK_S`$ regeneration and multiple scattering. The position resolutions are $`\sigma _{xy}150\mu \text{m}`$ and $`\sigma _z2\mathrm{mm}.`$ The momentum resolution is $`\sigma (p_{})/p_{}0.4\%.`$ Vertices are reconstructed with a spatial resolution of $`3\mathrm{mm}.`$

The calorimeter  is divided into a barrel and two endcaps containing a total of 88 modules, and covers 98% of the solid angle. The modules are read out at both ends by photomultipliers. The arrival times of particles and the three-dimensional positions of the energy deposits are determined from the signals at the two ends. The readout granularity is $`4.4\times 4.4\mathrm{cm}^2`$; the 2440 “cells” are arranged in layers five-deep. Cells close in time and space are grouped into a “calorimeter cluster.” For each cluster, the energy $`E_{\mathrm{CL}}`$ is the sum of the cell energies, and the time $`t_{\mathrm{CL}}`$ and position $`\stackrel{}{r}_{\mathrm{CL}}`$ are calculated as energy-weighted averages over the fired cells. The energy and time resolutions are $`\sigma _E/E=5.7\%/\sqrt{E(\mathrm{GeV})}`$ and $`\sigma _t=54\mathrm{ps}/\sqrt{E(\mathrm{GeV})}50\mathrm{ps},`$ respectively.

While the trigger  uses calorimeter and chamber information, only the calorimeter trigger, which requires two local energy deposits above threshold ($`50`$ MeV on the barrel and $`150`$ MeV on the endcaps), is used for the present measurement. Recognition and rejection of cosmic-ray events are also performed at the trigger level. Events with two energy deposits above a 30 MeV threshold in the outermost calorimeter plane are identified as cosmic-ray events and rejected. The trigger has a large time spread with respect to the beam crossing time. However, it is synchronized with the machine RF divided by 4, $`T_{\mathrm{sync}}=10.8\mathrm{ns},`$ with an accuracy of 50 ps. The time $`T_0`$ of the bunch crossing producing an event is determined after event reconstruction.

The mean decay lengths of the $`K_S`$ and $`K_L`$ are $`\lambda _S0.6\mathrm{cm}`$ and $`\lambda _L350\mathrm{cm}.`$ About 50% of $`K_L`$’s therefore reach the calorimeter before decaying. The $`K_L`$ interaction in the calorimeter ($`K_L`$-crash) is identified by requiring a cluster of energy above 200 MeV, which is not associated with any track, and whose time corresponds to a velocity $`\beta =r_{\mathrm{CL}}/ct_{\mathrm{CL}}`$ compatible with that of the $`K_L`$, $`\beta (K_L)0.214.`$ These $`K_L`$-crash events are used for $`K_S`$ tagging. The tagging efficiency depends slightly on the $`K_S`$ decay mode: the ratio $`R_{\mathrm{tag}}`$ of the tagging efficiencies for $`K_S\pi ^+\pi ^{}`$ and $`K_S\pi ^0\pi ^0`$ is determined from data.

The interaction time $`T_0`$, which must be known for the measurement of the cluster times, is obtained from the first particle reaching the calorimeter (pions or photons for the events of interest) assuming a velocity $`\beta =1.`$ In order to reduce the probability that $`T_0`$ is accidentally determined from a particle arising from machine background, we require the $`T_0`$ to be fixed by a cluster with energy $`E_{\mathrm{CL}}>50`$ MeV and distance to the beam line $`\rho _{\mathrm{CL}}>60\mathrm{cm}.`$ Charged pions from $`K_S\pi ^+\pi ^{}`$ arrive at the EMC $`3\mathrm{ns}`$ later than $`\gamma `$’s from $`\pi ^0`$ decays. The computed time $`T_0`$ is therefore delayed by one RF period, $`2.7\mathrm{ns}.`$ This results in a mismeasurement of the velocity for long-lived kaons accompanied by $`K_S\pi ^+\pi ^{}`$ decays. The distribution of the observed $`K_L`$ velocity transformed to the $`\varphi `$-rest frame, $`\beta ^{},`$ is shown in Fig. 1. The effect of the incorrect $`T_0`$ assignment for the case of $`K_S\pi ^+\pi ^{}`$ is evident.

Events are selected by requiring $`0.195\beta ^{}0.2475.`$ This window in $`\beta ^{}`$ encompasses both peaks and gives a tagging efficiency of $`20\%,`$ which differs slightly for each of the two $`K_S`$ decay modes. The ratio of tagging efficiencies $`R_{\mathrm{tag}}`$ is evaluated from a subsample of $`K_S\pi ^+\pi ^{}`$ events in which both pion tracks are associated to a calorimeter cluster. In this case, the event $`T_0`$ can be correctly evaluated making use of the measured pion track lengths and velocities, resulting in the $`\beta ^{}`$ spectrum shown in the inset of Fig. 1, which is very similar to that for $`K_S\pi ^0\pi ^0`$ events. This corrected $`\beta ^{}`$ distribution and that obtained from the standard $`T_0`$ estimate are normalized to the same number of $`K_L`$ interactions. $`R_{\mathrm{tag}}`$ is computed from the ratio of the total number of events inside the $`\beta ^{}`$ window in each of the two cases. The effects of residual differences between the $`\beta ^{}`$ distributions for $`K_S\pi ^0\pi ^0`$ and $`T_0`$-corrected $`K_S\pi ^+\pi ^{}`$ events are included in the determination of the systematic error.

Other sources of bias in the evaluation of $`R_{\mathrm{tag}}`$ have been taken into account as well. The probability that the event contains a cluster which determines $`T_0`$ is slightly different for $`\pi ^+\pi ^{}`$ and $`\pi ^0\pi ^0`$ decays. These probabilities are estimated using $`K_S`$ events tagged by $`K_L\pi ^+\pi ^{}\pi ^0`$ decays. For the identification of such events, knowledge of $`T_0`$ is not necessary. These events are identified by searching for two tracks forming a vertex in the drift chamber volume. The missing momentum of the $`\pi ^0`$ is determined by kinematic closure at the $`K_L`$ vertex. Photons from the $`\pi ^0`$ are identified using the difference between their times of flight and by requiring that the $`\gamma \gamma `$ pair be back-to-back in the $`\pi ^0`$ rest frame.

Accidental overlap of a $`\pi ^\pm `$ track with the $`K_L`$ cluster produces a loss of $`K_S\pi ^+\pi ^{}`$ events. The magnitude of this effect has been estimated using $`K_L`$-crash events selected without requiring that the cluster not be associated with any track.

Finally, a small fraction of $`\varphi K_SK_L`$ events is rejected by the cosmic-ray veto, since some $`K_L`$-crash interactions can deposit energy in the outermost EMC plane. $`R_{\mathrm{tag}}`$ is therefore corrected by the ratio of the inefficiencies induced by the cosmic-ray veto for $`K_S\pi ^+\pi ^{}`$ and $`K_S\pi ^0\pi ^0`$ events. This ratio is obtained using a subsample of events for which the cosmic-ray veto was recorded but not enforced.

All of the above effects result in corrections to $`R_{\mathrm{tag}}`$, which are listed in Tab. 1. The largest uncertainty is due to the correction for the $`\beta ^{}`$ cut.

$`K_S`$$``$$`\pi ^+\pi ^{}`$($`\gamma `$) events are selected by requiring the presence of two tracks of opposite charge that intersect a cylinder 4 cm in radius and 10 cm in length along the $`z`$-axis, centered at the interaction region (IR). Pion momenta and polar angles must satisfy the fiducial cuts $`120p300\mathrm{MeV}/c`$ and $`30^{}\theta 150^{}.`$ Pions must also reach the EMC, as opposed to spiralling in the drift chamber. The requirement that both tracks reach the EMC results in an acceptance which depends on the energy of the photon in radiative decays. Fig. 2 shows the dependence of the efficiency $`ϵ_{\pi \pi \gamma }`$ on the $`\gamma `$ energy in the $`K_S`$ rest frame $`E_\gamma ^{}`$, obtained by Monte Carlo (MC) simulation for $`K_S`$$``$$`\pi ^+\pi ^{}`$$`\gamma `$ events with $`E_\gamma ^{}>20`$ MeV and $`K_S\pi ^+\pi ^{}`$ events without radiation ($`E_\gamma ^{}=0`$). The efficiency for $`0<E_\gamma ^{}<20\mathrm{MeV}`$ is obtained by linear interpolation. The overall acceptance is obtained by folding the above efficiency with the photon spectrum from Ref. . The result differs from that computed assuming no radiation by $`0.3\%.`$

The single-track reconstruction efficiency is obtained in bins of $`(p,\theta )`$ using $`K_S\pi ^+\pi ^{}`$ events identified by a single reconstructed pion track. This selection is applied to both data and MC. The ratio of data and MC reconstruction efficiencies is found to be constant ($`0.990\pm 0.001`$) for pions with momenta and polar angles the satisfy the fiducial cuts used in the analysis. The overall MC acceptance is therefore scaled accordingly and used to correct the data. Systematic errors due to possible misalignments of the DC are evaluated by studying the dependence of the acceptance on the angular cut applied and are found to be negligible. The acceptance is listed in Tab. 1.

The trigger efficiency is also obtained from data. $`K_L`$ interactions always fire at least one trigger sector. About 40$`\%`$ of the time, however, the $`K_L`$-crash fires two nearby sectors, thus producing a valid trigger by itself. We use these events to find the probability that at least one $`K_S`$ decay product complements the $`K_L`$-crash cluster to satisfy the trigger condition in the remaining 60% of the events. The trigger efficiency is given in Tab. 1.

$`K_S\pi ^0\pi ^0`$ events are identified by the prompt photon clusters from $`\pi ^0`$ decays. A prompt photon cluster must satisfy $`|t_{\mathrm{CL}}r_{\mathrm{CL}}/c|5\sigma _t,`$ where $`\sigma _t`$ is the energy-dependent time resolution, and must not be associable to any track. Machine background is reduced by cuts on energy and polar angle; we require $`E_{\mathrm{CL}}>20\mathrm{MeV}`$ and $`\mathrm{cos}\theta <0.9.`$ Three or more prompt photons are required to accept a $`K_S\pi ^0\pi ^0`$ event.

The photon-detection efficiency in bins of energy and polar angle has been obtained from data using a sample of $`\varphi \pi ^+\pi ^{}\pi ^0`$ events. The acceptance for $`K_S\pi ^0\pi ^0`$ events is obtained from MC simulation and is listed in Tab. 1. The simulation incorporates the experimentally determined photon-detection efficiencies; we rely on it essentially for the averaging over event geometries and kinematics. The 0.2$`\%`$ error on the $`K_S\pi ^0\pi ^0`$ selection efficiency is dominated by the uncertainty in the photon-detection efficiency below 50 MeV due to non-linearity in the energy response of the calorimeter and/or the photon reconstruction algorithm. The trigger efficiency, also listed in Tab. 1, is evaluated by the same method used for $`K_S\pi ^+\pi ^{}`$ events.

Background in the sample of $`K_S\pi ^0\pi ^0`$ events is mainly due to $`K_S\pi ^+\pi ^{}`$ events in which the track-to-cluster association failed, and in which there were residual split and accidental clusters. This contamination induces an error on the counting of events with three or more prompt clusters that is lower than 0.1%, as estimated directly from data.

The present measurement of $`R_S^\pi `$ is based on an integrated luminosity of $`17\mathrm{pb}^1`$ collected during the year 2000. The data set contains $`N^\pm =\mathrm{1\hspace{0.17em}060\hspace{0.17em}821}`$ selected $`K_S\pi ^+\pi ^{}`$ events and $`N^{00}=\mathrm{766\hspace{0.17em}308}`$ selected $`K_S\pi ^0\pi ^0`$ events. The value for $`R_S^\pi `$ is obtained by correcting the ratio of these numbers for the ratio of tagging and selection (including trigger) efficiencies:

$$R_S^\pi =\frac{N^\pm }{N^{00}}\times \frac{1}{R_{\mathrm{tag}}}\times \frac{ϵ_{\mathrm{sele}}^{00}}{ϵ_{\mathrm{sele}}^\pm }$$

The result is

$$2.236\pm 0.003_{\mathrm{stat}}\pm 0.015_{\mathrm{syst}}$$

The contributions to the estimated errors are listed in Tab. 2. Additional details are available in Ref. . Our result is consistent with but slightly higher than the current world-average value $`R_S^\pi =2.197\pm 0.026`$ . This might be expected, since our measurement was designed to include the entire radiative spectrum.

We thank the DAΦNE team for their efforts in maintaining low background running conditions and their collaboration during all data-taking. We also thank F. Fortugno for his efforts in ensuring good operations of the KLOE computing facilities. This work was supported in part by DOE grant DE-FG-02-97ER41027; by EURODAPHNE contract FMRX-CT98-0169; by the German Federal Ministry of Education and Research (BMBF) contract 06-KA-957; by Graduiertenkolleg ’H.E. Phys. and Part. Astrophys.’ of Deutsche Forschungsgemeinschaft contract No. GK 742; by INTAS contracts 96-624 and 99-37; and by TARI contract HPRI-CT-1999-00088.


text:
# Focusing Revisited: an MN-dynamics Approach

## Abstract

The nonlinear Schrödinger (NLS) equation is a ubiquitous example of an envelope wave equation for conservative, dispersive systems. We revisit here the problem of self-similar focusing of waves in the case of the focusing NLS equation through the prism of a dynamic renormalization technique (MN dynamics) that factors out self-similarity and yields a bifurcation view of the onset of focusing. As a result, identifying the focusing self-similar solution becomes a steady state problem. The discretized steady states are subsequently obtained and their linear stability is numerically examined. The calculations are performed in the setting of variable index of refraction, in which the onset of focusing appears as a supercritical bifurcation of a novel type of mixed Hamiltonian-dissipative dynamical system (reminiscent, to some extent, of a pitchfork bifurcation).

Self-focusing and wave collapse have received a considerable amount of attention in many diverse branches of physics, ranging from optics  to hydrodynamics  and from gravity  to plasma physics . The interested reader can also consult .

The most frequently studied model equation in the context of self-focusing and the ensuing blow-up of the relevant solutions is the two-dimensional (2D) and three-dimensional (3D) nonlinear Schrödinger (NLS) equation. This benchmark example is well-known (see e.g.,  and references therein) to exhibit unstable solitary wave solutions which collapse in finite time. Numerous publications have addressed the focusing properties of the solutions (Refs.  can only claim to give a number of representative examples of this large volume of work). Recently, a revived interest in the subject has led to numerous new insights .

Here we take a slightly different approach to the same problem, viewing it from the perspective of dynamical systems theory and making the connections to both earlier work on the NLS problem  as well as more recent work on self-similar blowup problems . Many of the results that are presented here on NLS were previously known from the works of . Such results include the effective factoring out of the self-similarity to view focusing as a steady state problem, as well as the continuation of the branch of focusing solutions. What we believe is new in our exposition is the intuition it provides by identifying the onset of focusing as a steady state bifurcation of a (transformed) partial differential-algebraic equation (PDAE). This appears to us to be a new type of bifurcation reminiscent of a supercritical pitchfork. The identification of focusing solutions as steady states of this system also allows us to linearize around them and examine (as well as justify) their stability and the ensuing structure of the bifurcation diagram. Here these computer-assisted tasks are performed for a finite discretization of the problem.

Finally, it is worth noting that this approach has been implemented in a somewhat different but equivalent form of the regular blow-up NLS problem. In particular, in previous studies see e.g., , the radial form of the NLS equation was analyzed,

$`R_{rr}+{\displaystyle \frac{d1}{r}}R_rR+|R|^2R=0,`$ (1)
where $`R=R(r)`$ denotes the radial part of the solution. The analysis (when performed parametrically) would then involve a slightly unphysical “bifurcation parameter”, namely the spatial dimensionality $`d`$. In contrast to that, but equivalently for the purposes of the analysis, we use a different (but more physically realistic) bifurcation parameter. The latter is the power of the nonlinearity $`\sigma `$ in the model:

$`iu_t=\mathrm{\Delta }u|u|^{2\sigma }u;`$ (2)
$`u`$ is the complex envelope of the electric field (in optical applications) and the present setting arises when the index of refraction has a power law dependence ($`\sigma `$ being the relevant power) on the intensity of the field. It is well known that the model of Eq. (2) becomes critical (i.e., unstable to wave collapse) for $`\sigma 2/d`$, where $`d`$ is the spatial dimensionality of the Laplacian in (2) (see e.g., ). In view of that, we choose to study the model in one spatial dimension, using $`\sigma `$ as a bifurcation parameter, with $`\sigma _{cr}=2`$ in that case. Notice that for $`\sigma =1`$, the model becomes the celebrated integrable example of the 1D cubic NLS .

When considering self-focusing solutions of Eq. (2) it is customary to make a straightforward appropriate scaling. In particular, using $`u=L^{1/\sigma }v(\xi ,\tau )`$, with $`\xi =x/L`$, $`\tau =_0^tL^2(t)𝑑t`$, one obtains

$`iv_\tau =v_{\xi \xi }|v|^{2\sigma }via(\tau )({\displaystyle \frac{1}{\sigma }}v+\xi v_\xi ),`$ (3)
where $`a(\tau )=L_\tau /L`$ (see e.g., ). A direct pointer to the criticality of $`\sigma =2`$ is the fact that the norm is preserved under the rescaling only for that value of $`\sigma `$. Hence, dynamical (i.e., permitted by the time evolution) rescaling for the original NLS model of Eq. (2) is only permissible for $`\sigma =2`$. On the other hand, Eq. (3) has appeared in various forms in a number of works on the focusing problem such as  (see  for a review of the relevant methods and results), in the context of what was called “dynamic rescaling”. The idea behind such a rescaling is that the scaling is chosen in such a way that the rescaled time goes to infinity as the (finite) time of the singularity $`t^{}`$ is approached. The closure of the system (with an equation for $`a(\tau )`$) was given by imposing an arbitrary, but physically or numerically motivated, constraint, such as the preservation of a certain norm of the solution.

Here, we will follow a slightly different path, motivated from a general approach to self-similar problems recently presented in . In particular, if we consider the equation

$`iu_t=F(u),`$ (4)
typically the potential for self-similar solutions arises when a scaling is inherent in the nonlinear operator of the right hand side (RHS) of (4). In general:

$`F\left(Bu({\displaystyle \frac{x}{L}})\right)=L^cB^dF\left(u(\xi )\right),`$ (5)
where $`\xi =x/L`$ (keeping in line with the notation used in earlier works such as ). The dispersive operators with on-site nonlinearity considered here are however different from the ones of  in that the presence of the onsite potential imposes a certain restriction on the scaling through Eq. (5), namely for the $`F`$ of the RHS of 4)

$`B`$ $`=`$ $`L^{\frac{1}{\sigma }}`$ (6)
$`F\left(Bu({\displaystyle \frac{x}{L}})\right)`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}F\left(u(\xi )\right).`$ (7)
Using in the equation the ansatz $`u=B(\tau )v(x/L(\tau ),\tau )`$, where $`\tau `$ is an, as of yet, undetermined function of the original time $`t`$, we obtain:

$`i\left({\displaystyle \frac{B_\tau }{B}}v{\displaystyle \frac{L_\tau }{L}}\xi v_\xi +v_\tau \right)\tau _t={\displaystyle \frac{1}{L^2}}\left(v_{\xi \xi }+|v|^{2\sigma }v\right),`$ (8)
where the subscripts denote respective partial derivatives. Notice also that using Eq. (6), $`B_\tau /B=L_\tau /(\sigma L)`$ and hence the self-similarity has a single scaling factor, $`L`$. Time parametrization can be chosen on the basis of convenience and hence we use here $`\tau `$ such that $`\tau _t=1/L^2`$. If we look, as is customary (see e.g.,  and references therein) for a standing wave solution of Eq. (8), we use the so-called phase invariance to look for solutions of the form $`v\mathrm{exp}(i\mathrm{\Lambda }\tau )v`$, which gives the same equation as before but with an extra $`\mathrm{\Lambda }v/L^2`$ in the RHS. Without loss of generality, we set $`\mathrm{\Lambda }=1`$. Finally, one has the one-parameter family ($`L`$) freedom of the group orbit of self-similarity (see  for the case of traveling solutions; the case of self-similarity is treated in ). To determine the scaling factor, we lift the symmetry-induced degeneracy through a pinning condition  of the form

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}Re(v(\xi ,\tau ))T(\xi )𝑑\xi =C,`$ (9)
where $`C`$ is a constant. $`T(\xi )`$ is an arbitrary, so-called “template” function. Alternatively, one can consider maximazing the inner product with the template (such an approach has been considered in  for settings with travelling wave solutions, as well as ones with self-similar solutions). Differentiating Eq. (9) and using eq. (8), we arrive to the following equations for the real and imaginary parts of the solution $`v=U+iW`$

$`U_\tau `$ $`=`$ $`W_{\xi \xi }(U^2+W^2)^\sigma W+WG({\displaystyle \frac{1}{\sigma }}U+\xi U_\xi )`$ (10)
$`W_\tau `$ $`=`$ $`U_{\xi \xi }+(U^2+W^2)^\sigma UUG({\displaystyle \frac{1}{\sigma }}W+\xi W_\xi ),`$ (11)
where

$`G={\displaystyle \frac{L_\tau }{L}}={\displaystyle \frac{_{\mathrm{}}^{\mathrm{}}\left[W_{\xi \xi }+(U^2+W^2)^\sigma WW\right]T(\xi )𝑑\xi }{\frac{C}{\sigma }+_{\mathrm{}}^{\mathrm{}}\xi U_\xi T(\xi )𝑑\xi }}.`$ (12)
Eqs. (10)-(11) supplemented with Eq. (12) constitute the MN-dynamics formulation of the focusing problem for the NLS equation.

The regular solitary wave solutions of the equation (2) exist as stationary solutions of the equation ensuing from (2) upon the substitution $`u=\mathrm{exp}(it)v`$, or equivalently as “trivial” focusing solutions of the MN equations (10)-(11) with a zero blowup rate, $`G=0`$. These solutions are stable for $`\sigma <2`$. Their mechanism of instability involves a pair of linearization eigenvalues that bifurcate from the band edge of the continuous spectrum, for $`\sigma >1`$. The dispersion relation describing the continuous spectrum is $`\lambda =\pm i(1+k^2)`$, where $`\lambda `$ is the eigenvalue corresponding to wavenumber $`k`$, for our case of $`\mathrm{\Lambda }=1`$. This pair of eigenvalues moves towards the origin of the spectral plane, where 2 more pairs of eigenvalues lie due to the 2 invariances (translational and phase) of Eq. (2). At $`\sigma =2`$, the additional pair “arrives” at the origin, thereafter exiting along the real axis, rendering the soliton branch unstable as is shown in Fig. 1.

Beyond the instability threshold, the soliton is well-known to be unstable towards focusing, resulting in a self-similar blowup for $`\sigma >2`$. Such solutions are most commonly obtained through integration of the dynamic rescaling equations (or even possibly of the original equations in an adaptively refined mesh), see e.g.,  and references therein. Here we use an alternative perspective***Our work has a somewhat similar flavor to Phys. Rev. A 38, 3837 (1988), but is posed in a different, we believe slightly more general framework that allows us to perform a systematic bifurcation analysis of the results.. We identify the soliton branch as a special solution (with $`G=0`$) of Eqs. (10)-(11). We then use a “soliton-motivated” pinning condition. In particular, we use as our arbitrary template a Dirac mass, centered at a given point $`x_0`$, forcing the real part of the solution to have at $`x_0`$ the same value as the soliton does for $`\sigma =2`$. Other pinnings (for the imaginary part, or for combination of the real and the imaginary part) have also been used with the same results. The pinning condition in conjunction with Eqs. (10)-(11) allows us to construct for $`\sigma >2`$ the branch of self-focusing solutions as steady states through the use of Newton iteration and pseudo-arclength continuation .

The algebraic problem is solved using the boundary value solver of the package gPROMSsee e.g., www.psenterprise.com, a commercial simulation tool for solving systems of ordinary, and/or partial differential equations and/or of partial differential-algebraic equations of index 1. The pseudo-arclength continuation method was incorporated within gPROMS. The computational domain was chosen to be the right half-axis $`\xi =0`$ (due to symmetry of the solutions sought). The length and the tessellation of the domain were chosen appropriately to ensure reliability of the solutions, i.e., robustness against further discretization refinement and domain enlargement. The solutions were approximated using a centered finite difference method of order two in a domain of length equal to 14 and a total of 700 discretization intervals, while the absolute accuracy was set to $`1E06`$. The solutions we obtained are shown in Fig. 2.

An additional advantage of effectively “factoring out” the self-similarity through probing the solutions as steady states of the MN problem, stems from our ability to perform linear stability analysis computations for these solutions. It should be noted that linear stability within a renormalization group approach to self-similar solutions was also considered in gravitational problems, see e.g., . Such numerical computations for the truncated problem are shown in Fig. 2. These computations are for $`G>0`$, and clearly indicate that the dynamical system of Eqs. (10)-(11) is a dissipative (a non-Hamiltonian) one. This is also structurally clear in the equations, as terms of the form $`iv`$ are well-known dissipative perturbations of NLS. Analyzing the structure of the “umbrella-shaped” spectrum in comparison with the linearization spectrum of the Hamiltonian problem for $`G=0`$, we recognize that the genuinely complex eigenvalues constitute the continuous spectrum. In addition, we (typically) obtain $`3`$ eigenvalues along the real line. In the Hamiltonian untransformed problem, there are 6 eigenvalues that are near the origin for $`\sigma `$ close to 2, 4 of which are at the origin (two with a corresponding symmetric, zero node eigenvector, related to phase invariance and two with an antisymmetric, one node eigenvector, corresponding to translational invariance). The additional pair that eventually leads to the instability has a corresponding symmetric eigenvector with two nodes. In the dissipative MN system, we are left with only 3 eigenvalues, since 2 out of the 6 are absent due to considering only half the domain (where only symmetric eigenfunctions will be present). An additional eigenmode is absent due to imposing the pinning condition (decreasing by 1 the degrees of freedom of the system). One of the three eigenvalues is at the origin, as the MN dynamics also possess a conservation law, namely Eq. (9). The phase invariance is, however, explicitly broken here due to the presence of the term proportional to $`G`$. In fact, one can show that unless one is “exactly on” a steady state of Eqs. (10)-(11), neither the invariance $`vv\mathrm{exp}(i\varphi )`$ (for $`\varphi `$ an arbitrary constant), nor the corresponding conservation law $`du_2/dt=0`$ (where $`u_2`$ is the $`L^2`$ norm of the solution) are respected, for the pinning of Eq. (9). In fact, it can be straighforwardly predicted within the realm of (leading order) perturbation theory that the bifurcation of the phase eigenvalue for $`\sigma 2^+`$ can be approximated by

$`\lambda G{\displaystyle \frac{\sigma 2}{2\sigma }}.`$ (13)
This shows that, as is indeed verified in numerical computations, the bifurcation of the phase eigenvalue is in the real positive semi-axis of the spectral plane, for $`\sigma >2`$. However, such an instability (which is a genuine one for the non-phase invariant MN dynamics) will not be relevant for the original norm preserving (Hamiltonian) dynamics and hence the self-focusing solution will be a stable blowup solution of Eq. (2). I.e., the presence of a conservation law in the original NLS model does not permit the dynamical development of this instability, or, equivalently, does not permit the system to “explore” this unstable eigendirection.

As $`\sigma `$ is increased from the critical value, the generation of the blowup branches can be seen to occur in a form reminiscent of a supercritical pitchfork-like bifurcation in Fig. 1. Notice the steepness of the variation of $`G`$ in the neighborhood of the critical point. This very abrupt increase of the value of $`G`$, as well as earlier asymptotic arguments seem to support a dependence beyond all algebraic orders of $`G`$ on $`\sigma 2`$ (i.e., an exponential dependence). However, the rigorous proof of such an estimate is beyond the scope of the present work and will be addressed in a future publication. However, it should be highlighted that this is a rather unusual, and most probably novel form of a dynamical system, possessing solutions with as well as without focusing, the corresponding dynamics however being significantly modified (from dissipative to Hamiltonian) between the two cases. Identifying the normal form behavior of such dynamical systems close to criticality (close to the onset of focusing) is an interesting open problem. Furthermore, it should be noted that a natural setting for the study of the bifurcations of such dynamical systems appears to be the framework of the relevant differential-algebraic equations (due to the presence of the constraint).

Notice that for $`\sigma 2^+`$, $`G0`$ very abruptly. For $`\sigma =2`$, there is no solution with $`G0`$ . Dynamically, this is manifested as $`lim_\tau \mathrm{}G(\tau )=0`$, the approach to the limit being logarithmic . It would be interesting to explore whether travelling problems can exhibit similar critical behavior.

In summary, in this work we have presented a different, bifurcation-based perspective to the focusing problem in the ubiquitous NLS equation. We have “factored out” the self-similarity group action through the MN dynamics approach and set the problem up as a steady state problem in the latter framework. This allowed us to identify stable self-focusing solutions for the original problem (which actually are unstable for MN dynamics). We were able continue the relevant branches, constructing the supercritical pitchfork-like bifurcation diagram that constitutes the focusing instability for this novel (mixed Hamiltonian-dissipative) dynamical system. We were able to identify the eigenvalues of the focusing steady states of the MN dynamics and to explain their spectrum on the basis of the knowledge of the “soliton branch” spectrum. The theoretical study of the normal form of the MN-dynamics that could provide analytical justification of the scaling behavior near the transition is currently under study and will be reported in future publications.

We thank D. Aronson, S. Betelu, K. Lust, B. Malomed, D. McLaughlin and C. Rowley for stimulating discussions. The support of of AFOSR and NSF (IGK) and of the University of Massachusetts (PGK) is also gratefully acknowledged.