text:
# Mission: Impossible (Escape from the Lyman Limit) Based on observations taken with the NASA/ESA Hubble Space Telescope, which is operated by AURA under NASA contract NAS5-26555

## 1 Introduction

The diffuse ultraviolet background is one of the key ingredients in the recipe leading to galaxy formation (Miralda-Escudé & Ostriker 1990; Giroux & Shapiro 1996). The intensity of the diffuse ultraviolet background has been estimated at moderate and high redshifts using the proximity effect on QSO absorbers (Liske & Williger 2001; Scott et al. 2001; Fernández-Soto et al. 1995; Kulkarni & Fall 1993), and results obtained so far suggest that it is too large to be explained by QSOs alone (Madau & Shull 1996, but see also e.g. Giallongo et al. 1996). Although systematic errors may make these measurements uncertain (Loeb & Eisenstein 1995, Pascarelle et al. 2001), several suggestions have been put forward to explain the observed extra intensity. One of the most popular suggestions calls for the escape of a large fraction of ionising photons from high-redshift galaxies. If correct, this would imply that the interstellar media of high-redshift galaxies must be at least partially transparent to ionising photons.

But the best measurements available at low and moderate redshifts show that galaxies are (at least nearly) opaque to ionising photons. The large column densities of neutral hydrogen that surround most galaxies should be enough, depending on how they are distributed, to very effectively quench any outgoing flux of ionising photons. In fact, Leitherer et al. (1995) find, in a sample of four starburst galaxies at an average redshift $`z0.02`$, that the limits to the escape fraction of Lyman limit photons range from $`f_{\mathrm{esc}}<0.0095`$ to $`f_{\mathrm{esc}}<0.15`$ (with an observed flux ratio at 900 Å and 1500 Å of $`F_{9/15}<0.10-0.20`$). Hurwitz et al. (1997) analyse the same data and obtain slightly less stringent limits, ranging from $`f_{\mathrm{esc}}<0.032`$ to $`f_{\mathrm{esc}}<0.57`$. Deharveng et al. (2001) find a more stringent limit $`f_{\mathrm{esc}}<0.064`$ for a starburst galaxy of redshift $`z=0.0448`$. At higher redshifts ($`z1`$), Ferguson (2001) finds similar upper limits $`f_{\mathrm{esc}}<0.20`$ to the escape fraction.

It is much less clear whether galaxies at still higher redshifts are as opaque to ionising photons as galaxies at low and moderate redshifts. Accurate measurements become progressively more difficult at higher redshifts, because galaxies are much fainter and the presence of Lyman $`\alpha `$ forest absorption introduces further uncertainties. The only two measurements available so far are from Steidel, Pettini & Adelberger (2001) and Giallongo et al. (2002). Adopting a correction factor for intervening Lyman $`\alpha `$ absorption determined from a composite QSO spectrum, Steidel et al. reported $`F_{9/15}=0.22\pm 0.05`$ using 29 galaxies at $`z=3.40`$, while Giallongo et al. reported a $`1\sigma `$ upper limit $`F_{9/15}0.05`$ using two galaxies at $`z3`$. The former measurement indicates that galaxies at $`z3`$ are much more transparant to ionising photons, contributing an equal amount of ionising flux to the ultraviolet background radiation as the QSOs. The latter measurement indicates otherwise.

In this article, we present a new measurement of $`f_{\mathrm{esc}}`$ using 27 galaxies at redshifts $`1.9<z<3.5`$ in the Hubble Deep Field (HDF). All of these galaxies have secure spectroscopic redshifts (see Cohen et al. 2000–C00–and Dawson et al. 2001–D01) and accurate broad-band photometry from HST/WFPC2 observations in the F300W, F450W, F606W, and F814W bandpasses. Our measurement is based on deep space-based images, which provides two important advantages over the methods used by previous groups: First, our sky subtraction is more accurate, because two-dimensional sky subtraction in faint-galaxy images is more robust than one-dimensional sky subtraction in faint-galaxy spectra. Second, the wavelength interval below the rest-frame Lyman limit over which we integrate is larger, because the space-based broad-band images are sensitive to observed-frame wavelengths as short as $`\lambda 2800`$ Å, allowing us to integrate over rest-frame wavelength intervals as large as $`100`$ Å in some cases. Our method also improves upon previous work by accounting for the effect of the Lyman $`\alpha `$ forest absorption by performing a large numer of random realizations of a parameterized distribution, to account for variations of the absorbing clouds across different lines of sight.

Our results show that the escape fraction estimated from the ensemble of 27 galaxies has a (statistical) $`3\sigma `$ upper limit of $`f_{\mathrm{esc}}0.04`$. Given that our galaxy sample spans a range of properties in redshift, intrinsic colours, and luminosities, we also explore possible systematic effects by dividing the galaxies into different subsamples. We do not find any evidence to support that bluer galaxies have higher escape fractions of ionizing photons. Instead, we find that more luminous galaxies appear to have shallower ultraviolet spectral slopes and higher escape fractions (at the 2$`\sigma `$ level of significance). We do also search for other possible systematic effects that could be inherent to our method, and find that other effects cannot in any realistic way account for a large increase in the escape fraction.

## 2 Data

We consider all galaxies in the HDF with known spectroscopic redshifts in the range $`1.9<z<3.5`$. The lower limit is chosen so that the rest-frame Lyman limit begins to affect the photometry in the F300W bandpass. The upper limit is chosen so that we can still reasonably distinguish absorption due to the Lyman $`\alpha `$ forest from the intrinsic Lyman limit absorption. We previously measured high-quality photometry of these galaxies using HST/WFPC2 through the F300W, F450W, F606W, and F814W filters (Fernández-Soto, Lanzetta & Yahil 1999–FLY99 hereafter; Lanzetta et al. in preparation). The sample includes 27 galaxies, as given in Table 1. The reader is referred to FLY99, C00, D01, and Fernández-Soto et al. 2001 for more details on individual galaxies.

We determine the absolute $`AB`$ magnitude of each galaxy using a fiducial spectral energy distribution determined from a template fitting technique that we will describe in the next section. We adopt a cosmological model with vacuum energy density $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$, matter density $`\mathrm{\Omega }_M=0.35`$, and a Hubble constant $`H_0=65`$ km s-1 Mpc-1. Absolute $`AB`$ magnitudes of the galaxies at rest-frame wavelengths $`\lambda =1500`$ Å range from $`21.8`$ to $`25.4`$. Apart from these absolute magnitudes, all other results in this paper are independent of the chosen cosmological model.

## 3 Method

Our method is based on maximum-likelihood fitting of model galaxy spectral energy distributions, including the effects of intrinsic Lyman-limit absorption and random realizations of intervening Lyman-series and Lyman-limit absorption, to space- and ground-based broad-band images. In practice, our method proceeds in several steps as follows:

First, we determine a fiducial spectral energy distribution of each galaxy by fitting model spectral energy distributions to the broad-band photometric measurements that fall between rest-frame wavelengths $`\lambda =1250`$ and 2800 Å. The wavelength range is chosen to avoid the spectral discontinuity produced by the Lyman series (rest-frame wavelength $`\lambda <1215`$ Å) and the Balmer break (rest-frame wavelengths $`\lambda \begin{array}{c}>\\ \end{array}`$ 3000 Å). The fiducial spectral energy distributions are based on the six spectrophotometric templates (of E, Sbc, Scd, Irr, SB1, and SB2 galaxies) of our previous photometric redshift measurements (Yahata et al. 2000; see also Benítez 2000) together with a power-law form $`f_\nu \lambda ^\alpha `$. (The power-law form is required because a few of the galaxies are found to be bluer–in the ultraviolet rest-frame–than any of the starburst spectrophotometric templates. In each of these cases, $`\alpha <0.5`$.) The spectral energy distributions are selected according to a $`\chi ^2`$ analysis similar to the analysis used for our previous photometric redshift measurements (Lanzetta et al. 1996; FLY99).

Next, we modify the fiducial spectral energy distribution of each galaxy by adding the effect of intrinsic Lyman-limit absorption, which is characterized by an assumed Lyman-limit optical depth $`\tau _\mathrm{g}`$ that is taken to be universal across all galaxies. For a Lyman-limit optical depth $`\tau _\mathrm{g}`$, the output spectrum $`f(\lambda ,\tau _\mathrm{g})`$ in terms of the input spectrum $`f(\lambda ,0)`$ at rest-frame wavelength $`\lambda `$ is

$$f(\lambda ,\tau _\mathrm{g})=f(\lambda ,0)\mathrm{exp}[\tau _\mathrm{g}(\lambda /\lambda _{\mathrm{LL}})^3]$$
(1)
at rest-frame wavelengths $`\lambda <\lambda _{\mathrm{LL}}`$, where $`\lambda _{\mathrm{LL}}`$ is the rest-frame wavelength of the Lyman limit.

Next, we modify the fiducial spectral energy distribution of each galaxy by adding the effect of random realizations of intervening Lyman series and Lyman limit absorption. Here we characterize the distribution of Lyman $`\alpha `$ forest absorption systems by a distribution function $`F(N,b,z)`$, which is defined in such a way that $`F(N,b,z)dNdbdz`$ is the expected number of absorption systems in the neutral hydrogen column density interval $`dN`$ around $`N`$, Doppler parameter interval $`db`$ around $`b`$, and redshift interval $`dz`$ around $`z`$. We adopt as the functional form of the distribution function

$$F(N,b,z)(1+z)^\gamma N^\beta G^{}(\widehat{b},\sigma _b,b_{\mathrm{min}}),$$
(2)
where $`G^{}(\widehat{b},\sigma _b,b_{\mathrm{min}})`$ is a “truncated” Gaussian distribution of mean $`\widehat{b}`$, dispersion $`\sigma _b`$, and minimum truncation value $`b_{\mathrm{min}}`$. We adopt $`\gamma =2.5`$, $`\beta =1.6`$, $`\widehat{b}=23`$ km s-1, $`\sigma _b=8`$ km s-1, and $`b_{\mathrm{min}}=15`$ km s-1.

We have extensively checked that simulated spectra generated according to these parameters reproduce to within observed scatter the average Lyman $`\alpha `$ absorption $`D_A`$ and average Lyman series absorption $`D_B`$. In Figure 1 we show the average values of $`D_A`$ and $`D_B`$, compared to values found in the literature. The agreement is perfect within the scatter present in the observational data.

The incidence of Lyman limit systems (those with hydrogen opacity $`\tau >1`$) is an important factor in our simulation, as they are amongst the main drivers of the absorption bluewards of the galactic Lyman edge. We plot in Figure 2 the density of Lyman limit systems generated by our Lyman $`\alpha `$ forest model, and compare it with the numbers measured by different authors (as described in the Figure). As can be seen, our model produces slightly less Lyman limit systems at the lowest redshifts than have been observed. Whereas this effect could produce a bias favouring lower escape fraction values, we have checked that this effect is not important (see SS4.2 and 4.3 for details). This slight deffect is the result of our decission to adopt a single parameterization to generate all of the Lyman $`\alpha `$ forest–a change in the normalisation to improve the aspect of Figure 2 would conflict with the data shown in Figure 1. Globally, we consider that our model reproduces within reasonable limits the properties of the intergalactic medium.

