text:
# The First Detection of Cobalt In A Damped Lyman Alpha System

## 1 Introduction

Our understanding of the processes associated with galactic chemical enrichment, such as stellar nucleosynthesis, dust depletion and mixing mechanisms, can benefit greatly from the comparison of local abundances with those measured at high redshift. At low $`z`$, abundances can be determined for the Milky Way and its near neighbours from stars, H II regions and the interstellar medium (e.g. Chen et al. 1999; Kobulnicky and Zaritsky 1999; Shetrone et al. 2000; Izotov and Thuan 1999; Savage and Sembach 1996). At high redshift, the most accurate and extensive abundance measurements are those determined from Damped Lyman Alpha systems (DLAs). These systems are the highest H I column density members of the quasar absorption line menagerie and are traditionally defined as having N(H I) $`2\times 10^{20}`$ atoms cm-2. Coupling measurements of metal abundances in high redshift DLAs with our knowledge of local enrichment patterns plays a crucial role in our understanding of the star formation histories of these distant objects and the early stages in galaxy evolution.

In the context of investigating galactic chemical evolution, several key elements have previously been targeted as particularly informative diagnostics. The utility of these key metals is often based on our knowledge of their relative production and release timescales in stars (and subsequent supernovae) of different masses. For example, the $`\alpha `$ elements (Mg, Si, Ca etc.) have been extensively used to investigate the amount and rate of massive star formation. An overabundance of $`\alpha `$ relative to Fe in metal-poor Galactic halo stars chronicles the epoch of early star formation (Tinsley 1979; Ryan, Norris and Beers 1996) and in the bulge the same overabundances in solar metallicity stars are indicative of a ‘get rich quick’ population (McWilliam and Rich 1994; Rich and McWilliam 2000) i.e. one which is as old as the Galactic halo but as metal-rich as the Galactic disk. In DLAs, there is on-going discussion about the abundances of $`\alpha `$ elements, for which nucleosynthetic origins of abundance patterns must be disentangled from dust depletion effects. Although there is some evidence that Si may be overabundant in DLAs (Lu, Sargent and Barlow 1998), \[Si/Zn\] is solar in at least some systems (Pettini et al 2000). Moreover, since Si is highly refractory, the interpretation of \[Si/Zn\] measurements relies on a dust correction factor based on \[Zn/Cr\], which is quite uncertain. Conversely, the abundance of S, an $`\alpha `$ element not depleted onto dust, has been found by Centurion et al. (2000) to be underabundant relative to Fe, a result at odds with the Type II supernova pattern observed in the Milky Way. These results have been interpreted as possible evidence for star formation having occurred at a low rate or in sporadic bursts. However, the Centurion et al. (2000) data were obtained with resolution R = 5000 and need to be followed up with higher resolution data in order reliabley distinguish the S lines from the Lyman $`\alpha `$ forest. Further clues to the star formation history of DLAs may be gleaned from N/O ratios. The large observed spread in N/O in DLAs can be explained by the delayed release of primary N (relative to O) in metal-poor intermediate mass stars (Pettini et al. 2001; Lattanzio et al. 2001).

Here we focus on cobalt, an element which although relatively widely studied in Galactic stars, has not yet been featured in DLA abundance studies. The intriguing trends of \[Co/Fe\] as a function of \[Fe/H\] in the various Galactic stellar populations indicate that cobalt may be another element that can provide important clues to chemical evolution. In sub-solar metallicity thin disk stars, Co is found to be slightly underabundant with respect to Fe (\[Co/Fe\] $`0.1`$, Gratton and Sneden 1991), although this may be related to their solar log $`gf`$ scale being on average 0.06$`\pm `$0.04 higher than the laboratory scale adopted by many other observers (see Section 5.1 for more details). This trend continues in moderately metal-poor halo stars down to metallicities of \[Fe/H\] $`2.5`$. For more metal-poor halo stars, a considerable overabundance of Co with respect to Fe is observed, $`0<`$ \[Co/Fe\] $`<0.8`$, (McWilliam et al. 1995; Ryan, Norris and Beers 1996). Interestingly, this coincides with a downturn in relative Cr to Fe abundances, such that Co and Cr are closely anti-correlated in the metal poor population. Similar \[Co/Fe\] overabundances are observed in the bulge (\[Co/Fe\] $`+0.3`$, McWilliam & Rich 1994), which, we recall, is as old as the halo, and to a lesser degree in the thick disk (\[Co/Fe\] $`+0.1`$, Prochaska et al. 2000).

The detection of Co II transitions is a challenging possibility for ISM absorption line spectroscopy, both locally and at high $`z`$. Although there are numerous transitions that lie at convenient rest wavelengths for both local and high redshift studies, only two interstellar measurements currently exist. These detections have been achieved in the local ISM in the sightlines towards $`\zeta `$ Oph (Federman et al. 1993) and $`\rho `$ Oph A (Mullman 1998b) using the Co II $`\lambda `$1466 transition which has an equivalent width in these interstellar clouds of only 0.5 and 4.7 mÅ respectively. One of the advantages of moving to higher redshifts is that the relatively strong UV lines become shifted into the optical where one has access to more efficient detectors and larger telescopes. In addition, the observed equivalent width is increased by a factor of ($`z+1`$). However, the fundamental problem remains that the Co II lines are all intrinsically weak and that the solar abundance of Co is $`8\times 10^8`$ that of hydrogen. In addition, Mullman et al. (1998a,b) have provided new laboratory measurements for Co II $`f`$-values (oscillator strengths), typically revising the old values downwards by a factor of approximately two. This means that the spectral detection limit has effectively been increased by the same factor. Finally, the two local detections have brought to light an additional problem, namely that cobalt appears to be a highly refractory element that is depleted onto dust to a similarly high degree as iron (Mullman et al. 1998b).

Despite the technical challenge of observing Co in the high redshift ISM of DLAs, it is now feasible with high resolution spectrographs on 8 – 10 m class telescopes to detect this element, in reasonable exposure times, down to log $`N`$(Co) $``$ 12.5. Here we present the first ever detection of Co in a DLA, at $`z_{\mathrm{abs}}=1.92`$, determined from a spectrum obtained with the Ultra-violet and Visual Echelle Spectrograph (UVES) on UT2 of the VLT. In addition, we present a combined UVES and HIRES (High Resolution Echelle Spectrograph) spectrum obtained at the VLT and Keck telescopes respectively, which can also offer an interesting upper limit to \[Co/H\] in a second DLA. We discuss the implications for the star formation history of these high redshift galaxies and compare our observations with recent nucleosynthetic yield models.

## 2 Target Selection

In order to at least place interesting limits on the abundance of Co in DLAs, it is desirable to reach a sensitivity level of \[Co/Fe\] = 0. Such a detection limit is sufficient at least to distinguish between the marked overabundance of Co relative to Fe seen in the halo, bulge, and to a lesser degree in the thick disk, and the mild underabundance in the Galactic thin disk. Therefore, in order to maximise the likelihood of a Co II detection, we have selected bright QSOs with foreground DLAs known to contain high column densities of Fe. This criterion can either be achieved in moderate metallicity DLAs with high $`N`$(H I) or in lower column density systems with relatively high levels of enrichment. The former type of system is relatively rare due to the power law distribution of column densities in H I absorbers which shows that systems with $`N`$(H I) $`\stackrel{>}{}`$ 2 $`\times 10^{21}`$ cm-2 are very scarce. Similarly, few metal rich ($`Z>\frac{1}{10}Z_{}`$) DLAs have been identified at any redshift (Pettini et al. 1999), indicating that extinction due to dust may be biasing current samples (Ellison et al. 2000). Finally, as dust extinction may be greater in metal-rich and high-$`N`$(H) DLAs (e.g. Prantzos and Boissier 2000), it is particularly challenging to identify QSOs absorbed by such systems that also are bright enough to provide the requisite S/N.

Having searched the literature for candidate targets, we initially identified two QSOs for study, Q2206$``$199 ($`B=17.8,z_{em}=2.559`$) and Q1223+17 ($`B=18.5,z_{em}=2.936`$). The former of these targets actually exhibits 3 separate DLAs in its absorption spectrum (Prochaska and Wolfe 1997); here we focus on the system at $`z_{abs}=1.92`$. Although this system only has a moderate column density compared with other DLAs ($`N`$(H I) = 4.8$`\times 10^{20}`$ cm-2, see section 4), it is amongst the most metal-rich systems known with $`Z\frac{1}{3}Z_{}`$. In fact, Prochaska and Wolfe (1997) were able to detect some elements rarely seen in DLAs in this system, such as Ti. Conversely, Q1223+17 is a very high column density absorber, $`N`$(H I) = 3$`\times 10^{21}`$ cm-2 and therefore has a relatively high column density of Fe despite its low metallicity (Pettini et al. 1994; Prochaska et al. 2001).

## 3 Observations

We present high resolution echelle spectra for these two high redshift QSOs. The data were obtained with UVES on the VLT (Q2206$``$199 and Q1223+17) and with HIRES on the Keck telescope (Q1223+17). A summary of the observations is given in Table 1.

### 3.1 Details of Observations and Data Reduction of Q2206$``$199

Spectra of Q2206$``$199 were obtained over 3 nights from May 28-30 2000 with UVES on the VLT in conditions of good seeing, typically 0.5 – 0.7 arcseconds, although for brief periods extremes of 0.4 and 1.3 arcsecs were reached. A slit width of 1.0 arcsec was fixed throughout the observations, oriented along the parallactic angle. The data presented here were obtained with the dichroic 1 $`\lambda _{cen}=390+564`$ configuration and read out in 2 $`\times `$ 2 pixel-binned, high-gain mode. The blue arm CCD is a single thinned EEV chip whereas the red arm consists of a mosaic of a second EEV chip and an MIT chip, separated by approximately 0.96 mm. A single dichroic setting therefore has non-contiguous wavelength coverage and in the $`\lambda _{cen}=390+564`$ setting there are gaps between 4525 – 4625 and 5600 – 5680 Å.

The extraction of the spectra was achieved using a customized version of the UVES pipeline which is based on ECHELLE routines in the data reduction package MIDAS. A detailed description of this process can be found in Ballester et al. (2000). The spectra were optimally extracted separately for the three CCDs and the variation in the continuum flux level, due primarily to the blaze function of the echelle grating, was removed. The resolution of the spectra, as determined from ThAr arc lines, is $`R\mathrm{43\hspace{0.17em}000}`$, or approximately 7 km s-1 FWHM. 1-D spectra were obtained by joining the orders and converting the wavelengths to a vacuum heliocentric scale. The individual spectra were then co-added with a weighting proportional to their S/N. The final step was to normalise the spectrum by dividing through by a spline function fitted through absorption-free regions of data. Since most of the metal lines which we study here are found redward of the Lyman $`\alpha `$ emission, the continuum can be accurately fitted since there are relatively few absorption features and the S/N is quite high.

