text:
# Explosive Nucleosynthesis Associated with Formation of Jet-induced GRBs in Massive Stars

## 1 Introduction

A revolution in our understanding of gamma-ray bursts (hereafter GRBs) occurred in 1997, following accurate localization of GRBs by BeppoSAX Boella et al. (1997) and the discovery of the afterglow phenomenon. Some GRBs and afterglows are known to have extragalactic origins, which are nicely accounted for by fireball models M$`\stackrel{´}{\mathrm{e}}`$sz$`\stackrel{´}{\mathrm{a}}`$ros & Rees (1997), although the origin of the fireball is not yet clear.

However, there has been growing evidence linking GRBs to massive stars. The host galaxies of GRBs are star-forming galaxies and the position of GRBs appear to trace the blue light of young stars Tsvetkov et al. (2001); Bloom et al. (2002); Floc’h et al. (2003). Some of the host galaxies appear to be dusty with star-formation rates comparable to ultra-luminous infrared galaxies (Berger et al., 2001; Frail et al., 2002). On smaller scales, there is mounting evidence of association with regions of high ambient density (Galama & Wijers, 2001; Harrison et al., 2001) and the so-called dark GRBs arise in or behind regions of high extinction Piro et al. (2002).

Moreover, there has been tantalizing evidence for the existence of an underlying supernova (SN). The first association of a cosmologically distant GRB with the death of a massive star was found for GRB 980326, where a clear excess of emission was observed over and above the rapidly decaying afterglow component. This late-time bump was interpreted as a contribution of an underlying SN Bloom et al. (1999) since, unlike the afterglow, the bump was very red. GRB 970228, also with an intermediate-time bump and characteristic SN spectral rollover, is another good candidate (Reichart, 1999; Galama et al., 2000).

It was also reported that there seems to be a physical connection between GRB 980425 and SN 1998bw Galama et al. (1998). They discovered an optical transient within the error box of BeppoSAX Wide Field Camera for GRB 980425. It should be noted that this SN is categorized as a new type of SN (i.e. hypernova) with large kinetic energy ($`10^{52}`$ ergs), nickel mass ($`0.5M_{}`$), and luminosity (Iwamoto et al., 1998; Woosley et al., 1999). In the analysis of an underlying supernova mentioned above, the light curve of SN1998bw is usually used for the SN template.

Generally speaking, it is considered that it is too difficult to realize a GRB from a death of a massive star, because a fire-ball is required to be composed of enormous radiation energy ($`10^{51}`$ erg) with very small baryon number ($`10^6M_{}/m_p`$, where $`m_p`$ is the proton mass) M$`\stackrel{´}{\mathrm{e}}`$sz$`\stackrel{´}{\mathrm{a}}`$ros & Rees (1997). One of the most famous models to realize a GRB and large explosion energy from a death of a massive star is the collapsar model Woosley (1993); MacFadyen & Woosley (1999). Collapsar is defined as a massive star whose iron core has collapsed to a black hole that is continuing to accrete envelope matter at a very high rate. Woosley also pointed out that there will be two types for collapsars. One (type I collapsar) is that the central core immediately forms a black hole with an accretion disk. In the other case (type II collapsar), the central core forms a neutron star at first, but the neutron star collapses to a black hole with an accretion disk due to the continuous fall back. In both types, it is pointed out that an accretion disk should be formed around the central compact object and neutrinos should be emitted from the innermost region of the accretion disk through the electron (positron) capture and electron-positron pair annihilation Berezinskii & Prilutskii (1987); MacFadyen & Woosley (1999). Finally, a strong jet with large radiation energy and small baryon number, which is required to produce a GRB, is generated around the polar region due to the pair-annihilation of neutrinos that come from the accretion disk and/or MHD processes. In fact, it is reported that the narrow line of oxygen and broad lines of iron in SN 1998bw can be well reproduced when jet-like explosion is assumed Maeda et al. (2002).

One may have thought that the system of a GRB has been well known and there is no question about that. That is, a GRB is born from a death of a massive star which collapses and makes a black hole at the center forming a strong jet with a large photon to baryon ratio, large explosion energy of order $`10^{52}`$ ergs, and a great amount of $`{}_{}{}^{56}\mathrm{Ni}`$ ($`0.5M_{}`$) that makes the underlying SN so luminous that it is categorized as a hypernova. However, there are some important problems with GRBs.

First, it is reported that no coincident SN is detected at 99.7$`\%`$ confidence by Hubble Space Telescope in GRB 010921 to a limit 1.34 mag fainter than SN 1998bw. The host galaxy’s redshift is relatively low ($`z`$=0.451), which made intensive multi-color observations possible. It is also pointed out that the observed luminosity of SN 1998bw may be reproduced with a smaller amount of $`{}_{}{}^{56}\mathrm{Ni}`$ than that required in the spherically-symmetric models by invoking the angle-dependent luminosity expected from asymmetric explosion Wang et al. (2002). Thus it seems to be still uncertain whether a massive star, that makes a GRB, always produces much of $`{}_{}{}^{56}\mathrm{Ni}`$.

Second, there are some reports on the strong iron K$`\alpha `$ emission lines (Piro et al., 1998, 2000; Antonelli et al., 2000; Yoshida et al., 2001) in X-ray afterglows. Moreover, there is a report of emission lines of Mg, Si, S, Ar, and Ca with an outflow velocity of order 0.1$`c`$ in the X-ray afterglow of GRB 011211 Reeves et al. (2002), although it seems to be still controversial Rutledge & Sako (2002). At present, the origin of these lines is still open to arguments. We should investigate whether these features can be explained in the context of collapsar model.

Third, it is an open question whether $`{}_{}{}^{56}\mathrm{Ni}`$ is synthesized enough to explain the light curve of SN 1998bw in the collapsar model. As explained below, there is an inconsistent treatment between the collapsar model and the hypernova model. In a simulation of the type I collapsar, the mass accretion rate is estimated to be $`0.07M_{}`$ s-1 in the phase of jet formation MacFadyen & Woosley (1999), and the timescale of neutrino emission from the accretion disk is also estimated to be about 10 seconds since the temperature of the inner-most region of the accretion disk becomes lower than 1 MeV after that and electron-positron pair creation becomes impossible. As a result, the neutrino emission from the accretion disk ceases (MacFadyen & Woosley, 1999; Nagataki et al., 2003). In fact, the thermal energy is usually injected around the polar region at a rate $`\dot{E}10^{51}`$ erg s-1 for about 10 seconds in numerical simulations, which makes the total explosion energy $`10^{52}`$ ergs (MacFadyen & Woosley, 1999; Aloy et al., 2000). On the other hand, in a simulation of explosive nucleosynthesis in the context of a hypernova, total explosion energy is put at once below the mass cut that divides the ejecta from the collapsing core as an initial condition Maeda et al. (2002). In their simulations, considerable fraction ($`0.5`$) of the explosion energy is given as kinetic energy, as opposed to the simulation of collapsars. Also, the progenitor of a hypernova is assumed to be spherical. Thus it should be necessary whether these simulations are consistent with each other or not. That is, we should investigate whether a collapsar model can explain the phenomenon of SN 1998bw as in the hypernova model.

It should be noted here that the accretion disk around the central black hole may be the site where $`{}_{}{}^{56}\mathrm{Ni}`$ is synthesized (Woosley et al., 1999; MacFadyen & Woosley, 1999; Woosley, 2002; Pruet et al., 2003). Pruet et al. (2003) calculated electron fractions in the accretion disk and found the proper condition to produce $`{}_{}{}^{56}\mathrm{Ni}`$ in the disk. They also pointed out the possibility that a considerable amount of $`{}_{}{}^{56}\mathrm{Ni}`$ can be driven by the disk wind. Although their scenario is widely supported to be very promising, the quantative evaluation of the yields has yet to be done. In this study, we consider another possibility that $`{}_{}{}^{56}\mathrm{Ni}`$ is synthesized in the jet, as investigated in Maeda et al. (2002). Thus reader should be reminded that this is our assumption in this study. We also have to note that the opening angle of the jet launched by the neutrino heating has to be large in our scenario, as shown in the following sections, to produce a large enough amount of $`{}_{}{}^{56}\mathrm{Ni}`$ to explain the luminosity of the hypernova such as SN 1998bw. Such a wide jet will not make a GRB because a bulk Lorentz factor can not be so large as to realize a fireball. Thus, in our scenario, another narrow jet produced possibly by MHD effect ir as a small part of the wider jet will be needed to realize a GRB. In this study, however, we do not consider it explicitly, since, if any, it should have little influence on the product of explosive nucleosynthesis because the baryon mass in the narrow jet should quite small.

In this study, explosive nucleosynthesis is calculated in the context of type I collapsar model. At first, a progenitor collapses and forms a black hole at the center. Next, thermal energy is injected around the polar region to make a jet. Finally, explosive nucleosynthesis is calculated by the post-processing method. We investigate the influence of the progenitor structure and energy deposition rate on the results of explosive nucleosynthesis. We show the amount of $`{}_{}{}^{56}\mathrm{Ni}`$ is very sensitive to the energy deposition rate. A large amount of $`{}_{}{}^{56}\mathrm{Ni}`$ is produced when the energy deposition rate is high and vice versa. Thus we conclude that it is quite natural not to detect an underlying SN in some cases. In such cases, we consider that the energy deposition rate was relatively low and a small amount of $`{}_{}{}^{56}\mathrm{Ni}`$ was synthesized. We also give a discussion on the chemical composition of the ejecta of a collapsar and observed line features in some X-ray afterglows. We point out the possibility that the variation of chemical composition in the X-ray afterglow may be naturally explained by the variation of the energy deposition rate at the central engine. By this discussion, we can predict that there should be correlation between the line features in the X-ray afterglow and duration of the GRB. We also compare the calculated amount of $`{}_{}{}^{56}\mathrm{Ni}`$ with the observed value in SN 1998bw and investigate which model is best. We show that the type I collapsar model in which the energy deposition rate is relatively low ($`\dot{E}10^{51}`$ erg s-1) may have difficulty in reproducing the observed amount of $`{}_{}{}^{56}\mathrm{Ni}`$ of SN 1998bw and GRB 980425 and some improvements should be required to explain these observations in the context of the type I collapsar model. This means that the mechanism of the central engine of a hypernova and a GRB is constrained by the discussion of explosive nucleosynthesis.