Next, we measure the simulated spectrum of each galaxy—constructed from the fiducial spectral energy distribution of each galaxy modified by the effects of intrinsic Lyman-limit absorption with an assumed Lyman-limit optical depth $`\tau _\mathrm{g}`$ and of random realizations of intervening Lyman series and Lyman limit absorption—at short-wavelengths, i.e. over the F300W bandpass which is sensitive at wavelengths below the rest-frame Lyman limit, and we compare observed energy fluxes $`f_{\mathrm{obs}}^{(i)}`$ and uncertainties $`\sigma _{\mathrm{obs}}^{(i)}`$ and simulated energy fluxes $`f_{\mathrm{sim}}^{(i,j)}`$ of the $`i`$th galaxy in the $`j`$th simulation. To randomize the effects of intervening Lyman series and Lyman limit absorption, we repeat this procedure a large number of times, and we form the likelihood as a function of Lyman-limit optical depth $`\tau _\mathrm{g}`$ as

$$(\tau _\mathrm{g})=\underset{i=1}{\overset{N_{\mathrm{gal}}}{}}\left(\frac{1}{N_{\mathrm{sim}}}\underset{j=1}{\overset{N_{\mathrm{sim}}}{}}\frac{\mathrm{exp}\left[\frac{(f_{\mathrm{obs}}^{(i)}f_{\mathrm{sim}}^{(i,j)}(\tau _\mathrm{g}))^2}{2\sigma _{\mathrm{obs}}^{(i)}{}_{}{}^{2}}\right]}{\sqrt{2\pi \sigma _{\mathrm{obs}}^{(i)}{}_{}{}^{2}}}\right),$$
(3)
where the product extends over the $`N_{\mathrm{gal}}`$ galaxies of the sample and where the sum extends over the $`N_{\mathrm{sim}}`$ simulations of the analysis. We find that $`N_{\mathrm{sim}}=100`$ is enough to represent the variation of intervening Lyman series and Lyman limit absorption.

## 4 Results and Discussion

Figure 3 shows results of the method described in § 3 applied to five examples of the galaxies. The examples are chosen to span a range in spectral type, redshift, and luminosity. For clarity, only 10 of the 100 realizations are presented for each galaxy. It is obvious from inspection of Figure 3 that in many cases (especially at the lowest redshifts) the fluxes observed through the F300W filter are less than the fluxes expected, if the galaxies were completely transparent to ionising photons. This trend is also present in the other galaxies of the sample. Our quantitative maximum-likelihood analysis of this trend across the entire sample of galaxies establishes a lower limit to the optical depth at the Lyman limit or an upper limit to the escape fraction for these galaxies. In this section, we discuss these limits.

### 4.1 Galaxy Escape Fraction

Figure 4 and Table 2 show results of the method described in § 3 applied to the data described in § 2. The first panel of Figure 4 shows the relative likelihood as a function of $`\mathrm{log}(\tau _\mathrm{g})`$ determined from the entire sample of galaxies. The best-fit value of the galaxy Lyman-limit optical depth determined from the entire sample of galaxies is

$$\tau _g=5.2_{0.8}^{+1.4},$$
(4)
and the corresponding best-fit value of the galaxy escape fraction determined from the entire sample of galaxies is

$$f_{\mathrm{esc}}=0.008\pm 0.006.$$
(5)
These results apply at the average redshift of the entire sample of galaxies, which is

$$z=2.77$$
(6)
The $`3\sigma `$ upper limit to the galaxy escape fraction determined from the entire sample of galaxies is

$$f_{\mathrm{esc}}<0.039.$$
(7)
Apparently, our analysis rules out the possibility of a large galaxy escape fraction at high significance but cannot rule out (at the $`3\sigma `$ level) the possibility of a galaxy escape fraction identically equal to zero. Based on these results, we conclude that galaxies of redshift $`z3`$ are highly opaque to ionising photons.

### 4.2 Possible selection Effects

Results of our analysis are apparently at odds with results of Steidel et al. (2001) that galaxies of $`z3.4`$ are not highly opaque to ionizing photons. These authors observed in the average spectrum of 29 galaxies a flux ratio $`F_{9/15}=0.056`$. After correcting for the effects of the Lyman $`\alpha `$ forest absorption using a composite QSO spectrum at the same redshift, they found an intrinsic flux ratio of $`F_{9/15}=0.22\pm 0.05`$–where the effects of intrinsic galactic absorption and attenuation by dust or gas have not been eliminated. The authors further incorporated to this measurement a plausible estimate of the flux discontinuity at the Lyman limit due to absorption in stellar atmospheres to derive an escape fraction $`f_{\mathrm{esc}}0.10`$.111Steidel and collaborators used a definition of $`f_{\mathrm{esc}}`$ which is different from the usual one–they normalised the escape fraction at 900 Å to the escape fraction at 1500 Å. They reported $`f_{\mathrm{esc}}0.50`$, which is converted here to the usual definition assuming the escape fraction at 1500 Å is $`0.15-0.20`$ as quoted by the authors.

Steidel and collaborators suggested that the high escape fraction found in their work might result from some selection bias of their galaxy sample. The galaxies incorporated into their analysis include only 29 of 875 galaxies for which they obtained spectroscopic observations with the Keck telescope, selected on the basis of quality of the observations at blue wavelengths. This selection condition quite naturally selects galaxies that are either very blue, very luminous, or both. QSOs are known to be transparent to ionizing photons at ultraviolet wavelengths (Zheng et al. 1997), so it is not unreasonable to speculate that very blue or very luminous (or both) galaxies might exhibit larger escape fractions than avarage galaxies. To examine this possibilty, we repeated the analysis described in § 3 for various subsamples of the galaxy sample.

Color: First, we repeated the analysis for two subsamples of the galaxy sample divided according to color. Results are shown in the second panel of Figure 4 and the second group of entries in Table 2. Here “bluer” galaxies are galaxies best described by SB1 or power-law spectrophotometric templates and “redder” galaxies are galaxies best described by Scd, Irr, and SB2 spectrophotometric templates. Results obtained from both subsamples are consistent with each other and with results obtained from the entire sample of galaxies, although there is a (statistically insignificant) trend for the bluer galaxies to indicate a smaller galaxy escape fraction than is indicated by the redder galaxies. This result runs contrary to the idea described above.

Luminosity: Next, we repeated the analysis for two subsamples of the galaxy sample divided according to luminosity. Results are shown in the third panel of Figure 4 and the third group of entries in Table 2. Here “low-luminosity” galaxies are galaxies of absolute magnitude $`AB>23.75`$ and “high-luminosity” galaxies are galaxies of absolute magnitude $`AB<23.75`$. Results obtained from both subsamples are consistent with each other and with results obtained from the entire sample of galaxies, although there is a (statistically insignificant) trend for the lower-luminosity galaxies to indicate a smaller galaxy escape fraction than is indicated by the higher-luminosity galaxies.

Redshift: Finally, we repeated the analysis for two subsamples of the galaxy sample divided according to redshift. Results are shown in the third panel of Figure 4 and the third group of entries in Table 2. Here “low-redshift” galaxies are galaxies of redshift $`1.95<z<2.85`$ and “high-redshift” galaxies are galaxies of redshift $`2.85<z<3.50`$. Results obtained from both subsamples are consistent with each other and with results obtained from the entire sample of galaxies, although there is a (statistically insignificant) trend for the lower-redshift galaxies to indicate a smaller galaxy escape fraction than is indicated by the higher-redshift galaxies. To some extent, this statistically insignificant difference could be due to the already mentioned relative lack of Lyman limit systems in our model at low redshift. Given that the results from the high-redshift sample (where our model reproduces exactly the observed number of Lyman limit systems) are not significatively different from the results obtained using the low redshift sample, we infer that the difference in the model does not induce a large change.

It should be remarked that none of the subsamples analysed is different from the values given by the complete sample at the 3$`\sigma `$ level. Neither are any of the subsamples incompatible with their complementary subsamples at the same level. With this in mind, it is nevertheless interesting that the apparent effects of luminosity and ultraviolet blueness on $`f_{\mathrm{esc}}`$ seem to pull in opposite directions, suggesting that the higher escape fraction reported by Steidel and collaborators is unlikely due to selection effects.

### 4.3 Systematic uncertainties

All the confidence limits reported above refer only to the statistical properties of the sample of galaxies we are analysing. It is very likely that the method we are using to analyse the data may introduce errors of a systematic nature, that could potentially be of the same order of the statistical ones. We have performed some checks in order to estimate which could be the contribution from several sources of uncertainty.

#### 4.3.1 Uncertainties induced by the Lyman $`\alpha `$ forest model

As was shown in §3 above, our Lyman $`\alpha `$ forest model has been checked against the observed properties, and agrees with them within observational limits. In order to estimate the uncertainty in $`f_{\mathrm{esc}}`$ that could originate from the uncertainties in the model, we have performed new measurements of $`f_{\mathrm{esc}}`$ using forest models that differ from our “standard” one by applying a 10% increase/decrease in the normalisation of the line density. We must remark that this 10% change does largely overestimate the possible uncertainties in our model–the models with 10% more or less lines represent bad fits to the observed forest properties.

As expected, an increase in the forest density implies an increase in the measured value of the escape fraction: more of the UV photons can escape each galaxy while keeping the observed photometric properties, as their mean free path is reduced by the increased density of absorbers. Mutatis mutandis, when the forest density is reduced, the value of $`f_{\mathrm{esc}}`$ moves down. The change is, however, small: at the $`3\sigma `$ confidence level, the original limit ($`f_{\mathrm{esc}}<0.039`$) becomes $`f_{\mathrm{esc}}<0.030`$ ($`f_{\mathrm{esc}}<0.046`$) when the line density is decreased (increased) by 10%.

Another potential problem in our forest model could be the lack of clustering in our simulations, as opposed to what is shown by the observations. Our measurement is not affected by the small-scale clustering that has been observed in the Lyman $`\alpha `$ forest lines (Fernández-Soto et al. 1996) because we are integrating over much larger scales. Of potential importance is the observed clustering of high-column density absorbers (Sargent, Boksenberg & Steidel 1988). It is difficult to quantify its effect, but we expect it should be dilluted by the size of our sample, and estimate that it produces a smaller variance than the one induced by the $`\pm `$10% change in the number of lines presented above.

#### 4.3.2 Uncertainties induced by the choice of the spectral templates

In the case of some of the galaxies, the photometric data available could allow for the assignation of two different fiducial spectral templates with similar goodness-of-fit values. The choice of one of them over the other (which we perform algorithmically via the minimum chi-squared criterium) could produce an extra uncertainty in the measurement of $`f_{\mathrm{esc}}`$. The same effect can be produced if our selection of templates is not dense enough to cover the possible UV shapes, or if the photometric uncertainties are too big.

We have evaluated the uncertainty induced by this effect by performing the following exercise: we have forced all spectral templates to be one step bluer (always along the sequence formed by Ell, Sbc, Scd, Irr, SB2, SB1) than the one which is the best-fit choice333Exception made of the ones that are already fitted by a SB1 template or a power-law, and repeated the analysis. This does in fact result in very bad fits for many of the galaxies, and must thus be considered a very large overestimate of the possible uncertainties. Under this scenario, with the galaxies emitting far more ionising photons than what is indicated by their actual UV continua, the escape fraction diminishes to very low values: $`f_{\mathrm{esc}}<0.016`$ at the $`3\sigma `$ level.

