On ultra-high energy cosmic ray acceleration at the termination shock of young pulsar winds

# On ultra-high energy cosmic ray acceleration at the termination shock of young pulsar winds

## Abstract

Pulsar wind nebulae (PWNe) are outstanding accelerators in Nature, in the sense that they accelerate electrons up to the radiation reaction limit. Motivated by this observation, this paper examines the possibility that young pulsar wind nebulae can accelerate ions to ultra-high energies at the termination shock of the pulsar wind. We consider here powerful PWNe, fed by pulsars born with millisecond periods. Assuming that such pulsars exist, at least during a few years after the birth of the neutron star, and that they inject ions into the wind, we find that protons could be accelerated up to energies of the order of the Greisen-Zatsepin-Kuzmin cut-off, for a fiducial rotation period msec and a pulsar magnetic field G, implying a fiducial wind luminosity erg/s and a spin-down time s. The main limiting factor is set by synchrotron losses in the nebula and by the size of the termination shock; ions with may therefore be accelerated to even higher energies. We derive an associated neutrino flux produced by interactions in the source region. For a proton-dominated composition, our maximum flux lies slightly below the 5-year sensitivity of IceCube-86 and above the 3-year sensitivity of the projected Askaryan Radio Array. It might thus become detectable in the next decade, depending on the exact level of contribution of these millisecond pulsar wind nebulae to the ultra-high energy cosmic ray flux.

a]Martin Lemoine a]Kumiko Kotera b]Jérôme Pétri \affiliation[a]Institut d’Astrophysique de Paris,
CNRS, Université Pierre & Marie Curie,
98 bis boulevard Arago, F-75014 Paris, France \affiliation[b]Observatoire Astronomique de Strasbourg,
Université de Strasbourg, CNRS, UMR 7550,
11 rue de l’Université, F-67000 Strasbourg, France \emailAddlemoine@iap.fr \emailAddkotera@iap.fr \emailAddjerome.petri@astro.unistra.fr \keywordsCosmic rays – Pulsar wind nebulae

## 1 Introduction

The origin of the highest energy cosmic rays is a long-standing enigma of astroparticle physics, which has withstood some fifty years of intense experimental activity (see reviews by, e.g., [1, 2]). The existing data have brought in very significant results, such as the detection of a cut-off at the expected location of the Greisen-Zatsepin-Kuzmin (GZK) suppression [3, 4]. Such a spectral feature, combined with the absence of striking anisotropy in the arrival direction of the highest energy particles, could indicate that they originate from extragalactic sources. So far, however, no conclusive experimental evidence points towards one or the other of the many possible scenarios of ultra-high energy cosmic ray (UHECR) origin.

The central question in this field of research is how to accelerate particles to these extreme energies eV. Among the known particle acceleration scenarios, Fermi-type shock acceleration plays a special role. It is rather ubiquitous, since collisionless shock waves emerge as direct consequences of powerful outflows. Furthermore, shock acceleration, when operative, is known to dissipate into the supra-thermal particle population a substantial fraction of the kinetic energy that is inflowing into the shock, of the order of , see e.g. [5]. Shock acceleration also produces rather generically a spectrum with nearly constant energy per decade, which allows to transfer a sizable fraction of the energy in the ultra-high energy domain. Those are noticeable features in the context of the origin of ultra-high energy cosmic rays, because one indeed needs to extract a large fraction of the source energy in order to match the cosmic ray flux above eV, e.g. erg per transient source of ocurrence rate Mpcyr [6]1.

Particle acceleration to ultra-high energies in collisionless shock waves has been proposed in a number of scenarios [7], e.g. in gamma-ray bursts [8, 9, 10, 11, 12], in blazars [13, 12] or in radio-galaxy jets [13]. The possibility of accelerating ultra-high energy cosmic rays at the termination shock of pulsar winds has received so far little attention, except for Ref. [14], which has stressed the large energy gain associated to the first shock crossings. The present paper thus proposes a critical discussion of this issue. Let us recall here that the ultra-relativistic collisionless shock front separates the (inner) cold fast magnetized pulsar wind from the pulsar wind nebula, which itself is bounded by a shock propagating into the supernova remnant; this nebula thus contains the hot shocked wind material and hot shocked supernova remnant material, e.g. [15, 16, 17].

The main motivation of the present study comes from the realization that the Crab nebula represents so far the most efficient particle accelerator known to us, since the observation of a synchrotron spectrum extending up to ( the fine structure constant) attests of the capacity of the termination shock to accelerate electrons and positrons up to the radiation reaction limit at the Bohm rate, meaning an acceleration timescale in terms of the gyro-time , with , expressed here in the comoving blast frame, e.g. [18]. Furthermore, it is generally admitted that the highest energy pairs have been shock accelerated through a Fermi process, because their spectral index is remarkably similar to the predictions of relativistic shock acceleration for isotropic scattering, see e.g.  [19, 20, 21, 22, 23]. In this sense, it is natural to try to extend this result to the acceleration of ions. Of course, pulsar winds are usually modeled as pair winds, hence one central assumption of the present work is that such winds may also inject ions, see the discussion in Ref. [24], and see also  [25, 26, 27] which propose to interpret the morphological features of the Crab Nebula through the coupling between ion and pair dynamics at the termination shock.