In section 2, the method of calculation is explained in detail. Results are shown in section 3. Discussions are presented in section 4. Conclusion is given in section 5.

## 2 Method of Calculation

In this section, we present the method of the simulations in this study. It is true that there are some assumptions and simplifications in our simulations. In particular, the effects of gravitation and rotation are not included in this study. In order to realize the jet-induced explosion, we prepare an asymmetric progenitor as an initial condition following the radial infall with the infall timescale depending on the zenith angle. In this way, we can mimic an asymmetric progenitor models of MacFadyen and Woosley (1999). After that, we realize the jet-induced explosion by injecting the thermal energy around the polar region in the same way as MacFadyen and Woosley (1999) did. Of course, this needs improvement in the future. We explain the assumption and simplification in more detail in the following subsections.

### 2.1 Hydrodynamics

#### 2.1.1 The Scheme

We performed 2-dimensional hydrodynamic calculations using a spherical coordinate taking account of special relativity. The calculated region corresponds to a quarter of the meridian plane under the assumption of axial-symmetry and equatorial symmetry. The number of mesh points is 500 in the radial direction and 50 in the angular one for the calculation of explosive nucleosynthesis. The radial mesh size is set to be constant in space, 2$`\times 10^7`$ cm. With this resolution, the structure of the central accretion disk shown by MacFadyen and Woosley (1999) can be marginally resolved. We did not perform 3-dimensional hydrodynamic calculations simply for saving CPU time and memory. We employed an approximate Riemann solver, Maquina’s flux formula Donat & Marquina (1996); Donat et al. (1998), for special relativistic hydrodynamic equations. The code is of first order accuracy both spatially and temporally. In this study, we assume that the gas is ideal with the equation of state $`p=\rho ϵ(\gamma 1)`$, where $`p`$, $`\rho `$, $`ϵ`$, and $`\gamma `$ are pressure, rest mass density, specific internal energy and adiabatic exponent (=4/3 constant in this paper), respectively. Since the hydrodynamics code is Eulerian, we use the test particle method Nagataki et al. (1997) in order to obtain the informations on the time evolution of the physical quantities along the fluid motion, which are then used for the calculations of the explosive nucleosynthesis. Test particles are scattered in the progenitor and are set at rest at first. They move with the local fluid velocity at their own positions after the passage of the shock wave. The temperature and density that each test particle experiences at each time step are preserved. This is the test particle method we use. See Nagataki et al. (1997) and Nagataki et al. (1998) for more details.

Calculations of hydrodynamics and explosive nucleosynthesis are performed separately, since the entropy produced during the explosive nucleosynthesis is much smaller ($``$ a few$`\%`$) than that generated by the shock wave. In calculating the total yields of elements, we assume that each test particle has its own mass determined from their initial distribution so that their sum becomes the mass of the layers where these are scattered. It is also assumed that the nucleosynthesis occurs uniformly in each mass element. These assumptions will be justified since the movement of the test particles is not chaotic (i.e. the distribution of test particles at the final time still reflects the given initial condition) and the intervals of test particles are sufficiently narrow to give a smooth distribution of the chemical composition in the ejecta. The number of the test particles are 42000. The test particles are distributed uniformly within the radius of 7$`\times 10^9`$ cm.

#### 2.1.2 Initial and Boundary Conditions

As for the progenitor of SN 1998bw, it is thought to have had the mass $`40M_{}`$ in the main-sequence stage (Iwamoto et al., 1998; Woosley et al., 1999) and had $``$ 16$`M_{}`$ helium core Iwamoto et al. (1998). In this study, the presupernova model obtained from the evolution of 16 $`M_{}`$ helium core Nomoto & Hashimoto (1988) is used for the calculation of explosive nucleosynthesis.

As explained above, we make an asymmetric progenitor from the spherical presupernova model by assuming that matter free-falls radially during the times which depend on the zenith angle. The infall time is determined so that the mass of the central black hole becomes 3$`M_{}`$ at appropriate time and the asymmetric progenitor resembles to the result of MacFadyen and Woosley (1999) in which the effects of gravitation and rotation are taken into account (see their figures 7 and/or 9). The adopted function for the free-fall timescale is

$`t_{\mathrm{fall}}=c_1\mathrm{exp}((\theta /\sigma )^2)[\mathrm{sec}]`$ (1)
where $`\theta `$ is the zenith angle (in units of radian), and $`c_1`$ and $`\sigma `$ are set to be 1.0 and $`\pi /4.6`$, respectively. After the radial infall, we can determine the Schwarzschild radius and the baryon mass of the black hole. As mentioned above, the parameters in Eq. (1) are so shosen to give the baryon mass of the central black hole 3$`M_{}`$. As will be explained in section 4, the neutrino heating is effective as long as the mass of the central black hole is less than$`4M_{}`$. Hence, if the mass accretion rate is $`0.1M_{}`$ s-1 and the duration of the neutrino heating is $`10`$ s as inferred from the result of MacFadyen and Woosley (1999), the mass of the central black hole will be $`3M_{}`$ at the time of the jet launch. In fact, MacFadyen and Woosley (1999) launched the jet when the mass of the central black hole becomes 3.5$`M_{}`$ in their calculation, which is consistent with our treatment. Aloy et al. (2000) also launched the jet when the mass of the central black hole becomes 3.762$`M_{}`$ in their calculation. The density, on the other hand, at the innermost region around the polar axis is about $`10^6`$ g cm-3 and that around the equatorial plane is about $`10^8`$ g cm-3 in our model, which is similar to the result of MacFadyen and Woosley (1999). This indicates that our prescription Eq. (1) approximately reproduces the more realistic models of MacFadyen and Woosley (1999). It is noted that the accretion disk obtained numerically by MacFadyen and Woosley (1999) well matches the steady state disk (slim disk) model for 3$`M_{}`$ Schwarzschild black hole (kerr parameter $`a=0`$), viscosity parameter $`\alpha =0.1`$, and accretion rate of 0.1$`M_{}`$ s-1 Popham et al. (1999). Since we here mimiced the numerical result of MacFadyen and Woosley (1999), the situation we consider corresponds to the accretion disk around the black hole with the parameters mentioned above. The dependence of the explosive nucleosynthesis on different initial conditions will be reported in the forth-coming paper. For comparison, we also prepare a spherical model in which the free-fall time does not depend on the zenith angle. The mass of the central black hole is again set to be 3$`M_{}`$. The density contours of the progenitors in this study are shown in Figs. 1 and 2, with the former corresponding to the spherical model and the latter to the asymmetric model. In the asymmetric model, the upper hemisphere is shown.

It will become useful to evaluate the free-fall timescale in the central region of the progenitor and mass accretion rate assuming the free-fall. The free-fall timescale can be estimated as Woosley (1986)

$`\tau _{\mathrm{ff}}={\displaystyle \frac{1}{\sqrt{24\pi G\rho }}}{\displaystyle \frac{446}{\rho }}[\mathrm{sec}].`$ (2)
The mass accretion rate in the case of free-fall can be also estimated using the mass coordinate and $`\tau _{\mathrm{ff}}`$ as

$`\dot{M}(r)={\displaystyle \frac{M_r}{\tau _{\mathrm{ff}}}}[M_{}\mathrm{sec}^1],`$ (3)
where $`M_r`$ is the mass coordinate which means the enclosed mass within the radius $`r`$ from the center of the progenitor. In Fig. 3, we show the estimated free-fall timescale and mass accretion rate in the case of free-fall as a function of the mass coordinate. We can check from the figure that the mass of the central black hole becomes $`3M_{}`$ in 1 second in the case of free-fall. As briefly explained in section 1, the neutrino emission from the accretion disk around the black hole ceases when the temperature of the inner-most region of the accretion disk becomes lower than 1 MeV. This corresponds to the time when the mass of the black hole becomes $`4M_{}`$ (see section 4 and Fig. 12). Thus we can find from the figure that the energy deposition from the accretion disk ceases within 2 seconds in the case of free-fall. Further, we can estimate the thermal energy deposition rate due to the pair annihilation of (anti-)neutrinos as a function of the mass accretion rate, which we will discuss in section 4.

After making the asymmetric progenitor, the thermal energy is injected at the innermost grid around the polar axis at a rate $`\dot{E}=10^{51}`$ ergs s-1, which is a common treatment for the calculations of collapsars (MacFadyen & Woosley, 1999; Aloy et al., 2000). We do not consider explicitly how infalling matter are inverted. Instead, we assume that the momentum balance is achieved when the mass of the black hole becomes 3$`M_{}`$. After the launch of the jet, the thermal pressure dominates the ram pressure since the accretion energy has declined to about a few times 1050 erg s-1 at that time MacFadyen & Woosley (1999). For comparison, we also perform the computations in which the total explosion energy (= $`10^{52}`$ ergs) is put instantaneously at the innermost grid around the polar axis. The models investigated in this study are summarized in Table 1. In Models Sa/Aa, the spherical/asymmetric structure is adopted as a progenitor into which thermal energy is injected at a rate $`\dot{E}=10^{51}`$ ergs s-1. In Models Sb/Ab, the total explosion energy is deposited instantaneously to the spherical/asymmetric structure.