We have also performed the opposite exercise, forcing all galaxies to be one step redder than their best-fit template444Exception made of those galaxies fitted by a power-law. Under this assumption the galaxies emit far less ionising photons than what their UV continua actually suggests, and this results in an increase in the acceptable escape fraction, which reaches $`f_{\mathrm{esc}}=0.065`$ , with a $`3\sigma `$ upper limit $`f_{\mathrm{esc}}<0.156`$.

As we explained before, these limits must be taken as very large overestimates of the possible systematic effects induced by our spectroscopic templates, as we do not realistically expect to mistake all galaxies in the same direction. However, even under this extreme assumption, and incorporating also the effects described in §4.3.1, the upper limit to the measured escape fraction is still far below the values measured by Steidel et al. (2002).

## 5 Conclusions

We have presented a new method to measure the escape fraction of ionising photons from galaxies at $`z3`$. Our method provides a cleaner measurement compared to the classical spectroscopic technique because the problem of the sky subtraction is reduced almost to insignificance and we can use a wider bandpass.

The value we obtain ($`f_{\mathrm{esc}}=0.008`$, with $`f_{\mathrm{esc}}<0.039`$ at a $`3\sigma `$ confidence level) is in agreement with most previous measurements of the escape fraction of ionising photons at high and low redshift, being more stringent than most previous upper limits. It is largely in contradiction with the value presented by Steidel et al. (2001) for high-redshift galaxies555Steidel et al. (2001) remark that their results “should be treated as preliminary until high quality observations of individual galaxies exhibiting clear evidence for Lyman continuum photon leakage become available”. We have studied possible selection effects that could cause this apparent difference and conclude that, within the limits set by the size of our sample, they cannot explain the large discordance between our results. Neither can the difference be explained by possible systematic effects induced by our method, as we also have shown. A (very conservative) estimate of the confidence limits when all systematic effects are included leads to a $`3\sigma `$ upper limit to the escape fraction that cannot be higher than $`f_{\mathrm{esc}}\begin{array}{c}<\\ \end{array}0.15`$.

If the escape fraction of ionising photons reported here represents the general situation at high redshift, then normal galaxies cannot be responsible for any significant fraction of the high-redshift ionising background. On the other hand, if the proximity effect measurements of the background flux are confirmed, then the problem persists to find the objects responsible for the ionisation state of the high-redshift universe. Deep, narrow-band imaging of local and high-redshift galaxies at wavelengths slightly above and below their intrinsic Lyman limits could settle this argument. Such imaging campaigns would be less expensive, in terms of observing time, than the spectroscopic observations done to date.

## 6 Acknowledgments

We thank our referee, Emanuele Giallongo, for his useful comments that have improved the clarity of our paper. AFS gratefully acknowledges support by a Marie Curie Fellowship of the European Community programme “Improving the human research potential and the socio-economic knowledge base” under contracts number MCFI-1999-00494 and MCFI-2002-00472. HWC acknowledges the hospitality of the Osservatorio di Brera-Merate where this work was completed.


format: text

text:
# Dobiński-type relations and the Log-normal distribution

## References


text:
# Reaction mechanisms involved in the production of neutron-rich isotopes

## 1 Introduction

The production of nuclei far from stability has provided a new ground to investigate the structure and the dynamics of nuclei. The present experimental facilities are able to produce a large variety of neutron-deficient isotopes, however, the neutron-rich side of the chart of the nuclides is only accessible up to the drip lines for the lightest nuclei. Consequently, the challenge for most of the future planned rare-beam facilities is the production of neutron-rich isotopes.

Together with the technological improvements concerning high-current accelerators, large-acceptance separators and larger extraction and charge-breeding efficiencies, the optimal choice of the reaction mechanisms will play a major role on the final production rates of neutron-rich nuclei. Different reaction mechanisms can be used to produce these nuclei, however, heavy-ion collisions at low energies like fusion, deep inelastic or multinucleon transfer can only be applied with thin targets limiting the final production rates. Better suited seem to be fragmentation or spallation at high energies and fission. In addition, these two reaction mechanisms allow to produce a large variety of neutron-rich nuclei.

In this paper we discuss the main issues of these two reaction mechanisms. The large collection of data obtained during the last years has brought new model descriptions of these reactions providing reliable predictions on the expected production rates of neutron-rich isotopes in future planned rare-beams facilities.

## 2 Production of heavy neutron-rich nuclei

Heavy exotic nuclei (Z$`>`$70) can be produced by means of fusion-evaporation reactions or in fragmentation(spallation) of heavy nuclei. Both reaction mechanisms lead mainly to the production of neutron-deficient residues. However, it has recently been shown that fragmentation reactions at relativistic energies present large fluctuations in the N/Z distribution of the final residues and in its excitation-energy distribution. In particular, the proton-removal channels have been investigated in cold-fragmentation reactions where only protons are abraded from the projectile, while the induced excitation energy is below the particle-emission threshold. These reactions can lead to the production of heavy neutron-rich nuclei beyond the present limit of the chart of the nuclides as shown in Fig. 1.

Cold-fragmentation reactions can be described in terms of the abrasion-ablation model as a two-step process. First, the interaction between projectile and target leads to a projectile-like residue with a given excitation energy which statistically de-excitates by particle evaporation or fission. A new analytical formulation of the abrasion-ablation model, the code COFRA, has been developed in order to calculate the expected low production cross sections of extremely neutron-rich nuclei which are not reachable with Monte Carlo codes.

The results of these calculations are shown in Fig. 1. The different lines correspond to calculations with different enhancement factors of the excitation energy induced in the collision due to particle-hole excitations (see reference for details). As can be seen, the production cross sections of very neutron-rich nuclei can nicely be described with the new analytical formulation of the abrasion-ablation model. In Fig. 2 we compare the predictions of the code with the production cross sections of proton-removal channels populated in different reactions .

Since the COFRA code describes correctly all the measured proton-removal channels in relativistic heavy-ion collisions, it was used to determine the expected production of heavy neutron-rich nuclei in future rare-beam facilities. The results of these calculations are shown in Fig. 3. In this figure we report the expected production cross sections of heavy neutron-rich nuclei that can be obtained in the fragmentation of 238U, 208Pb and 174W. According to these calculations, large progress is expected in this region of the chart of the nuclides, where the r-process path may even be reached around the end point N=126.

## 3 Production of medium-mass neutron-rich nuclei

### 3.1 Fission

Fission has largely been used to produced medium-mass neutron-rich nuclei up to the present limits of the chart of the nuclides . The isotopic distribution of residues produced in fission can be understood in terms of the potential governing this process. The Coulomb term of the nuclear potential is responsible for the neutron excess of the stable fissile nuclei leading to fission residues with an even larger neutron excess with respect to the valley of beta stability. However, the asymmetry term preserves the same N/Z of the fissioning nucleus in the fission residues. Shell effects and temperature induce a polarization effect which allow to produce even more neutron-rich residues.

In order to investigate the fluctuations in N/Z and mass asymmetry induced by the temperature of the fissioning system we have performed several simulations with the fission code of Ref.. In figure 4 we represent on top of a chart of the nuclides the distributions of residues after the fission of 238U at different excitation energies. As can be seen in the figure, when increasing the excitation energy, shell effects (double humped distribution) disappear, and the fluctuations in mass asymmetry and N/Z increase, populating a larger variety of neutron-rich residues. However, at high excitation energies neutron evaporation becomes more important, and the residue distribution moves to the neutron-deficient side. One can define an optimum excitation energy of the fissioning system around 50 MeV to produce the largest variety of neutron-rich nuclei.

Once we know the optimum conditions to produced neutron-rich nuclei in fission reactions, the final production rates will be defined by the evolution of the production cross sections. In Fig.5, we report the production cross sections of several neutron-rich tin and nickel isotopes produced in the fission of 238U projectiles on Be, Pd and H2. These data show a dramatic decrease of the production cross section of around one order of magnitude per additional neutron. Consequently, one can not expect a spectacular expansion of the chart of the nuclides in the region of medium-mass neutron-rich isotopes by using fission reactions with the future rare-beam facilities.

### 3.2 Two-step scenario: fission + cold fragmentation

Recently it has been proposed to use both fission and cold-fragmentation reactions, in a two-step reaction scheme, in order to overcome the limitations of fission to produce extremely neutron-rich nuclei in the medium-mass region. However, nowadays it is difficult to make reliable predictions of the final production rates using this idea. As shown in Fig. 6, present fragmentation codes show clear discrepancies in the predicted production rates when neutron-rich projectiles are used. In addition, the energy at which the fragmentation stage takes place plays a major role as shown in Fig. 7. In this figure we report the measured production cross sections of different tin isotopes in the reaction 129Xe on aluminum at 790 A MeV  and 50 A MeV . According to these data the higher energies enhance the production of neutron-rich residues.

Nevertheless, we used the cold-fragmentation code COFRA, which is one of the more reliable codes to predict the production of fragmentation residues from neutron-rich projectiles to estimate the production rates in a two-step reaction scheme at energies above 100 A MeV. In these calculations, the primary production cross sections in the fission step were taken from measured data in the reactions 238U(1 A GeV)+p  and 238U(750 A MeV)+Be . The most representative results of these calculations are shown in figure 8. In this figure we represent the production of different neutron-rich isotopes along the neutron shells N=50 and N=82. The thick line represents the direct production in the fission of 238U induced by 1 GeV protons, while the thin lines correspond to the two-step production after the cold fragmentation in a Be target of different Ga and Sn isotopes produced by fission. As can be observed in both pictures, the two-step scheme is competitive with the direct production only for the largest neutron excess.

However, if one considers extraction, ionization and re-acceleration efficiencies, the situation could change. In this realistic scenario the two-step scheme could be used to produce by fission an abundant long-lived neutron-rich nuclei like 132Sn and fragmenting it to produce those neutron-rich isotopes that have low extraction efficiencies.

## 4 Conclusion

In this paper we have reviewed the two reaction mechanisms best suited to extend the present limits of the chart of the nuclides in the neutron-rich side. Recent experiments have shown that the cold-fragmentation process constitutes the appropriate tool to produce heavy neutron-rich residues. In the intermediate-mass region it is well established that fission allows to produce abundantly moderately neutron-rich isotopes. However, it will be difficult to enlarge the present limits of the chart of the nuclides using this reaction mechanism. An alternative would be to use a two-step scenario where the cold fragmentation of neutron-rich isotopes produced by fission is foreseen. Although this new scenario is only competitive with fission for isotopes with a large neutron excess, it would be an optimum solution to produce neutron-rich isotopes with low extraction efficiencies.


text:
# 1 Starting point

## 1 Starting point

For the first part of this story we will define what a scattering cross section is. We will only consider in this work elastic scattering of pions. Two pions can scatter at a certain energy which we will call $`M_{\pi \pi }`$. The differential cross section $`\sigma `$ at a given $`M_{\pi \pi }`$ is

$$\frac{d\sigma }{d\varphi d\theta }=\frac{1}{K^2}\left|\underset{\mathrm{}}{}(2\mathrm{}+1)T_{\mathrm{}}P_{\mathrm{}}(cos(\theta ))\right|^2$$
(1)
where $`\varphi `$ and $`\theta `$ are the azimuthal and scattering angles, respectively. $`T_{\mathrm{}}`$ is a complex scattering amplitude and $`\mathrm{}`$ is the angular momentum. $`P_{\mathrm{}}`$ is the Legendre polynomial, which is a function of $`cos(\theta )`$. $`K`$ is the flux factor equal to the pion momentum in the center of mass. For this note we only consider $`\mathrm{}=`$ 0 and $`\mathrm{}=`$ 1. The $`T_0`$ and $`T_1`$ elastic scattering amplitudes are complex amplitudes described by one real number which is in units of angles. The form of the amplitude is