The confinement energy of cosmic rays in Crab-like nebulae is nevertheless quite low, being of the order of eV; here, and indicate respectively the magnetic field and the size of the blast. We will thus be interested in more powerful pulsar wind nebulae, able to confine particles up to higher energies. As we show in the following, young pulsars born with periods s do fulfill this criterion; moreover, their huge rotational energy reservoir erg is also highly beneficial for producing a substantial flux of UHECRs [28, 29, 30, 31].

The release of this tremendous rotational energy into the surrounding supernova ejecta should modify its radiative properties. In particular, such objects could lead to ultra-luminous supernovae lasting for months to years, with distinctive bright gamma-ray and X-ray counterparts [32, 33, 34, 35, 36]. These scenarios could provide an explanation to some of the observed ultra-luminous supernovae [37], that would then constitute an indirect probe of the existence of pulsars born with millisecond periods. Interestingly, it has been suggested that the Crab pulsar itself was born with a msec period on the basis of the observation of the large number of radio-emitting pairs in the nebula [38].

Engine-driven supernovae, or trans-relativistic supernovae, whose characteristic high speed ejecta is believed to be powered by some internal source such as a magnetized neutron star, have also been considered as potential sources of ultra-high energy cosmic rays [39, 40, 41]. However, in those scenarios, acceleration is argued to take place in the outer fastest parts of the mildly relativistic external shock which propagates in the wind of the progenitor star. Our present discussion proposes another view of these objects, in which particle acceleration to ultra-high energies takes place well inside the remnant, at the ultra-relativistic shock that is running up the pulsar wind.

Finally, young pulsar winds themselves have also been considered as potential sites of UHECR acceleration, mainly through the electric field associated to the rotating magnetic dipole, e.g.  [42, 28, 29, 43, 24, 30, 31]. We stress that the present scenario is wholly different in terms of acceleration physics. In particular, pulsar magnetospheric and wind physics do not play any role in our scenario, beyond controlling the spin-down time of the engine, while it directly sets the maximum energy that might be reached in those wind acceleration models. In the present work, particle acceleration takes place in two steps: in a first step, the initial Poynting flux of the pulsar wind is assumed to be dissipated down to near equipartition with ions and electrons before or around the termination shock; this is a generic assumption, motivated by observational results on known pulsar wind nebula, e.g.  [44, 45]; the injected high energy ions are then accelerated at the termination shock, up to a maximal energy determined by energy losses and escape considerations.

The lay-out of this paper is as follows. In Sec. 2, we discuss the dynamics and the radiative properties of the nebula, which limit the acceleration through radiative losses and escape, then we discuss the maximal energy as a function of the various parameters. In Sec. 3, we discuss the neutrino signal associated to the acceleration of ultra-high energy cosmic rays as well as the dependence of our results on the model assumptions and parameters. We provide a summary of our results and conclusions in Sec. 4.

## 2 Particle acceleration in young powerful PWNe

The physics of young PWNe is controlled by the amount of energy injected by the pulsar wind, of luminosity , and the velocity-density profile of the supernova ejecta which confines the nebula. Detailed morphological studies of the Crab nebula, along with numerical MHD simulations [46, 47] indicate that the geometry is mostly axisymmetric. It is however possible to reproduce the main characteristics of a young pulsar wind nebula with a spherically symmetric picture, in particular the location of the termination shock radius and the size of the nebula, e.g., [42, 44]. Given the other astrophysical uncertainties described below, this suffices for our purposes and, in the following, we assume spherical symmetry.

Throughout this study, neutron stars have a moment of inertia with fiducial value gcm, instantaneous and initial rotation velocities and (corresponding initial period ), radius and dipole magnetic field ; these should not be confused with magnetic field strengths and spatial scales of the nebula.

The pulsar rotational energy reservoir amounts to . The wind luminosity decreases as

 Lw(t)=Lp/(1+t/tsd)(n+1)/(n−1) , (1)

in terms of the braking index (defined by ) and spin-down time , with initial luminosity erg/s. For magneto-dipole losses in the vacuum, , while observations rather indicate . Nevertheless, we will be interested in the structure of the nebula at time , at which a substantial fraction of the rotational energy has been output into the nebula; the braking index controls the later evolutionary stages, therefore it will not impact significantly our results. We thus adopt for simplicity in what follows.

The spin-down timescale is then given by

 tsd≃9I⋆c38B2⋆R6⋆Ω2i∼3.1×107sI⋆,45B−2⋆,13R−6⋆,6P2i,−3 . (2)

For convenience, we indicate our results in terms of , and , as well as in terms of and .