### 3.2 Details of Observations and Data Reduction of Q1223+17

The spectrum presented here of the $`z_{abs}=2.466`$ DLA towards Q1223+17 is a composite of HIRES plus UVES data. The UVES data were obtained during the same run as the Q2206$``$199 spectra and reduced in an identical way. The only difference for the data of Q1223+17 is that they were obtained with the dichroic 2 $`\lambda _{cen}=437+860`$ setting. This mode offers a larger, redder overall wavelength coverage at the expense of slightly larger inter-arm gaps (between 5000 – 6750 Å and 8500 – 8750 Å). In addition to the UVES data, we have obtained $``$ 5 hours of Keck HIRES data, which are described in more detail by the comprehensive abundance analysis of Prochaska et al. (2001). In brief, the C5 decker was implemented providing FWHM $`8`$ km s-1 resolution with the setup covering $`\lambda 47507200`$ Å. The data were reduced with the MAKEE software package developed by T. Barlow for the extraction of HIRES spectra. The coadded spectrum has S/N $`30`$ per 2 km s-1 pixel.

To obtain the strictest upper limit on the Co II column density for the damped Lyman $`\alpha `$ system towards Q1223+17, we coadded the HIRES and UVES spectra in the region spanning the Co II $`\lambda `$2012 transition. Each spectrum was rebinned to a common wavelength scale and then optimally coadded resulting in a final spectrum with S/N$`36`$ per 2 km s-1 pixel.

## 4 Neutral Hydrogen And Cobalt Column Densities

The H I column density of the two QSOs was determined by fitting fully damped profiles to the normalised spectra using the Starlink package Dipso. The UVES spectrum of Q2206$``$199 was found to have a best fit $`N`$(H I) = 4.8$`\pm 0.3\times 10^{20}`$ cm-2 and a redshift of $`z_{abs}=1.9207`$, Figure 1. Note that this fit is very well constrained by the shape of the central trough and damping wings and the $`N`$(H I) is in good agreement with previously determined estimates (e.g. Pettini et al. 1994) for this DLA. The HIRES spectrum of Q1223+17 exhibits a relatively large column density of H I which is found to have a best fit of $`N`$(H I) = 3$`\pm 0.3\times 10^{21}`$ cm-2 and a redshift of $`z_{abs}=2.4661`$, also in good agreement with the Pettini et al. (1994) measurement.

A detailed analysis of other metal abundances in the DLAs towards Q2206$``$199 and Q1223+17 can be found in Prochaska and Wolfe (1997) and Prochaska et al. (2001); here we simply concentrate on the detection of Co II. However, the analysis of other metal lines plays an important role in determining whether lines are blended and distinguishing weak features from the noise. One of the standard procedures is to fit Voigt profiles to the absorption system which is described as a complex of ‘clouds’, each defined by a Doppler parameter ($`b`$-value), column density ($`N`$) and redshift. It has been unanimously found in other high resolution DLA studies that singly ionized species such as Fe II, Si II and Ni II can be fitted using identical cloud models (e.g. Pettini et al 1999). Therefore, line-fitting is often performed first for relatively strong lines with several observed transitions, such as Fe II, in order to fix the $`b`$-values and redshifts of the constituent components. For the other species, excellent profile fits can be achieved by allowing only the column density to vary as a free parameter, indicating that the singly ionized metal lines have similar kinematical structure. Fixing the cloud model therefore offers a significant advantage when fitting weak lines, especially if there is a risk of some contaminating absorption from other species at a different redshift, e.g. contamination in the Ly$`\alpha `$ forest.

### 4.1 Detection Of Cobalt in the DLA Towards Q2206$``$199

In the case of Q2206$``$199, the absorption system is quite complex in comparison with other DLA metal line systems. The fit was initially constrained using Fe II $`\lambda `$2249 and Fe II $`\lambda `$2260 and refined by performing a simultaneous fit with Ni II $`\lambda `$1750, Ni II $`\lambda `$1741 and Ni II $`\lambda `$1709. Together, these five transitions can provide excellent constraints on the cloud model and, in addition, are two elements that when considered along with Co may provide interesting clues to nucleosynthetic processes (see section 5). The $`b`$-values and redshifts of each cloud were tied within the fitting procedure such that these parameters could vary, but remained identical for the two elements. The best fit results for this procedure yielded a 6 component cloud model which extends over almost 150 km s-1, see Table 2. Having determined the cloud model for Fe II and Ni II, we use this template to fit Co II (and the other metal lines), again allowing only the column densities to vary.

The revised $`f`$-values of Co (Mullman et al. 1998a,b) show that the two strongest ground state resonance lines covered by this spectrum are at rest wavelengths of 1466.21 Å ($`f`$=0.031) and 2012.16 Å ($`f`$=0.037). A third line of similar strength is located at $`\lambda `$ = 1941.29 Å ($`f`$=0.034) but falls in the gap between the two red arm CCDs. The line at 1466 Å for a system at $`z_{abs}=1.92`$ lies in the Lyman $`\alpha `$ forest of the $`z_{em}=2.559`$ QSO. The Co II $`\lambda `$2012 line lies redward of the QSO’s Lyman $`\alpha `$ emission where the risk of mis-identification or contamination is much lower.

The $`\lambda `$2012 line was used to determine the $`N`$(Co) of individual components given in Table 2 which also lists the fits to the other Fe-peak elements covered by our spectrum. Although the column density error for any individual component can be quite significant, particularly for closely-spaced, blended lines, the error on the total N(X) is usually substantially smaller. Also, the error clearly depends on the wavelength of the transition (since the S/N varies along the spectrum) and is much reduced if several transitions are used in the fit. All of these factors are taken into account by VPFIT when determining the total column density errors from the chi-squared reduction. In addition, we consider the effect of the continuum fit on our results, which may be particularly important for Co since a large fraction of the absorption is associated with relatively broad, shallow features. We firstly note that the continuum fit was done with a low order polynomial across the flat continuum around the CoII $`\lambda `$2012 line and that any significant departure from this would require quite a contrived placement. Secondly, VPFIT may adjust the continuum in order to improve the fit, although no change was made in this case. As can be seen in Table 3, the relative amounts of Co in the red and blue components are in good agreement with the other elements, a sign that at least the shape of the continuum is correct. In terms of a systematic offset, in Figure 2 it can be seen that had the continuum been placed too low, there would be clear evidence of flux above the unity level. If the continuum had been placed too high, this would only increase $`N`$(Co), thus reinforcing the observed overabundance of \[Co/Fe\]. Nonetheless, we have investigated the possible systematic uncertainty induced by continuum placement around Co II and other elements and determine the effect on derived column densities to be not more than 0.10 dex, and typically significantly less than this.

Over-plotting the fit from Co II $`\lambda `$2012 onto the Co II $`\lambda `$1466 line (Figure 2) shows that although there is evidence for weak Co II absorption there is considerable blending, probably due to Lyman $`\alpha `$ forest lines. There is also some possibility of blending in the $`\lambda `$2012 transition, since it can be seen in Figure 2 that an additional component is present at the blue edge of the absorption complex that is not fitted by the cloud model we adopt (see the top panel of Figure 2 for a comparison with the Fe II $`\lambda `$2250 line). In order to investigate the possible impact of this component on our column density determination, we consider the relative contributions of the individual clouds to the total absorption seen in the other metal lines. Table 3 shows the column density of metal enriched gas in the blue ($`N(X)_b`$) and red ($`N(X)_r`$) components, where the blue component is defined to be made up of clouds 1 and 2 from Table 2, whereas the red component is taken to consist of clouds 3 – 6. Also listed in Table 3 is the fraction of gas in the blue component ($`f_b`$) for all of the metals measured in this DLA. From these statistics it can be seen that, excluding Co, on average 27% of the enriched gas resides in the blue component, a value which remains remarkably constant for all the ions considered here. For Co, however, we see that a slightly higher fraction of gas, 29%, resides in the blue component. This small increase may be caused by the wing of the contaminating line. We attempted to identify the possible cause of contamination, but were unable to associate the extra component with any metal line from either of the other 2 DLAs. To show that the effect of the blending is probably negligible, we assume that the fractional distribution of gas is the same for Co as for the other metals and that therefore the $`N_r`$(Co)=12.94 represents 73% of the total, implying a total $`N`$(Co) = 13.08, almost exactly what is determined from the original fit.

### 4.2 The Upper Limit of Cobalt in the DLA Towards Q1223+17

There is no Co II $`\lambda `$2012 absorption visible in the combined UVES+HIRES spectrum of Q1223+17, see Figure 2. To place an upper limit on the $`N`$(Co II) value and in turn the Co/Fe ratio, we integrated the apparent column density111$`N_a(v)=m_ec\tau _a(v)/\pi e^2f\lambda `$ where $`\tau _a(v)=\mathrm{ln}[I_i(v)/I_a(v)]`$, $`I_i`$ and $`I_a`$ are the incident and observed intensity, $`f`$ is the oscillator strength, and $`\lambda `$ is the rest wavelength of the transition (Savage & Sembach 1991). over a small region covering the Co II $`\lambda `$2012 transition. For high S/N data (i.e. in the limit that $`\mathrm{ln}(1x)x`$), this is equivalent to deriving the equivalent width over the same velocity range and then calculating the column density assuming the linear curve of growth. Although the stronger metal line transitions span over 200 km s-1 (e.g. Fe II 1608), to place the tightest constraints on the Co/Fe ratio we integrated from $`20`$ km/s to $`+20`$ km/s which covers the strongest Fe II feature in the velocity profile. Restricting our analysis to this velocity range, we find $`N`$(Co II)$`<12.63`$, $`N`$(Fe II)$`=15.05`$ and \[Co/Fe\] $`<+0.18`$ where the limits are at the $`3\sigma `$ level.

### 4.3 Corrections Due to Ionization and Dust?

It is usually assumed in DLA abundance studies that the column densities of the singly ionized metal species are a good approximation to the total column density of this element, based on observations and theoretical modeling (e.g. Viegas 1995). However, there has been some suggestion (Howk and Sembach 1999; Izotov et al. 2000) that ionization corrections may in fact be important when deriving DLA abundances. Although there is some temptation to invoke ionization corrections to explain a few unusual abundance trends in DLAs, the predictions of these models are often not borne out by observations (Pettini et al. 2000; Lattanzio et al 2001). In addition, Levshakov et al. (2000) have recently shown that the physical conditions adopted by these models are unrealistic when compared with those observed in local DLA analogs. Here, we are specifically concerned with the relative abundances of Co and Fe. Even in the unlikely event that some correction is relevant to the absolute abundances as suggested by Izotov et al. (2000), the first and second ionization stages of cobalt and iron have very similar potentials (e.g. 7.88 eV and 7.90 eV for Co I and Fe I) so that it is highly unlikely that an ionization correction is required when converting N(Co II) and N(Fe II) to \[Co/Fe\]. The same is true for nickel, since the ionization potential for Ni I is 7.64 eV.