$$T_{\mathrm{}}=e^{i\delta _{\mathrm{}}}sin(\delta _{\mathrm{}})$$
(2)
We note that $`\delta `$ depends on the value of $`\mathrm{}`$ and $`M_{\pi \pi }`$. We will use the $`\mathrm{}=`$ 0 and $`\mathrm{}=`$ 1 $`\delta `$’s given in .

## 2 Defining A

Let us consider two pions scattering in the final state of the heavy ion collision. The scattering will be either $`\mathrm{}=`$ 0 or $`\mathrm{}=`$ 1 partial waves. The $`M_{\pi \pi }`$ of the scattering di-pion system will depend on the probability of the phase space of the overlapping pions. The pions emerge from a close encounter in a defined quantum state with a random phase. We will call this amplitude $`A`$ and note that the absolute value squared of the amplitude is proportional to the phase space overlap. The emerging pions can re-interact or re-scatter through the quantum state of the pions, which is a partial wave or a phase shift. We have amplitude $`A`$ plus $`A`$ times the re-scattering of pions through the phase shift consistent quantum state of $`A`$. The correct unitary way to describe this process is given by  equation (4.5)

$$T=\frac{V_1U_1}{D_1}+\frac{\left(V_2+\frac{D_{12}V_1}{D_1}\right)\left(U_2+\frac{D_{12}U_1}{D_1}\right)}{D_2\frac{D_{12}^2}{D_1}}$$
(3)
In the above equation we have two terms, 1 and 2. The first term denoted by 1 is the $`\pi \pi `$ scattering through p-wave which will become the amplitude $`A`$ mentioned above, where $`V`$ is the incoming and $`U`$ is the outgoing $`\pi \pi `$ system. The second term denoted by 2 is the direct production of the $`\rho `$ meson with $`V`$ being the production, the propagation being $`D`$ and the decay being $`U`$. We see that there are terms $`D_{12}`$ which involves a loop of pions between scattering pions and the formation of a $`\rho `$ by the pions.

## 3 Final equation

The complete derivation is in the appendix. From the appendix we get two terms, one being the direct production of the $`\rho `$ or $`\pi \pi `$ p-wave phase shift and the second being the $`\rho `$ from re-scattering. The final equation number 21 has two important factors, one is two-body phase space and the other is a coefficient $`\alpha `$. This coefficient is related to the real part of the $`\pi \pi `$ rescattering loop and is given by equation 14. When the pions rescatter or interact at a close distance or a point $`\alpha `$ has its maximum value of one. While if the pions rescatter or interact at a distance determined by the diffractive limit the value of $`\alpha `$ is zero.

In the equation 21 $`|T|^2`$ is the cross section for $`\rho `$ production, where $`D`$ is the direct production amplitude and $`A`$ is the amplitude introduced above for the re-scattering pions into the $`\rho `$ meson. $`\delta `$ is the $`\pi \pi `$ phase shift . $`q`$ is the $`\pi \pi `$ center of mass. At a given $`p_T`$ and $`y`$ bin, $`D`$ will be a constant as a function of $`M_{\pi \pi }`$. The $`\alpha `$ factor is one for rescattering coming from a point source, but goes to zero when rescattering is diffractive. The dependence of $`A`$ is calculated by the phase space overlap of di-pions added as four vectors and corrected for proper time, with the sum having the correct $`p_T`$ and $`y`$ for a given $`M_{\pi \pi }`$.

Finally we most use the correct two body phase space. For a two body system of pions, phase space goes to an constant as $`M_{\pi \pi }`$ goes to infinity. Let us choose this constant to be unity. Phase space which is denoted by PS is equal to

$$PS=\frac{2qB_{\mathrm{}}(q/q_0)}{M_{\pi \pi }}$$
(4)
where $`B_{\mathrm{}}`$ is a Blatt-Weisskopf-barrier factor for $`\mathrm{}`$ angular momentum quantum number. The $`q_0`$ is the momentum related to the range of interaction of the $`\pi \pi `$ scattering. 1 fm is the usual interaction distance which implies that $`q_0`$ is 200 MeV/c. For the $`\rho `$ meson $`\mathrm{}`$=1 and $`B_1`$ = $`\frac{(q/q_0)^2}{(1+(q/q_0)^2)}`$.

The phase space factor PS as a function of $`q`$ near the $`\pi \pi `$ threshold is given by $`q^{2\mathrm{}+1}`$. Thus in the appendix we use $`q^{2\mathrm{}+1}`$ for the factor PS except for equation 21 which is the final equation.

## Appendix A Appendix

Starting with equation (4.5) from

$$T=\frac{V_1U_1^{}}{D_1}+\frac{\left(V_2+\frac{D_{12}V_1}{D_1}\right)\left(U_2^{}+\frac{D_{12}U_1^{}}{D_1}\right)}{D_2\frac{D_{12}^2}{D_1}}$$
(5)
In order to have the correct threshold kinematics, we define

$$U_1^{}=U_1\sqrt{q^{2\mathrm{}+1}}$$
(6)
$$U_2^{}=U_2\sqrt{q^{2\mathrm{}+1}}$$
(7)
where $`q`$ is the $`\pi \pi `$ center of mass momentum and $`\mathrm{}`$ is the value of the angular momentum. The amplitude $`A`$ of the text is given by

$$\frac{V_1U_1}{D_1}=A$$
(8)
Thus we have

$$\frac{V_1U_1^{}}{D_1}=A\sqrt{q^{2\mathrm{}+1}}$$
(9)
The phase shift for the $`\mathrm{}^{th}`$ partial wave will be given by $`\delta _{\mathrm{}}`$, where

$$\frac{U_2^{}U_2^{}}{D_2}=e^{i\delta _{\mathrm{}}}sin(\delta _{\mathrm{}})$$
(10)
The above equality is true if the $`D_1`$ mode plays no role in the $`\pi \pi `$ scattering in the $`\mathrm{}^{th}`$ partial wave. But in the initial state there is a large production of $`D_1`$. The $`U^{}`$s are the basic coupling of the $`D^{}`$s to the $`\pi \pi `$ system. In order to decouple $`D_1`$ from the $`\pi \pi `$ system $`U_1`$ must go to zero. We can maintain a finite production of $`D_1`$ if we define

$$V_1=\frac{1}{U_1}$$
(11)
Thus the first term in the equation becomes

$$\frac{V_1U_1^{}}{D_1}=\frac{\frac{1}{U_1}U_1\sqrt{q^{2\mathrm{}+1}}}{D_1}=\frac{\sqrt{q^{2\mathrm{}+1}}}{D_1}=A\sqrt{q^{2\mathrm{}+1}}$$
(12)
The form of $`D_{12}`$ is given by

$$D_{12}=\alpha U_1U_2+iq^{2\mathrm{}+1}U_1U_2$$
(13)
$`D_{12}`$ is the real and imaginary part of the two pion loop from state 1 to state 2. The $`U^{}`$s are the $`\pi \pi `$ couplings and the imaginary part goes to zero at the $`\pi \pi `$ threshold. The $`\alpha `$ factor is one for rescattering coming from a point source, but goes to zero when rescattering is diffractive. A simple form for $`\alpha `$ is given by

$$\alpha =(1.0\frac{r^2}{r_0^2})$$
(14)
where $`r`$ is the radius of rescattering in fermis and $`r_0`$ is 1.0 fermi or the limiting range of the strong interaction. The second term of the first equation is

$$\frac{\left(V_2+\frac{D_{12}V_1}{D_1}\right)\left(U_2^{}+\frac{D_{12}U_1^{}}{D_1}\right)}{D_2\frac{D_{12}^2}{D_1}}$$
(15)
Rewriting

$$\frac{\left(V_2+\frac{\alpha U_1U_2V_1}{D_1}+iq^{2\mathrm{}+1}\frac{U_1U_2V_1}{D_1}\right)\left(U_2^{}+\frac{D_{12}U_1^{}}{D_1}\right)}{D_2\frac{D_{12}^2}{D_1}}$$
(16)
Let us make substitutions

$$V_1=\frac{1}{U_1},U_2=\frac{U_2^{}}{\sqrt{q^{2\mathrm{}+1}}},\frac{1}{D_1}=A,D_{12}=0$$
(17)
The second term becomes

$$\frac{\left(V_2+\frac{A\alpha U_2^{}}{\sqrt{q^{2\mathrm{}+1}}}+i\sqrt{q^{2\mathrm{}+1}}AU_2^{}\right)U_2^{}}{D_2}$$
(18)
The first term is

$$\frac{V_1U_1^{}}{D}=A\sqrt{q^{2\mathrm{}+1}}$$
(19)
Adding the first and the second terms and substituting the phase shift,

$$T=\frac{V_2}{U_2}\frac{e^{i\delta _{\mathrm{}}}sin(\delta _{\mathrm{}})}{\sqrt{q^{2\mathrm{}+1}}}+A\left(\frac{e^{i\delta _{\mathrm{}}}\alpha sin(\delta _{\mathrm{}})}{\sqrt{q^{2\mathrm{}+1}}}+\sqrt{q^{2\mathrm{}+1}}e^{i\delta _{\mathrm{}}}cos(\delta _{\mathrm{}})\right)$$
(20)
The term with the factor $`\frac{V_2}{U_2}`$ is the direct production of the di-pion system. We shall call this amplitude $`D`$. The re-scattered amplitude is $`A`$ and is modified by the di-pion phase shift. These two amplitudes have some random phase and are not coherent. Thus the cross section is

$$|T|^2=|D|^2\frac{sin^2(\delta _{\mathrm{}})}{PS}+\frac{|A|^2}{PS}\left|\alpha sin(\delta _{\mathrm{}})+PScos(\delta _{\mathrm{}})\right|^2$$
(21)


text:
# Clumps in large scale relativistic jets.

## 1 Introduction

Thanks to Chandra observations, we now know that the jets of extragalactic radio sources have X–ray knots which are tens to hundreds of kpc away from the nucleus. This seems to be true both for powerful flat radio spectrum sources (FSRQ) with a jet probably aligned with the line of sight (Chartas et al. 2000; Pesce et al. 2001; Sambruna et al. 2002; Sambruna et al. 2001; Schwartz 2002; Schwartz et al. 2000; Siemiginowska et al. 2002; Tavecchio et al. 2000), and also for closer radio–galaxies whose jets are observed at large viewing angles (e.g. Wilson & Yang, 2002; Kraft et al. 2002; Hardcastle et al., 2002; Hardcastle et al. 2001; Wilson et al. 2001).

In the powerful FSRQ the X–ray luminosity often exceeds the optical and radio power emitted in the same knots. Thermal and synchrotron self Compton models have difficulties in explaining the data (see e.g. Chartas et al. 2000; Schwartz et al. 2000; Siemiginowska et al. 2002; Harris & Krawczynski 2002), while interpreting the X–rays as due to inverse Compton radiation off the cosmic microwave background (CMB), and the radio–optical emission as due to synchrotron, can well explain the observations and at the same time minimize the energy requirements (Celotti, Ghisellini & Chiaberge 2001; Tavecchio et al. 2000; Ghisellini & Celotti 2001). This however requires the jet to be highly relativistic even at the largest scales: in the case of PKS 0637–712 the bulk Lorentz factor is of the order of 10–15 at hundreds of kpc away from the core.