We concentrate mainly on pulsars with magnetic fields G and not on magnetars (with G). As Eq. (2) indicates, magnetars spin down on a timescale much shorter than a year, and at the early times when the highest energy particles are accelerated, the density of the surrounding supernova ejecta does not allow their escape [30]. A magnetar scenario thus require the disruption of the supernova envelope by the wind [24] or that particles escape through a region punctured by a jet, like in a strongly magnetized proto-magnetar scenario discussed by [48]. Gravitational wave losses are negligible for pulsars with G, hence we neglect them, see e.g.  [24].

### 2.1 General input from the Crab nebula

The extrapolation of the phenomenology of the Crab pulsar wind nebula to the young PWNe that we are interested in is by no means trivial, because it involves an increase by some six orders of magnitude in luminosity, and it is hampered by three unsolved issues: the so-called problem, the physics of particle acceleration at the termination shock of pulsar winds and the origin of radio emitting electrons in PWNe. Let us recall briefly these issues in order to motivate our model of the nebula.

The magnetization parameter relates the Poynting flux to the matter energy flux; in the comoving wind frame, for a cold plasma of rest mass energy density , it is defined by . For a mixed pair-ion composition of respective densities and , , defining the multiplicity factor for pairs achieved through pair cascade in the magnetosphere.

Observationally, the problem results from the difficulty in reconciling the large value of the magnetization parameter at the pulsar light cylinder with that inferred downstream of the termination shock () through a leptonic model of the emission [42, 44, 45, 18, 49, 46]. A generic estimate for in the Crab nebula is , assuming that pairs are injected at the Goldreich-Julian rate (recalled further below), with a low energy, at the base of the wind, with multiplicity . In contrast, models of pulsar wind nebulae and their comparison to observations rather suggest  [44, 45, 50, 51, 52, 53], although more recent three-dimensional MHD simulations suggest that values as large as could reproduce the morphological data for the Crab nebula [47, 54]. From a more theoretical perspective, this problem characterizes the difficulty of pushing cold MHD winds to large Lorentz factors through the Poynting flux, e.g. [55, 56, 57, 58, 59]. Indeed, in a radially expanding MHD wind, should remain constant from close to the light-cylinder up to the termination shock.

How particle acceleration takes place in the Crab nebula is another puzzle. Although the termination shock offers an obvious site for particle acceleration, and the high energy spectral index of the reconstructed electron distribution conforms well to the expectations of a relativistic Fermi process with isotropic scattering [19, 20, 21, 22, 23], so far our understanding of particle acceleration rather suggests that the Fermi process should be inefficient in mildly-magnetized – namely for a magnetization – ultra-relativistic shocks, see  [60, 61] for an analytical discussion and  [5] for simulations. To summarize such discussions briefly, Fermi acceleration can take place at ideal – meaning planar and steady – relativistic shock waves only if intense small-scale turbulence has been excited in the shock vicinity2 [62]. Such small-scale turbulence may in principle be excited by streaming instabilities between the supra-thermal particles and the background unshocked plasma in the shock precursor. At the termination shock of pulsar winds, this may come through a current-driven instability if  [63], or through the synchrotron maser instability, if  [25, 26]. Nevertheless, the finite magnetization, even if is as small as , should prevent acceleration to very high energies, because the efficiency of scattering in small-scale turbulence relatively to the gyration in the background field decreases in inverse proportion to the particle energy [64, 60]. As discussed in these references, this implies a maximum energy beyond which scattering (hence acceleration) becomes ineffective. This maximum energy has been observed in recent PIC simulations  [5]. Its exact value is not of importance for the present discussion; it suffices to say that it scales as so that, at mildly magnetized shock waves with , Fermi acceleration should not be able to accelerate the particles to the energies observed, in ideal conditions.

Efficient dissipation of the magnetic field in the nebula, is probably the key to resolving both the problem and the issue of particle acceleration. Regarding the former, recent 3D MHD simulations, which account for dissipative effects inside the nebula, alleviate somewhat the question of conversion of Poynting flux in the wind by allowing for values  [47, 54]. Although the exact mechanism of dissipation remains open to debate, reconnection in the current sheet separating the “stripes” of opposite magnetic polarity (the “striped wind” [65, 66, 55]), upstream of the termination shock, has long been discussed as a possible way of converting part of the Poynting flux, e.g. [55, 67]. Dissipation around the termination shock may also support efficient particle acceleration in two ways: by seeding large scale turbulence, causally disconnected from the upstream magnetic topology, in which case particle acceleration could proceed unimpeded [62]; or, by providing an extra mechanism of particle acceleration, which would feed into the Fermi process at higher energies, e.g. through the dissipation of MHD waves [68], or reconnection [69, 70, 71]. Furthermore, the MHD simulations of Ref. [72] reveal that the termination shock is unsteady and corrugated, leading to the excitation of mildly relativistic turbulence immediately downstream. Finally, we also note that the high energy spectral index is typical of a relativistic Fermi process in a mildly or weakly magnetized shock, whereas magnetized shock waves lead to a weaker compression ratio, hence a softer spectrum [20]. This, again, argues in the favor of substantial dissipation around the termination shock of the wind.