As for the boundary condition, the reflective boundary is adopted on the equatorial plane and symmetry axis for simplicity. It is true that some fraction of explosion energy and matter falls into the black hole. However, we do not observe these energy and matter and are not interested in their chemical composition, either. Hence we consider that the deposited energy in this study corresponds to the observed explosion energy, $`10^{52}`$ erg.

### 2.2 Nuclear Reaction Network

Since the chemical composition behind the shock wave is not in nuclear statistical equilibrium, the explosive nucleosynthesis has to be calculated using the time evolution of $`(\rho ,T)`$ and a nuclear reaction network. We obtain $`(\rho ,T)`$ comoving with the matter by means of the test particle method mentioned in subsection 2.1.1. The nuclear reaction network contains 250 species (see Table 2). We add some species around $`{}_{}{}^{44}\mathrm{Ti}`$ to Hashimoto’s network that contains 242 nuclei Hashimoto et al. (1989), although it turned out that the result was not changed essentially by the addition.

## 3 Results

We first show the results of hydrodynamical calculations. In Figs. 45, and 6, we draw the contours of temperature for the spherical explosion models at $`t=1.0`$ sec, 3.0 sec, and 5.0 sec, respectively. Here we assumed the axial symmetry and equatorial symmetry. The left panel corresponds to Model Sa and the right panel to Model Sb. It is clearly found that the velocity of the shock wave is faster in Model Sb than in Model Sa. This is because the energy density behind the shock is greater in Model Sb, which means that the pressure behind the shock is greater in Model Sb.

In Figs. 78, and 9, we show the contours of temperature for the asymmetric explosion models. In these figures, only the upper hemisphere is shown. The same tendency as in the spherical explosion model, that the velocity of the shock wave is faster in Model Ab than in Aa, is confirmed. It is also clearly shown that the velocity of the shock wave depends on the zenith angle and jetlike explosion is induced, although (as opposed to Maeda et al. 2002) only the thermal energy is injected. This reflects the initial density structure of the progenitor and is consistent with MacFadyen and Woosley (1999).

In Figs. 10 and 11, we show the positions at $`t=0`$ of the test particles with the mass fraction of $`{}_{}{}^{56}\mathrm{Ni}`$ greater than 0.3. It is clearly seen that there is almost no region where $`{}_{}{}^{56}\mathrm{Ni}`$ is synthesized in Models Sa and Aa, whereas $`{}_{}{}^{56}\mathrm{Ni}`$ is synthesized in large region for Models Sb and Ab. In Table 3, we summarize the abundance of heavy elements in the ejecta for each model, assuming that all unstable nuclei produced in the ejecta decay to the corresponding stable nuclei. The amount of $`{}_{}{}^{56}\mathrm{Ni}`$ is also shown in the last line. From the table, we can easily confirm that $`{}_{}{}^{56}\mathrm{Ni}`$, which decays to $`{}_{}{}^{56}\mathrm{Fe}`$, is hardly synthesized in Models Sa and Aa (Note that almost all of 56Fe is generated from the decay of $`{}_{}{}^{56}\mathrm{Ni}`$). It is emphasized that the amount of $`{}_{}{}^{56}\mathrm{Ni}`$ synthesized in Models Sa and Aa is less than the required amount to explain the luminosity of SN 1998bw ($`0.7M_{}`$; Iwamoto et al. 1998). On the other hand, the amount of $`{}_{}{}^{56}\mathrm{Ni}`$ synthesized in Models Sb and Ab is almost sufficient to explain the light curve of SN 1998bw. It is noted that the progenitor before the explosion is mainly composed of O, Ne, and Mg Nomoto & Hashimoto (1988), while the elements heavier than Si are synthesized through the explosive nucleosynthesis Hashimoto (1995). It is evident in Table 3 that explosive nucleosynthesis occurs and the initial chemical composition of the progenitor is changed to the heavier elements in Models Sb and Ab, while nearly no explosive nucleosynthesis occurs in Models Sa and Aa.

Here we consider why iron elements are hardly synthesized in Models Sa and Aa while much of them are synthesized in Models Sb and Ab. The most significant parameter in explosive nucleosynthesis in a massive star is the temperature. The criterion for the complete silicon burning is $`T_{\mathrm{max}}5\times 10^9`$ \[K\] Thielemann et al. (1996). It is well known that the matter behind the shock wave is radiation dominated and $`T_{\mathrm{max}}`$ can be well estimated by equating the supernova (hypernova) energy with the radiation energy inside the radius $`r`$ of the shock front

$`E_{\mathrm{HN}}=10^{52}\left({\displaystyle \frac{E_{\mathrm{HN}}}{10^{52}\mathrm{erg}}}\right)={\displaystyle \frac{11\pi ^3}{45}}{\displaystyle \frac{k^4}{\mathrm{}^3c^3}}r^3T_{\mathrm{max}}^4[\mathrm{erg}],`$ (4)
where $`E_{\mathrm{HN}}`$ is the total explosion energy of hypernova, $`k`$ is the Boltzmann constant, $`\mathrm{}`$ is the Planck constant divided by $`2\pi `$, $`c`$ is the velocity of light. Here a spherical explosion is assumed. This equation gives

$`r=5.7\times 10^8\left({\displaystyle \frac{5\times 10^9\mathrm{K}}{T_{\mathrm{max}}}}\right)^{4/3}\left({\displaystyle \frac{E_{\mathrm{HN}}}{10^{52}\mathrm{erg}}}\right)^{1/3}[\mathrm{cm}].`$ (5)
In the case of Model Sb, $`{}_{}{}^{56}\mathrm{Ni}`$ is synthesized within the edge for the complete silicon burning ($`5.7\times 10^8`$ cm). In the case of Model Sa, on the other hand, matter start to move outwards after the passage of the shock wave, and almost all of the matter move away ($`r6\times 10^8`$ cm) before the injection of all the thermal energy (= $`10^{52}`$ erg). This is the reason why almost no complete silicon burning occurs in Model Sa. The situation should be the same in the asymmetric explosion models. It is apparent from Table 3 that the chemical composition of the ejecta depends more sensitively on the timescale of thermal energy injection than on the asymmetry of the initial density structure. It should be emphasized that the discussion presented here can also be adopted to the pulsar-powered models for GRBs Woosley (2002), in which, as Woosley (2002) pointed out, very little nickel will be produced unless a pulsar deposit at least $`10^{51}`$ erg in a few tenths of a second. In the next section, we discuss the implication of these results.

## 4 Discussion

In this study, we assume that $`{}_{}{}^{56}\mathrm{Ni}`$ is synthesized in the jet as done in Maeda et al. (2002). In this section, we give further discussion under this assumption.

First, as stated in section 1, it is reported that no coincident SN is detected in GRB 010921 to a limit of 1.34 mag fainter than SN 1998bw at 99.7$`\%`$ confidence is detected by Hubble Space Telescope. If we simply interpret this result as a lack of $`{}_{}{}^{56}\mathrm{Ni}`$ in GRB 010921, the amount of $`{}_{}{}^{56}\mathrm{Ni}`$ synthesized in GRB 010921 is less than 0.7$`M_{}`$/(2.5121.34)$``$0.2$`M_{}`$, which can be naturally explained by Models Sa and Aa. We want to emphasize that the amount of $`{}_{}{}^{56}\mathrm{Ni}`$ in the ejecta can be small even if the total explosion energy is much larger than that of the normal collapse-driven supernova, as long as the energy deposition rate is relatively low. We insist that observations like GRB 010921 can be naturally explained in the context of the type I collapsar model.

Second, we discuss the features of line emission in some X-ray afterglows. As stated in section 1, strong iron K$`\alpha `$ emission lines are found in X-ray afterglows (Piro et al., 1998, 2000; Antonelli et al., 2000; Yoshida et al., 2001). There is also a report of emission lines of Mg, Si, S, Ar, and Ca with an outflow velocity of order 0.1$`c`$ in the X-ray afterglow of GRB 011211 Reeves et al. (2002). According to Reeves et al. (2002), the relative abundance of the light metals (Mg, Si, S, Ar, Ca) was found to be $`9`$ times the solar value, whereas the relative abundance of iron was found to be $`1.4`$ times the solar value. Although the mechanism of such emission lines is still controversial (Paczyński, 1998; Lazzati et al., 1999; Rees & M$`\stackrel{´}{\mathrm{e}}`$sz$`\stackrel{´}{\mathrm{a}}`$ros, ; Vietri et al., 2001; Kotake & Nagataki, 2001; Kosenko et al., 2002), we want to point out that the chemical composition in the X-ray afterglow may be naturally explained if the effects of the energy deposition rate are taken into account, as shown in this study. Here the origin of the emission lines is assumed to be the ejecta from the underlying SN (HN). One may think that the relative abundance of iron becomes small when the effect of fallback is taken into account. However, if the jetlike explosion occurs to generate a GRB, iron elements are synthesized around the jet axis as shown in Fig. 11. In this case, it will be hard to consider that the relative abundance of iron becomes small due to the effects of the fallback. According to Reeves et al. (2002), the velocity of the outflow which emits the lines of heavy elements is of order 0.1$`c`$, which is about ten times faster than the velocity of the ejecta in the normal collapse-driven supernova (Nagataki et al., 1998; Nagataki, 2000). The required kinetic energy is easily estimated as

$`E_{\mathrm{kin}}=9\times 10^{51}\left({\displaystyle \frac{M}{1M_{}}}\right)\left({\displaystyle \frac{v}{0.1c}}\right)^2\mathrm{erg},`$ (6)
which means that the underlying SN is a hypernova. We emphasize again that the abundance of the elements with intermediate mass number such as Si and S becomes greater than that of iron even if a hypernova explosion is adopted as long as the timescale of the energy deposition is long. If this discussion is true, there should be correlation between the chemical composition of the ejecta and the timescale of the energy deposition. This should mean that there should be correlation between the line features in the X-ray afterglow and the duration of the GRB. If the emission lines are mainly composed of light elements such as Si and S, the duration of the GRB should be relatively long and vice versa. In fact, it is reported that the duration of GRB 011211 is 270 seconds, making it the longest burst ever observed by Beppo-SAX Reeves et al. (2002). Although the interpretation of the line feature of X-ray afterglow of GRB 011211 seems to be still controversial Rutledge & Sako (2002), we hope the increase of observations in near future.