In this case the CMB energy density is seen boosted in the comoving frame of the jet, by a factor $`\mathrm{\Gamma }^2`$, causing it to dominate over the synchrotron and the magnetic field energy densities existing locally.

Another consequence of this scenario is that it is possible to constrain the low energy cut off of the emitting particle distribution, $`\gamma _{\mathrm{min}}m_ec^2`$. Taking PKS 0637–752 as an example, we require $`\gamma _{\mathrm{min}}<30`$ in order that the inverse Compton spectrum starts before the X–ray band (assuming $`\mathrm{\Gamma }=14`$). We also require $`\gamma _{\mathrm{min}}>10`$ to avoid the inverse Compton process overproducing the observed optical emission. This is true independent of whether the optical emission belongs to the initial part of the IC spectrum or to the high energy tail of the synchrotron flux. However, a synchrotron origin of the optical flux corresponds to emission of ultrarelativistic electrons of relatively short cooling time, which nicely corresponds (albeit within the large uncertainties) to the dimensions of the optical knots as observed by HST (Schwartz et al. 2000; Celotti, Ghisellini & Chiaberge 2001).

In this paper we consider sources observed at small viewing angles, that can emit beamed X–rays by their plasma moving at bulk relativistic speeds. Since we see bright knots, we know that there is some process accelerating particles at very high energies: to produce optical radiation by synchrotron emission with reasonable values of the magnetic field, we require electrons of TeV energies. At these scales, the cooling time of the electrons is relatively long for all but the highest energy particles, and it is therefore of interest to study what is the predicted spectral energy distribution (SED) outside the bright knots, taking into account that the energy density (in the comoving frame) of the CMB is constant for a non–decelerating jet, and therefore the inverse Compton scattering of these photons is always important. On the other hand it is reasonable to assume that the magnetic field becomes smaller for increasing distances from the nucleus, causing the synchrotron luminosity to decrease even if the emitting particles do not cool. In this respect, it is intriguing that we may be able to detect inverse-Compton radiation produced by particles that are only mildy relativistic or even sub-relativistic in the comoving frame but have substantial bulk motions. Is this “bulk Compton” luminosity observable by existing or future instruments, such as ALMA and NGST?

Furthermore, we will discuss observations which already exist, and pose an interesting problem: the X–ray flux, outside the bright knots, is dimming as fast as the optical and the radio fluxes, despite the fact that it should be produced by low energy electrons, which do not cool radiatively. Adiabatic losses are required, but this of course implies that the emitting region expands, and by calculating how much expansion is required we will reach interesting conclusions about the geometry of the acceleration sites.

In order to attack these problems we will first study how the emitting particle distribution evolves in time. In the presence of bulk motion, the time evolution translates to a radial profile from the acceleration site. However, this correspondence is not unique, since it depends on the nature of the bright knots we see: we do not know yet if these correspond to standing or moving features. In the former case relativistic particles originating in a shock move out of it forming a trail; in the latter case there is no connection between the bright knots (which are solidly moving together with the emitting particles) and the interknot regions.

The main result of our analysis is that if the knot has a homogeneous structure, i.e. it is homogeneously filled by a single relativistic particle distribution, radiative and adiabatic losses alone cannot account for the fast decrease of the emission out of the knot. We discuss the possible solutions to this problem and finally we describe a scenario that seems to account for the main properties of the knot emission.

## 2 Adiabatic and radiative losses

In this section we analyze the evolution of particles subject to both radiative and adiabatic cooling and we apply the results to the overall knot emission.

### 2.1 The evolution of the electron distribution

Let us call $`N(\gamma ,r)`$ the particle distribution at some distance $`r`$ from some initial site $`r_0`$, where $`N`$ describes the total number of electrons (not the density). The plasma is assumed to have bulk motion with velocity $`\beta c`$ and Lorentz factor $`\mathrm{\Gamma }`$. In the absence of (re)–acceleration, $`N(\gamma ,r)`$ is described by (Sikora et al. 2001)

$$\frac{N}{r}=\frac{}{\gamma }\left(N\frac{d\gamma }{dr}\right)$$
(1)
This is the usual continuity equation in which time has been substituted by the radial coordinate $`r`$ through $`r=\beta c\mathrm{\Gamma }t^{}`$, and where $`t^{}`$ is the time measured in the comoving frame. Assuming that inverse Compton losses (off CMB photons) dominate over synchrotron losses we have

$$\frac{d\gamma }{dr}=\frac{1}{\beta c\mathrm{\Gamma }}\frac{4\sigma _TcU_{\mathrm{rad}}}{3m_ec^2}\gamma ^2A\frac{\gamma }{r}$$
(2)
where $`U_{\mathrm{rad}}=\mathrm{\hspace{0.17em}4}\times 10^{12}(1+z)^4\mathrm{\Gamma }^2`$ erg cm<sup>-3</sup>, and where the latter term is the adiabatic term: $`A=1`$ for 3D expansion and $`A=2/3`$ for 2D expansion. The solution of Eq. (2) is

$$\gamma =\frac{(1A)}{\alpha r+(r_0/r)^A[(1A)/\gamma _0\alpha r_0]}$$
(3)
where $`\alpha \mathrm{\hspace{0.17em}4}\sigma _TU_{\mathrm{rad}}/(3\beta \mathrm{\Gamma }m_ec^2)`$.

All particles initially between $`\gamma _0`$ and $`\gamma _0+d\gamma _0`$, in a time $`t`$, (or after some distance $`r`$), will form the distribution between the corresponding interval $`\gamma `$ and $`\gamma +d\gamma `$, where $`\gamma `$ is given by Eq. (18). We can therefore set $`N(\gamma ,r)d\gamma =N_0(\gamma _0,r_0)d\gamma _0`$. From Eq. (3) we have:

$$\frac{d\gamma }{d\gamma _0}=\frac{(1A)^2}{\gamma _0^2(r_0/r)^A}\left[\alpha r+\left(\frac{r}{r_0}\right)^A\left(\frac{1A}{\gamma _0}\alpha r_0\right)\right]^2$$
(4)
Therefore the particle distribution at a given distance $`r`$ and energy $`\gamma `$ is given by calculating the corresponding initial $`\gamma _0`$ (which is transformed in $`\gamma `$ in the time to go from $`r_0`$ to $`r`$) and taking into account of the corresponding differentials:

$$N(\gamma ,r)=\frac{N_0(\gamma _0,r_0)}{d\gamma /d\gamma _0}$$
(5)
(This technique has been already applied in the context of gamma–ray bursts by Malesani 2002). In Fig. 1 we show an example of particle evolution calculated according to Eq. 5. Note that total number is conserved: since the low energy cutoff moves to lower energies (because of adiabatic expansion losses), the normalization decreases with $`r`$.

### 2.2 Results

Using Eq. (5) we have computed the evolution of the spectral energy distributions (SED) as a function of $`r`$, assuming that $`N(\gamma _0,r_0)`$ is a power law between $`\gamma _{\mathrm{min},0}`$ and $`\gamma _{\mathrm{max},0}`$. The initial parameters, listed in Tab. 1, have been chosen to roughly match those found by fitting the SED of a knot of a large scale jet visible in radio, optical and X–rays. Note that the magnetic field does not affect the particle evolution (since we have assumed that the radiation energy density of the CMB dominates the cooling), but it determines the amount of the emitted synchrotron radiation. We have assumed that the magnetic field scales as $`r^1`$ (as in the case in which the magnetic field has a prevalent toroidal component or in the case of constant magnetic flux along the jet). We have then assumed that the total number of electrons in a “slice” of the jet of width $`\mathrm{\Delta }r`$ is conserved, and to calculate the emitted flux we have assumed a fixed width $`\mathrm{\Delta }r`$. We have then consistently set the adiabatic constant $`A=2/3`$, corresponding to 2D expansion.

Having computed the SEDs in this way, we have constructed the flux profiles at different frequencies. Jets with a bulk Lorentz factor of $`\mathrm{\Gamma }=4`$ and $`\mathrm{\Gamma }=15`$ have been considered, since this is the range of values found when fitting quasars with large scale X–ray jets (Tavecchio et al. 2000; Sambruna et al. 2002; Celotti et al. 2001). Furthemore we have performed our calculations assuming that the optical flux is produced by the initial tail of the IC process (SEDs in Fig. 2 and Fig. 5; profiles in Fig. 3 and Fig. 6.) or by the high energy tail of the synchrotron process (profiles shown in Fig. 4 and Fig. 7). The only input parameter changing between the two sets of figures is the value of $`\gamma _{\mathrm{max},0}`$.

In Fig. 2 we show the computed SEDs (for the case of $`\mathrm{\Gamma }=4`$) assuming that the optical emission, at the starting point $`r_0`$, is due to initial part of the inverse Compton scattering off the CMB photons. Therefore the electrons producing this radiation have a small energy, and their cooling time is long. In this case the inverse Compton spectra are dominating the bolometric luminosity (this is also true at early times, although this may not be immediately clear from Fig. 2 which cuts the Compton spectrum before its peak), but not by a large factor at the beginning. At late times (or at some distance from $`r_0`$) the Compton dominance increases, since the magnetic field is assumed to decrease as $`r^1`$, while the energy density of the CMB is constant. Fig. 3 shows the expected profiles of the flux in different bands using the same parameters used to construct the SEDs as shown in Fig. 2 and listed in Table 1. The x–axis of this figure can be expressed as a distance from the starting point (bottom axis) or equivalently as the time elapsed from the beginning of our calculations (top axis). In the latter case the light curves shown do not take into account light travel time effects: since the source is very extended, we will receive light from different evolutionary phases of the source (evolved from the near part, younger from the far part of the source). This is important of course for $`t<r_0/c`$: at these times the observed light curve is not the one shown in Fig. 3 (and the analogous Figs. 4, 6 and 7), but should be calculated (as e.g. in Chiaberge & Ghisellini 1999) considering the different travel paths of light rays produced in different parts of the source. For the purposes of the present paper these effects are not crucial, but they will be properly considered in our future studies. Note that if the profiles (i.e. bottom x–axis) shown in Fig. 3 correspond to emission from a standing shock (see below), and so a quasi–stationary case, then light travel time effects are not important, and Fig. 3 (and the other analogous figures) already fully describes the expected behavior.

Fig. 4 shows the profiles assuming that the optical emission is due to the high–energy tail of the synchrotron component. The only parameter different from those adopted for Fig. 2 and Fig. 3 is $`\gamma _{\mathrm{max},0}`$, which is now slightly smaller (see Tab. 1).