One can argue further in the favor of dissipation inside the pulsar wind nebula, as follows: ad absurdum, one could not construct a stationary model with a super-fast magnetosonic wind in the absence of dissipation [44], as the shock crossing conditions at the termination shock would then lead to a ultra-relativistic bulk velocity for the shocked wind material; however, any unsteady solution is bound to populate the nebula with magnetized turbulence, with a fast magnetosonic velocity close to , which would lead to particle acceleration on a fast timescale, hence to efficient dissipation.

To summarize, both observational and theoretical arguments indicate that efficient dissipation of the magnetic field and the tapping of this energy into kinetic energy are relevant processes at the pulsar wind termination shock.

The last issue concerns the origin of radio-emitting electrons in the Crab nebula, which are about more numerous than the optical to -ray emitting electrons [58]. In terms of multiplicity, the radio emission suggests a pair multiplicity , well above the theoretical expectations \citepHibschman01,Timokhin10, which are in much better agreement with the multiplicity associated to higher energy electrons. Two generic, diverging interpretations are usually given: either the multiplicity indeed reaches values , in which case the pulsar spin-down power divided by the total kinetic energy leads to a rather low value for the Lorentz factor of the wind at the termination shock [69], or and, for one interpretation, the radio emitting low energy electrons were injected at an earlier stage of the nebula [45, 18, 38]. Interestingly, the latter interpretation requires that the Crab pulsar was born with a period of order msec, i.e. quite close to the range of values that we are interested in [38].

Let us note that the present discussion implicitly assumes values , as if were as high as , the ions could carry only a tiny fraction of the wind energy and it would become difficult to match the cosmic-ray flux at the highest energies; this issue is discussed in Sec. 3.

In spite of the above unknowns, the fact is that the Crab accelerates electron-positron pairs efficiently, up to the radiation reaction limit, with a high energy index very similar to that expected in a relativistic Fermi process, and this remains our main motivation to discuss the possibility of pushing ions to ultra-high energies. In what follows, we build a one-zone model of the synchrotron nebula in order to quantify the various losses that limit the acceleration of particles to ultra-high energies.

Following the above discussion, we assume that the wind energy is efficiently dissipated into random particle energies inside the nebula, i.e. ; hereafter, is understood as the average magnetization parameter inside the nebula, after dissipation has taken place. This dissipation can either take place upstream of the termination shock, in which case the shock itself transforms the ordered kinetic energy into random particle energies; or, it can take place at or downstream of the termination shock, through one of the various processes discussed above. We also assume that the termination shock is strong, which implies super-fast magnetosonic velocities of the wind; this, however, is not a stringent requirement on wind physics, since it only requires a wind Lorentz factor at the termination shock for our fiducial parameters, see the definition of in Eq. (13) below, and corresponds to the initial magnetization of the wind [58].

### 2.2 The millisecond nebula structure

The structure of the nebula can be approximated by analytical solutions at times , when is approximately constant [75]. The pulsar wind nebula radius can then be written

 RPWN = (12599β3PWNc3L0M0)1/5t6/5(tc≳t) (3) = (815L0M0)1/2t3/2(tsd≳t≳tc) , (4)

in terms of the time at which the external shock of the PWN has swept up the mass of the supernova ejecta, , which we assume of constant density (with a rapidly declining density profile beyond) and in terms of the wind power ; represents the velocity of the pulsar wind nebula in the source rest frame. For parameters of interest, one finds that , since for a core mass and an ejecta velocity km/s [75, 35], hence we will use mostly Eq. (4) in the following.

At times , the pulsar input into the nebula decreases rapidly. For this phase, we then assume that the blast evolves in free expansion, meaning  [35]

 RPWN(t)=RPWM(tsd)ttsd (5)

Of course, we are mostly interested in the PWN structure at the time , since it corresponds to the time of maximum energy injection into the nebula. The temporal scalings of at times and thus do not play a crucial role in the forthcoming analysis, but they help in understanding how the various quantities evolve in time.

The above estimates neglect the interaction of the blast with the outer shocked region of the supernova, i.e. the forward and reverse shocks associated with the interaction with the circumstellar medium. In more typical, less powerful (erg/s) pulsar wind nebulae, the interaction with the reverse shock takes place on timescales of thousands of years and this leads to the compression of the nebula  [50]. In the present case, the interaction takes place shortly after , i.e. shortly after the PWN external shock into the ejecta has swept up the inner core of the remnant. However, the pulsar energy output ergs dominates the kinetic energy of the outer blast (ergs), therefore the nebula will dominate the dynamics. We neglect this interaction phase; again, it should not modify appreciably the values of that we derive at time , which is our prime objective here.