Third, we consider that the observations of SN 1998bw can be well explained by Models Sb and/or Ab, but not by Models Sa and Aa. This means that the energy deposition timescale should be relatively short ($`10`$ sec) in order to explain the luminosity of SN 1998bw by the type I collapsar model. One may think that the amounts of $`{}_{}{}^{56}\mathrm{Ni}`$ produced in Models Sb and Ab (0.14$`M_{}`$ and 0.16$`M_{}`$, respectively) are still less than the observed value ($`0.7M_{}`$; Iwamoto et al. 1998). However, the total explosion energy is assumed to be $`3\times 10^{52}`$ ergs in Iwamoto et al. (1998) while it is decreased to $`10^{52}`$ ergs in this study. Thus, we consider that the observed amount can be explained with no problem when the total explosion energy is set to be $`3\times 10^{52}`$ ergs in Models Sb and/or Ab. If we believe the conclusion derived by Maeda et al. (2002) that the line features in SN 1998bw can be well reproduced by the jet-induced explosion model, Model Ab seems to be the best one. However, we argue below that it seems to be difficult even for the models Sb and Ab to generate a GRB. At present, we are led to the conclusion that there seems to be no model to generate a GRB associated with a hypernova synthesizing much iron.

Now we discuss whether the results obtained in this study can constrain the model of the central engine of a hypernova accompanied by a GRB. According to the scenario of the type I collapsar model, the explosion energy larger than that of the normal collapse-driven supernova may be given by the pair annihilation of neutrinos from the accretion disk around the central black hole. It takes about several seconds from the beginning of the core collapse to form an accretion disk. In the meantime, the matter around the polar region falls directly into the central black hole. This is because the matter around the polar region possess smaller angular momentum. As a result, the polar region may become an adequate environment for generating a fireball. The timescale to form a fireball is determined by the timescale of neutrino emission from the accretion disk, which is about 10 seconds MacFadyen & Woosley (1999). Although we pointed out in this study that the type I collapsar model can naturally explain the observations of GRB 010921 and GRB 011211 in which little amount of iron is required, this scenario may face difficulty in producing a large amount of $`{}_{}{}^{56}\mathrm{Ni}`$ required in SN 1998bw, since the energy deposition rate is too low. If the energy deposition timescale is forced to be shorter ($`10`$ sec), it seems to become doubtful whether a GRB can be really generated, as explained below. To make the energy deposition rate greater means the larger the mass accretion rate. This corresponds to the case in which the total angular momentum in the central region of the progenitor is small. Then, it is doubted that whether an asymmetric explosion occurs and an adequate environment for generating a fireball is achieved around the polar region. That is, we suspect that the density around the polar region can not become low enough to realize a fireball when the mass accretion around the equatorial plane is large. This is because there is little time for the mass around the polar region to fall into the black hole before the explosion due to the pair annihilation of neutrinos occurs.

Let us investigate the above discussion quantitatively. Popham et al. (1999) proposed an accretion-disk model including the effect of general relativity as a collapsar model Popham et al. (1999). This model and the numerically computed collapsar model MacFadyen & Woosley (1999) are reproduced well by the analytical model of the accretion disk by Nagataki et al. (2002). The density, the temperature and the disk thickness are fitted as follows (see Nagataki et al. 2002 for details).

$`\rho \left[\mathrm{g}\mathrm{cm}^3\right]=8.23\times 10^8\left({\displaystyle \frac{M}{3M_{}}}\right)^{1.7}\left({\displaystyle \frac{\dot{M}}{0.1M_{}\mathrm{s}^1}}\right)^{1.03}{\displaystyle \frac{1}{\left({\displaystyle \frac{r}{r_s}}\right)^{1.07}\left(1+\left({\displaystyle \frac{r}{r_s}}\right)\right)^{0.76}}},`$ (7)
$`T\left[\mathrm{MeV}\right]=2.3\times \left({\displaystyle \frac{M}{3M_{}}}\right)^{0.2}\left({\displaystyle \frac{\dot{M}}{0.1M_{}\mathrm{s}^1}}\right)^{0.108}{\displaystyle \frac{1}{\left({\displaystyle \frac{r}{r_s}}\right)^{0.425}\left(1+\left({\displaystyle \frac{r}{r_s}}\right)\right)^{0.21}}},`$ (8)
$`H\left[\mathrm{cm}\right]=5.8\times 10^6\left({\displaystyle \frac{M}{3M_{}}}\right)^{0.9}\left({\displaystyle \frac{\dot{M}}{0.1M_{}\mathrm{s}^1}}\right)^{0.0183}{\displaystyle \frac{1}{\left({\displaystyle \frac{r}{r_s}}\right)^{1.66}\left(1+\left({\displaystyle \frac{r}{r_s}}\right)\right)^{0.3867}}},`$ (9)
where $`M`$ is the mass of the central black hole, $`r`$ is the radial coordinate, and $`r_s`$ (= 107cm) is the core radius, respectively. Note that the Schwarzschild radius is $`8.862\mathrm{km}(M/3M_{})`$. We see from Eq. (8) that as the mass of the central black hole grows, the temperature becomes lower, and the flux of neutrinos decreases. In Fig. 12, we show the duration of neutrino emission. This (almost) corresponds to the timescale in which the temperature at the innermost region of the accretion disk becomes lower than $``$ 1 MeV. This will approximately reflect the timescale of GRB. In this figure the initial mass of the central black hole was set to be 3$`M_{}`$. The time when the temperature of the inner-most region of the accretion disk becomes lower than 1 MeV almost corresponds to the time when the mass of the black hole becomes $`4M_{}`$. From this figure, we find that the mass accretion rate has to be larger than several $`M_{}`$ s-1 so that the energy deposition timescale is shorter than 1 second. This mass accretion timescale is shorter than that presented in MacFadyen and Woosley (1999). On the other hand, if the mass accretion is so rapid, there will be no time to make an adequate environment to produce a fireball around the polar region as mentioned above. Of course, we can not conclude that the type I collapsar model is ruled out by the discussion of the explosive nucleosynthesis as a model of hypernova associated with GRB. We guess that there will be some ways to solve the difficulty presented in this study. However, we believe that it is very important to point out that the central engine of a GRB accompanied by SN (HN) explosion can be constrained by the discussion of the explosive nucleosynthesis. At present, we have to say that the Type I collapsar model is required to be improved so as to explain the observed luminosity in SN 1998bw. It is noted that MacFadyen and Woosley (1999) pointed out the possibility that a substantial amount of $`{}_{}{}^{56}\mathrm{Ni}`$ is produced in the accretion disk and a part of it is conveyed outwards by the viscosity-driven wind (see their figure 16; see also Pruet et al. (2003)). Their scenario will be important and it should be investigated whether the behavior of the light curve and the line features of iron and oxygen can be explained in such a framework as well as the jetlike explosion model (Maeda et al. 2002). In the study of Maeda et al. (2002), the narrow line of oxygen and broad lines of iron can be explained when a jet-induced explosion is assumed.

We consider that the variation of the energy deposition rate is naturally caused by the variation of the progenitor’s initial angular momentum. Thus it may be still too early to consider that a GRB is always accompanied by a hypernova that generates a large amount of iron elements. This is what we want to stress in this paper and we expect the increase of observations in near future.

As mentioned in section 1, the accretion disk around the central black hole may be the site of $`{}_{}{}^{56}\mathrm{Ni}`$ synthesis Pruet et al. (2003). This scenario is widly supported to be most promising, although the quantitative evaluation of the yield is remaining to be done. If true, there may not be any correlation between the duration of GRB and the amount of synthesized $`{}_{}{}^{56}\mathrm{Ni}`$ (the collapsar without the disk wind may make too little $`{}_{}{}^{56}\mathrm{Ni}`$ in any case). Thus it might be possible to distinguish two scenarios by future observations.

In this study, we used the nuclear reaction network which contains 250 species (see table 2). Thus we could not discuss the nucleosynthesis of heavier elements such as in the r-process nucleosynthesis. We consider that the polar region within the jet may be an adequate site for r-process nucleosynthesis (MacFadyen & Woosley, 1999; Nagataki, 2000, 2001). We are planning in near future to perform the calculation of r-process nucleosynthesis in the context of jetlike explosion presented in this study.

In this study, we modeled the jet launched by the neutrino heating and a large opening angle of the jet was necessary for the jet to produce a large enough amount of $`{}_{}{}^{56}\mathrm{Ni}`$ to explain the luminosity of the hypernova such as SN 1998bw. As stated in section 1, such a wide jet will not make a GRB because a bulk Lorentz factor can not be sufficiently large to realize a fireball. Therefore, to realize a GRB in our scenario, another narrow jet might be needed supposedly launched by MHD effect. If the jet is launched by the MHD effect (Blandford & Znajek, 1977), a considerable fraction of the explosion energy should be in the form of kinetic energy and/or magnetic field energy Zhang et al. (2003). As a result the temperature should be lower compared with the neutrino heating model until the reverse shock converts the kinetic energy into the thermal energy Zhang et al. (2003). Thus the explosive nucleosynthesis will be unlikely to occur in this scenario unless the reverse shock propagates inwards immediately and makes the temperature high enough to cause the explosive nucleosynthesis.