Fig. 5 shows, by contrast, the SEDs predicted in the case of $`\mathrm{\Gamma }=15`$. We have assumed that in this case the Compton emission is more dominant than in the previous case, due to the increased energy density of the CMB as seen in the comoving frame ($`\mathrm{\Gamma }^2`$). The other input parameters are listed in Table 1. Note that since the comoving CMB energy density is larger than in the previous case, the corresponding radiative cooling times are shorter, making the evolution of the high energy electrons faster, as can be seen comparing the radio and far IR spectra in Fig. 5 and Fig. 2 for late times (or for large distances). As in Fig. 2, the SEDs shown in Fig. 5 correspond to the case of the optical emission produced by the initial part of the inverse Compton radiation. Fig. 6 shows the corresponding profiles of the flux at different frequencies, while Fig. 7 shows the profiles assuming that, at the start, the optical flux is due to the tail of the synchrotron emission. Fig. 5 shows also that the synchrotron and IC emission components becomes “detached” at late times, when only the low energy electrons survive.

On inspecting the figures showing the flux profiles, a general trend emerges: the synchrotron produced fluxes are dimming faster than the IC/CMB ones, due to i) the different electron energies involved and ii) the decreasing magnetic field. In all cases the X–ray fluxes are roughly linearly decreasing with the distance $`r`$, while the initial decrease of the fluxes produced by the synchrotron process (away from the cutoffs) scale roughly as $`r^3`$. As discussed below, this different behavior will be crucial when comparing our results with already existing data.

As mentioned in the Introduction, there are two possible cases for interpreting the results shown in Figs. 3, 4 and in Figs. 6, 7. Let us call them:

case A: streaming of particles from a stationary (in the observer frame) shock (bottom x–axis, distance from the acceleration site);

case B: rigidly moving blob whose internal particles have been accelerated at $`r_0`$ (top x–axis, time from the acceleration event in the blob reference frame).

If we are in case A (standing shock), then Figs. 3, 4, 6 and 7 directly show the profiles that would be observed in deep images of the jet for $`r>r_0`$ (assuming that the standing shock lives for a sufficiently long time to allow particles to move from $`r_0`$ to $`r`$). In the case of Fig. 3 and Fig. 6 for which the optical flux is due to the low–energy tail of the IC component, the optical flux increases with $`r`$. This happens despite the decrease in the normalization of $`N(\gamma )`$, and it is due to the fact that the optical is produced by the initial tail of the Compton emission (steeply rising with frequency), which increases its flux because of the decrease of $`\gamma _{\mathrm{min}}`$. In other words, in this case the overall spectrum shifts to lower frequencies and to lower fluxes: in the power law part of the spectrum the flux density decreases, while at the low energy end the flux density increases. In fact, when the optical flux is due to the power law part of the external Compton spectrum (for $`r/r_0>10`$), it decreases as the X–ray flux.

In case B (evolving blob), the particle distribution at a given $`r`$ is not directly related to the bright knots we are detecting, but is instead the “fossil” distribution of other unknown, acceleration phases occurred in the past (at $`t^{}=t_0^{}`$ on the top x–axis). Assuming then that there have been such phases in the past, with the optical flux initially produced through IC (Fig. 3 and Fig. 6), we should observe some bright optical knots with relatively weak X–ray and radio fluxes (with respect to the HST/Chandra bright knots). Knots which are bright only in the optical band have not been detected, and this could lead to the conclusion that the optical emission is always due to the high energy tail of the synchrotron spectra, and therefore to very high energy electrons. However, as will be discussed below, this conclusion is premature, since there may be an alternative possibility, in which the knot emission is considered as the result of the contribution of many sub–knot emission sites.

## 3 Comparison with observations

The results of the analysis described above can be compared with real profiles of nearby jets. One of the sources with the best published data is the luminous jet of 3C 273 (Sambruna et al. 2001; Marshall et al. 2001)

The jet presents numerous knots visible in radio, optical and X–rays. Sambruna et al. (2001) and Marshall et al. (2001) reach different conclusions about the viability of a synchrotron interpretation for the X–ray emission in knots A and B1 (the two most luminous portions of the jet). This is due to the different slopes of the optical spectrum derived by the two groups: Sambruna et al. (2001) present an optical spectrum with a slope systematically steeper than that of Marshall et al. (2001). For this reason Sambruna et al. can easily rule out the possibility that a unique spectral component from radio to X–rays can account from the overall spectrum, while the opposite conclusion is preferred by Marshall et al.

Despite this disagreement, both groups agree on the fact that the optical flux is due to synchrotron, and on the fact that maps at different frequencies yield knots of the same dimensions, and similar flux profiles. In other words, there appears to be no frequency–dependence of the size of knots.

However, according to the analysis given above, we would expect to find the radio knots to be larger than the optical ones, due to the longer lifetime of the radio-emitting electrons. An even larger region is expected in the case of X–ray knots produced by the IC process (i.e. by very low–energy electrons). Contrary to this expectation, it appears that knots have almost the same size at different frequencies.

In the following we examine what at first sight can be the possible and most plausible solutions:

* Consider a small and single standing shock. Suppose that the acceleration site is very compact, much smaller than the cross sectional radius of the jet, i.e. smaller than one arcsec. The angular resolution of Chandra in this case is not enough to resolve it in X–rays. In this way the adiabatic losses are much more effective. They corresponds to take, as $`r_0`$ in our figure, a region which is much less than a kpc (say a factor 10 less). In this case there are enough “doubling radii” to let adiabatic losses to decrease the flux (as in our figures). But the angular resolution of HST allows to resolve the emission region in the optical, and show that it cannot be smaller than about 1 arcsec. This implies that the acceleration site is extended (at least 1 arcsec), or made by several sub–units, contrary to our assumption.
* Consider moving blobs. In this case we see a snapshot of the blobs which are now active, and moving. The emission at different frequencies can in this case be cospatial, with no emission between the knots. The knots are observed to be active simultaneously according to our (observing) time, implying that they have been lighted up at different intrinsic times (because of light travel time effects). If these knots were the only existing ones, this requires a strong fine tuning. On the other hand, if the jet has many knots that randomly (in time and in location) light up, we always have the possibility to observe a few of them “on” (thus avoiding the fine tuning problem), but automatically this requires the existence of many “fossils” of previous active phases, emitting some X–rays, some radio, and no optical \[i.e. with only the low energy part of $`N(\gamma )`$ surviving\]. Since we observe no “fossil” knot, we disfavor this possibility. Yet another alternative is to assume that the jet really has only a few knots, which are always active. In this case we do not observe fossils, but we require that the particles in each knot are continuously reaccelerated. However this leads to another problem: in fact if the same particles are continuously reaccelerated, they should reach a typical energy (at which the cooling and acceleration timescales are equal), causing a pile–up at some energy, which is not observed. A possible way out is to assume that particles are only episodically reaccelerated, but this would require a strong fine–tuning between the duty cycle of the acceleration episodes and the cooling time of the optical electrons (the two times must be comparable).

We conclude that simple alternatives to the basic scenario (basic in the sense of a single homogeneous region emitting by synchrotron and external Compton) do not work or require fine–tuning. We are then forced to explore a slightly more complex alternative. In the next section we will then consider the possibility of knots composed of several sub–units (which we may identify as the acceleration sites), in order to make expansion losses more effective.

## 4 The clumping scenario

Assume that the knot is not homogenous, but is instead formed by several smaller substructures, i.e. several clumps. The key point of this idea is that small clumps can expand much more than a unique knot and therefore adiabatic losses can play a more important role than in the homogeneous case discussed above.

Assume then several acceleration sites which are much smaller than the transverse dimension of the jet. Let us call them “mini–knots”. Assume a number $`N`$ of mini–knots, and for simplicity assume that all of them are the same, with the same density $`n_i`$, radius $`R_i`$, volume $`V_i`$ and bulk Lorentz factor $`\mathrm{\Gamma }_i`$. The filling factor of the mini–knots is $`fNV_i/V`$ where $`V`$ is the volume of the large knot. $`R_i`$ is limited to be smaller than – say – one tenth of the total size of the knot (to have enough expansion losses within the big knot). With these assumptions also the magnetic field $`B_i`$ must be the same, whatever the radius and number of mini–knots. This follows from the fact that in each mini–knot the magnetic field to CMB energy density ratio must be equal to the one-zone case, since the resulting spectrum must have the same synchrotron to IC luminosity ratio. Since $`U_{\mathrm{rad}}`$ in the mini-knot scenario is equal to the the one-zone scenario (as long as $`\mathrm{\Gamma }_i`$ is the same), also the magnetic field must assume the same value in both scenarios.

Just for illustrative purposes, let us compare the case of a single homogeneous model in which the size of the emitting region is the entire big knot (we may call it the one–zone case: this is the case discussed up to now in the literature) with the clumping case. Assume for simplicity that the total number of emitting particles is the same in the two cases, as well as the magnetic field and the bulk Lorentz factor. What varies is the density and the size of the mini–knots: $`n_i`$ must be a factor $`f`$ larger than in the one–zone case. With these assumptions, the average bulk kinetic power (dominated by the cold proton component) in the case of the mini–knots is the same as in the one–zone case: the total number of particles is the same, and the time needed for them to cross a section of the jet is the same in both scenarios.

There is no obvious lower limit to the size $`R_i`$ of the mini–knots, but it must be noted that (within our assumptions) for a fixed number $`N`$ of mini–knots, the smaller their size, the larger their density, and the larger the deviation from equipartition between particle and magnetic energy densities. This deviation scales as the ratio of the particle density in each mini–knot and the density in the one–zone case, and therefore it scales as $`f`$. Therefore it is possible to have many (and small) mini–knots with a small filling factor, but this implies a large deviation from equipartition.

With respect to the one–zone case, each mini–knot has an enhanced particle pressure, and it is thus likely that it is overpressured with respect to its surroundings, implying expansion.

In conclusion, we propose that there are several mini–knots within the larger blob, which are the sites of particle acceleration, which are denser and hotter than the surroundings, leading automatically to the expansion of the entire mini–knot. What we see is then the convolution of the radiation produced by several expanding mini–knots.

The formation mechanisms of these “mini–knots” is outside the main goal of this paper; here we would only like to mention possible instabilities, clouds (filaments) crossing the jet (as proposed by, e.g., Blandford & Konigl 1979, although they considered this process at much smaller scales), or instabilities triggered by entrained gas that can act as an obstacle for the flow, forming mini–shocks. Another possibility is that particle acceleration is induced by reconnection events in the jet (e.g. Drenkhan & Spruit 2002). In this case the mini–knots could be identified as the reconnection sites.

If such mini–knots exist, then the problem to let low energy electrons adiabatically cool before they reach the borders of the bigger knot as detected by Chandra is solved. To illustrate this, we show in Fig. 8 the SEDs corresponding to a single mini–knot and in Fig. 9 the corresponding flux profiles at different frequencies. Since $`r_0`$ is much smaller than in the previous cases, all the timescales are shorter, changing the relative importance of the radiative and adiabatic losses. Note that, in Fig. 9, the entire range of the x–axis is now 1 arcsec, as opposed to 10 arcsec of the previous ones (for a source at $`z1`$).

While cooling for radiative and adiabatic losses, and while the mini–knot expands, the electrons emit some radiation even outside the acceleration sites. High energy particles, producing the optical synchrotron radiation, are the ones cooling the most, making the optical flux to decrease at the faster rate initially. Then, when $`\gamma _{\mathrm{min}}`$ is low enough to produce optical radiation by the inverse Compton process with the CMB photons, we have an increasing optical flux (until it is produced by the rising initial tail of the IC spectrum). The same situation occurs (but later) for the infrared flux. The X–ray flux instead steadily decreases, being always produced by low energy electrons. In conclusion, having enough angular resolution, one should resolve each mini–knot, especially in the optical and infrared bands, where the contrast between the emission inside and outside the mini–knots is the largest. However, the flux between the knots never vanishes, and remains at a considerable level especially in the X–ray band.