The presence of dust in DLAs and its tendency to remove a fraction of the total metals from the gas phase has now been well documented (e.g. Pettini et al. 1994; Pettini et al. 1997). This renders the task of determining the true abundances in these systems somewhat uncertain, although several attempts have been made to recover intrinsic abundances by correcting for Galactic depletion patterns (Vladilo 1998; Savaglio 2000). In the case of Co, some of the more interesting Galactic (stellar) abundance patterns are seen not in absolute (e.g. \[Co/H\]) measurements, but in the values of relative abundances like \[Co/Fe\] and \[Co/Cr\]. The important consideration is therefore the relative depletion fractions of these elements, which can be determined from local interstellar observations. From Savage and Sembach (1996) and Mullman et al. (1998b) we can see that Co, Cr, Ni and Fe are all depleted to a very similar level. In addition, at relatively low levels of depletion, as is the case for Q2206$``$199, the small differences between the dust correction required for these elements becomes negligible (Pettini et al. 2000).

Therefore, we assume that we can safely investigate relative abundance ratios such as \[Co/Fe\] without the need for either corrections for dust or ionization.

## 5 Discussion Of Cobalt Abundances

In Table 4, we list the abundances of the iron peak elements (Fe, Ni, Cr, Zn, Co) measured for the two DLAs analysed herein and two other systems found in the literature for which useful \[Co/H\] upper limits can be determined.

### 5.1 Homogeneity of Stellar Cobalt Abundances

The cobalt abundances measured in the four DLAs will be compared in Section 5.2 with those for stars spanning a range of Galactic populations and metallicities (Sneden & Gratton 1991; McWilliam & Rich 1994; McWilliam et al. 1995; Ryan, Norris, & Beers 1996; Prochaska et al. 2000). The stellar measurements use solely Co I lines.

All stellar measurements utilised laboratory $`gf`$ values from Cardon et al. (1982), with the exception of Sneden & Gratton who derived solar $`gf`$ values. These two scales appear to be very similar; five lines in common give a mean difference log $`gf_{\mathrm{GS90}}`$ log $`gf_{\mathrm{CSSTW82}}`$ = +0.06 dex with a standard error in the mean of 0.04 dex. As this is significant only at the 1.5$`\sigma `$ level, it would be premature to adjust for it. Furthermore, by using solar $`gf`$ values, Sneden and Gratton have effectively analysed their stars differentially with respect to the Sun, so changing the $`gf`$ values would alter the inferred solar abundance as well. We note, nevertheless, that even in their most metal-rich stars, \[Co/Fe\] tends towards $`0.05`$ at \[Fe/H\] = 0, which supports the view that their Co abundances are marginally low. The problem (assuming there is one) may perhaps be traced to differences in the assumptions in the solar model used for their solar analysis and those employed in the other stellar calculations. We did, however, correct the three stars of Sneden & Gratton affected by the probable 0.47 dex error in the solar log $`gf`$ value of the 4118 Å line; see Norris, Ryan & Beers (2001) for details. A comparison with the improved $`gf`$ values of Nitz et al. (1999) confirms the accuracy of the Cardon et al. scale: log $`gf_{\mathrm{CSSTW82}}`$ log $`gf_{\mathrm{NKWL99}}`$ = $`0.02`$ dex and $`\sigma _{\mathrm{CSSTW82}}`$ $``$ 0.08 dex. We are thus confident that all stars are on a uniform, accurate $`gf`$ scale. This is quite an achievement since the stars span four orders of magnitude in abundance, and have been analysed by several independent groups.

As the DLA Co abundances are derived from ionised lines whereas the stellar transitions are neutral, there is (in principle) a possibility that the stellar and DLA abundance scales differ. However, as the $`gf`$ scales for both ionisation states are tied to modern lifetime measurements, we doubt that differences exceed 5-10%, or 0.02 – 0.04 dex.

Co I lines have extensive hyperfine structure which can affect abundance measurements by up to several $`\times 0.1`$ dex (e.g. Ryan et al. 1996, Fig. 1). The stellar studies cited have all accounted for hyperfine structure, which was not always the case in earlier analyses. Moreover, work by Pickering and colleagues (1998 and priv.comm.) have confirmed that the Co II lines we have observed in the DLAs will not be broadened sufficiently to affect our results.

Both Co and Fe are primarily singly ionised in the temperature range of the stars investigated, so the Co I measurements are of a minority state. However, because these elements have almost identical first ionisation potentials, uncertainties in their temperatures or surface gravities will affect the ionisation degree of Co I just as much as for Fe I. Representative errors in \[Co/Fe\] — see Ryan et al. 1996, Table 3 — appear to be $``$0.05-0.10 dex (1$`\sigma `$) for good S/N ($``$100) data, but can be $`0.2`$ dex for poorer data with S/N $``$ 30. It is believed that the lower S/N levels achieved in the first studies of the most metal-poor stars may explain most of the spread in \[Co/Fe\] seen in those stars. (See Norris et al. 2001 for much improved \[Co/Fe\] measurements in such objects.)

### 5.2 Comparison with Relative Stellar Abundances

The stellar data discussed in the previous section are plotted in Figure 3, as a function of \[Fe/H\]. As discussed in the Introduction, Galactic stellar \[Co/Fe\] ratios show distinct — but non-monotonic — trends as a function of metallicity, which indicate key differences between the various stellar populations. The halo and bulge data, for example, at high and low metallicity, share great age and the highest \[Co/Fe\] ratios, whereas later-forming stars in the disk have progressively lower \[Co/Fe\] ratios. In addition, the Co detection in Q2206$``$199 and three upper limits are also plotted, although for the DLAs \[Zn/H\] is used as a metallicity indicator rather than \[Fe/H\]. The first determination of \[Co/Fe\] = $`+0.31\pm 0.05`$ in a DLA is consistent with the overabundance seen in Galactic halo stars and the transition regime between the thick disk and bulge at a metallicity consistent with the latter populations. The point plotted for Q2206$``$199 is the value determined from our Voigt profile fit although, as discussed in section 4, this may be a small overestimate if there is some blending of the blue components. In addition, Co is also overabundant relative to Cr in Q2206$``$199, \[Co/Cr\] = +0.23, again consistent with thick disk, moderately metal poor halo and metal-rich bulge stars.

The trend of high \[Co/Fe\] in metal-poor stars was first pointed out by McWilliam et al. (1995) and confirmed by Ryan et al. (1996). The observation that stellar \[Co/Fe\] abundances are highest in the oldest populations might suggest that elevated abundances arise from the fastest evolving stars, i.e. the most massive supernova (SN) progenitors. However, halo star \[Co/Fe\] values are not uniformly high, and clearly vary as a function of metallicity. McWilliam et al. (1995) suspected that the trends were due to a metallicity dependence of the yields of elements produced in supernova (SN) nucleosynthesis, while Ryan et al. (1996) noted that differences in the explosion energy could affect the yields in this way. It is unclear whether the metallicity-dependence is secondary, being driven by lower Co/Fe yields from later-evolving, lower mass SN, or whether it reflects an explicit dependence on progenitor metallicity. More than one factor may be (and probably is) at work. Umeda et al. (2000) have shown that the abundances of these elements do not depend strongly on the initial stellar metallicity. Nakamura et al. (1999) have shown that the abundance trends of Cr, Co and Mn observed in the most metal-poor stars can be reproduced when changing the position of the mass cut, but at the expense of overproducing Ni, which is only rarely observed in the stars (Ryan et al. 1996). More recently, Nakamura et al. (2000) have also been able to reproduce the trends, though not the absolute values, in Co and Cr by using nucleosynthesis models of hypernovae, i.e. supernovae with very high explosion energies. These authors remark that as well as re-producing Galactic abundances, they can also explain the unusual patterns seen in some starburst galaxies, such as M82. Indeed, several other observational clues indicate that star formation in DLAs proceeds in short bursts separated by extended quiescent periods, such as low N/O and \[S,Si/Zn\] values (Pettini et al. 2001; Pettini et al 2000) and the lack of detected star formation through H$`\alpha `$ imaging down to levels of $`1`$ M yr-1 (Bouche et al. 2000; Kulkarni et al. 2001). However, both the models of Nakamura et al. (1999; 2000) also predict an overabundance of Ni relative to Fe, something that is rarely seen in stars, but quite commonly seen in DLAs (Figure 4). In fact all 4 of the DLAs plotted in Figure 3 exhibit \[Ni/Fe\] $``$ 0.2 222We note that Prochaska & Wolfe (1997) use the old NiII f-values, requiring a correction by a factor of 1.9, which brings their N(Ni) into very good agreement with ours.. However, since at least HE1104$``$1805 shows no strong overabundance of Co, it appears that neither stars nor DLAs exhibit consistent or comprehended relative abundances of Co and Ni. Finally, we note that since the revision of Ni II $`f`$-values, a systematic overabundance is now seen for \[Ni/Fe\] in DLAs, perhaps indicating that further refinement of the atomic data may still be required.

As detailed supernova nucleosynthesis models have yet to provide a convincing explanation of the \[Co/Fe\] abundances of stars, we would not presume to be able to solve that mystery with the first ever DLA observations. Nevertheless, by comparing the DLA data with the stellar measurements, useful insights can be obtained. As noted, the \[Co/Fe\] ratio in the $`z_{\mathrm{abs}}=1.92`$ DLA of Q2202$``$199 falls near the boundary between bulge and thick disk stars. It is seen at a lookback time of $``$8.5 Gyr, around the time of (or a few Gyr after) the formation of the oldest thick disk and bulge objects. However, the possible absence of an elevated \[Co/Fe\] ratio in Q1223+17, with a lookback time slightly ($`0.5`$ Gyr) longer, emphasises that age alone does not explain the elevated \[Co/Fe\] values. If massive stars do produce the elevated \[Co/Fe\] ratios, then possibly star formation started later in the $`z=1.92`$ DLA than in the others. This postulate could be further explored by considering the abundance of Si, an $`\alpha `$ element. We find \[Si/Zn\] =$`0.07`$, although this value is likely be be an underestimate of the true Si abundance due to dust depletion. It is possible to make a crude adjustment for the depletion fraction by considering Zn and Cr ratios, which are found to be solar in stars. Therefore, any significant departure from \[Zn/Cr\] = 0 can be explained by the incorporation of Cr into grains. For Q2206$``$199, we determine \[Zn/Cr\] = 0.35 dex, applying this correction to Si gives \[Si/Zn\] $``$ 0.3, although this is likely to be an upper limit since Si is less depleted than Cr in the local ISM (Savage and Sembach 1996). In any case, it is likely that Si exhibits at least a mild overabundance relative to Zn, consistent with the expectation if the enrichment of this galaxy’s ISM has been dominated by the products of Type II SN.