We made in this study an asymmetric initial progenitor model by a simple prescriptism. This treatment should be, of course, improved. We are planning to perform numerical calculations of core collapse in a massive star that leads to a black hole formation including the effects of gravity, rotation and neutrino transfer (Yamada & Sato, 1994; Yamada, 1997; Kotake et al., 2003). We will then investigate in detail the dependence of explosive nucleosynthesis on the density structure at the time of the jet launch as a function of initial angular momentum in a progenitor based on their results in the near future.

## 5 Conclusion

We have performed numerical calculations of explosive nucleosynthesis in the context of the collapsar model using 2D relativistic hydrodynamical simulations. We have investigated the influence of the structure of the progenitor and energy deposition rate on the explosive nucleosynthesis, assuming that $`{}_{}{}^{56}\mathrm{Ni}`$ is mainly synthesized in the jet. We have shown the amount of $`{}_{}{}^{56}\mathrm{Ni}`$ is very sensitive to the energy deposition rate rather than the structure of the progenitor. We conclude that it is quite natural to detect no underlying SN in some cases, such as GRB 010921. We have pointed out that the relative abundance of the elements with intermediate mass number such as Si and S in the X-ray afterglow of GRB 011211 is naturally explained if the energy deposition rate at the central engine is relatively small. We also predicted that there should be correlation between the line features in the X-ray afterglow and the duration of the GRB. Interestingly, the duration of GRB 011211 is 270 seconds, the longest burst ever observed by Beppo-SAX although it suffers from the effect of red-shift ($`z_{\mathrm{host}}=2.14`$). Our results also imply that the type I collapsar model in which the energy deposition rate is relatively low ($`\dot{E}10^{51}`$ erg s-1) as shown in MacFadyen and Woosley (1999), may have difficulty in reproducing the observed amount of $`{}_{}{}^{56}\mathrm{Ni}`$ in SN 1998bw. This means that the mechanism of the central engine of hypernova and GRB is constrained by the discussion of explosive nucleosynthesis. We hope that the increase of observations will lead to the improvement of theories on GRB in near future.

The authors thank to Dr. Blinnikov for useful comments. S.N. thanks the Yukawa Institute for Theoretical Physics at Kyoto University, where this work was initiated during the YITP-W-01-06 on ’GRB2001’ and completed during the YITP-W-99-99 on ’Blackholes, Gravitational Lens, and Gamma-Ray Bursts’. Hydrodynamic calculations were carried out on NEC SX-5 system at Cybermedia Center, Osaka University. This research has been supported in part by the Grant-in-Aid by the Ministry of Education, Culture, Sports, Science and Technology of Japan (No.14079202, No.S14102004, No.S14740166).


format: text

text:
# Possible Fulde-Ferrell-Larkin-Ovchinnikov State in CeCoIn5

## Abstract

Angular dependent magnetothermal conductivity experiments on CeCoIn5 indicate that this compound is a d$`_{x^2y^2}`$-wave superconductor. In this study, the low-temperature behavior of the upper critical field is measured in a single crystal of CeCoIn5 along the directions $`\stackrel{}{H}\stackrel{}{a}`$ and $`\stackrel{}{H}\stackrel{}{c}`$. The data is compared with model calculations of the upper critical field in a d$`_{x^2y^2}`$-wave superconductor. It is found that the observed $`H_{c2}(T)`$ along $`\stackrel{}{H}\stackrel{}{a}`$ is consistent with a Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state at low temperatures, $`T<`$ 0.7K, whereas for $`\stackrel{}{H}\stackrel{}{c}`$ the FFLO state appears to be absent in CeCoIn5. Furthermore, it is predicted that the quasiparticle density of states in the FFLO state exhibits a complex peak structure which should be observable by scanning tunneling microscopy.

I. Introduction

Recent measurements on CeCoIn5 have led to a renewed discussion of a possible high-field Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state in unconventional superconductors. In this state, the coupling of the magnetic field to the quasiparticle spins dominates over the orbital coupling, leading to pairing between exchange-split Fermi surfaces, and hence to a spatially non-uniform superconducting order parameter. For conventional superconductors, its realization appears to be practically impossible because of two reasons. First, the sample quality has to be in the superclean limit, i.e. the quasiparticle mean path needs to be much larger than the coherence length. Second, the Ginzburg-Landau parameter $`\kappa `$, measuring the ratio of the magnetic penetration depth vs. the superconducting coherence length, should be very large, i.e. $`\kappa 10`$.

The recent synthesis of quasi-two-dimensional nodal superconductors, such as the high-$`T_c`$ cuprates, the $`\kappa `$-(ET)2 salts, and CeCoIn5 has changed this situation dramatically. It appears that the above two conditions can be met in high-quality single crystal samples of these compounds.  These systems are quasi-two-dimensional, leading to a large Ginzburg-Landau parameter in a planar magnetic field. Furthermore, unlike in the conventional s-wave superconductors, the stability region of the FFLO state is much more extended in d$`_{x^2y^2}`$-wave superconductors compared to conventional ones.

Indications for possible FFLO states in organic superconductors were already reported in $`\lambda `$-(BEDTS)2GaCl4 and $`\kappa `$-(BEDT-TFF)2Cu(NCS)2. In the first compound, a kink in the thermal conductivity points to a transition from a FFLO state to a vortex lattice. In the second material a similar feature in the magnetization was identified. Moreover, recent evidence for d$`_{x^2y^2}`$-wave order parameter symmetry was found in $`\kappa `$-(BEDT-TFF)2Cu(NCS)2 by angular dependent magnetothermal conductivity measurements in a rotating magnetic field within the conducting crystal plane. Moreover, it was observed that the upper critical field $`H_{c2}(T)`$ in the FFLO regimes decreases quasi-linearly with temperature, in contrast to the rather weak temperature dependence of $`H_{c2}`$ in the absence of a FFLO state as $`T0`$.

More recently, a new heavy fermion compound, CeCoIn5, was discovered. This material superconducts below a critical temperature $`T_c`$ = 2.3K,  and it has a layered structure similar to the high-$`T_c`$ cuprates. Angular dependent magnetothermal conductivity experiments indicate d$`_{x^2y^2}`$-wave superconductivity in this material. Furthermore, the temperature dependence of the upper critical field $`H_{c2}(T)`$ for both $`\stackrel{}{H}\stackrel{}{a}`$ and $`\stackrel{}{H}\stackrel{}{c}`$ was measured in single crystals of CeCoIn5. It was observed that $`H_{c2}(T)`$ for $`\stackrel{}{H}\stackrel{}{a}`$ exhibits a quasi-linear temperature dependence in the proposed FFLO regime, T $`<`$ 0.7K.

In this paper, the upper critical field $`H_{c2}(T)`$ in a single crystal sample of CeCoIn5 is determined from thermal expansion and magnetorestriction measurements with field orientations $`\stackrel{}{H}\stackrel{}{a}`$ and $`\stackrel{}{H}\stackrel{}{c}`$. A d$`_{x^2y^2}`$-wave model calculation is used to explain the temperature dependence of $`H_{c2}(T)`$ along both directions, considering the orbital effect and the Pauli term. In particular, the possibility of a FFLO state is addressed. For $`\stackrel{}{H}\stackrel{}{c}`$ it is found that the temperature dependence of $`H_{c2}(T)`$ can be fitted consistently to the experimental data without invoking a FFLO state. On the other hand, for $`\stackrel{}{H}\stackrel{}{a}`$ we observe that the inclusion of a $`\stackrel{}{v}\stackrel{}{q}`$ term arising from a FFLO state is crucial for a consistent description of the observed $`H_{c2}(T)`$ below T = 0.7K. In particular, the quasi-linear T-dependence of $`H_{c2}(T)`$ at low temperatures can be understood within this framework, which provides compelling evidence for a FFLO state in CeCoIn5. In order to help to further test and scrutinize this model, we also determine the corresponding quasiparticle density of states which should be accessible to scanning tunneling microscopy experiments.

II. $`H_{c2}(T)`$ for $`\stackrel{}{H}\stackrel{}{c}`$

The temperature dependence of the upper critical field along the crystal c-direction can be obtained from two coupled integral equations

$`\mathrm{ln}t`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{du}{\mathrm{sinh}u}}\left[1\mathrm{exp}(\rho u^2)\mathrm{cos}(hu)(1+2\rho ^2u^4C)\right],`$ (1)
$`C\mathrm{ln}t`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{du}{\mathrm{sinh}u}}\left\{C\mathrm{exp}(\rho u^2)\mathrm{cos}(hu)\left[{\displaystyle \frac{\rho u^2}{12}}+C\left(18\rho u^2+12\rho ^2u^4{\displaystyle \frac{\rho ^3u^6}{3}}\right)\right]\right\},`$ (2)
where $`tT/T_c`$, $`\rho (v^2eH)/(8\pi ^2T^2)`$, and $`h(g\mu _BH)/(2\pi T)`$. Here it is assumed that the wave function of the Abrikosov state for a d$`_{x^2y^2}`$-wave superconductor can be written as

$`|\mathrm{\Psi }=\left(1+C\left(a^{}\right)^4\right)|0,`$ (3)
where the “vacuum” $`|0`$ is the Abrikosov state of an s-wave superconductor, and $`a^{}`$ is the raising operator of the Landau level. In other words, $`|0`$ is a combination of the N = 0 Landau states. For d$`_{x^2y^2}`$-wave superconductors an admixture of higher Landau states that are allowed by symmetry needs to be included in order to account for structural changes in the vortex lattice.