The above analytical solutions also fail when radiative cooling of the blast becomes important. As we show in the following, the latter possiblity is to be considered, because the electrons cool through synchrotron faster than a dynamical timescale, contrary to what happens in PWNe such as the Crab. Therefore, if dissipation of the Poynting flux is efficient, and if ions represent a modest part of the energy budget, most of the wind luminosity input into the nebula is actually lost into radiation. In order to account for this effect, we use an improved version of Eqs. (3),(4), in which represents the actual power deposited into the nebula, representing the fraction of luminosity converted into radiation through pair cooling (: fraction of energy in the magnetic field, : fraction of energy in ions, in the nebula). These approximations are used to provide analytical estimates of the various quantities characterizing the nebula at time .

We complement these estimates with a detailed numerical integration of the following system:

 ˙RPWN = βPWNc (6) ˙Res = βesc (7) ˙Usw = (βw−βts)Lw−Pem−4πR2PWNβPWNcpPWN, (8) ˙Mse = (βes−βej)4πr2esρejc, (9) ˙Use = (βes−βej)4πr2esρejc3+4πR2PWNβPWNcpPWN. (10)

All quantities are defined in the source rest frame; they are as follows: corresponds to the radius of the contact discontinuity, interpreted as the size of the nebula; represents the location of the outer shock of the nebula, propagating in the supernova remnant; to a very good approximation, (thin-shell approximation); consequently represents the velocity of this outer shock while denotes the velocity of the termination shock; represents the energy contained in the shocked wind region, beneath the contact discontinuity; represents the power lost through radiation; since the electrons cool faster than an expansion timescale (see below), one can write , in terms of the fraction of power injected into the nebula in magnetic field () and ions (). represents the pressure inside the nebula, which can be well approximated by  [51]; the term associated to consequently represents adiabatic losses for ; denotes the mass accumulated in the shocked ejecta region, between the contact discontinuity and the outer shock; corresponds to the velocity of the supernova remnant ejecta, in the source frame; finally, denotes the energy contained in the shocked ejecta region. To a good approximation, since the ejecta is non-relativistic.

These equations can be obtained by integrating the equations of particle current density and energy-momentum conservation over the spatial variables, between the boundaries of interest. This procedure introduces the brackets for and for , which correspond to the fact that the boundaries of the shocked wind and shocked ejecta are delimited by the moving shock waves. The velocity of the termination shock in the source frame depends non-trivially on the degree of magnetization of the shock; for , however,  [44] and , therefore we approximate . For the outer shock, assuming it is strong, non-radiative and non-relativistic, one has . This closes the system.

In order to evaluate the dynamics of the nebula, we assume that the supernova ejecta consists of a core mass of constant density. An analytical estimate of the size of the pulsar wind nebula is then given by Eqs. (4) and (5):

 RPWN ≃ 4.1×1016L1/2p,45t3/2sd,7.5ˇt3/2^tcm (11) ≃ 3.2×1016P−3B−2⋆,13R−6⋆,6I3/2⋆,45ˇt3/2^tcm .

The quantities and indicate the scaling of these values at times respectively short and long of , obtained respectively through Eq. (4) and Eq. (5). Figure 1 presents the evolution in time of pulsar wind dynamical quantities ( and ) calculated analytically and by numerical integration. The prefactors match the numerical evaluation shown in Fig. 1 for . The scaling departs slightly from the (resp. ) behaviour at short (resp. late) times compared to , but we neglect this difference in the following. Radiative nebulae, in which is closer to unity, tend to be more compact; this difference can be read off Fig. 1 and inserted in the relations that follow. For the sake of simplicity, we assume an adiabatic case in the following, i.e. .

At this stage, it may be useful to make contact with known pulsar wind nebulae, such as the Crab: its radius is about pc, i.e. about a hundred times larger than the above. In the Crab nebula, the radius of the termination shock is estimated to be pc [44], while in the present case, the termination shock is located close to the contact discontinuity, see Sec. 2.4.1, making the above nebula not only more compact, but also much thinner.

The mean magnetic field in the nebula is then obtained as follows. Recall that corresponds to the magnetic fraction of the energy actually injected into the nebula after a proper account of dissipation, i.e. ; one thus has

 BPWN = ⎛⎝6ηB∫t0Lw(t′)dt′R3PWN⎞⎠1/2 (12) ≃ 14η1/2B,−1L−1/4p,45t−7/4sd,7.5ˇt−7/4^t−3/2G ≃ 12P−5/2−3η1/2B,−1I−7/4⋆,45B3⋆,13R9⋆,6ˇt−7/4^t−3/2G .

The numerical values are obtained by plugging into the first equation the temporal scalings of and obtained previously. Analytical and numerical estimates agree for the adiabatic case at the spin-down time, see Fig. 1. The magnetic field strength is of course much larger than that seen in more standard PWNe [50, 76, 51], as a result of the larger input energy and of the younger age, which implies a more compact nebula.

The right panel of Fig. 1 depicts the evolution of the distribution of the energy fractions , and , discussed below.