Nakamura et al. (1999, 2000) have achieved partial success in attempting to explain the iron-group element ratios (Ni excepted) in their nucleosynthesis models, firstly by moving the mass-cut and secondly by changing the explosion energy. The latter, as they note, partially simulates the former effect because the zone of complete Si-burning is shifted outward in high-energy models compared to classical 1051 erg models. However, increasing the energy has the additional beneficial outcome that it provides a solution, for the first time, to the long-standing problem of explaining why Ti is overabundant in halo stars. It is found that the greater radial extent over which Si-burning occurs in the high-energy models necessarily incorporates lower density zones, which allows for increased $`\alpha `$-rich freeze-out-production of 48Ti. Previous classical SN nucleosynthesis models always yielded \[Ti/Fe\] $``$ 0. If this finding may be considered independent evidence that higher-energy explosions are the norm in the more massive SN, i.e. those that enriched the extremely metal-poor halo and the ‘get-rich-quick’ bulge, then we should understand the higher \[Co/Fe\] and simultaneously the higher \[Ti/Fe\] ratios in these populations that were enriched by massive-stars. However, Prochaska & Wolfe (1997) find that \[Ti/Fe\] = $`0.07`$ in Q2206$``$199, although considering that Ti is slightly more depleted than Fe (Savage and Sembach 1996), this DLA is likely to have a roughly solar Ti/Fe ratio.

Once high-energy SN models are explored more widely and included into galactic chemical evolution calculations, and a larger range of elements including Co and possibly Ti have been measured in more DLAs, the complementary constraints presented by DLA and stellar systems can be considered jointly. Nevertheless, it is important to appreciate through all of this that the similarities of DLA \[Co/Fe\] and \[Co/Cr\] ratios to those of a given Galactic population do not necessarily imply that this DLA, or the DLA population in general, represents the precursor to modern day halos/bulges. Rather, it is indicative of similar physical properties such as star formation rate and mix of SN types between the populations.

## 6 Summary And Conclusions

The distinct trends exhibited by Co in the various Galactic stellar populations makes this an interesting element to target in DLA abundance studies. However, the weakness of the Co II lines renders this a challenging prospect, even for an 8-m class telescope such as the VLT. Nonetheless, we have shown that by selecting DLAs with high column densities of Fe towards bright QSOs, the detection of Co is feasible with instruments such as UVES and HIRES down to \[Co/Fe\] $``$ 0. We have achieved the first detection of \[Co/H\] = $`0.50`$ in one of the DLAs studied herein and a sensitive upper limit for a second system. In addition, we have determined useful upper limits for \[Co/H\] for two other DLAs from the literature. The detection of Co towards Q2206$``$199 shows an overabundance relative to Fe, \[Co/Fe\] = +0.31$`\pm 0.05`$, at a metallicity of approximately $`\frac{1}{3}Z_{}`$, consistent with both the ratios observed in Galactic bulge and thick disk stars. In fact, Prochaska et al (2000) noted that there is evidence from Galactic bulge/thick disk stellar abundances that these two populations may have shared a similar formation history. One current limitation of the comparison with bulge stars, however, is that the only measurements presently available in the literature were obtained with a 4-m telescope with relatively low ($`R17,000`$) resolution (McWilliam and Rich 1994). Repeating these observations with Keck/HIRES has shown that these abundances may undergo significant revision, for example \[Fe/H\] has already been found to be higher by $``$ 0.2 compared with the 4-m data of lower resolution and S/N (Rich and McWilliam 2000). The revised Co abundances from the HIRES data are not yet available, but will eventually provide an interesting comparison for our DLA detection. Conversely, the upper limit of \[Co/Fe\] $``$ +0.18 that we determine for Q1223+17 is not consistent with the bulge star pattern and is more in line with Galactic disk abundances, as is the result for HE1104$``$1805.

However, the interpretation of the data presented here are limited by statistics and we clearly can not presume to understand the star formation history of DLAs based on a single detection. In addition, Q2206$``$199 is already known to exhibit large abundances of other elements which are not normally seen in DLAs, such as Ti (Prochaska and Wolfe 1997). We can therefore not exclude the possibility that Q2206$``$199 exhibits peculiar abundances that may not be seen in other such systems. Finally, we have deliberately chosen DLAs with high $`N`$(Fe) for this initial study in order to maximise the probability of a Co detection. In order to extend this analysis, such a selection bias will eventually have to be overcome if we are to truly understand the trends of relative abundances. Nonetheless, this is an interesting new result and the study of Co in DLAs clearly invites further study.

Our understanding of the processes that drive abundance trends in stars has only truly become possible thanks to large datasets, since intrinsic scatter and errors will dominate small samples. Therefore, more DLA data points are clearly required in order to draw firm conclusions about the star formation histories implied by the observed abundance trends. Perhaps most importantly, we have shown that with current 8-m class telescopes and high resolution echelle spectrographs that such measurements are not only possible, but also that Co may be a rewarding element to study in the context of testing nuclesynthetic models and understanding galactic star formation histories.

## Acknowledgments

We are very grateful to our collaborators Max Pettini, Jacqueline Bergeron and Patrick Petitjean for consenting to release the data for Q2206$``$199, originally obtained for a different project. It is a pleasure to acknowledge the staff of the VLT and Keck for their assistance during the acquisition of these data. In particular we are grateful to Andreas Kaufer for his support during both the VLT observations and data reduction. We would like to thank Chris Howk and Emmanuel Jehin for helpful discussions and to Max Pettini and Norbert Christlieb for useful comments on an earlier draft. We are grateful to Sebastian Lopez for providing his Co upper limit in the DLA towards HE1104$``$1805.


format: text

text:
# A simple one-dimensional model of heat conduction which obeys Fourier’s law

## Abstract

We present the computer simulation results of a chain of hard point particles with alternating masses interacting on its extremes with two thermal baths at different temperatures. We found that the system obeys Fourier’s law at the thermodynamic limit. This result is against the actual belief that one dimensional systems with momentum conservative dynamics and nonzero pressure have infinite thermal conductivity. It seems that thermal resistivity occurs in our system due to a cooperative behavior in which light particles tend to absorb much more energy than the heavier ones.

J.B.J. Fourier developed his theory of heat conduction at the beginning of 19th century. It states (in actual notation) that the temperature profile of an isolated system will evolve following the conservation equation

$$c_v(T)\frac{}{t}T(\underset{¯}{r},t)=\underset{¯}{}\left[\kappa \underset{¯}{}T\right]$$
(1)
where $`T(\underset{¯}{r},t)`$ is the temperature measured by a probe at position $`\underset{¯}{r}`$ at time $`t`$, $`c_v(T)`$ is the specific heat per unit volume and $`\kappa `$ is the thermal conductivity. Fourier’s law may be applied, in particular, to a system in contact with two heat reservoirs at different temperatures placed at $`x=0`$ and $`x=L`$. In this case, the stationary state has the property of

$$J=\kappa (T)\frac{dT}{dx}=const.$$
(2)
where $`J`$ is the stationary heat flux through the system. Notice that one should assume that there is not mass transport and/or other mechanisms different than heat conduction. This law has been extensively tested in experiments in fluids and crystals. However, we do not understand yet many of its fundamental aspects (see for instance the review in ). In particular, its derivation from a microscopic Hamiltonian dynamics and the existence of a local equilibrium that gives meaning to the local temperatures are open problems.

Actually, heat transport in one dimensional systems is an interesting problem in the context of both nonlinear dynamics and statistical physics. Long ago, Peierls proposed a successful perturbative theory, based on a phonon scattering mechanism, in order to explain the thermal conductivity in solids. In particular, Peierls theory predicts that we do not expect a finite thermal conductivity in one-dimensional monoatomic lattices with interactions between nearest neighbors. Accordingly, the temperature profile of a chain of $`N`$ harmonic oscillators is flat and its thermal conductivity goes like $`\kappa (N)N`$ for large enough $`N`$. More generally, any integrable Hamiltonian system is expected to have such divergent conductivity because its associated normal modes behave like a gas of non-interacting particles carrying energy from the hot source to the cold one with no loss. On the other hand, there are non-integrable systems, to which the Peierls theory does not apply directly, whose behavior is known to agree with its prediction. For example, the thermal conductivity of the Fermi-Pasta-Ulam-$`\beta `$ model goes like $`\kappa (N)N^\alpha `$, with $`\alpha =0.43`$. The value of $`\alpha `$ can be predicted by using a mode coupling approximation to the corresponding interacting gas of phonons . An exception to these results are systems with translation invariant Hamiltonian that have zero pressure. For instance, a one dimensional chain of rotators show normal heat conduction .

In general, as an extension of the prediction from Peierls argument, it is presently believed that one shouldn’t expect in general a finite thermal conductivity in one-dimensional systems with momentum conserving interactions and nonzero pressure. The goal of this letter is to show a counterexample to the above belief. We introduce a system that, although its particle interaction conserves momentum and the pressure is non zero, the energy behavior has a diffusive character and Fourier’s law holds. Therefore, we think that in one dimensional systems with nonzero pressure, the conservation of momentum does not seem to be a key factor to find anomalous heat transport.

Let us introduce our model. In a line of length $`L`$, there are $`N`$ point particles of different masses interacting exclusively via elastic collisions. In order to minimize the finite size effects, the particles have only two different masses and they alternate along the line, i.e. $`m_{2l1}=1`$ and $`m_{2l}=(1+\sqrt{5})/2`$ with $`l=1,\mathrm{},N/2`$. We have chosen the masses of the even particles to be the most irrational number in order to minimize possible periodicities, resonances or non-ergodic behaviors. At the extremes of the line there are thermal reservoirs at fixed temperatures $`T_1=1`$ and $`T_2=2`$ at $`x=0`$ and $`x=L`$, respectively. We simulate the reservoirs by using the following process: each time particle $`1`$ ($`L`$) hits the boundary at $`x=0`$ ($`x=L`$) with velocity $`v`$, the particle is reflected with the velocity modulus

$$v^{}=\left[\frac{2}{m_{1(L)}\beta _{1(2)}}\mathrm{ln}\left(1e^{\frac{\beta _{1(2)}}{2}m_{1(L)}v^2}\right)\right]^{1/2}$$
(3)
where $`\beta _{1(2)}=1/T_{1(2)}`$. This reversible and deterministic map is due to H. van Beijeren (private communication). In order to check the influence of the type of reservoir into the system properties we have also used more conventional stochastic boundary conditions, but only different finite size effects and no other relevant behavior has been observed. For $`T_1T_2`$ there is a flow of energy from the high temperature reservoir to the low one, and the system then evolves to a non-equilibrium stationary state. A version of this model in which the masses are randomly placed was already studied by us . In this work, the system thermodynamic limit behavior was not considered but the local equilibrium property was demonstrated. Our goal in this letter is to check if the system has a finite thermal conductivity in the thermodynamic limit $`N,L\mathrm{}`$ with $`N/L=1`$. With this aim we performed a detailed numerical analysis along the following lines:

(1) The existence of a nontrivial thermal profile: We define the local temperature by measuring the mean kinetic energy of each particle and its mean position at the stationary state. We computed the profiles for $`N=50`$, $`100`$, $`500`$, $`1000`$, $`2000`$, with fixed $`N/L=1`$, $`T_1=1`$ and $`T_2=2`$. Fig.1 shows the local temperature as a function of $`x/N`$ (by seeking clearer figures, we have performed local averaging of the temperatures and positions to draw only $`100`$ points; no difference is found by drawing all the points). We see in Fig. 1 that the temperatures follow linear profiles in the interval $`x/N[0.4,0.6]`$ with slopes depending on the system size. This slope apparently tend to converge to unity but we find that the convergence is so slow that we cannot (numerically) exclude the possibility that the limiting temperature profile is nonlinear (which is not against Fourier’s law). In any case, a non-flat profile is clearly expected in the thermodynamic limit.

(2) The averaged heat current: If Fourier law holds and the heat conductivity is finite, the mean heat current, $`J=N^1_{i=1}^Nm_iv_i^3/2`$ should go to zero as $`1/N`$ whenever $`T_{1,2}`$ and $`N/L`$ are kept fix. The data does not give us a conclusive answer. In fact we fitted the experimental points ($`J`$ corresponding up to seven different $`N`$’s) to behaviors like $`J=aN^{0.71}`$, $`J=aN^1(1+bN^1)`$, $`J=aN^1(1+b\mathrm{ln}N)`$ and $`J=aN^1(1+b/\mathrm{ln}N)`$ all of them with regression parameters of order $`0.999`$. This reflects that the corrections to the leading order are dominant and that we are far from the asymptotic regime for the observable heat current. Therefore, the direct use of the Fourier Law $`\kappa =JN/(T_2T_1)`$ does not clarifies (from the numerical point of view) the existence of a finite heat conductivity in the thermodynamic limit.

(3) The current-current self correlation function: The heat conductivity is connected to the current-current self correlation function evaluated at equilibrium via its time integral (Green-Kubo formula). The integral has some meaning whenever the correlation function decays as $`t^{1\mathrm{\Delta }}`$ with $`\mathrm{\Delta }>0`$. We measured the current current correlation function with periodic boundary conditions and total momentum equal to zero. We find that the simulation has strong finite size effects. Nevertheless we can see a clear tail of order $`t^{1.3}`$ (see (a) in Fig. 2). In order to confirm such result we computed the local current-current correlation function that it has much better averaging properties. We show in Fig. 2 (b) the behavior of the logarithm of $`c(t)=j_i(0)j_i(t)`$ with $`j_i(t)=m_iv_i(t)^3/2`$, where the average is evaluated at the equilibrium state with $`T_1=T_2=1.5`$ and we average over all the particles and different initial states. In order to check that we are doing right we first computed $`c(t)`$ when the system has equal masses. In this case we know from Jepsen that the exact solution  behaves like $`c(t)t^3`$ for $`t`$ large enough. In Fig. 2 (b) we see how this behavior is obtained numerically and we also see that for the different masses case $`c(t)`$ decays as $`t^{1\mathrm{\Delta }}`$ where $`\mathrm{\Delta }`$ is again close to $`0.3`$. This implies that we can define a finite thermal conductivity via Green-Kubo. We think that the decay of correlations is so slow that it explains the strong finite size effects observed in the temperature profile and in the mean heat current. In fact we can argue that $`JN/(T_2T_1)=\kappa AN^\mathrm{\Delta }`$ which explains why we do not see a clear behavior of $`J`$ with $`N`$ with system sizes of order $`10^3`$ (the corrections are of order unity for those sizes).

(4) The energy diffusion: We also wanted to check if the dynamical version (see Eq.1) of Fourier’s law holds. With this aim, we prepared the system with zero energy (all particles at rest) and positions $`x(i)=i1/2`$, $`i=1,\mathrm{},N`$. Then, we give to the light particle $`i=N/2+1`$ a velocity chosen from a Maxwellian distribution with temperature $`T=1.5`$. We monitored how the energy flows through the system until any boundary particle moves. Finally, we average over many initial conditions. If the system follows the Fourier’s law we should see a diffusive type of behavior (if the thermal conductivity is constant). Figure 3 shows the energy distribution for $`N=100`$ and different times measured in units $`t_0=0.032`$. Let us remark here that to apply Eq.1 the temperature should have a smooth variation in the microscopic scale to guarantee that local equilibrium holds. In Figure 3 we see that, for times larger than $`t=200t_0`$, the average variation in the local temperature is of order $`0.001`$. Therefore, we may assume that we are in a regime where Eq.1 holds. Initially, the energy of the light particle is transfered to the neighbors very fast and then the particle stays very cold, much colder than its neighbors. In fact, in this initial regime, the energy maxima are moving outwards at constant velocity. This behavior ends at around $`t=100t_0`$. The system then begins to slow down and, at $`t300t_0`$, the structure of the energy distribution changes, and one can then differentiate the behavior corresponding to light particles and heavy ones at least around the maxima of the distribution. We measured the mean square displacement of the energy distribution at each time: $`s(t)=_n(n51)^2e(n,t)`$. We found that we can fit $`\mathrm{ln}s(t)=6.39(0.04)+2.05(0.01)\mathrm{ln}t`$ for $`t/t_0(30,100)`$, thus, a ballistic behavior that changes smoothly until for $`t>400t_0`$ were we find a diffusive behavior $`\mathrm{ln}s(t)=1.00(0.01)+1.005(0.002)\mathrm{ln}t`$. This last result confirms that our system follows even the dynamical aspects of Fourier’s law.

As we noticed above, in Fig. 3 we see that the light and heavy particles seem to follow different energy distributions, at least for times longer than $`t=300t_0`$. In order to get some more insight about such behavior, we computed the evolution of the total energy stored in the light (heavy) particles. The result is shown in Fig. 4 where we can detect five different time regions: (I) $`t/t_0(0,16)`$; only the light particle and the two heavy nearest neighbors have a nonzero velocity. (II) $`t/t_0(16,23)`$; the five central particles (three light and two heavy ones) are moving. The total energy stored in the light particles reaches a minimum. (III) $`t/t_0(23,233)`$; the heavy particles begin to release energy (on the average) until, at $`t=233t_0`$, both types of particles have the same amount of energy. (IV) $`t/t_0(233,600)`$ Light particles keep getting energy until we reach region. (V) $`t>600t_0`$; where the total energy stored in the light particles reaches a constant value that exceeds to the one corresponding to the heavier ones. Let us remark that, in the asymptotic regime $`t>600t_0`$, the energy distribution is still evolving and, therefore, this partition of energy between both degrees of freedom is an asymptotic dynamical property of the system.

In order to discard any non-ergodic behavior of our system we included reflecting boundary conditions at the extremes of the chain and we did much longer simulations. We saw that the isolated system tends to the equilibrium in which equipartition of energy between all degrees of freedom holds. That is, the total energy stored in the light particles is equal to the one stored in the heavy ones. Moreover, we have checked that the system at any stationary state (equilibrium or non-equilibrium) does not present the property of non-equipartition of the energy. We think that this non-equipartition of the energy between degrees of freedom is responsible for the normal thermal conductivity. In fact, we see that, around the distribution maxima, the particles arrange in the form that hot light particles are surrounded by cold heavy ones. The energy is then trapped and released in a diffusive way. But we also see that the release is diffusive when an enough large number of those hot-cold structures develop. Therefore we think that the mechanism for the thermal resistance is somehow cooperative.

In conclusion, Peierls arguments have successfully explained the observed thermal conductivity in solids by applying a perturbative scheme around the lattice harmonic interaction. The actual belief is that strong anharmonicity is not enough to guarantee a normal thermal conduction in one dimensional systems. Moreover, it has been proposed that the key lacking ingredient is that the dynamics of the system should not conserve linear momentum via the existence of local potentials through the line (think about particles attached to the one dimensional substrate through some kind of nonlinear springs). In such a way, local potentials should act as local energy reservoirs that slow down the energy flow. These properties, anharmonicity and non-conservation of momentum, are in some way the ones used on the original Peierls argument. We have shown a model that does not follow such clean picture. Although our one dimensional model is nonlinear and it conserves linear momentum (with non-zero pressure), we find that it follows Fourier’s law. We think that there are other cooperative mechanisms that can do the job of the local potentials. Maybe, systems having degrees of freedom that acquires easily energy but releases it in a very long times scale have, in general, normal thermal conductivity. In any case, we think that it is worth to explore such possibility.

ACKNOWLEDGMENTS- It is a pleasure to acknowledge J.L. Lebowitz, G. Gallavotti, R. Livi, F. Bonetto, V. Ricci, L. Rey-Bellet, J. Marro and E. Presutti for useful discussions and criticisms. This work has been partially supported by Junta de Andalucía, project FQM-165 and by the Ministerio de Educación: projects DGESEIC PB97-0842 and PB97-1080.


text:
# A semi-analytical approach to non-linear shock acceleration

## 1 Introduction

Diffusive shock acceleration is thought to be responsible for acceleration of cosmic rays in several astrophysical environments. Most of the observational evidence for this mechanism, also known as first order Fermi acceleration, has been provided by studies of heliospheric shocks, but there are indirect lines of evidence that acceleration occurs at other shocks. A particularly impressive example was provided a few years ago by the observation of gamma ray emission from the supernova remnant SN1006 . These observations could be interpreted as inverse Compton emission of very high energy electrons, accelerated at the shock on the rim of SN1006, though other radiation processes may contribute .

Shock acceleration has been studied carefully and a vast literature exists on the topic. Some recent excellent reviews have been written . Some of the more problematic aspects of the theory of particle acceleration at astrophysical shocks have been understood, while others are still subject of investigation.

One of the problems that are harder to face is the problem of the injection of particles in the acceleration region. Only particles with a Larmor radius larger than the thickness of the shock are actually able to feel the discontinuity at the shock. The shock thickness is of the order of the Larmor radius of thermal protons, so that only a small fraction of the particles can be accelerated.

The calculation of the injection efficiency is quite problematic, for several reasons: first, the distribution function of thermal particles is steeply decreasing with momentum, so that the number of accelerated particles changes wildly with changing injection momentum. Moreover the distribution of the particles in the shock frame is strongly anisotropic at these low momenta, which adds to the difficulty of obtaining a straight analytical answer. Second, the injection of particles from the thermal distribution and their subsequent diffusive transport is thought to be due to the scattering against plasma waves, which are likely to be excited by the particles themselves, which makes the problem intrinsically non-linear. This non-linearity is exacerbated by the backreaction of the particles on the structure of the shocked fluid . The only way to have a complete quantitative picture of the problem of shock acceleration is to use numerical simulations .