In Fig. 1(a), experimental data of $`H_{c2}(T)`$ along the c-axis in CeCoIn5 is compared with the numerical solution of Eqs. 1 and 2, describing a d$`_{x^2y^2}`$-wave vortex lattice. The temperature dependence of the upper critical field $`H_{c2}(T)`$ is extracted from low-temperature thermal expansion, $`\mathrm{\Delta }l(T,B=\text{const})`$, and magnetostriction $`\mathrm{\Delta }l(B,T=\text{const})`$ measurements utilizing a high-resolution capacitive dilatometer at temperatures down to 15 mK and in magnetic fields up to 18 T.  For temperatures above $`T_00.7\text{K}`$, the superconducting-to-normal phase transition is of second order and the $`H_{c2}(T)`$ is determined from the midpoint of idealized jumps in $`\mathrm{\Delta }l/T`$ and $`\mathrm{\Delta }l/B`$. Below $`T_0`$ sharp jumps in $`\mathrm{\Delta }l`$ are observed upon crossing $`H_{c2}(T)`$, indicative of a first order phase transition. The best fit with a numerical solution of Eqs. 1 and 2 is obtained with $`v=3.2738\times 10^8`$ cm/sec, and $`g`$=1.5. For comparison, solutions with the same Fermi velocity $`v`$, but $`g`$= 0 and 2 are also shown. The fit to the experiment appears to be very good in the low-temperature regime $`T<`$0.7K. Therefore, there are presently no obvious indications for a FFLO state in this material by measurements of $`H_{c2}(T)`$ along the crystal c-direction.

In Fig. 1(b) we show the numerical solution for the admixture parameter $`C`$. Interestingly, for g$`1.2`$, $`C`$ changes its sign as the temperature is decreased. Consequently, for g$`=1.5`$ one finds that the conventional hexagonal vortex lattice, which is stable at high temperatures, may change into a square vortex for $`T/T_c<0.3`$. This transition should be observable by small angle neutron scattering (SANS) with $`\stackrel{}{H}\stackrel{}{c}`$. A similar instability to a square vortex lattice was previously predicted for the high-$`T_c`$ cuprates in the vicinity of $`H_{c2}`$. Both, SANS and STM on the vortex state of YBCO single crystals indeed indicate significant deviations from a hexagonal towards a square vortex lattice in this compound.  Furthermore, the above theory was recently extended to lower magnetic fields. The predicted square vortex lattice was observed by SANS in a single crystal of LSCO at $`H`$=2 Tesla.

III. $`H_{c2}(T)`$ for $`\stackrel{}{H}\stackrel{}{a}`$

In order to match the experimental data for $`H_{c2}(T)`$ along the crystal a-axis, we explore the effect of a $`\stackrel{}{v}\stackrel{}{q}`$ term arising from the formation of a FFLO state.  This leads to a new set of coupled integral equations,

$`\mathrm{ln}t`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{du}{\mathrm{sinh}u}}\left[1\mathrm{exp}(\rho u^2|s|^2)\mathrm{cos}\left[h(1p\mathrm{cos}\varphi )u\right]\left(1+\mathrm{cos}(4\varphi )\right)\left(12\rho u^2s^2C\right)\right]`$ (4)
$`C\mathrm{ln}t`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{du}{\mathrm{sinh}u}}\left[C\mathrm{exp}(\rho u^2|s|^2)\mathrm{cos}\left[h(1p\mathrm{cos}\varphi )u\right]\left(1+\mathrm{cos}(4\varphi )\right)\left(\rho u^2s^2+C\left(14\rho u^2|s|^2+2\rho ^2u^4|s|^4\right)\right)\right],`$ (5)
where $`s\mathrm{sin}\chi +i\mathrm{sin}\varphi `$, $`\rho (vv_ceH)/(8\pi ^2T^2)`$, $`p\mathrm{cos}\varphi (\stackrel{}{v}\stackrel{}{q})/(2h)`$, $`\chi ck_z`$, and $`\mathrm{}`$ is the angular average over $`\varphi `$ and $`\chi `$. Here $`\sqrt{vv_c}=1.63\times 10^8`$cm/s is used, and following Gruenberg and Gunther, we chose $`\stackrel{}{q}\stackrel{}{H}\stackrel{}{a}`$. In this configuration, the vortex state is represented by

$`|\mathrm{\Psi }=\left(1+C\left(a^{}\right)^2\right)|0,`$ (6)
where the vacuum $`|0`$ is again the Abrikosov state of a simple s-wave superconductor, but mixing now occurs with the N=2 Landau level.

Furthermore, let us note that for the purpose of In Fig. 2(a) $`H_{c2}(T)`$ is shown for g= 0.64 and 2 along with the measurements of the upper critical field along the a-axis of CeCoIn5. We find that without the FFLO state ($`p=0`$) one obtains a good fit to the experiment down to $`T`$ = 0.7K with g=0.64. However, for $`T<`$ 0.7K the measured upper critical field is approximately linear in temperature. This feature can be reproduced by including a $`\stackrel{}{v}\stackrel{}{q}`$ term due to the FFLO state with $`p=0.9`$. For comparison, we also show results for g = 2 which yield a zero-temperature critical field that is less than half the value detected in the experiment. In Fig. 2(b) $`p(t)`$ is shown for g = 0.64 and g = 2. From this plot it is clear that the FFLO region expands as g is increased. In Fig. 2(c) the coefficient $`C(t)`$ is shown. Here we observe that $`C`$ exhibits a significant temperature dependence for g = 2, whereas for g = 0.64 the admixture is almost negligible.

Furthermore, let us note that for the purpose of the present discussion it was assumed that the transition at $`H=H_{c2}`$ is of second order. However, a number of experiments on CeCoIn5 indicate a possible first order transition, and onset of magnetic order at $`T`$ 0.8K.  At the moment the nature of this magnetic order is unknown. In case it is a spin density wave, the condensation energy is expected to be relatively small, and consequently its effect on $`H_{c2}(T)`$ should be small as well.

IV. Quasiparticle Density of States

Let us conclude this discussion of a possible FFLO state in CeCoIn5 by calculating the shape of the associated quasiparticle density of states. In the vicinity of $`H=H_{c2}`$ and for $`\stackrel{}{H}\stackrel{}{a}`$ this observable is well approximated by

$`{\displaystyle \frac{N(E)}{N_0}}1={\displaystyle \frac{\mathrm{\Delta }^2}{4\sqrt{\pi }}}{\displaystyle \underset{\pm }{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑u{\displaystyle \frac{\mathrm{exp}(u^2)\mathrm{cos}^2(2\varphi )}{\left[E\pm \stackrel{~}{H}(1p\mathrm{cos}\varphi )ϵ|s|u\right]^2}},`$ (7)
where $`\stackrel{~}{H}(\mu _BgH)/2`$, $`|s|=\sqrt{\mathrm{sin}^2(\varphi )+\mathrm{sin}^2(\chi )}`$, and $`ϵ\sqrt{vv_ceH}`$. Again, a finite $`p`$ indicates the presence of a FFLO state. This quasiparticle density of states is plotted in Fig. 3. For the parameters, we have chosen $`p=0.9`$ and $`ϵ/h=0.2`$, appropriate for CeCoIn5 in the temperature regime $`T0.1K`$. In the absence of the FFLO state (Fig. 3(a)), there are two sharp resonances close to $`E=\pm H`$. In the presence of the FFLO state, more structure appears in the spectral response, as shown in Fig. 3 (b), with resonances at $`E=\pm H(1\pm p)`$. Therefore, precision measurements of the quasiparticle density of states can provide a clear signal for the presence of FFLO states and the symmetry of the underlying superconducting order parameter. It should thus be of interest to conduct a scanning tunneling microscope study of the quasiparticle density of states in CeCoIn5 at $`T<`$ 0.7K in order to further scrutinize the proposed FFLO state.

V. Conclusions

In summary, the model calculation in this study incoporates consistently (i) the d$`_{x^2y^2}`$-wave symmetry of the superconducting order parameter, (ii) the orbital effect, and (iii) a $`\stackrel{}{v}\stackrel{}{q}`$-term due to the formation of a FFLO state. The model appears to describe well the observed temperature dependence of the upper critical field in CeCoIn5. Furthermore, it indicates a significant renormalization of the g-factor in this compound, as well as the presence of a FFLO state at low temperatures if the applied field has an in-plane component. In this phase, the quasiparticle density of state is predicted to have a more complex structure. In order to further scrutinize the proposed model, it will be interesting to determine further relevant properties, such as the specific heat and the thermal conductivity.

VI. Acknowledgements

We thank J.L. Sarrao for providing a high-quality CeCoIn5 single crystal. Furthermore we are grateful to Y. Matsuda, A. Ardavan, T. Ishiguro, R. Movshovich, M.-S. Nam, T. Roscilde, and J. Sarrao for useful discussions on the FFLO state. S.H. acknowledges financial support by the Petroleum Research Foundation and the National Science Foundation, DMR-0089882.


text:
# Stochastic ‘fuzzy confinement’ of intrinsic localised modes

## Abstract

The long time diffusive behaviour of intrinsic localised modes (discrete breathers) in the discrete damped-driven sine-Gordon chain under Gaussian white noise (to simulate temperature) is studied. We present a theoretical model for an approximate description of the diffusion, derive an expression for the diffusion constant and compare with results from simulations. It turns out that an increase of the temperature inhibits the diffusive motion in such a way that the breather, propagating with a well-defined velocity in the noise-free case, is, with increasing probability, confined to a finite region of the system. As all physical processes in the real world occur at temperatures $`T>0`$, these results also have a bearing on the experimental detection of mobile breathers, e.g. in parallel arrays of Josephson junctions.

Intrinsic localised modes (ILMs) are localised excitations that can exist even in perfectly homogeneous discrete nonlinear systems due to an interplay of discreteness and nonlinearity. These modes, also referred to as discrete breathers (DBs), have received much attention over the past 10-15 years. At least in part this is due to the fact that little is needed for their existence, as shown by existence proofs  for this type of localised excitation. Therefore DBs are rather generic and can occur for example in magnetic lattices , nonlinear optical waveguide systems , Josephson junction arrays  or nonlinear oscillator networks, e.g. : thus DBs are of broad interest. More extensive references to the literature, covering both theoretical and experimental aspects, can be found in the cited publications.
We study the discrete sine-Gordon equation under damping and driving,