Another point that must be considered is the nature of the medium between the “mini-knots”. If particle acceleration is concentrated in the “mini-knots”, the possibility that appears the most natural is that the medium is mainly composed by cold particles, perhaps embedded in a very low magnetic field. These particles would emit through IC/CMB: the emission is concentrated in the IR-optical regime, at $`10^{14}`$ Hz. In this case there would be the possibility that this material is detectable in a near future with NGST, apearing as a “diluting” emission between mini-knots. The possibility to detect the emission mainly depends on the density of the medium: even an upper limit could help to better characterize the density contrast between the clumps and the external medium.

## 5 Summary

We have discussed how the presence of low energy electrons in jets, predicted by the IC/CMB model, poses important problems related to their extremely long radiative cooling time. In fact we know that the X–ray emission (as well as the optical and radio emission) comes in localized regions of the jet (the knots), despite the fact that the X–ray emitting electrons cannot radiatively cool inside these small regions. This motivated the exploration of the possible effects of adiabatic losses that for low energy electrons are much more important than radiative losses. We have shown the predicted profiles at different frequencies under different assumptions.

In particular we have identified two possible cases: in case A we considered a standing (in the observer frame) shock, while in case B the emission we observe is produced in a moving blob. Some general consequences appropriate for future and more detailed observations will be briefly discussed below and in a forthcoming paper.

However, the existing observations already pose a problem. In fact comparing our results with the data of nearby aligned radio sources (i.e. nearby blazars, such as 3C 273), we found that in any case (A or B) radiative and adiabatic losses cannot explain what is observed. In case A we would expect radio (and perhaps X–ray) knots larger than the optical ones, contrary to what is observed. In case B we would expect a number of unobserved “fossil knots”, bright in the radio and especially in the X–ray bands, but much fainter in the optical.

A possible way out of this puzzle is to admit that the emitting region is clumped in several mini–knots, allowing particles to cool by adiabatic expansion much more than in the one–zone case. These clumps can originate from different processes. An important consequence of this scenario is that, since the magnetic field energy density has to be the same as in the case of the one-zone scenario, the emitting region is far from equipartition, with the particle energy density dominating over the magnetic one.

## 6 Discussion

In large scale jets, with moderate values of the energy densities involved (both magnetic and radiative), the radiative cooling times are very long for all but the TeV electrons responsible for the high energy tail of the synchrotron emission. In particular, in the IC/CMB models, the X–ray flux is produced by electrons with random Lorentz factors of the order of 10–100, implying that they never cool through radiation losses. Adiabatic losses must then play a crucial role to explain the knotty morphology of large scale X–ray/optical/radio jets.

What we see can be the result of particle acceleration occurring in two different scenarios: acceleration can occur either at the front of a standing shock or accelerated particles can be injected throughout a blob moving along with the jet.

In the first case we expect to observe an extension of the emission region downstream to the shock (at larger distances from the shock front), whose size depends on the observing frequency. In particular, since X–rays are produced by low-energy electrons, the size of the X–ray emitting region must be larger than the radio one. The optical behavior is more complex, being due initially to the highest energy and rapidly cooling electrons, and later being due to the lowest energy electrons through IC.

In the moving blob case, particles are confined and evolve within the blob, and the particle evolution translates into a time dependent appearance of the blob itself. In particular, also old, “fossil” blobs are expected to be observable, bright especially in the X–ray band, relatively less in the radio and optical. Even taking into account adiabatic losses, the lifetime of the X–ray emission is long (compared to $`R/c`$, where $`R`$ is the size of the blob), and therefore the number of these “fossil blobs” should be larger than the young ones.

These considerations point towards the possibility to discriminate between the two scenarios: observing ”X–ray trails” more elongated than the radio and the optical ones would favor the standing shock scenario, while observing isolated X–ray knots (with weaker radio and optical emission, but cospatial) would favor the moving blob scenario. In order to perform such discriminating observations one needs sub–arcsec angular resolution at flux levels of the order of $`10^{15}`$ erg cm<sup>-2</sup> s<sup>-1</sup> or better, at least in the optical/IR, possibly in the mm and radio domains. In general, NGST, LBT and ALMA are then required.

Nearby blazars can already give some indications, since for them HST already can probe the interesting spatial scales. We have compared our results with the existing data of 3C 273. What we found is at first sight puzzling, since the flux outside the knots decreases fast and no “fossils” are observed. We have then proposed that what appears a single knot is instead the convolution of several smaller emission sites. This allows particles to adiabatically cool within the observed knot, allowing the flux to decrease fast outside it.

A natural consequence of this clumping scenario is that the emission can be variable on short (months–years) timescales. We note that, very recently, the knot HST–1 in M87 (emitting by synchrotron from the radio to the X–ray band), has been claimed to vary on a timescale of months in the X–ray and optical bands, and not in the radio (Harris et al. 2002). These variations have been attributed to the radiative cooling of the synchrotron emitting electrons. However, the dimension of knot HST–1 is of the order of $`100`$ pc, implying a minimum observable variability timescale of $``$300 years, corresponding to the light crossing time. If the HST–1 blob is instead composed of several mini–knots, the minimum variability timescale becomes the light crossing time of each single mini–knot (as long as the total number of mini–knots is not large, otherwise the amplitude of variability is diluted and becomes unobservable). Also in M87’ – suggest rewriting this sentence as follows ‘M87, then, also shows some evidence for clumped substructure in the jet knots, as we have suggested in the case of 3C273 using a completely independent argument.

###### Acknowledgements.

We thank the Italian MIUR and ASI for financial support. We thank T. Belloni and D. Malesani for useful discussions.


text:
# Dynamic transition in driven vortices across the peak effect in superconductors

## Abstract

We study the zero-temperature dynamic transition from the disordered flow to an ordered flow state in driven vortices in type-II superconductors. The transition current $`I_p`$ is marked by a sharp kink in the $`V(I)`$ characteristic with a concomitant large increase in the defect concentration. On increasing magnetic field $`B`$, the $`I_p(B)`$ follows the behaviour of the critical current $`I_c(B)`$. Specifically, in the peak effect regime $`I_p(B)`$ increases rapidly along with $`I_c`$. We also discuss the effect of varying disorder strength on $`I_p`$.

The model of driven interacting particles in presence of quenched disorder captures essential features of the dynamics occurring in many condensed matter system. Of specific interest is the dynamics of vortices above the critical (or depinning) current $`I_c`$ in type II superconductors. Early numerical simulation of the transport characteristics in 2D showed that the topological defects are generated close to $`I_c`$, and are annealed at high driving force shi . This suggested a dynamic transition from the disordered flow to an ordered flow at a current $`I_p>I_c`$. Theoretically, Koshelev and Vinokur (KV) considered this transition at $`I_p`$ similar to an equilibrium “melting” transition of a clean system in the moving reference frame kosh . The effect of quenched disorder on the moving vortices is included in the theory through an effective shaking temperature $`T_{sh}I^1`$. In Koshelev-Vinokur theory, the re-crystallization current $`I_p`$ increases as $`(T_mT)^1`$ on approaching the equilibrium melting temperature $`T_m`$ of the static system. Later calculations scheidl ; balents ; giamarchi  and simulations moon ; olson ; kolton ; fangohr  identified the free flowing state above $`I_p`$ as the smectic phase. In the smectic phase, the particles move in static channels with quasi-long range order perpendicular to the flow compared with liquid like short range order within the channels comment1 .

The dynamic transition in driven vortices was first studied experimentally in the transport measurements by Bhattacharya and Higgins (BH) shobo , and was later observed in neutron scattering yaron  and Bitter decoration duarte  experiments. BH studied the behaviour of transition current $`I_p`$ as a function of the magnetic field $`B`$, particularly in the peak effect (PE) regime. The PE is marked by a rapid increase in $`I_c(B)`$, generally close to the upper critical field $`B_{c2}`$. BH found that the $`I_p(B)`$ increases rapidly in the field range in which the PE occurs. This behaviour of $`I_p(B)`$ in the PE regime is similar to the behaviour of $`I_p(T)`$ around $`T_m`$ as observed in simulation by KV in Ref. kosh . The similarity can be traced to the softening of the shear modulus $`c_{66}`$ of the vortex lattice. But an important difference separates the two : unlike BH, KV do not see any corresponding increase in $`I_c(T)`$, i.e., no peak effect in the $`I_c(T)`$ is observed in the simulation. This makes it difficult to correlate the increase in $`I_p(T)`$ close to $`T_m`$ to the enhanced coupling of the vortex lattice to the quenched disorder. The increase in $`I_p(T)`$ could very well be due to increased thermal fluctuations on approaching $`T_m`$, thus requiring larger currents to anneal the lattice defects.

Recently, we showed that the PE occurs in a system of 2D vortices close to $`B_{c2}`$ at zero temperature chandran1 . In this paper, we revisit the dynamic transition, particularly across the PE. We find that $`I_p(B)`$ increases rapidly in the field range in which the $`I_c(B)`$ shows the PE. The dynamic transition at $`I_p`$ is characterized by hysteresis in $`V(I)`$ and a sharp peak in the dynamic resistance $`R_d(I)`$. The topological defect concentration also shows a jump at $`I_p`$. We present the dynamic phase $`I`$-$`B`$ diagram, and discuss the effect of increasing the disorder strength $`\mathrm{\Delta }`$ on $`I_p`$.

Consider a 2D cross-section of a bulk type-II superconductor perpendicular to the magnetic field $`𝐁=B\widehat{𝐳}`$. Within London’s approximation, the dynamics of vortices are governed by the overdamped equation of motion

$`\eta {\displaystyle \frac{d𝐫_i}{dt}}={\displaystyle \underset{ji}{}}U^v(𝐫_i𝐫_j){\displaystyle \underset{k}{}}U^p(𝐫_i𝐑_k)+𝐅_{ext}.`$
Here $`\eta `$ is the flux-flow viscosity, and the first term represents the inter-vortex interaction given by the potential $`U^v(r)=\frac{\varphi _0^2}{8\pi ^2\lambda ^2}K_0(\stackrel{~}{r}/\lambda )`$, where $`K_0`$ is the zeroth-order Bessel function, and $`\stackrel{~}{r}=(r^2+2\xi ^2)^{1/2}`$. $`\lambda `$ and $`\xi `$ are the penetration depth and coherence length of the superconductor, respectively. $`\varphi _0`$ represents the flux quantum. The vortex pinning is added through the second term which is an attractive interaction with parabolic potential wells $`U^p(r)=U_0(\frac{r^2}{r_p^2}1)`$ for $`r<r_p`$, and 0 otherwise, centered at the random $`𝐑_k`$ locations. The third term $`𝐅_{ext}=\frac{1}{c}𝐉\times \varphi _0\widehat{𝐳}`$ is the Lorentz force due to transport current density $`𝐉`$. The length is in units of $`\lambda (B=0)=\lambda _0`$, and $`J`$ is in units of $`cf_0/\varphi _0`$ where $`f_0=\frac{\varphi _0^2}{8\pi ^2\lambda _0^3}`$. The time $`t`$ is in units of $`\eta \lambda /f_0`$ whereas the velocity $`v`$ of vortices is in units of $`f_0/\eta `$.