All these complicated effects, which seem to be important in several astrophysical situations, are nevertheless often neglected, mainly because of the lack of an approach that allows to take them into account without the use of complicated numerical simulations which are usually of restricted use. As a consequence, in most of the applications of the diffusive shock acceleration to astrophysical situations, the assumption of test particles is adopted, even in those cases where this approximation works poorly. Numerical simulations show however that even when the fraction of particles injected from the plasma is relatively small, the energy channelled into these few particles can be close to the kinetic energy of the unshocked fluid, making the test particle approach unsuitable. The most visible effect is on the spectrum of the accelerated particles, which shows a peculiar flattening at the highest energies, due to the backreaction of accelerated particles on the fluid. The consequences on the spectra of secondary particles and radiation processes are clear.

The need to have a theoretical understanding of the non-linear effects in particle acceleration fueled many efforts in finding some effective though simplified picture of the problem. The structure of shocked fluids with a backreaction of accelerated particles was investigated in  in a fluid approach. The thermodynamic quantities were calculated including the effects of cosmic rays, but the approach did not provide information on the spectral shape of the accelerated particles.

In Ref.  a perturbative approach was adopted, in which the pressure of accelerated particles was treated as a small perturbation. By construction this method provides an answer only for weakly modified shocks.

An alternative approach was proposed in . This approach is based on the assumption that the diffusion of the particles is sufficiently energy dependent that different parts of the fluid are affected by particles with different average energies. The way the calculations are carried out implies a sort of separate solution of the transport equation for subrelativistic and relativistic particles, so that the two spectra must be connected at $`pmc`$ a posteriori.

Recently, in , the effects of the non-linear backreaction of accelerated particles on the maximum achievable energy were investigated, together with the effects of geometry. The solution of the transport equation was written in  in an implicit form, and then expanded in terms of the unperturbed (linear) solution.

Recently, some analytical solutions were also presented for the non-linear shock acceleration, in the particular case of Böhm diffusion coefficient .

The need for a practical solution of the acceleration problem in the non-linear regime was recognized in , where a simple analytical approximation of the non-linear spectra was presented. In this model the spectrum of the accelerated particles was assumed to consist of a broken power law, with three slopes characterizing the low, intermediate and high energy regimes. The basic features of the spectra derived from numerical simulations were reproduced with this method.

In the present paper we propose an approach that puts together some of the elements introduced in ,  and  and provides a semi-analytical solution for the spectrum of accelerated particles and for the structure of the shocked fluid. The method proposed is of simple use, can be adapted to several situations and provides results in very good agreement with numerical simulations, and with simplified models as that in .

The paper is structured as follows: in §2 we describe the general problem of linear and non-linear shock acceleration; in §3 we explain in detail our approach to non-linear effects in shock acceleration. In §4 we discuss the results of our model and compare them with the predictions of the model in Ref.  and with the results of some numerical simulations. Our discussion and conclusions are presented in §5.

## 2 Shock acceleration: linear and non-linear

In this section we discuss the basic elements of shock acceleration and introduce our approach to the description of the non-linear effects due to the backreaction of the accelerated particles on the shocked fluid.

For simplicity we limit ourselves to the case of one-dimensional shocks, but the introduction of different geometrical effects is relatively simple, and in fact many of our conclusions will be not affected by geometry. The non-linear effects are restricted to the mutual action of the shocked fluid and the accelerated particles. In other words, the present work does not include self consistently the production and the absorption of plasma waves by the accelerated particles. This simplification is common to the approaches in  and .

The equation that describes the diffusive transport of particles in one dimension is

$$\frac{f}{t}=\frac{}{x}\left[D\frac{}{x}f\right]+u\frac{f}{x}+\frac{1}{3}\frac{u}{x}p\frac{f}{p}+Q,$$
(1)
where $`f(x,p)`$ is the distribution function, $`u`$ is the fluid velocity and $`D`$ is the diffusion coefficient. The injection of particles is assumed to occur only immediately upstream of the shock, so we write the source function as $`Q=Q_0(p)\delta (x)`$, where $`x=0`$ corresponds to the position of the shock front. For monoenergetic injection, the function $`Q_0(p)`$ has the following form:

$$Q_0(p)=\frac{N_{inj}u_1}{4\pi p_{inj}^2}\delta (pp_{inj}),$$
(2)
where $`p_{inj}`$ is the injection momentum and $`u_1`$ is the fluid velocity immediately upstream ($`u_1=u(0^+)`$). $`N_{inj}`$ is the number density of particles injected at the shock, parametrized here as $`N_{inj}=\eta N_{gas,1}`$, where $`N_{gas,1}`$ is the gas density at $`x=0^+`$. The boundary condition at the shock reads:

$$\frac{u_1u_2}{3}p\frac{f_0}{p}=\left(D\frac{f}{x}\right)_1\left(D\frac{f}{x}\right)_2+Q_0(p),$$
(3)
where $`u_2`$ is the fluid velocity downstream. Here we called $`f_0=f(0,p)`$ the distribution function at the shock position.

A useful way of handling eq. (1) was suggested in Ref.  (a similar approach was also adopted in Ref. ), and consists of integrating this equation in the variable $`x`$ from $`x=0^+`$ (upstream) to $`x=+\mathrm{}`$ (far upstream). After some simple algebraic steps, in which we make use of eq. (3), we obtain the following equation:

$$p\left(\frac{f_0}{p}\right)=\frac{3}{u_pu_2}\left\{f_0\left[u_p+\frac{1}{3}\frac{du_p}{d\mathrm{ln}p}\right]Q_0(p)\right\},$$
(4)
where the stationarity assumption was adopted and we assumed $`(Df/x)_2=0`$ . We have introduced the quantity:

$$u_p=u_1+\frac{1}{f_0(p)}_0^{\mathrm{}}𝑑x\left(\frac{du}{dx}\right)f(x,p),$$
(5)
and again we called $`u_1`$ and $`u_2`$ the fluid velocities at $`x=0^+`$ and $`x=0^{}`$ respectively. With this formalism the compression factor at the shock is $`R_{sub}=u_1/u_2`$. The function $`u_p`$, at each momentum $`p`$ has the meaning of average fluid velocity felt by a particle with momentum $`p`$ while diffusing upstream. Since the diffusion is in general $`p`$-dependent, particles with different energies will feel a different compression coefficient, and the correspondent local slope of the spectrum will be $`p`$-dependent. Note that, according to eq. (5), the velocity $`u_p`$ must be a monotonically increasing function of $`p`$.

The function $`u_p`$ describes the mutual interaction between the accelerated particles and the fluid. In other words, if we find the way of determining the function $`u_p`$, as we show later, we also determine the spectrum of the accelerated particles.

Eq. (5) is clearly an implicit definition of $`u_p`$, meaning that $`u_p`$ depends on the unknown function $`f`$. However, eq. (5) allows us to extract important physical information. For the monoenergetic injection in eq. (2) the solution can be implicitly written as

$$f_0(p)=\frac{N_{inj}q_s}{4\pi p_{inj}^3}\mathrm{exp}\left\{_{p_{inj}}^p\frac{dp}{p}\left[\frac{3u_p}{u_pu_2}+\frac{1}{u_pu_2}\frac{du_p}{d\mathrm{ln}p}\right]\right\},$$
(6)
where we put (for $`\gamma _g=5/3`$):

$$q_s=\frac{3R_{sub}}{R_{sub}1}.$$
(7)
Eq. (6) tells us that the spectrum of accelerated particles has a local slope given by

$$Q(p)=\frac{3u_p}{u_pu_2}\frac{1}{u_pu_2}\frac{du_p}{d\mathrm{ln}p}.$$
(8)
The problem of determining the spectrum of accelerated particles would then be solved if the relation between $`u_p`$ and $`p`$ is found 111In order to determine the spatial distribution of fluid velocity it is necessary to specify the exact diffusion coefficient as function of $`x`$ and $`p`$. In this paper we do need this information, therefore the choice of $`D`$ does not affect our conclusions.. This is the scope of the next section.

## 3 The gas dynamics of modified shocks

The velocity, density and thermodynamic properties of the fluid can be determined by the usual conservation equations, including now the pressure of the accelerated particles. We write these equations between a point far upstream ($`x=+\mathrm{}`$), where the fluid velocity is $`u_0`$ and the density is $`\rho _0=mN_{gas,0}`$, and a generic point where the fluid upstream velocity is $`u_p`$ (density $`\rho _p`$). The index $`p`$ will denote quantities measured at the point where the fluid velocity is $`u_p`$. We call this generic point $`x_p`$.

The mass conservation implies:

$$\rho _0u_0=\rho _pu_p.$$
(9)
Conservation of momentum reads:

$$\rho _0u_0^2+P_{g,0}=\rho _pu_p^2+P_{g,p}+P_{CR,p},$$
(10)
where $`P_{g,0}`$ and $`P_{g,1}`$ are the gas pressures at the point $`x=+\mathrm{}`$ and $`x=x_p`$ respectively, and $`P_{CR,p}`$ is the pressure in accelerated particles at the point $`x_p`$ (we used the symbol $`CR`$ to mean cosmic rays, to be interpreted here in a very broad sense). In writing eqs. (9) and (10) we implicitly assumed that the average velocity $`u_p`$ as defined in eq. (5) coincides with the fluid velocity at the point where the particles with momentum $`p`$ turn around to recross the shock.

Our basic assumption, similar to that used in , is that the diffusion is $`p`$-dependent and that therefore particles with larger momenta move farther away from the shock than lower momentum particles. This assumption is expected to describe what actually happens in the case of diffusion dependent on $`p`$. As a consequence, at each fixed $`x_p`$ only particles with momentum larger than $`p`$ are able to affect the fluid. Strictly speaking the validity of the assumption depends on how strongly the diffusion coefficient depends on the momentum $`p`$, but the results should not be critically affected by this assumption. Moreover, in case of strong shocks there are arguments that suggest that the strong turbulence excited by the shock should produce a Böhm diffusion coefficient, so that the dependence $`D(p)`$ on $`p`$ should be at least linear.

According to this assumption, only particles with momentum $`\stackrel{>}{}p`$ can reach the point $`x=x_p`$, therefore

$$P_{CR,p}=\frac{4\pi }{3}_p^{p_{max}}𝑑pp^3v(p)f(p),$$
(11)
where $`v(p)`$ is the velocity of particles whose momentum is $`p`$, and $`p_{max}`$ is the maximum momentum achievable in the specific situation under investigation. In realistic cases, $`p_{max}`$ is determined from geometry or from the duration of the shocked phase, or from the comparison between the time scales of acceleration and losses. Here we consider it as a parameter to be fixed a priori. From eq. (10) we can see that there is a maximum distance, corresponding to the propagation of particles with momentum $`p_{max}`$ such that at larger distances the fluid is unaffected by the accelerated particles and $`u_p=u_0`$.