$$\frac{\text{d}^2u_n}{\text{d}t^2}=\frac{1}{2\pi }\mathrm{sin}\left(2\pi u_n\right)\alpha \frac{\text{d}u_n}{\text{d}t}+F\mathrm{sin}(\omega _0t)+C\left(u_{n+1}+u_{n1}2u_n\right)+\sigma \xi _n$$
(1)
which is already written in terms of dimensionless quantities. This model describes a chain of pendula in a homogeneous gravitational field, neighbouring pendula being harmonically coupled. Here $`u_n`$ is the angle of deviation from the equilibrium position in units of $`2\pi `$, $`\alpha `$ measures the strength of the damping, $`F\mathrm{sin}(\omega _0t)`$ is a harmonic external driving, $`C`$ gives the strength of the coupling between nearest-neighbour sites and the $`\xi _n`$ are Gaussian white noise terms of vanishing mean and $`\xi _n(t)\xi _m(t^{})=\delta _{mn}\delta (tt^{})`$, $`\sigma `$ scaling the strength of the noise. The noise represents thermal fluctuations, the relation between $`\sigma `$ and the temperature $`T`$ given by $`\sigma ^2=2\alpha k_BT`$. A more appropriate treatment of the white noise could be achieved in terms of increments of Wiener processes , which also would avoid the $`\delta `$-singularity in its correlation. We do, however, not use the formalism of stochastic differential equations in this article, therefore we do not deem it necessary to introduce these expressions here. Equation (1) can for instance also model a parallel array of Josephson junctions , in which case $`u_n`$ is the phase difference between the two superconductors in the $`n`$-th junction, given in units of $`2\pi `$ and the driving force corresponds to a bias current. In  mobile intrinsic localised modes have been found to exist in model (1) (without the noise term) as attractors of the dynamics in certain regions of parameter space. Correspondingly, the velocity of the mobile ILMs depends on the system parameters ($`C,\alpha ,\omega ,F`$), but, after a transient, not on the precise initial conditions. At some values of the parameters more than one basin of attraction exists, each corresponding to a different value of the velocity. In all cases there is a degeneracy of the model due to symmetry, as for a given velocity $`v`$ there are always basins of attraction corresponding to $`+v`$ and $`v`$. The diffusive behaviour of the localised excitation under the influence of the white noise $`\sigma \xi `$ can be described, as we will show in this article, by a hopping of the system between these two degenerate basins of attraction. At the parameter values we are using in the simulations for this article, $`\alpha =0.02,F=0.02,\omega _0=0.2\pi ,C=0.890`$ (a choice which to our knowledge is generic within the parameter regime where mobile ILMs exist), there are known three different basins of attraction with different (absolute) values of the velocities. We solve equation (1) numerically using the Heun-algorithm, with a time step $`\mathrm{\Delta }t=0.01`$, which corresponds to $`0.001\times 2\pi /\omega _0`$. It turns out in simulations that transitions between basins corresponding to different velocities do not play a role in the process, so we focus only on one value of the velocity, in our case $`v0.0186`$.
A noise induced flip of the velocity, i.e. a transition from one attracting configuration to a different one, involves a probably complex dynamical reconfiguration of the system. Unable at present to exactly describe this process, we aim at approximately capturing the main features of the transition in terms of a simple model. We assume that for fixed parameters there is a constant probability $`\mathrm{\Theta }`$ per unit time for a transition to be excited. At least in the ranges of parameter space we have explored numerically, we have found the plausible result that the transition probability increases with the noise strength (temperature), see figure 1.

For the sake of completeness we add here the comment that breathers also can be destroyed or created by thermal fluctuations, or form out of generic energy distributions in a system; in the temperature range and on the time scales we have used here, such processes have not been observed, and therefore will not play a role in this work. Related discussions can be found for instance in . We have to take into account that the reconfiguration accompanying a velocity flip is not an instantaneous process. Rather than considering details of this process we introduce a time scale $`\tau _A`$, which can be considered the (average) duration of the reconfiguration. We model the transition process in the following way: With probability $`\mathrm{\Theta }`$ per time unit a transition is initiated. If this occurs, in our simplified description the velocity of the ILM immediately assumes the value $`0`$ and the breather remains immobile for a time span $`\tau _A`$; this drop of the velocity to $`v=0`$ suggested itself as the most tractable approximation to the behaviour of $`v`$ in a transition that flips the sign of the velocity. After the delay $`\tau _A`$ the velocity assumes the value opposite to its initial value. As during the ‘waiting time’ $`\tau _A`$ the system is already involved in reconfiguration, we do not allow any further jumps to be initiated during the waiting time. If therefore a jump has been initiated at time $`t_i`$ the probability for no further jump to occur at least up to time $`t`$ is

$$1H(tt_i\tau _A)+H(tt_i\tau _A)\mathrm{exp}\left[\mathrm{\Theta }(tt_i\tau _A)\right]$$
(2)
where $`H`$ is Heaviside’s step function. The probabilities for $`n`$ jumps in time $`t`$ can also be derived and enter in the calculation of the expectation values and variances of the position, which involve rather tedious algebra and can be found in . Here we are only interested in the long time behaviour, which can be obtained in a much simpler way. Long time in this context means observation times $`t\frac{1}{\mathrm{\Theta }}`$, the latter being the mean time before a jump of the velocity once the delay period $`\tau _A`$ has passed; note that for the parameter values we have been looking at our numerical results have shown that $`\tau _A<\mathrm{\Theta }^1`$. Thus we will be considering the case of stochastic trajectories with many jumps in the observation time $`t`$. Suppose that a jump has been initiated at $`t=0`$ and a second jump will be initiated at some time $`\tau _A<t^{}<t`$. Then the mean distance the ILM travels between the jumps is

$$\mathrm{\Delta }X=v\mathrm{\Theta }\underset{\tau _A}{\overset{t}{}}(t^{}\tau _A)e^{\mathrm{\Theta }(t^{}\tau _A)}\text{d}t^{}\frac{v}{\mathrm{\Theta }}$$
(3)
Here we have taken into account that the velocity is zero and jumps are forbidden for the time interval $`\tau _A`$ after the jump at $`t=0`$; also we have made use of the long time condition $`t\mathrm{\Theta }^1`$ to obtain the approximate equality. Furthermore we find

$$\mathrm{\Delta }X^2=v^2\mathrm{\Theta }\underset{\tau _A}{\overset{t}{}}(t^{}\tau _A)^2e^{\mathrm{\Theta }(t^{}\tau _A)}\text{d}t^{}2\frac{v^2}{\mathrm{\Theta }^2}$$
(4)
so that we get

$$\text{Var}\left(\mathrm{\Delta }X\right)\frac{v^2}{\mathrm{\Theta }^2}$$
(5)
Introducing the parity $`\pi (n)`$ of a natural number $`n`$, i.e. $`\pi (n)=+1`$ if $`n`$ is even and $`\pi (n)=1`$ if $`n`$ is odd, we can write for the distance travelled after $`N`$ jumps:

$$X_N=\underset{i=1}{\overset{N}{}}\pi (i)\mathrm{\Delta }X_i$$
(6)
where $`\mathrm{\Delta }X_i`$ is the length of the $`i`$-th segment of the path of the discrete breather, i.e. the segment terminated by the $`i`$-th jump. In (6) we obviously start with an ILM of positive velocity. We straightforwardly obtain

$$X_N=\frac{1\pi (N)}{2}\frac{v}{\mathrm{\Theta }}$$
(7)
and

$$\text{Var}(X_N)=\left[N\frac{1\pi (N)}{2}\right]\frac{v^2}{\mathrm{\Theta }^2}$$
(8)
We are more interested in the dependence of the mean value and variance on time than on the number of jumps. The mean number of jumps after time $`t`$ is

$$N=\frac{t}{\tau _A+\frac{1}{\mathrm{\Theta }}}$$
(9)
because the average time between subsequent jumps is the sum of the waiting time $`\tau _A`$ and the mean time $`\mathrm{\Theta }^1`$ the system stays in one attractor after the passage of the waiting time. As we consider the long time limit, we have $`N1`$, and therefore the relative weight of trajectories with an even number of jumps and those with an odd number of jumps will be approximately equal. In this way we arrive at the following results:

$$X(t)=\frac{v}{2\mathrm{\Theta }}$$
(10)
and

$$\text{Var}\left[X(t)\right]=\left[\frac{t}{\tau _A+\frac{1}{\mathrm{\Theta }}}\frac{1}{2}\right]\frac{v^2}{\mathrm{\Theta }^2}$$
(11)
the latter giving a diffusion constant $`D=v^2/[\mathrm{\Theta }^2(\tau _A+\mathrm{\Theta }^1)]`$. The nonvanishing mean in (10) is due to the asymmetry introduced by the choice of the initial condition in (6). When the average numbers of realisations in the positive and negative velocity attractors have equilibrated, the mean position does not change any more.

In figure 2 these results are compared with simulations. The numerical results for the mean value and the variance of the position have been obtained as averages over 2000 realisations. Of course there still occur fluctuations in these numerical averages. In order to obtain a definitive value for the mean of the position to be used in (10) we average the numerically obtained mean value over the final 10000 time units (which can be considered to be in the ‘long time regime’) and use this averaged value in equation (10) to calculate $`\mathrm{\Theta }`$. Likewise, the mean slope of the numerically obtained variance over the final 10000 time units and the calculated value of $`\mathrm{\Theta }`$ are introduced in equation (11) to obtain $`\tau _A`$. The two parameters are then used in the full model  to generate a prediction for the mean and the variance of the position. For the mean value we see from figure 2 that there is a distinct deviation between the numerical and analytical results in the transition range between short time and long time behaviour. This deviation occurs because of the crudeness of our model and would disappear if the details of the reconfiguration during a velocity flip were taken into account more fully. Numerical results show that this overshooting increases, if the ratio $`\tau _A/\mathrm{\Theta }^1`$ increases, i.e. if the time the configuration spends in transition increases relative to the time the system is in one of the attracting states. This can easily be understood: Our model does not take into account the details of the transitions between the attractors, so here the description is only approximate. Consequently, if the system spends more time in transition, our description gets less accurate. A closer look at the transition region reveals that the observed deviation is also responsible for the shift between the numerical result and the prediction for the variance in the