### 2.3 Energy injected into particles in the nebula

As argued in Section 2.1, we assume efficient dissipation of the initial Poynting flux, into random particle energy. In the absence of energy losses, this conversion implies that particles (electron-positron pairs or ions) acquire a typical Lorentz factor

 γdiss. ≃ 11+σPWNLw˙Nmc2 (13) ≃ 2.2×1091−ηB1+xiκ−14L1/2p,45^t−1 ≃ 1.8×1091−ηB1+xiκ−14P−2−3B⋆,13R3⋆,6^t−1 .

This Lorentz factor can also be written as: in terms of , the magnetization of the flow short of the termination shock. The particle rest mass power injected into the nebula is written here: , with the Goldreich-Julian rate [77]; is the ratio of the power injected into ions, relatively to that injected into pairs: if ions of charge are injected at a rate , , so that for .

Up to radiation reaction effects, heating through dissipation proceeds equally for pairs and ions, meaning that both acquire a same Lorentz factor. This is guaranteed for all dissipation processes mentioned earlier. Among others, this implies that, notwithstanding radiation reaction effects, the ratio of the energy injected into pairs to that injected into ions is conserved in the dissipative processes.

At the present time, one cannot predict what kind of ions the pulsar would output. We therefore remain general and consider ions of mass number , charge . Note that the composition of the highest energy cosmic rays is not very well-known either: while the Pierre Auger Collaboration reports a light composition at eV, transiting to an intermediate composition at higher energies [78], the results of the Telescope Array experiment seemingly point towards a light composition [79], even though the depths of shower maximum of both experiments appear compatible, see  [80].

Actually, energy losses may limit the typical Lorentz factor of electron-positron pairs to the minimum of and the radiation reaction limiting Lorentz factor . Of course, cannot be lower than the Lorentz factor of the wind at the termination shock, , which is unknown. We assume that . This is not a strong assumption, since the latter two Lorentz factors are quite large.

Equating the synchrotron cooling time with the gyration time of the particle, the radiation reaction limiting Lorentz factor is given by the usual formula

 γe−loss = 32mec2e3/2B1/2PWN (14) ≃ 3×107η−1/4B,−1L1/8p,45t7/8sd,7.5ˇt7/8^t9/8 ≃ 3×107η−1/4B,−1P5/4−3I7/8⋆,45B−3/2⋆,13R−9/2⋆,6ˇt7/8^t9/8 .

The synchrotron power of ions of Lorentz factor , atomic number and mass number scales as , therefore the radiation reaction Lorentz factor for ions is a factor larger than . For ions, radiation reaction therefore does not limit the efficiency of dissipation.

Accounting for dissipation and radiation reaction limitations, one can thus write the fractions of energy and carried by the electrons and the ions inside the nebula as:

 ηe=1−ηB1+ximin(1,γe−lossγdiss.) ,ηi=(1−ηB)xi1+xi . (15)

The quantity is understood as characterizing the energy injected in pairs in the nebula, after dissipation/acceleration processes, but before synchrotron cooling has taken place. Clearly, for the above fiducial parameters, at time , meaning that most of the energy dissipated into the electrons has been radiated at the radiation reaction limit, producing photons of energy MeV. This radiation does not contribute to the radiation losses of ultra-high energy ions but it may lead to a specific signature of dissipation processes in such young PWNe. In contrast, in the Crab nebula, so that the electrons can take away most of the dissipated energy without losing it to radiation.

In such compact PWNe, the electrons cool through synchrotron radiation on a timescale that is much shorter than a dynamical timescale, down to non-relativistic velocities, since the cooling Lorentz factor is given by

 γc ≃ 6πmec2βPWNσTB2PWNRPWN (16) ∼ t2sd,7.5βPWNη−1B,−1ˇt2^t2 ∼ η−1B,−1P4−3βPWNB−4⋆,13R−12⋆,6I245ˇt2^t2 .

This represents a major difference with respect to the case of the Crab nebula, for which , so that most electrons do not cool on an expansion timescale, due to the smaller amount of energy injected into the wind and to the larger size of the nebula.

The above allows us to characterize the spectral energy distribution (SED) of the nebula; in particular, the low-frequency spectral luminosity is represented by

 Lν,syn≃ηeLw(ϵϵe)1/2(ϵc<ϵ<ϵe) (17)

with and in terms of the synchrotron peak frequencies associated to Lorentz factors and . In Fig. 1, we plot the time evolution of the various quantities that characterize the nebula, i.e. the mean magnetic field, the mean nebular radius and the fractions of energy , and .

### 2.4 Ion acceleration

After their injection through the termination shock, the ions are energized through dissipative processes up to , then shock accelerated to the maximal Lorentz factor that we seek to determine here. This latter Lorentz factor is given, as usual, by the competition between shock acceleration, escape from the PWN and energy losses in the synchrotron nebula.