We will show later that for strongly modified shocks the integral in eq. (11) is dominated by the region $`pp_{max}`$. This improves even more the validity of our approximation $`P_{CR,p}=P_{CR}(>p)`$. This also suggests that different choices for the diffusion coefficient $`D(p)`$ may affect the value of $`p_{max}`$, but at fixed $`p_{max}`$ the spectra of the accelerated particles should not be appreciably changed.

Assuming an adiabatic compression of the gas in the upstream region, we can write

$$P_{g,p}=P_{g,0}\left(\frac{\rho _p}{\rho _0}\right)^{\gamma _g}=P_{g,0}\left(\frac{u_0}{u_p}\right)^{\gamma _g},$$
(12)
where we used the conservation of mass, eq. (9). The gas pressure far upstream is $`P_{g,0}=\rho _0u_0^2/(\gamma _gM_0^2)`$, where $`\gamma _g`$ is the ratio of specific heats ($`\gamma _g=5/3`$ for an ideal gas) and $`M_0`$ is the fluid Mach number far upstream. Note that eq. (12) cannot be applied at the shock jump, where the adiabaticity condition is clearly violated.

We can rewrite eq. (10) in a convenient way, by dividing it by $`\rho _0u_0^2`$ and using the mass conservation. We then obtain:

$$1+\frac{1}{M_0^2\gamma _g}=U+\frac{1}{M_0^2\gamma _g}U^{\gamma _g}+\frac{4\pi }{3}\frac{1}{\rho _0u_0^2}_p^{p_{max}}𝑑p^{}p^3v(p^{})f_0(p^{})=$$

$$U+\frac{1}{M_0^2\gamma _g}U^{\gamma _g}+\frac{N_{inj}q_s}{\rho _0u_0^2p_{inj}^3}_p^{p_{max}}𝑑p^{}p^3v(p^{})\mathrm{exp}\left\{_{p_{inj}}^p^{}\frac{dp^{\prime \prime }}{p^{\prime \prime }}\left[\frac{3u_{p^{\prime \prime }}}{u_{p^{\prime \prime }}u_2}+\frac{1}{u_{p^{\prime \prime }}u_2}\frac{du_{p^{\prime \prime }}}{d\mathrm{ln}p^{\prime \prime }}\right]\right\},$$
(13)
where we have explicitely written the distribution function $`f(p)`$ as a power law with local slope $`Q(U)`$, and we put $`U=u_p/u_0`$.

Differentiating the previous equation with respect to $`p`$ we obtain:

$$\frac{dU}{dp}\left[1\frac{1}{M_0^2}U^{(\gamma _g+1)}\right]=\frac{1}{3}\frac{N_{inj}q_s}{\rho _0u_0^2}\left(\frac{p}{p_{inj}}\right)^3v(p)\mathrm{exp}\left\{_{p_{inj}}^p\frac{dp}{p}\left[\frac{3u_p}{u_pu_2}+\frac{1}{u_pu_2}\frac{du_p}{d\mathrm{ln}p}\right]\right\}$$
(14)
or

$$\frac{dU}{d\mathrm{ln}p}\left[1\frac{1}{M_0^2}U^{(\gamma _g+1)}\right]=$$

$$\frac{1}{3}\frac{N_{inj}q_s}{\rho _0u_0^2}v(p)p_{inj}\mathrm{exp}\left\{4\mathrm{ln}\left(\frac{p}{p_{inj}}\right)_{p_{inj}}^p\frac{dp}{p}\left[\frac{3u_p}{u_pu_2}+\frac{1}{u_pu_2}\frac{du_p}{d\mathrm{ln}p}\right]\right\}.$$
(15)
Note that the velocity $`u_p`$ changes as a consequence of the pressure added by non-thermal particles, therefore the function $`U(p)`$ must be a monotonically increasing function of the particle momentum. Since $`U(p)`$ and $`dU/d\mathrm{ln}p`$ are always non-zero, we can calculate their logarithm, so that eq. (15) becomes:

$$\mathrm{ln}𝒟U+\mathrm{ln}\left[1\frac{1}{M_0^2}U^{(\gamma _g+1)}\right]$$

$$=\mathrm{ln}\left[\frac{1}{3}\frac{N_{inj}q_s}{\rho _0u_0^2}v(p)p_{inj}\right]+4\mathrm{ln}\left(\frac{p}{p_{inj}}\right)_{p_{inj}}^p\frac{dp}{p}\left[\frac{3R_{tot}U}{R_{tot}U1}+\frac{R_{tot}}{R_{tot}U1}𝒟U\right],$$
(16)
where we put $`𝒟U=dU/d\mathrm{ln}p`$. The following equality is easily demonstrated:

$$_{p_{inj}}^p\frac{dp}{p}\frac{R_{tot}}{R_{tot}U1}𝒟U=\mathrm{ln}\left(\frac{R_{tot}U1}{R_{sub}1}\right),$$
(17)
so that the equation for $`𝒟U`$ becomes:

$$\mathrm{ln}𝒟U+\mathrm{ln}\left[1\frac{1}{M_0^2}U^{(\gamma _g+1)}\right]$$

$$=\mathrm{ln}\left[\frac{1}{3}\frac{N_{inj}q_s}{\rho _0u_0^2}v(p)p_{inj}\right]+4\mathrm{ln}\left(\frac{p}{p_{inj}}\right)\mathrm{ln}\left(\frac{R_{tot}U1}{R_{sub}1}\right)_{p_{inj}}^p\frac{dp}{p}\frac{3R_{tot}U}{R_{tot}U1}$$
(18)
Solving this differential equation provides $`U(p)`$ and therefore the spectrum of accelerated particles, through eq. (6).

The operative procedure for the calculation of the spectrum of accelerated particles is simple: we fix the boundary condition at $`p=p_{inj}`$ such that $`U(p_{inj})=u_1/u_0`$ for some value of $`u_1`$ (fluid velocity at $`x=0^+`$). The evolution of $`U`$ as a function of $`p`$ is determined by eq. (18). The physical solution must have $`U(p_{max})=1`$ because at $`p\stackrel{>}{}p_{max}`$ there are no accelerated particles to contribute any pressure. There is a unique value of $`u_1`$ for which the fluid velocity at the prescribed maximum momentum $`p_{max}`$ is $`u_{p_{max}}=u_0`$ (or equivalently $`U(p_{max})=1`$). Finding this value of $`u_1`$ completely solves the problem, since eq. (18) provides $`U(p)`$ and therefore the spectrum of accelerated particles, calculated according to eq. (4). Conservation of energy can be easily checked.

An illustration of this procedure is presented in fig. 1 where we considered the following special set of parameters: the Mach number far upstream is $`M_0=43`$, and the gas temperature is $`T_0=10^6`$ K, corresponding to $`u_0=5\times 10^8`$ cm/s. The parameter $`\eta `$ is taken to be $`4\times 10^3`$, and the injection of particles occurs at $`p_{inj}=10^2mc`$, where $`m`$ is the mass of the accelerated particles. In this paper we assume that the accelerated particles are protons. The maximum momentum is $`p_{max}=10^5mc`$. The three curves which are plotted correspond to the function $`U(p)`$ for $`R_{sub}=2`$ (dotted line), $`R_{sub}=1.5`$ (dashed line) and $`R_{sub}=1.817`$ (solid line). It is immediately evident that only for $`R_{sub}=1.817`$ the function $`U(p)`$ reaches unity at $`p=p_{max}`$. Our method then provides the value of $`R_{sub}`$ and consequently the values of the other parameters.

Fig. 1 is very useful for understanding the physical meaning of the local velocity $`u_p`$. Particles with large momenta feel compression factors $`u_p/u_2`$ which are larger than those felt by low momentum particles. Large compression factors correspond to locally flatter spectra, so that the spectrum of accelerated particles is expected to become flatter at large momenta. We define the total compression factor $`R_{tot}=u_0/u_2`$. Note that the compression factor at the gas subshock is now $`R_{sub}4`$, the value expected for a strong shock in the test particle regime.

## 4 Results

The calculations illustrated in the previous section are here tested versus previous models and numerical simulations. In particular we compare our results with the predictions of the simple approach presented in , that we briefly summarize. In  the spectrum of accelerated particles has a prescribed shape, made of three power laws, in the low ($`p_{inj}pmc`$), intermediate ($`mcp10^2p_{max}`$) and high ($`10^2p_{max}pp_{max}`$) energy regimes. The slope in the three regions is then calculated by requiring mass, momentum and energy conservation.

We first check that our model reproduces the results in the linear regime, where the test particle approximation can be adopted, and then we study the transition to the non-linear regime. For this purpose, we consider a shock with Mach number $`M_0=5`$, with a gas temperature $`T_0=10^8`$ K. We choose $`p_{inj}=10^2mc`$ and $`p_{max}=10^5mc`$, and we study the result for different values of $`\eta `$, as plotted in Fig. 2. The solid line is obtained for $`\eta =10^5`$; our model gives for this case $`R_{sub}=R_{tot}=3.57`$, which is exactly the value obtained from test particle theory: $`R_{sub}=(8M_0/3)/((2/3)M_0+2)`$, for $`\gamma _g=5/3`$. The test particle approximation would provide the same compression factor for any value of $`\eta `$. Our results for $`\eta =10^3`$ and $`\eta =10^2`$ are plotted in Fig. 2 as dotted and dashed lines respectively. These two cases result in $`(R_{sub},R_{tot})=(3.45,3.83)`$ for $`\eta =10^3`$ and $`(R_{sub},R_{tot})=(2.85,4.49)`$ for $`\eta =10^2`$. The transition from unmodified shocks to strongly modified shocks is evident.

In all the cases the common feature of the modified shocks is a steepening of the spectrum (in comparison with the linear result) at low energy and a flattening at high energies. From the phenomenological point of view this can be of paramount importance since it may change the spectral features of the secondary radiation produced by the interactions of the accelerated particles.

The more interesting case is that of strongly modified shocks, where we expect the pressure in accelerated particles to become comparable with the kinetic energy of the upstream fluid. Intuitively this is mainly the case for high Mach number shocks, though we also show that there are exceptions. We can compare our results with those of the model in Ref. .

Let us consider the case $`M_0=43`$, $`u_0=5\times 10^8`$ cm/s, $`p_{inj}=10^2mc`$, and let us study the resulting spectra for $`\eta =10^3`$ and different values of $`p_{max}`$. The results are plotted in Fig. 3.

The results of  are plotted in the form of crosses ($`p_{max}=10^3mc`$), stars ($`p_{max}=10^5mc`$) and diamonds ($`p_{max}=10^7mc`$). The corresponding continuous lines are the results of our model.