long time regime. The slope of the variance in the long time regime is well reproduced by the full model, as is the long time behaviour of the mean value.
We see from equation (11) that an increase in the temperature (meaning also an increase in $`\mathrm{\Theta }`$) reduces the slope of the variance as a function of time, which has been observed in simulations as well, see figure 3. Furthermore this result is also plausible, because an increase in $`\mathrm{\Theta }`$ implies an increase in the mean number of jumps that occur along a trajectory within a given time. Thus the distance a breather can travel away from its average position before the next flip of the velocity is reduced. In average for high $`\mathrm{\Theta }`$, the ILM is therefore confined to a small region of the system, where it is erratically moving back and forth. Of course, this confinement is not strict, but an issue of probabilities. With increasing time the probability that an ILM exits from a fixed part of the chain grows, but the higher the temperature the slower this growth is. We would like to point out that small noise is sufficient to effect this behaviour; all noise strengths $`\sigma `$ for which we have shown results in this article are below one percent of the driving amplitude $`F`$. This said, it may look as if by reducing the noise strength $`\sigma `$ to arbitrarily low values the slope of the variance of $`X`$, as appearing for example in figure 3, could be increased without limit. We have not observed this, instead we found behaviour like in figure 4. However, our observation times necessarily have been of limited extension, whereas our results here apply to the long time case, i.e. a large average number of jumps should occur within the observation time. By reducing $`\sigma `$ the transition probability between the attracting configurations of the system becomes smaller and smaller and we reach the point where, on average, only a few velocity flips occur within the limited observation time. In this case the long time condition is violated, and our results no longer apply. What would happen if we were to increase the observation time (and necessarily also the number of realisations) by several orders of magnitude, while lowering the noise strength to significantly smaller values than the ones we have used, is a question we deem unsuitable for our computational resources. Will there eventually appear a long time behaviour as for higher noise values or will discreteness effects like in figure 4 dominate the dynamics?
In summary we have shown that the presence even of small thermal noise can change the behaviour of a mobile intrinsic localised mode in the damped-driven sine-Gordon chain from deterministic translation with a well defined velocity to propagation determined by random flips of the sign of the velocity. With increasing temperature the number of flips increases and the mode is stochastically confined to a part of the system of decreasing size. These results clearly have to be taken into account in any experiments that aim at the detection of propagating ILMs in a system which can be described by the damped-driven sine-Gordon model, like for instance a parallel array of Josephson junctions under an ac bias current. As our description of the diffusion and qualitatively also the observed phenomenon itself depend only on the existence of attracting configurations corresponding to propagation velocities $`+v,v`$ of a localised mode, we expect analogous behaviour of propagating discrete breathers in other systems which show these characteristics.
MM and LV acknowledge support from the European Commission within the Research Training Network (RTN) LOCNET, contract HPRN-CT-1999-00163. LV also was partially supported by the Ministry of Science and Technology of Spain through grant BFM2002-02359.


text:
# Nucleon Momentum Distributions From a Modified Scaling Analysis of Inclusive Electron-Nucleus Scattering.

## Abstract

Inclusive electron scattering from nuclei at low momentum transfer (corresponding to $`x1`$) and moderate $`Q^2`$ is dominated by quasifree scattering from nucleons. In the impulse approximation, the cross section can be directly connected to the nucleon momentum distribution via the scaling function $`F(y)`$. The breakdown of the $`y`$-scaling assumptions in certain kinematic regions have prevented extraction of nucleon momentum distributions from such a scaling analysis. With a slight modification to the $`y`$-scaling assumptions, it is found that scaling functions can be extracted which are consistent with the expectations for the nucleon momentum distributions.

Quasielastic (QE) electron scattering can provide important information about the distribution of nucleons in nuclei. With simple assumptions about the reaction mechanism, functions can be deduced that should scale (i.e. become independent of momentum transfer), and which are directly related to the nucleon momentum distribution. The concept of $`y`$-scaling of the quasielastic response was first proposed in 1975. It was shown that in the plane wave impulse approximation (PWIA) a scaling function, $`F(y)`$, could be extracted from the inclusive cross section which was related to the nucleon momentum distribution. In the simplest approximation, the scaling variable $`y`$ is the initial momentum of the struck nucleon along the direction of the virtual photon. We determine $`y`$ from energy conservation assuming a spectator model of the interaction and neglecting the transverse momentum of the struck nucleon:

$$\nu +M_A=\sqrt{M_N^2+(y+q)^2}+\sqrt{M_{A1}^2+y^2},$$
(1)
where $`M_A`$ is the mass of the target nucleus and $`M_{A1}`$ is the ground state mass of the $`A1`$ nucleus (assumed to be in an unexcited state).

Measurements of inclusive electron-nucleus scattering from deuterium and heavy nuclei at $`x>1`$ have been performed at JLab up to $`Q^27`$ GeV2. At low $`Q^2`$ values the scaling function depends strongly on $`Q^2`$ due to final state interactions (FSIs). As these FSIs become small the extracted scaling function becomes nearly independent of $`Q^2`$ and depends only on $`y`$, as predicted in the $`y`$-scaling picture. However, while the data show scaling in $`y`$, this by itself does not ensure that the scaling function is connected to the momentum distribution. We present here an attempt to test the assumptions of the scaling analysis and the extraction of the nucleon momentum distributions.

Figure 1 shows $`F(y)`$ for deuterium, as extracted from the cross sections measured in E89-008. As the momentum distribution is related to the derivative of $`F(y)`$, the lack of high precision data on deuterium at large $`|y|`$ makes it difficult to directly extract the momentum distribution from this data. We can, however, compare the scaling function to what we expect based on a calculation of the deuteron momentum distribution. The solid line is a calculation of $`F(y)`$ using a momentum distribution calculated from the Argonne-v14 N-N potential. The normalization of the scaling function extracted from the data is consistent with unity (as it must be if it is related to the momentum distribution) and the distribution is in generally good agreement with the calculation. In particular, they are in very good agreement at very large values of the nucleon momentum, ($`y`$). This region is especially important because these high momentum components are generated by short range interactions of the nucleons. It has been suggested that final state interactions in this region, where the nucleons are close together, may not disappear as $`Q^2`$ increases. If there were large final state interactions that were nearly independent of $`Q^2`$, one might see scaling but the scaling function would not yield the proper momentum distribution in the tails. The data from deuterium indicate that such $`Q^2`$-independent FSIs are small or absent, although higher precision data at high $`Q^2`$ and large $`|y|`$ would allow a much stronger limit to be set.

While the $`y`$-scaling analysis of the deuterium data appears to yield the correct deuteron momentum distribution, this is not the case for the heavier nuclei. The momentum distribution extracted from $`F(y)`$ for heavy nuclei falls off much more rapidly at large $`y`$, indicating that the high momentum components in heavy nuclei are much smaller than in deuterium, which is the opposite of what one might expect. In addition, the normalizations of the scaling functions for heavy nuclei are $``$20-30% lower than they should be if $`F(y)`$ is related to the nucleon momentum distribution. These problems indicate that there is a failure of some kind in the scaling analysis for heavy nuclei. The breakdown for $`A>2`$ nuclei comes from the assumption that the residual $`(A1)`$ nucleus remains in an unexcited state. This is a reasonable approximation when removing a single nucleon from a shell at low missing energy. However, the high momentum nucleons are predominantly generated by short range correlations, meaning that the momentum of the struck nucleon is mostly balanced by a single nucleon, leaving a high momentum nucleon in the residual nucleus. In the following analysis, we take this into account by assuming a simple three-body breakup of the nucleus, where the struck nucleon is assumed to be one of a correlated pair of nucleons moving within the residual $`(A2)`$ nucleus. The scaling variable in this case is $`y^{}=k+K_{2N}/2`$, where $`y^{}`$ is the total momentum of the struck nucleon, coming from the relative momentum of the two correlated nucleons, $`k`$, and the momentum of the pair within the residual nucleus, $`K_{2N}`$.

Figure 2 shows $`F(y^{})`$ from iron, along with the fit to the deuterium scaling function (note that for deuterium there is no $`(A2)`$ residual nucleus, so $`y^{}=y`$). The high momentum behavior is identical for deuterium and heavy nuclei (carbon, iron, and gold), indicating that the two-nucleon correlations that dominate in deuterium are the main source of high momentum components in heavy nuclei. Using the modified scaling variable, the normalization of $`F(y^{})`$ is also consistent with unity, as it should be if the scaling function is related to the nucleon momentum distribution.

While this data indicates that the modified scaling analysis is valid and allows extraction of the nuclear momentum distributions, the data at large nucleon momentum is somewhat limited, especially for few-body nuclei where the extracted distributions can be compared to essentially exact calculations of nuclear structure. Future measurements are planned with 6 GeV beam which will significantly increase the amount of data in the scaling region at large nucleon momenta. This data will significantly improve the data at large $`|y|`$ and $`Q^2`$, and will include measurement on 3He and 4He. We can then use this to extract information on the momentum distributions in heavy nuclei, and study in more detail the nature of their short range correlations.

This work is supported (in part) by the U.S. DOE, Nuclear Physics Division, under contract W-31-109-ENG-38.