We assume here that acceleration proceeds at the Bohm rate, an assumption that is motivated and supported by the two following remarks. The first is of a more empirical nature as it follows from the observation that the Crab nebula does accelerate electron-positron pairs up to the radiation reaction limit. If the acceleration timescale is written , then the comparison of with the synchrotron loss timescales leads to an upper bound on the maximum photon energy, . The fact that the synchrotron spectrum extends up to MeV or so in the Crab nebula indicates that .

The second line of argument originates from our theoretical understanding of particle acceleration at relativistic shock waves; to put it briefly, if some large scale turbulence, seeded downstream of the termination shock by dissipative processes, mediates the scattering process at the termination shock. In order to see this, one must recall that supra-thermal particles probe a short length scale of order behind a relativistic oblique shock, before returning to the shock or being advected away [62]. Therefore, if some turbulence is transmitted from upstream to downstream through the shock, the shock crossing conditions imply the continuity of the magnetic field lines through the shock, which then prevent repeated Fermi cycles, unless most of the turbulent power lie on short length scales [81, 62, 64, 60]. If, however, the turbulence is seeded downstream of the shock, which requires additional dissipative processes as in the present case, then this continuity is broken, hence Fermi cycles can take place, as modelled in test-particle Monte Carlo simulations, e.g. [19, 20, 21, 22]. The downstream and upstream residence timescales are then both of order in the shock rest frame, so that the acceleration timescale corresponds to .

#### Confinement at the shock and in the nebula

Confinement in the nebula itself leads to a maximal Lorentz factor:

 γconf ∼ ZieBPWNRPWNmic2 (18) ≃ 1.5×1011ZiAiη1/2B,−1L1/4p,45t−1/4sd,7.5ˇt−1/4^t−1/2 ≃ 1.4×1011ZiAiη1/2B,−1P−3/2−3I−1/4⋆,45B⋆,13R3⋆,6ˇt−1/4^t−1/2 . (19)

The large value of confirms that young fast pulsars input enough power into the nebula to confine particles up to ultra-high energies, a non-trivial result in itself.

As a matter of fact, the finite size of the termination shock limits the maximum acceleration energy to a factor below the above confinement energy, because the shock radius with . This value is derived from the numerical simulations discussed in Sec. 2.2, but it can be understood as follows. In the source rest frame, one can write the velocity of the termination shock to first order in as follows,

 βts≃1−3β2−3+β2 (20)

where denotes the velocity of the shocked wind, immediately downstream of the termination shock. This equation assumes hydrodynamic jump conditions at the shock to simplify the analysis. It is essentially a rewriting of the shock velocity in a frame in which the downstream is moving at velocity ; the expression in the downstream frame is obtained by , which implies , as expected for a strong ultra-relativistic hydrodynamic shock wave [82]. One commonly assumes that the flow velocity then evolves as a function of radius according to , see  [44, 45] for a detailed discussion; this implies that the blast velocity . Since , one can solve as a function of , using the above in conjunction with Eq. (20). One finds

 rts≃[√3βb+O(β2PWN)]RPWN . (21)

Therefore, for a typical velocity at (a value checked in numerical calculations).

At a gyroradius a factor of order unity to a few above , the particle feels the shock curvature and it gradually decouples from the flow as its size exceeds the size of the accelerator. Acceleration therefore stops at Lorentz factor

 γmax≃αtsγconf (22)

#### Synchrotron losses

Due to the strong magnetic field in the nebula, synchrotron losses represent a potential limitation to the maximum energy. In the course of acceleration, synchrotron losses limit to

 γi,syn−acc ≃ 6.2×1010η−1/4B,−1AiZ−3/2iL1/845t7/8sd,7.5ˇt7/8^t3/4 (23) ≃ 6.0×1010η−1/4B,−1AiZ−3/2iP5/4−3B−3/2⋆,13R−9/2⋆,6I7/8⋆,45ˇt7/8^t3/4

As discussed in Sec. 2.3, this maximum Lorentz factor is simply times . It lies a factor of a few below the confinement energy, but it nevertheless allows protons to be accelerated up to energies of the order of the GZK cut-off for our fiducial parameters.

Given that , particles also suffer synchrotron energy losses during their escape out of the nebula. The corresponding cooling Lorentz factor is obtained by matching and the particle energy, which leads to

 γi,syn−esc ≃ 2.6×1010η−1B,−1A3iZ−4it2sd,7.5ˇt2^t2 (24) ≃ 2.4×1010η−1B,−1A3iZ−4iP4−3B−4⋆,13R−12⋆,6I2⋆,45ˇt2^t2

As a result of its combination of powers of , and , the above equation can actually be rewritten as , with a numerical prefactor entirely controlled by the core mass of the supernova ejecta and by fundamental constants. However, we are interested in computing the limiting Lorentz factor at a time , at which most of the rotational energy has been output by the pulsar; for this reason, one should actually read Eq. (24) above with , which de facto introduces a dependence on .