We use the reduced magnetic field $`b=B/B_{c2}`$ with $`B_{c2}=\frac{\varphi _0}{2\pi \xi ^2}`$, and is calculated from the lattice constant $`\frac{a_0}{\lambda }=(\frac{4\pi }{\sqrt{3}})^{\frac{1}{2}}(\frac{1}{\kappa ^2b})^{\frac{1}{2}}`$. The $`B`$ dependence of $`\lambda `$ is follows $`\lambda (b)=\lambda _0/(1b^2)^{\frac{1}{2}}`$, with a similar expression for $`\xi `$. The simulation is for the Ginzburg-Landau parameter $`\kappa =\frac{\lambda }{\xi }=10`$ which is typical of the low-$`T_c`$ materials. The number of vortices $`N_v`$ was chosen in the range of 800-1200. The prefactor $`U_0`$ of the pinning potential is distributed randomly between $`\mathrm{\Delta }\pm 0.01`$. The results presented below are for the pin density $`n_p=2.315`$. We calculate the transport characteristics $`V(I)`$, where the current $`IJ_y`$, and the voltage $`Vv_x`$. The topological defect concentration $`f_d(I)`$ is obtained by Delaunay triangulation of the real space configuration. We choose the free flowing vortex lattice at high driving currents as the initial configuration, and decrease the current to 0 in small steps across $`I_c`$. This particular method of preparing the system minimizes the influence of non-equilibrium defects which are present in field cooled simulations. The equation of motion is time integrated by standard techniques. The $`V(I)`$ and $`f_d(I)`$ are time averaged in the steady state. Parallel algorithms were employed to speed up the simulation at high densities chandran2 .

Fig. 1(a) shows the $`V(I)`$ curve for $`b=0.7`$ and $`\mathrm{\Delta }=0.04`$. Broadly, three current regimes can be identified: (a) a pinned state for $`I<I_c`$, (b) a disordered flow state between $`I_c<I<I_p`$ where some vortices remain immobile and large transverse excursions of the active channels are present, and (c) an ordered flow state for $`I>I_p`$ with all vortices moving and interchannel hopping takes place with, at most, neighboring channels. The $`f_d(I)`$, and the dynamic resistance $`R_d(I)=\frac{dV}{dI}`$ is plotted in the inset of Fig. 1(a). The current $`I_p`$ is marked by a sharp peak in $`R_d(I)`$ and a large change in $`f_d`$. The response function $`V(I)`$ shows a kink at $`I_p`$, as evident from Fig. 1(a). We find that $`V(I)`$ and $`f_d(I)`$ are also hysteretic across $`I_p`$, as shown in Fig. 1(b). The kink and the hysteresis in the response function $`V(I)`$ suggests that the dynamic transition occurs at $`I_p`$. The transition is of first order nature, as was identified in Ref. kosh .

In Fig. 2, real space configurations are shown for currents above and below $`I_p`$ corresponding to parameters in Fig. 1. The instantaneous configuration shows a relatively ordered lattice above $`I_p`$. The vortex trajectories are well separated channels aligned parallel to $`v_x`$. The average transverse wandering of the channel $`\delta y_{chnl}<\alpha a_0^{}`$, where $`a_0^{}=\frac{\sqrt{3}}{2}a_0`$ and $`a_0`$ is the lattice constant, and $`\alpha 0.05`$. Just below $`I_p`$, $`\alpha `$ becomes $`0.3`$, as shown in the inset of Fig. 1(b). The slowing down of vortices effectively couples the transverse velocity component $`v_y`$ to the quenched disorder. This induces large scale dislocations, and the dynamic friction ($`R_d`$) increases. The increase in transverse wandering of vortices is also reflected in the dense braiding of the active channels. The behaviour of $`\delta y_{chnl}`$ across $`I_p`$ is similar to the Lindemann criterion for thermal melting. Note that some of the vortices appear as immobile for $`I=0.055`$ in Fig. 2(b) (rightside panel). The region of immobile vortices grow with decreasing $`I`$ until $`I_c`$ is reached when the whole system is pinned.

The $`R_d(I)`$ has been used previously to identify the phase boundary. In Ref. shobo , the transition from the disordered flow to the ordered flow state is identified as the current at which $`R_d(I)R_{BS}`$ where $`R_{BS}`$ is the asymptotic Bardeen-Stephen flux-flow resistance. This current is marked as $`I_{cr}`$ in Fig. 1(a). In Ref.kolton , the state between $`I_p`$ and $`I_{cr}`$ is identified as the smectic phase, whereas above $`I_{cr}`$, the state is defined as a transversely frozen phase. The later state is distinguished from the smectic phase by the absence of transverse jumps by vortices between adjacent channels.

We note that neither $`V(I)`$ nor $`f_d(I)`$ shows any feature at $`I_{cr}`$. Therefore, a more appropriate description of the state between $`I_p`$ and $`I_{cr}`$ is possible in terms of fluctuation which allows occasional transverse excursions of the vortices between adjacent channels. The rate of such transverse jumps can be calculated by employing Arrhenius relation which is conventionally used to describe thermal activation over energy barrier, but with temperature replaced by $`T_{sh}`$ chandran3 . Such a description can account for the suppression of transverse jumps above $`I_{cr}`$, and also for the experimentally observed $`R_d(I)`$ curves. To summarize, the dynamic transition in the moving lattice occurs at $`I_p`$ at which the dynamic resistance $`R_d(I)`$ shows a peak. The flowing state above $`I_p`$ is a smectic phase with small longitudinal correlation and quasi-long range order in the transverse direction.

We now consider the behaviour of $`I_p`$ as the magnetic field $`b=B/B_{c2}`$ is changed, particularly across the PE in $`I_c(b)`$. For $`\mathrm{\Delta }=0.04`$, the onset field for PE is $`b_{op}0.75`$ whereas the peak occurs at $`b_{peak}0.9`$. Fig. 3 shows the $`V(I)`$ curves for four values of $`b`$ between 0.5 and 0.8. The $`I_p`$ can be easily identified by the sharp kink in the $`V(I)`$ curves. With increasing $`b`$, the $`I_p`$ first decreases before increasing rapidly in the PE field regime. The $`f_d(I)`$ and $`R_d(I)`$ characteristics across $`I_p`$ are similar to that shown in Fig. 1 for all values of $`b`$. In the inset of Fig. 3, the effect of increasing disorder strength $`\mathrm{\Delta }`$ on $`I_p`$ is shown for $`b=0.6`$. Other than the expected increase in $`I_p`$ with increasing $`\mathrm{\Delta }`$, we find that the transition at $`I_p`$ is broadened for $`\mathrm{\Delta }0.05`$. Also, the hysteresis in $`V(I)`$ across $`I_p`$ is too small to be distinguished. This suggests that the first-order nature of the dynamic transition at $`I_p`$ does not survive for $`\mathrm{\Delta }0.05`$.

Fig. 4(a) shows the dynamic phase diagram $`I`$-$`B`$ for $`\mathrm{\Delta }=0.04`$ which summarizes the main result of the paper. The $`I_p(b)`$ and $`I_c(b)`$ shows similar behaviour as $`b`$ approaches the upper critical field value. Particularly, the rapid increase in $`I_p(b)`$ coincides with the peak effect in $`I_c(b)`$. The region between $`I_p(b)`$ and $`I_c(b)`$ constitutes the disordered flow regime, whereas above $`I_p`$ vortices flow in ordered channels. For $`\mathrm{\Delta }=0.04`$, the $`I_p(b)=k(b_0b)^1`$ with $`b_0=0.92`$ and $`k=0.015`$ gives a reasonable fit for $`b>b_{op}`$, as shown by the thick line in Fig. 4(a). This form of $`I_p(b)`$ was motivated by the behaviour of $`I_p(T)`$ on approaching the equilibrium melting transition in Ref. kosh . The $`I`$-$`B`$ plot in Fig. 4(a) thus supports the picture that both the static and the dynamic friction of the vortex system increases rapidly in the field range in which the peak effect occurs. As shown in Ref. chandran1 , the PE in the critical current $`I_c`$ is driven by the softening of the vortex interaction on approaching the upper critical field $`B_{c2}`$. Softening of the vortex interaction also reduces the shaking temperature $`T_{sh}I^1`$ required for the dynamic ordering (or recrystallization) of the vortices which explains the increase in $`I_p(b)`$ in the PE regime. Thus, an increase in $`I_c`$ (as a function of $`B`$ or $`T`$) implies an increase in $`I_p`$. We emphasize that the reverse does not hold, i.e., increase in $`I_p`$ does not imply an increase in $`I_c`$. This can be seen from the behaviour of $`I_p(T)`$ and $`I_c(T)`$ in Ref. kosh  where thermal depinning causes $`I_c(T)`$ to decrease monotonically even as $`I_p(T)`$ increases on approaching the equilibrium melting temperature $`T_m`$. Overall, we find the behaviour of $`I_p(B)`$ and $`I_c(B)`$ in good agreement with the experimental observation shobo .

Fig. 4(b) shows the $`I_c`$ and the $`I_p`$ as a function of $`\mathrm{\Delta }`$ for $`b=0.6`$. For the system size used in the simulation ($`N_v1000`$), the depinning transition is elastic for $`\mathrm{\Delta }0.03`$ (hence, $`I_p=I_c`$), whereas for $`\mathrm{\Delta }0.05`$, $`I_p\mathrm{\Delta }`$. The two disorder regimes can be identified as the weak pinning and the strong pinning regime, respectively. The crossover from weak to strong pinning occurs for $`\mathrm{\Delta }`$ between 0.03 and 0.05, and is also the range in which a sharp dynamic transition is observed. We find that for the intermediate pinning the real space configuration at $`I=0`$ show domains of ordered lattice separated by domain walls. The presence of ordered region ensures that the distribution of $`v_x\frac{1}{T_{sh}}`$ is narrow which consequently gives a sharp dynamic transition. On the other hand, increasing $`\mathrm{\Delta }`$ above 0.05 decreases the domain size to $``$ 2-3$`a_0`$. This leads to a broader distribution of $`v_x`$ which in turn broadens the dynamic transition. Thus, a sharp dynamic transition at $`I_p`$ implies the existence of ordered regions with a narrow distribution of the domain size. In Ref. paltiel , the coexistence of the ordered and the disordered regions was shown to underlie the shape of the $`V(I)`$ curves in the vicinity of the PE. The disordered state in Ref. paltiel  is thought to appear due to the injection of vortices across the surface barrier. Our simulation suggests that the shape of the $`V(I)`$ curve depends crucially on the distribution of the size of the ordered domains in equilibrium.

In conclusion, we have shown that for disordered type-II superconductors, the dynamic transition current $`I_p`$ and the critical current $`I_c`$ shows similar behaviour as the magnetic field $`B`$ is varied. $`I_p`$ decreases with $`B`$ for fields below the PE regime. In the peak effect regime, $`I_p(B)`$ increases rapidly along with the $`I_c(B)`$. The dynamic transition is sharp for the intermediate pinning strength due to the presence of large domains of ordered vortex lattice in equilibrium.

M.C. acknowledges useful discussions with E. Zeldov during the course of the work, and thanks the University of New Mexico for access to their Albuquerque High Performance Computing Center. This work was supported by NSF-DMR 9985978.