Some comments are in order: the slope of the spectrum of accelerated particles predicted in our model at $`pp_{max}`$ is approximately equal to that obtained in . This is not surprising since the model in  is based on a three power law approximation, and the slope of the spectrum at the highest energies is calculated using an asymptotic expression derived from eq. (6). At low values of $`p_{max}`$ our results are in very good agreement with the general features of the solutions in . At increasingly larger values of $`p_{max}`$ the agreement is not extremely good but still reasonable, if one takes into account that the following three assumptions were made in : 1) the velocity of the particles is assumed to be $`p/m`$ for $`pmc`$ and equal to $`c`$ for $`pmc`$; 2) the momentum at which there is a change in slope at low energy is forced to be at $`p=mc`$; 3) the position of the point where there is the change in slope at intermediate energies is forced to be at $`p=10^2p_{max}`$. In our model all the three assumptions are released and the spectra are smooth. This is the reason for the slightly different position of the dip in our spectra ($`p^4f(p)`$) when compared with those derived according to . Similar conclusions hold for different values of the parameters.

A more interesting comparison, to test the effectiveness of our model is that with the results of numerical simulations. Since a similar comparison was carried out in  to test that model, we consider here the same situation, so that a full cross-check is possible. The case we consider is that of a shock with Mach number $`M_0=128`$ and a fluid temperature of $`T=10^6`$ K. The results provided by a numerical simulation are plotted in Fig. 5 of Ref.  and reproduced in our Fig. 4 as a dashed line. From the simulation we can extract the values of some parameters: the injection momentum is $`p_{inj}=7\times 10^3mc`$, while the maximum momentum is $`p_{max}=10^5mc`$. Note that the simulation does not provide these parameters in a clear way, because there the thermal and non-thermal particles are treated in the same way, therefore some approximation is involved in deriving these numerical values. The best fit of the model in  to the results of the simulation implies $`R_{sub}=2.68`$ and $`R_{tot}=52`$ using $`\eta =5\times 10^3`$. Using the same value of the injection efficiency $`\eta `$, our model predicts $`R_{sub}=2.365`$ and $`R_{tot}=51.9`$. The spectra for the model in  and for our model are plotted in Fig. 4 as a solid light broken line, and as a thick solid line respectively. The improvement provided by our model is evident (note that we only used the values of $`\eta `$ and $`p_{inj}`$ derived in , though an even better agreement can be found by slightly changing these values, that, as we stressed above, are not clearly determined by the simulation).

A basic issue raised in  is that of the accurate description of the slope of the spectrum of accelerated particles at $`pp_{max}`$. The importance of this point is due to the flat spectral shape at high momenta, which gives the main contribution to the total energy budget in accelerated particles, for strongly modified shocks. We show our prediction for the local slope in Fig. 5, as a function of the momentum. The most noticeable feature of this figure is the slope $`3.5`$ at $`pp_{max}`$, very close to that predicted by numerical simulations.

One comment is in order concerning the injection: in all the cases considered above, the injection momentum and $`\eta `$, the injection efficiency, are chosen independently. In a realistic case (for instance in the simulation) the two parameters would be actually related to each other. For instance, if the particles are extracted from the thermal distribution, then one can write in a very general way, that the injection momentum is proportional to the sound speed downstream ($`c_{s,2}`$), times the particle mass, $`p_{inj}=\xi mc_{s,2}`$. The parameter $`\xi `$ must be large enough that the particles can feel the thickness of the shock, determined either by the interaction pathlength of the particles, or by the Larmor radius of the thermal protons in the local magnetic field. Hence, $`\xi `$ must be at least a few, and indeed it is usually assumed to be in the range $`\xi 410`$. If the expression $`p_{inj}=\xi mc_{s,2}`$ is adopted, and the injection is assumed to occur from the thermal distribution, then the parameter $`\eta `$ is no longer free, and can be determined self-consistently. However it is wise to keep in mind that in this case $`\eta `$ would change wildly for small changes in $`\xi `$, due to the exponential suppression in the distribution function at momenta larger than the thermal average.

We continue now the analysis of the predictions of our model and the comparison with the results of Ref. . An important issue is that of determining where the transition from unmodified to strongly modified shock occurs, as a function of the parameters of the calculation.

In Fig. 6 we plot the compression ratios (total and at the subshock) for different values of the maximum momentum of the accelerated particles. The solid lines represent the results of our model, compared with the values predicted by the simple approach in  (dashed lines). There is a good agreement between the two approaches. For $`\eta =10^3`$ the shock becomes strongly modified already for $`p_{max}`$ a few times larger than $`mc`$, confirming that even relatively low injection efficiencies cause the shock to be affected by the backreaction of the accelerated particles.

In Fig. 7 the same compression factors are plotted versus the minimum (injection) momentum $`p_{inj}`$. Here some differences between our model and that in Ref.  are visible. The results can be interpreted as an evidence that our model predicts the onset of the modified shock phase, at slightly larger values of $`p_{inj}`$ compared to the simple model.

The results of a similar study are plotted in Fig. 8, where the dependence of $`R_{tot}`$ and $`R_{sub}`$ is investigated versus the injection efficiency $`\eta `$. Again, our model predicts that the shock starts to be modified at values of $`\eta `$ slightly larger than for the simplified model. In particular for the values of the parameters reported in Fig. 8, we obtain that the shock starts to be modified at $`\eta \stackrel{>}{}10^4`$, while the simple model gives a modified shock for $`\eta \stackrel{>}{}3\times 10^5`$.

Finally in Fig. 9 we investigate the dependence of $`R_{tot}`$ and $`R_{sub}`$ on the Mach number of the fluid at infinity (upstream). We specialize our prediction to the case $`\eta =10^3`$, $`p_{inj}=10^2mc`$ and $`p_{max}=10^5mc`$, but clearly similar plots can be produced for different regions of the parameter space. The results of our model are in very good agreement with the results of . Our model however predicts a slightly lower value of the critical Mach number $`M_{cr}`$, above which the shock is no longer modified, so that $`R_{tot}`$ and $`R_{sub}`$ settle down at the usual (linear) value of $`4`$. As found in , the relation between $`R_{tot}`$ and $`M_0`$ is $`R_{tot}M_0^{3/4}`$, for $`M_0M_{cr}`$.

## 5 Discussion and Conclusions

We proposed a semi-analytical method to calculate the spectrum of particles accelerated at shocks, and the velocity profile of the shocked fluid, including the non-linear effects due to the backreaction of the non-thermal particles on the dynamics of the shocked fluid.

The approach presented in this paper provides results that are in good agreement with the output of numerical simulations on shock acceleration and also with the results of a previous model aimed to a simple description of the non-linear effects . Compared with the latter, our model provides a better fit to numerical results, and is not based on a priori assumptions on the spectrum of accelerated particles in some energy ranges.

We confirm the main results on shock acceleration, known from simple models and from simulations: i) the backreaction of the accelerated particles is important, even in those cases in which a small fraction of the particles injected at the shock is actually accelerated to suprathermal energies. For $`\eta \stackrel{>}{}10^4`$ the shock becomes modified by the non-thermal pressure and both the velocity field of the shocked fluid and the spectrum of the accelerated particles are affected. The minimum value of $`\eta `$ for which the shock is modified in our model is slightly larger than that predicted in . ii) For linear shocks the maximum compression factor that can be achieved is $`4`$ (for $`\gamma _g=5/3`$) which corresponds to spectra $`p^4`$. This result is obtained for $`M_0\mathrm{}`$. When the backreaction becomes relevant, the structure of the shocked fluid is changed into a smooth decrease of the fluid speed from $`u_0`$ at infinity (effectively at some distance $`D_{max}`$) to $`u_1`$ at the position of the gas (ordinary) subshock. While at the subshock it is still true that the maximum compression factor is $`4`$, the overall compression factor ($`R_{tot}`$) between downstream and $`D_{max}`$ can be arbitrarily large. iii) The large total compression factors for strongly modified shocks result in a flattening of the spectra of accelerated particles at high energy. The slope there tends to $`3.5`$ (flatter than $`4`$). The slope of the spectrum at low energy is determined by the compression factor at the subshock, and is usually steeper than the slope predicted in linear theory. iv) In general, increasing the Mach number corresponds to an increasingly more modified shock, up to a critical Mach number $`M_{cr}`$. At $`M>M_{cr}`$, the total and subshock compression factors both converge to the linear value of $`4`$ and the shock behaves as an ordinary strong shock (it is no longer an efficient accelerator).

Although the calculations reported in the paper refer to the case of a plane (one dimensional) shock and the effects of waves (Alfvèn heating) have not been included, these changes are almost straightforward. In fact, for spherical shocks, the expression for the distribution function $`f_0(p)`$ remains unchanged \[eq. (6)\], although the definition of $`u_p`$ is formally different (it reflects the geometry). All the rest of the calculation remains unaffected and most of the results still hold. The case of spherical symmetry would be relevant for shocks related to supernova remnants, as discussed in .

The introduction of the heating due to damping of Alfvén waves implies a change in the relation between $`P_{g,p}`$ and $`P_{g,0}`$ but this change can be easily written in a way which is useful for our purposes. For instance, generalizing a discussion in , in the approximation of large Alfvén Mach number ($`M_A=u_0/v_A`$, where $`v_A`$ is the Alfvén speed, assumed constant in all the fluid) we can write

$$\frac{P_{g,p}}{P_{g,0}}\left(\frac{\rho _p}{\rho _0}\right)^{\gamma _g}\left\{1+(\gamma _g1)\frac{M_0^2}{M_A}\left[1\left(\frac{\rho _0}{\rho _p}\right)^{\gamma _g}\right]\right\}.$$
(19)
Clearly the effects of Alfvén heating can be neglected as long as $`M_0^2M_A`$. Introducing eq. (19) \[instead of eq. (12)\] in eq. (10), it is easy to derive an equation similar to eq. (18), that can be solved for $`U(p)`$. The spectrum of accelerated particles is then obtained using the same procedure illustrated above. The basic effect of the Alfvén heating is to reduce the total compression factors, and make the spectra at high energy slightly steeper than found in §4.

We conclude by stressing that the non-linear effects discussed here appear to be relevant even for a small injection efficiency, and their phenomenological consequences can be critically important. For instance, as an application of the simple model presented in , in Ref.  the spectra of secondary radio and gamma radiation in supernova remnants were calculated. The differences with respect to the results obtained by adopting the test particle approximation are impressive.

In this prospective, the semi-analytical model presented here, being of easy use and having an immediate physical interpretation, provides the suitable tool to estimate the phenomenological consequences of non-linearity in shock acceleration, in the cases where it is unpractical to have access to numerical simulations.

Acknowledgments This work was supported by the DOE and the NASA grant NAG 5-7092 at Fermilab.