This limit appears more severe than the previous ones, but with a rather strong dependance on the spin-down time , or alternatively, on the dipole moment of the neutron star, see Eq. (2). Therefore, a spin-down time larger by a factor two, or a (modest) factor decrease in would push this maximum energy [as well as Eq. (23)] at or above eV for protons, while the maximal energy associated to the finite size of the termination shock at would remain of the order to eV. Note that the total injected energy does not depend on . Alternatively, if one assumes , a magnetic field G would still guarantee that the limiting synchrotron loss energy lies above the GZK energy (for protons), just as the maximum energy associated to the finite size of the termination shock.

Of course, synchrotron energy losses are much weaker for ions with ; the corresponding maximal energy is larger than that of protons by a factor .

#### Photohadronic losses

Assuming that the injected ions are protons, the impact of photopion interactions can be evaluated in the standard approximation, according to which interactions take place with photons of energy , with the proton Lorentz factor and cross-section mb. It is straightforward to check that all protons interact with photons in the low-frequency part of the synchrotron SED, which justifies our neglect of the high-frequency part of the synchrotron SED. The number density of such low-frequency synchrotron photons is3

 ϵγdnγdϵγ≃ηe(ϵγϵe)1/2Lw4πR2PWNcϵγ(ϵc<ϵγ<ϵe) (25)

so that the pion production optical depth, reads:

 τγπ ≃ 7.9×10−3ηeγ1/2p,11L1/2p,45t−3/2sd,7.5ˇt−3/2^t−3 (26) ≃ 6.6×10−3ηeγ1/2p,11P−5−3B4⋆,13R12⋆,6I−3/2⋆,45ˇt−3/2^t−3 ,

with . Note that depends on time, in particular at early times, see Eq. (15) and Fig. 1, while at late times.

The ratio of the acceleration timescale to pion production timescale,

 tacctγπ ≃ γpγconfτγπ , (27)

indicates that pion production is not a limiting factor for acceleration to the highest energies.

The left panel of figure 2 presents the time evolution of the proton confinement Lorentz factor , dissipation Lorentz factors , and limiting Lorentz factor from synchrotron losses and pion production interactions , for initial periods ms.

The evolution over time of the pion production optical depth and the ratio of the acceleration to pion production timescales, calculated at the at energy are shown in the right panel of Fig. 2.

The same type of calculations can be performed for heavier nuclei, considering the Giant Dipole Resonance (GDR) as the main channel for energy losses on the background photons. The parameters for such interactions read: cm for the cross-section, and [83]. The results for iron nuclei are presented in Fig. 3.

All in all, Equations (22), (23), (24) and (27) and the accompanying discussion thus indicate that proton (and heavier ion) acceleration to energies of the order of or above the GZK cut-off appears possible in PWNe with parameters close to those chosen in this paper, namely , ms, .

## 3 Discussion

Equations 19, (23), (24), 26 and 27 summarize the confinement and acceleration properties of young powerful pulsar wind nebulae at time , when most of the rotational energy is output into the nebula. A third constraint on potential ultra-high energy cosmic ray sources comes from the energy output into cosmic rays above eV. For a proton dominated composition, this energy output must match ergs/Mpc/yr once it is folded over the population of sources with ocurrence rate  [6].

Over , the pulsar injects into ions

 Ei=ηiLwtsd=(1−ηB)xi(1+xi)Lptsd . (28)

Assuming that these msec PWNe constitute a fraction of the total core collapse supernovae, with ocurrence rate Mpcyr, the normalization to the flux of ultra-high energy cosmic rays thus implies

 ηSN∼3×10−5qix−1iL−1p,45t−1sd,7.5 , (29)

assuming here , i.e. , and . Recall that , indicating that these young msec pulsars should constitute a fraction of the supernova rate. If one rather wishes to normalize the flux at an energy eV with , then this fraction would go up to . The prefactor for an injection spectral index , accounts for the difference in normalization induced by the lower cosmic-ray injection energy limit.

### 3.1 Neutrino signal for a proton-dominated composition

The detection of neutrinos associated with hadronic and photo-hadronic interactions of nuclei in the nebula would provide an unambiguous test of the present scenario. Let us first consider the yield of neutrinos through photo-hadronic interactions on the nebula SED. Since the neutrino yield for heavy nuclei is much smaller than that for protons, we assume in this section that . The neutrino spectrum is then shaped by the accelerated proton spectrum and by the conversion efficiency. A parent proton produces 3 neutrinos per conversion, which takes place with probability in each photopion interaction, thus on a timescale . At the detector, the neutrino carries a fraction of the parent proton energy, i.e. , with the redshift of the source.

The time-dependent neutrino spectrum emitted by one PWN at luminosity distance can be written:

 E2νΦν(Eν,t)=(1+z)4πD2LE2ν3tγπ+(Ei)dEidEνdNidEi , (30)

in terms of the parent proton spectrum in the source. The latter is formally given by

 dNi(t)dEi=ηi∫t0q