text:
# Characterizing the reflectance of periodic layered media

## I Introduction

Wave propagation in periodic layered media is very similar to the motion of electrons in crystallyne solids (hence, they are also known as photonic crystals). Thus, it is hardly surprising that some well-established concepts of solid-state physics (such as Bloch waves, Brillouin zones, and forbidden bands) have been incorporated in the doctrine of photonic crystals JO95 .

The standard approach YE88  identifies the periodic structure with a one-dimensional lattice that is invariant under lattice translations, which requires the system to be strictly infinite. Of course, this is unattainable in practice. In solid-state physics one goes around this difficulty by considering that the one-dimensional lattice closes into a ring; i.e., by imposing periodic boundary conditions. While this seems to be reasonable for a lattice of atoms, is obviously unfeasible for periodically stratified media.

In spite of this fact, only under the assumption of an endlessly repeating stack of dielectric slabs this standard treatment gets an eigenvalue equation for the Bloch factor and then one finds YE77 ; FI98 ; DO98 ; YA98 ; LE00  that strong reflection (stop bands) occurs when the trace of the $`2\times 2`$ matrix describing the basic period exceeds 2 in magnitude.

In practice not very many periods (say, ten) are needed to have a stop band. Therefore, one is led to consider stacks of $`N`$ periods, which are often appropriately called “finite” periodic structures LE94 . The matrix describing these structures can be computed as the $`N`$th power of the basic period in terms of Chebyshev polynomials BO99 ; LE87 . This gives an exact expression for the reflectance, but so complicated that it is difficult to get, at first glance, any feeling about the performance of the system. Moreover, the distinction between allowed and forbidden bands is obtained only numerically.

The aim of this paper is to provide an alternative setting for dealing with these finite periodic structures that stresses the different scaling laws appearing in the reflectance. To this end we resort to the Iwasawa decomposition to represent the action of any multilayer as the product of three fundamental matrices of simple interpretation. We show that these matrices present quite distinct reflectances. The trace of the basic period allows us to classify periodic structures in three types with properties closely related to one (and only one) of the three fundamental matrices. Finally, we study the reflectance associated to each one of these types.

The advantage of this approach is, in our view, twofold: it naturally emphasizes the different laws for the forbidden bands, allowed bands, and band edges and, at the same time, provides universal expressions remarkably simple for the corresponding reflectances.

## II The origin of the bandgap in infinite periodic structures

To keep the discussion as self-contained as possible, we first briefly summarize the essential ingredients of multilayer optics that we shall need for our purposes. The basic period of the structure is illustrated in Fig. 1 and consists of a stack of $`1,\mathrm{},j,\mathrm{},m`$ plane-parallel layers sandwiched between two semi-infinite ambient ($`a`$) and substrate ($`s`$) media, which we shall assume to be identical, since this is the common experimental case. Hereafter all the media are supposed to be lossless, homogeneous, and isotropic.

We consider an incident monochromatic linearly polarized plane wave from the ambient that makes an angle $`\theta _0`$ with the normal to the first interface and has amplitude $`E_a^{(+)}`$. The electric field is either in the plane of incidence ($`p`$ polarization) or perpendicular to the plane of incidence ($`s`$ polarization). We consider as well another plane wave of the same frequency and polarization, and with amplitude $`E_s^{()}`$, incident from the substrate at the same angle $`\theta _0`$.

As a result of multiple reflections in all the interfaces AZ87 , we have a backward-traveling plane wave in the ambient, denoted $`E_a^{()}`$, and a forward-traveling plane wave in the substrate, denoted $`E_s^{(+)}`$. If we take the field amplitudes as a vector of the form

$$\stackrel{}{E}=\left(\begin{array}{c}E^{(+)}\\ E^{()}\end{array}\right),$$
(1)
which applies to both ambient and substrate media, then the amplitudes at each side of the multilayer are related by the $`2\times 2`$ complex matrix $`𝖬_{as}`$, we shall call the transfer matrix, in the form

$$\stackrel{}{E}_a=𝖬_{as}\stackrel{}{E}_s.$$
(2)
The matrix $`𝖬_{as}`$ can be shown to be of the form YO02 ; MO02

$$𝖬_{as}=\left(\begin{array}{cc}1/T_{as}& R_{as}^{}/T_{as}^{}\\ R_{as}/T_{as}& 1/T_{as}^{}\end{array}\right)\left(\begin{array}{cc}\alpha & \beta \\ \beta ^{}& \alpha ^{}\end{array}\right),$$
(3)
where the complex numbers

$$R_{as}=|R_{as}|\mathrm{exp}(i\rho ),T_{as}=|T_{as}|\mathrm{exp}(i\tau ),$$
(4)
are, respectively, the overall reflection and transmission coefficients for a wave incident from the ambient. Note that we have $`det𝖬_{as}=+1`$, which is equivalent to $`|R_{as}|^2+|T_{as}|^2=1`$, and then the set of lossless multilayer matrices reduces to the group SU(1,1), whose elements depend on three independent real parameters MO99a ; MO99b .

The identity matrix corresponds to $`T_{as}=1`$ and $`R_{as}=0`$, so it represents an antireflection system without transmission phase shift. The matrix that describes the overall system obtained by putting two multilayers together is the product of the matrices representing each one of them, taken in the appropriate order. So, two multilayers, which are inverse, when composed give an antireflection system.

When we repeat the basic period $`𝖬_{as}`$ the fields at the $`n`$th cell are given by

$$\stackrel{}{E}_{n+1}=𝖬_{as}\stackrel{}{E}_n.$$
(5)
In an infinite structure the Floquet theorem ensures that these cells are equivalent and so the two vectors in Eq. (5) are proportional YE88 ; LE94

$$\stackrel{}{E}_{n+1}=\lambda \stackrel{}{E}_n,$$
(6)
and therefore the Bloch factor $`\lambda `$ is determined as the eigenvalue of the basic-period transfer matrix $`𝖬_{as}`$:

$$(𝖬_{as}\lambda )\stackrel{}{E}_n=0,$$
(7)
which gives

$$\lambda ^2[\mathrm{Tr}(𝖬_{as})]\lambda +1=0.$$
(8)
If $`[\mathrm{Tr}(𝖬_{as})]^2>4`$, the solution grows or decays exponentially and no propagating wave is possible (stop band). On the contrary, Bloch waves appear in this algebraic approach when $`[\mathrm{Tr}(𝖬_{as})]^2<4`$ and then the matrix $`𝖬_{as}`$ has two eigenvalues (allowed band). The band edges are obtained by setting $`[\mathrm{Tr}(𝖬_{as})]^2=4`$. Such a periodic structure exhibits resonance reflection very much like the diffraction of X rays by crystal lattice planes and, in consequence, it is usually called a Bragg reflector YE88 . By properly designing the basic period it is possible to achieve extremely high reflectance for some selected spectral region SO99 .

## III Finite periodic structures and Iwasawa decomposition

As anticipated in the Introduction, an infinite multilayer is unrealistic in practice, and one should consider $`N`$-period finite structures, for which the previous standard approach fails BE96 ; CO01 .

The overall transfer matrix $`𝖬_{as}^{(N)}`$ for a $`N`$-period structure is equal to $`(𝖬_{as})^N`$. Explicit but complicated formulas for the corresponding matrix elements exist BO99 ; LE87 , but they do not give analytic expressions for the bandgap and are difficult to analyze.

We have recently advocated using the Iwasawa decomposition to study layered media MO01b , since it provides a remarkable factorization of the matrix representing any multilayer (no matter how complicated it could be) as the product of three matrices of simple interpretation. The decomposition can be stated as

$$𝖬_{as}=𝖪(\varphi )𝖠(\xi )𝖭(\nu ),$$
(9)
where

$`𝖪(\varphi )`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{exp}(i\varphi /2)& 0\\ 0& \mathrm{exp}(i\varphi /2)\end{array}\right),`$ (12)
$`𝖠(\xi )`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{cosh}(\xi /2)& i\mathrm{sinh}(\xi /2)\\ i\mathrm{sinh}(\xi /2)& \mathrm{cosh}(\xi /2)\end{array}\right),`$ (15)
$`𝖭(\nu )`$ $`=`$ $`\left(\begin{array}{cc}1i\nu /2& \nu /2\\ \nu /2& 1+i\nu /2\end{array}\right).`$ (18)
The parameters $`\varphi ,\xi `$, and $`\nu `$ are given in terms of the elements of the multilayer matrix by

$`\varphi /2`$ $`=`$ $`\mathrm{arg}(\alpha +i\beta ),`$
$`\xi /2`$ $`=`$ $`\mathrm{ln}(1/|\alpha +i\beta |),`$ (19)
$`\nu /2`$ $`=`$ $`\mathrm{Re}(\alpha \beta ^{})/|\alpha +i\beta |^2,`$
and their ranges are $`\xi ,\nu `$ and $`2\pi \varphi 2\pi `$. Now, we can interpret the physical action of the matrices appearing in Eq. (9). $`𝖪(\varphi )`$ represents the free propagation of the fields $`\stackrel{}{E}`$ in the ambient medium through an optical phase thickness of $`\varphi /2`$, which reduces to a mere shift of the origin of phases. Alternatively, this can be seen as an antireflection system. The second matrix $`𝖠(\xi )`$ represents a symmetric system with real transmission coefficient $`T_𝖠=\mathrm{sech}(\xi /2)`$ and reflection phase shift $`\rho _𝖠=\pm \pi /2.`$ Finally, the third matrix, $`𝖭(\nu )`$, represents a system having $`T_𝖭=\mathrm{cos}(\tau _𝖭)\mathrm{exp}(i\tau _𝖭)`$ and $`R_𝖭=\mathrm{sin}(\tau _𝖭)\mathrm{exp}(i\tau _𝖭)`$, with $`\mathrm{tan}(\tau _𝖭)=\nu /2`$. In Ref. YO02  we have provided simple examples of realistic systems that accomplish these requirements.

For our purposes here, the essential point is that the reflectance $`=|R|^2`$ associated to each one of these matrices is

$`_𝖪`$ $`=`$ $`0,`$
$`_𝖠`$ $`=`$ $`\mathrm{tanh}^2(\xi /2),`$ (20)
$`_𝖭`$ $`=`$ $`(\nu /2)^2/[1+(\nu /2)^2].`$
While $`_𝖪`$ is identically zero, $`_𝖠`$ and $`_𝖭`$ tend to unity when $`\xi `$ and $`\nu `$, respectively, increase. However, they have distinct growth: $`_𝖠`$ goes to unity exponentially, while $`_𝖭`$ goes as $`O(\nu ^2)`$.

The Iwasawa decomposition proves also to be a powerful tool for the classification of layered media. Indeed, after Ref. SA01  one is led to introduce the following criterion: a matrix is of type $`𝖪`$ when $`[\mathrm{Tr}(𝖬_{as})]^2<4`$, is of type $`𝖠`$ when $`[\mathrm{Tr}(𝖬_{as})]^2>4`$, and finally is of type $`𝖭`$ when $`[\mathrm{Tr}(𝖬_{as})]^2=4`$. Although this trace criterion has an elegant geometrical interpretation MO02  and coincides with the one giving the stop bands in Eq. (8), let us remark that if a multilayer has a transfer matrix $`𝖬_{as}`$ of type $`𝖪`$, $`𝖠`$, or $`𝖭`$, one can always find a family of matrices $`𝖢`$ \[also in SU(1,1)\] such that

$$\widehat{𝖬}_{as}=𝖢𝖬_{as}𝖢^1$$
(21)
is just a matrix $`𝖪(\varphi )`$, $`𝖠(\xi )`$, or $`𝖭(\nu )`$, respectively. This conjugation by $`𝖢`$ preserves the trace: $`\mathrm{Tr}(\widehat{𝖬}_{as})=\mathrm{Tr}(𝖬_{as})`$. Moreover, Eq. (21) allows us to recast the multilayer action in Eq. (2) as

$$\widehat{\stackrel{}{E}}_a=\widehat{𝖬}_{as}\widehat{\stackrel{}{E}}_s,$$
(22)
where the new field vectors are

$$\widehat{\stackrel{}{E}}=𝖢\stackrel{}{E}.$$
(23)
In other words, the matrix $`𝖢`$ gives a new vector basis such that the basic period of the system, when viewed in such a basis, presents a reflectance exactly of the form in Eq. (III). This provides a full characterization of the reflectance of any periodic system, as we shall demonstrate in next Section.

## IV Characterizing the three fundamental behaviors of the reflectance

Since the matrix $`\widehat{𝖬}_{as}`$ associated to the basic period belongs to one of the subgroups $`𝖪(\varphi )`$, $`𝖠(\xi )`$, or $`𝖭(\nu )`$ of the Iwasawa decomposition, and all these subgroups are, in our special case, Abelian and uniparametric, we have that

$$\widehat{𝖬}_{as}(\mu _1)\widehat{𝖬}_{as}(\mu _2)=\widehat{𝖬}_{as}(\mu _1+\mu _2),$$
(24)
where $`\mu `$ represents the appropriate parameter $`\varphi `$, $`\xi `$, or $`\nu `$. For a $`N`$-period system the overall transfer matrix $`𝖬_{as}^{(N)}`$ is

$$𝖬_{as}^{(N)}=𝖢^1[\widehat{𝖬}_{as}(\mu )]^N𝖢=𝖢^1\widehat{𝖬}_{as}(N\mu )𝖢,$$
(25)
which does not depend on the explicit form of the basic period.

After Eq. (25), one must expect three universal behaviors of the reflectance according the basic-period transfer matrix is of the type $`𝖪`$, $`𝖠`$ or $`𝖭`$. We shall work in what follows the detailed structure of these three basic laws.

Because the stop bands are given by the condition $`[\mathrm{Tr}(𝖬_{as})]^2>4`$, we first consider the case when $`𝖬_{as}`$ is of type $`𝖠`$. Then, Eq. (21) gives

$$𝖬_{as}=𝖢^1𝖠(\chi )𝖢,$$
(26)
where we have denoted

$$\mathrm{Re}(\alpha )=\frac{1}{2}\mathrm{Tr}(𝖬_{as})=\mathrm{cosh}(\chi )>1,$$
(27)
because we are taking into account only positive values of $`\mathrm{Re}(\alpha )`$, since negative values can be treated much in same way. If we put

$$𝖢=\left(\begin{array}{cc}𝔠_1& 𝔠_2\\ 𝔠_2^{}& 𝔠_1^{}\end{array}\right),$$
(28)
one solution of Eq. (26) is

$$𝔠_1=F(\beta ^{}+i\mathrm{sinh}\chi ),𝔠_2=iF\mathrm{Im}(\alpha ),$$
(29)
where the factor $`F`$ is

$$F=\frac{1}{\sqrt{2\mathrm{sinh}\chi [\mathrm{sinh}\chi \mathrm{Im}(\beta )]}}.$$
(30)
Carrying out the matrix multiplications in Eq. (25) it is easy to obtain the reflectance of the $`N`$-period system as

$$_𝖠^{(N)}=\frac{|\beta |^2}{|\beta |^2+[\mathrm{sinh}(\chi )/\mathrm{sinh}(N\chi )]^2}.$$
(31)
This is an exact expression for any value of $`N`$. As $`N`$ grows, $`_𝖠^{(N)}`$ approaches unity exponentially, as one could expect from a band stop.

The band edges are determined by the limit condition $`[\mathrm{Tr}(𝖬_{as})]^2=4`$; that is, when $`𝖬_{as}`$ is of type $`𝖭`$. A calculation very similar to the previous one shows that now

$$_𝖭^{(N)}=\frac{|\beta |^2}{|\beta |^2+(1/N)^2},$$
(32)
with a typical behavior $`_𝖭^{(N)}1O(N^2)`$ that is universal in the physics of reflection. The general results (31) and (32) have been obtained in a different framework by Yeh YE88  and Lekner LE94 .

Finally, in the allowed bands we have $`[\mathrm{Tr}(𝖬_{as})]^2<4`$; $`𝖬_{as}`$ is of type $`𝖪`$, and

$$_𝖪^{(N)}=\frac{𝒬^22𝒬\mathrm{cos}(2N\mathrm{\Omega })}{1+𝒬^22𝒬\mathrm{cos}(2N\mathrm{\Omega })},$$
(33)
where

$$𝒬=\frac{|\beta |^2}{|\beta |^2|\alpha e^{i\mathrm{\Omega }}|^2},$$
(34)
and the phase $`\mathrm{\Omega }`$ is determined by

$$e^{i\mathrm{\Omega }}=(\mathrm{Re}\alpha )+i\sqrt{1(\mathrm{Re}\alpha )^2}.$$
(35)
Now the reflectance oscillates with $`N`$ between the values $`(𝒬^22𝒬)/(𝒬1)^2`$ and $`(𝒬^2+2𝒬)/(𝒬+1)^2`$.

Equations (31), (32) and (33) are the three fundamental reflectances we were looking for. In our opinion, they have the virtue of simplicity. In practice, the basic period usually consists of two homogeneous layers and the expressions become simpler BO99 . We emphasize that the scaling laws expressed by these equations are universal features of the reflection in finite periodic systems.

## V Conclusions

We have used a trace criterion to provide in a direct way three fundamental forms for the reflectance that are universal for any multilayer. When the system is finite periodic this result has allowed us to obtain in a direct way explicit expressions for the reflectance of the overall system.

This approach has physical meaning on its own, irrespective whether the system is periodic or not. Moreover, it leads naturally to a further analysis of finite periodic bandgap structures that does not need the strict requirement of infinite periodicity, as it is assumed in the standard theory.

We acknowledge illuminating discussions with José F. Cariñena.


format: text

text:
# On the Evolution and Appearance of a Surviving Companion after a Type Ia Supernova Explosion

## 1 Introduction

In recent years, Type Ia supernovae (SNe Ia) have been used successfully as cosmological probes of the Universe (Riess et al. 1998; Perlmutter et al. 1999). However, the nature of their progenitors has remained somewhat of a mystery. There is almost universal agreement that they represent the disruption of a degenerate object; but this is also where the agreement ends. There are numerous progenitor models (for detailed reviews see, e.g., Branch et al. 1995; Ruiz-Lapuente, Canal & Isern 1997), but most of these have serious theoretical/observational problems or do not appear to occur in sufficient numbers to explain the observed frequency of SNe Ia in our Galaxy ($`3\times 10^3`$yr-1; Cappellaro & Turatto 1997).

The progenitor models can roughly be divided into three classes: double-degenerate (DD) models, sub-Chandrasekhar models and single-degenerate Chandrasekhar models. The double-degenerate model (Iben & Tutukov 1984; Webbink 1984) involves the merger of two CO white dwarfs with a combined mass in excess of the Chandrasekhar mass ($`1.4\text{M}_{}`$). While this model has the advantage of being quite common (see, e.g., Iben & Tutukov 1986; Yungelson et al. 1994; Han et al. 1995; Iben, Tutukov & Yungelson 1997; Han 1998; Nelemans et al. 2001), it seems more likely that the disruption of the lighter white dwarf and the accretion of its debris onto the more massive one leads to the transformation of the surviving CO white dwarf into an ONeMg white dwarf which subsequently collapses to form a neutron star (i.e. undergoes accretion-induced collapse; AIC) rather than experience a thermonuclear explosion (Nomoto & Iben 1985; Saio & Nomoto 1985, 1998; Timmes, Woosley & Taam 1994; Mochkovitch, Guerrero & Segretain 1997). There may be a small parameter range where AIC can be avoided, but it is unlikely to account for more than a small number of SNe Ia.

Woosley & Weaver (1995) speculated that a sub-Chandrasekhar-mass white dwarf could produce a SN Ia if helium was ignited violently in a shell surrounding the CO core and triggered a detonation wave that propagated inwards and ignited the CO core (see also Livne & Arnett 1995). While it is not clear whether this can even work in principle, simulations of such explosions show that they do not reproduce the properties of observed SNe Ia (see, e.g., the discussion in Wheeler 1996).

The arguably most favoured class of models at the present time involves single-degenerate scenarios, where the white dwarf accretes from a non-degenerate companion star (Whelan & Iben 1973; Nomoto 1982). In these models, the companion star can be in different evolutionary phases and may either be a hydrogen-rich star or a helium star. One of the major problems with these models is that it is generally difficult to increase the mass of a white dwarf by accretion due to the occurrence of nova explosions and/or helium flashes (Nomoto 1982) which may eject most of the accreted mass. There is a narrow parameter range where a white dwarf can accrete hydrogen-rich material and burn it in a stable manner, but this requires rather special circumstances.

One promising channel that has been identified in recent years relates them to supersoft X-ray sources (van den Heuvel et al. 1992; Rappaport, Di Stefano & Smith 1994; Li & van den Heuvel 1997). In this channel, the companion star is a somewhat evolved main-sequence star or subgiant of 2 – 3 $`\text{M}_{}`$, transferring mass on a thermal timescale to a white dwarf. In Figure 1 we present a typical binary calculation that illustrates this channel (also see Langer et al. 2000; Han & Podsiadlowski 2003). The initial system consists of a 2.1$`\text{M}_{}`$ somewhat evolved main-sequence star and a 0.8$`\text{M}_{}`$ white dwarf. In this calculation we adopted the estimates for the accretion efficiency of Hachisu, Kato & Nomoto (1999) and assumed that any excess matter that could not be acrreted by the white dwarf was ejected from the system. As the calculation demonstrates, the white dwarf is able to grow very effectively. When it reaches the Chandrasekhar mass, it has parameters very similar to U Scorpii, a supersoft binary and recurrent nova where the white dwarf is already close to the Chandrasekhar mass (Hachisu et al. 2000; Thoroughgood et al. 2001).

While U Scorpii provides an excellent candidate for a SN Ia progenitor, consistent with theoretical expectations, it would be even better to have a more direct observational test for progenitor models (e.g. Ruiz-Lapuente 1997), in particular since it is quite possible – perhaps even likely – that there is more than one channel that leads to a SN Ia. Such direct tests could involve the detection of hydrogen or helium in the ejecta or the supernova environment, which could come from the outer layers of the exploding object, circumstellar material that was ejected from the progenitor system (e.g. Cumming et al. 1996), or matter that was stripped from the secondary by the supernova interaction and was mixed into the ejecta. A particularly conclusive test would be the detection of a companion star that has survived the supernova explosion in the supernova remnant. At present the detection of a surviving companion would only be feasible in our Galaxy, e.g. in the supernova remnants SN 1006 and Tycho (Ruiz-Lapuente et al. 2002).

A surviving companion can be identified by its expected unusual properties (Marietta, Burrows & Fryxell 2000; Canal, Méndez & Ruiz-Lapuente 2001): (1) it will be a runaway star moving directly away from the centre of the supernova remnant, where the runaway velocity will be similar to the orbital velocity of the companion star in the pre-supernova binary (e.g. $`170\text{ km\hspace{0.17em}s}\text{-1}`$ for a system like U Sco; $`20\text{ km\hspace{0.17em}s}\text{-1}`$ for a symbiotic binary with a giant companion); (2) the secondary will have been strongly affected by its interaction with the supernova ejecta and have an unusual appearance; and (3) the composition of the secondary may show chemical anomalies and be enriched with heavy elements produced in the thermonuclear explosion (similar to what is observed in Nova Sco, a black-hole binary where the companion star has been enriched with supernova material; Israelian et al. 1999; Podsiadlowski et al. 2002).

Marietta et al. (2000) have recently performed the most detailed and most systematic hydrodynamical investigation of the impact of the supernova shell in a SN Ia with hydrogen-rich secondaries of various types, representing different progenitor channels (for earlier studies see Wheeler, Lecar & McKee 1975; Fryxell & Arnett 1981; Taam & Fryxell 1984; Livne, Tuchman & Wheeler 1992). In particular, they showed that main-sequence and subgiant companions will lose $`15`$per cent of their mass due to the supernova impact, while giant companions will be stripped of almost all of their hydrogen-rich envelopes. Furthermore, the interior of the secondaries will be significantly shock-heated by the supernova impact; as a consequence the immediate post-impact luminosity will be increased dramatically.

It is the purpose of this paper to model the further evolution of a subgiant secondary to clarify the principles that govern their post-impact evolution and to determine their expected appearance which will be useful for searches of surviving companions in historical SN Ia remnants which are presently under way. In Section 2 we systematically examine the consequences of heating and mass stripping for the evolution and appearance of the secondary and in Section 3 we discuss the implications of these results.

## 2 Post-supernova heating and stripping

Immediately after the impact with the supernova ejecta, the secondary will be strongly out of thermal equilibrium and is likely to be significantly puffed up and overluminous (Marietta et al.  estimate that a subgiant may have a luminosity as high as $`5000\text{L}_{}`$). The subsequent evolution is governed by how the star re-establishes thermal equilibrium, which depends on the amount of heating and mass loss it has experienced. Because these effects have very different consequences for the evolution of the star, the secondary may appear overluminous or underluminous at different times during the equilibration phase, as we will demonstrate below.

For this study we employed a standard Henyey-type stellar evolution code, described in detail in Podsiadlowski, Rappaport & Pfahl (2002), which we have previously used to study the effects of external irradiation and tidal heating in stars (Podsiadlowski 1991, 1996). Since we cannot follow the dynamical interaction of the secondary during the impact phase, we adopt the following procedure to model the effects of mass stripping and supernova heating. We first perform a mass-loss calculation where we take off mass from the star at a very high rate (typically $`10^2\text{M}_{}\text{yr}\text{-1}`$), so that mass loss can be considered essentially adiabatic (except in the outermost layers). To model the residual heating of the envelope (i.e. the extra heating of the part of the envelope that was not stripped off in the supernova impact), we add a constant and uniform heating source to the remaining outer layers of the star with a total luminosity of $`10^5\text{L}_{}`$ for a certain of time and then follow the subsequent re-equilibration of the star. This procedure is of course only very approximate, but should nevertheless be sufficiently realistic to study the basic physical aspects of the problem, the main objective of this paper111We have tested that the subsequent evolution of the star is not very sensitive to the timescale on which the initial stripping and heating is simulated (as long as it is short compared to the thermal timescale of the outer layers of the star); while the initial appearance may differ substantially for different implementations, the evolutionary tracks and internal structures tend to converge rapidly in the subsequent re-equilibration phase..

In Figure 2 we present evolutionary tracks in the Hertzsprung-Russell (H-R) diagram illustrating the post-impact equilibration phase for a subgiant (with an initial pre-supernova mass of 1$`\text{M}_{}`$) that was stripped of the outermost 0.2$`\text{M}_{}`$ of its envelope by the supernova impact and was heated by various different amounts of energy (as indicated at the beginning of each track). The initial subgiant has a radius of 2.5$`\text{R}_{}`$ and is very similar to the model of the secondary shown in Figure 1 at the time of the supernova. In these simulations, the energy was deposited uniformly in the outermost 90 per cent of the radial extent of the star which contained 0.57$`\text{M}_{}`$. The energy in the calculation with the largest amount of energy deposited ($`4\times 10^{47}`$erg) is close to the energy needed to unbind these layers completely, which we estimate to be $`4.9\times 10^{47}`$erg. As a consequence, the star is initially very puffed up with a radius of 230 $`\text{R}_{}`$ and has a luminosity of 4300$`\text{L}_{}`$ (cf. Marietta et al. 2000). Because of the large radius and the large luminosity, the thermal timescale of the puffed-up envelope is very short initially, and the star contracts very rapidly along a Hayashi track. After $`10^3`$ yr its luminosity has dropped to 200$`\text{L}_{}`$ and after $`10^4`$ yr to 27$`\text{L}_{}`$. In this context we emphasize that in all the models shown in Figure 2 the evolution is initially much faster than the Kelvin-Helmholtz timescale, $`t_{\mathrm{KH}}`$, of the pre-supernova subgiant ($`10^6`$yr in our case), since the timescale on which the star evolves in the H-R diagram is determined by the thermal timescale of the rapidly evolving outer layers of the star, which can be many orders of magnitude shorter than $`t_{\mathrm{KH}}`$. (Of course, it will take a full Kelvin-Helmholtz time before the star has re-established thermal equilibrium completely.)

In the calculations with $`10^{47}`$ and $`2\times 10^{47}`$erg, the evolution is similar although less dramatic. On the other hand, in the case where the star is only heated by $`5\times 10^{46}`$erg or is not heated at all, the evolution is dramatically different, and the secondary becomes significantly underluminous during the early thermal re-equilibration phase (with a luminosity of $`10^1\text{L}_{}`$ after $`10^3`$yr and $`10^2\text{L}_{}`$ after $`10^4`$yr). This behaviour is similar to the behaviour commonly found in binary mass-transfer calculations where a secondary strongly out of thermal equilibrium can be very underluminous (e.g. Langer et al. 2000; King et al. 2001; Podsiadlowski et al. 2002).

The principal reason for this behaviour is that, in a subgiant whose outer layers have been suddenly stripped off, the layers in the deep interior of the envelope have to expand significantly since the final thermal equilibrium radius is similar to the initial radius. The energy needed to drive this expansion comes from the nuclear luminosity produced in the core of the star, reducing the radiative luminosity that reaches the outer parts of the star. In the case where the star is dramatically out of thermal equilibrium, as is the case here, essentially all the nuclear luminosity is used up in the expansion of the envelope. As a consequence, the very outer layers of the star are no longer re-supplied with energy from below; but because of the continued radiative losses from the surface, the outer layers start to cool dramatically.

To illustrate this behaviour, Figure 3 shows the early evolution of the thermodynamic structure in the outer parts of the envelope for the stripped subgiant in Figure 2 that has not been heated (the region shown contains about 1 per cent of the mass of the star). As the outer layers cool (particularly clearly seen in the drop of the entropy), this outer layer contracts, and the star in fact becomes somewhat undersized relative to its final equilibrium radius (it takes a full Kelvin-Helmholtz time, $`10^6`$yr in the case shown, for the star to re-establish thermal equilibrium). During this cooling phase, the star evolves roughly along a line of constant radius in the H-R diagram and would have the appearance of an anomalously red subgiant. Since the surface luminosity becomes much smaller than the luminosity inside the star, only a small fluctuation of the luminosity profile (e.g. due the generally non-linear behaviour of the equilibration process which can drive luminosity waves through the star) may lead to a sudden change of the surface luminosity, and we therefore suspect that such a star would be highly variable222In our calculations we regularly observe sudden changes of the luminosity profile, where the star suddenly brightens by an order of magnitude or more in a single time step and subsequently continues to cool. However, at present this behaviour is not well resolved numerically..

Because our models were constructed in a very different fashion from the models of Marietta et al. (2000), it is not entirely clear which of our sequences most closely resembles their subgiant case. The fact that the subgiant is stripped of a significant fraction ($`15`$per cent) of its mass certainly implies that the remaining envelope will also be strongly heated. Comparing the radii of our models to the radius ($`9\text{R}_{}`$) found by Marietta et al. (2000) then suggests that the sequences with $`10^{47}`$ and $`2\times 10^{47}`$erg of heating are probably the most comparable. In these sequences, the secondary resembles a typical subgiant with a luminosity of 10 –25$`\text{L}_{}`$ after $`10^3`$yr and 6 – 10$`\text{L}_{}`$ after $`10^4`$yr, still significantly overluminous relative the initial model which had a luminosity of $`3\text{L}_{}`$.

The timescale on which the surface luminosity evolves after the supernova impact is determined mainly by the thermal timescale of the envelope which depends, among other factors, on the evolutionary phase of the initial subgiant, the amount of mass that has been transferred to the companion before the supernova and the amount of mass stripped by the supernova impact itself. To illustrate this sensitivity we performed a similar series of thermal equilibration calculations for the same initial subgiant but where we stripped off the outermost 0.8$`\text{M}_{}`$, and heated the remaining outer 83 per cent in radius containing 0.03 $`\text{M}_{}`$, by various amounts of energy from $`5\times 10^{45}`$erg to $`3.5\times 10^{46}`$erg; the latter is close to the energy needed to unbind the envelope in this case ($`3.7\times 10^{46}`$erg). Because of the lower remaining mass, the star evolves much more rapidly. In Figure 2 the shaded (yellow) region shows the location of the secondary in the H-R diagram between $`10^3`$ and $`10^4`$ yr in these calculations. Note that, in all cases, the secondary is significantly underluminous (by up to a factor of 10) – even in the strongly heated case – and may have the appearance of a main-sequence star similar to the Sun (but with much lower gravity).

## 3 Discussion and Conclusions

In this paper we have examined the evolution of a typical subgiant secondary (as expected in one of the most popular progenitor models for SNe Ia) that has been heated and stripped by the impact of a supernova shell. We have demonstrated that the secondary may be significantly underluminous or overluminous at different times during the post-impact phase when it tries to re-establish thermal equilibrium. This depends on the amount of heating it has experienced and the amount of envelope mass that has been stripped off by the impact or has been lost by mass transfer before the explosion. Based on our calculations, we typically expect its luminosity $`10^3`$yr after the impact to be within an order magnitude of the pre-supernova luminosity (either above or below). If the outer envelope cools dramatically (which happens if most of the internal luminosity goes into driving internal expansion), the secondary may become extremely underluminous. While it would be quite faint, it would be a rather unusual, extremely red and probably highly variable subgiant that should stand out because of its anomalous appearance. Even if it has the appearance of a more typical subgiant, where it may in fact be somewhat hotter (more yellow) than a standard subgiant, one may be able to see changes in its appearance on a human timescale because of the short thermal timescale of the outermost layers (which may be many orders of magnitude shorter than the Kelvin-Helmholtz timescale of the pre-supernova star).

The main implication for searches of surviving companion stars in historical supernova remnants is that, in order to be able to rule out a subgiant companion with any degree of confidence, one has to go to much fainter magnitudes than previously assumed (down to a luminosity of $`10^1\text{L}_{}`$ for a $`10^3`$yr-old remnant, a luminosity of $`10^2\text{L}_{}`$ for a $`10^4`$yr-old one). We still expect that it will have an appearance and/or show signs of variability that makes it easily identifiable as the supernova companion (in addition to its runaway space velocity and possible chemical anomalies).

If no plausible candidate for a subgiant companion is found at this level or for any other type of pre-supernova companion, this would give some weight to other, presently less favoured progenitor models where no companion is expected (e.g. the double-degenerate merger model) or where the companion may already have become a very cool white dwarf (provided that there is a long time delay between the white-dwarf accretion phase and the supernova, as expected in some models).

## Acknowledgements

This work was initiated during a workshop on Type Ia supernovae in July 2002 at the Lorentz Center in Leiden organized by B. Schmidt and N. Langer. The author thanks the organizers and many other participants for numerous fruitful discussions and thanks the Center for its support. This work was also in part supported by a European Research and Training Network Grant on Type Ia Supernovae (HPRN-CT-20002-00303).


text:
# Threshold Electrodisintegration of 3He

## I Introduction

During the past decade much experimental and theoretical progress has been made in the study of the trinucleon system—the first non-trivial test of the adequacy of phenomenological $`NN`$ potentials. Especially instructive tests are provided by photo- and electro-disintegration reactions. For example, Carlson et al.Carlson02  recently reviewed the status of the 3He quasielastic response functions measured in inclusive electron scattering far from the elastic scattering peak. In this complementary work we report upon the status of the kinematic region near the break-up threshold of $`5.5`$ MeV.

In the early seventies the electrodisintegration of the deuteron near break-up threshold was recognizedHockert73  as a decisive test of the understanding of meson exchange in the traditional picture of the NN force. This reaction is unusually informative because the wave functions of the initial and final states are relatively simple and well-known and, if the electron is deflected to far-backward angles, the break-up is dominated by a pure M1, $`\mathrm{\Delta }T=1`$ transition. The contribution of meson exchange currents (MEC) generally grows with increasing three-momentum transfer $`q`$: at $`q=2.5`$ fm-1 MEC raise the threshold cross section by about a factor of $`3`$; near $`q=3.5`$ fm-1 MEC account for nearly 100% of the transverse cross section due to destructive interference between the one-body transition amplitudes.

Although the threshold electrodisintegration of 3He aroused similar interest, it took another three decades before Viviani et al.Viviani00  were finally able to confirm the importance of MEC in the trinucleon break-up. Earlier, Hadjimichael et al.Hadjimichael83  had established the need for MEC in the elastic cross sections, but due to the requirement of knowing not just the ground state wave functions, but also those of the continuum, the break-up poses a more challenging test. Following quickly on the paper by Viviani et al., additional evidence for MEC was given in two lettersXu00 ; Xiong01  reporting asymmetry measurements for longitudinally-polarized electrons scattered from a polarized 3He target. The firstXu00  of these measurements was performed near the quasielastic peak, where MEC effects are small. Stronger evidence for MEC was given in a subsequent paperXiong01  on the threshold region, where measurements at $`q=1.60`$ and $`2.27`$ fm-1 were presented. Near threshold the effect of MEC on the spin-dependent asymmetry is calculated to be large, and although the measurements strongly support this prediction, the agreement is not exact. On the other hand, spin-dependent asymmetries represent an especially demanding test of nuclear theory.

The delayed confirmation of significant MEC effects in the trinucleon break-up stems from the recent parallel developments in precise empirical NN potentials and powerful theoretical methods for calculating exactly the $`A=3`$ wave function. These requirements have now been met, with such success that we can now claim a detailed understanding of most of the basic properties of the trinucleon, at least at low-to-moderate energies and momenta.

The advances in NN potentials were allowed by high-quality measurements and analyses of pp and np scattering. Precise non-relativistic potentialsStoks94 ; Wiringa95 ; Machleidt96  were constructed that fit the vast databases with $`\chi ^2`$-per-datum values close to unity. In addition to the usual charge-independent parts, charge-dependence and asymmetry terms were introduced to account for pp and np scattering simultaneously. The electromagnetic parts of these potentials contain Coulomb, Darwin-Foldy, vacuum polarization, and magnetic moment terms with finite-size properties. Although the calculations shown in this paper rely upon just oneWiringa95  of these potentials, the Argonne AV18, for the properties investigated here little sensitivity would be expected to the differences between these modern potentials.

The theoretical techniques devised to solve the three-body Schrödinger equation are described in two comprehensive reviewsGlockle96 ; Carlson98 . Monte Carlo methods, Faddeev techniques, and variational procedures that utilize correlated hyperspherical harmonics have all been successfully employed. Because theoretical predictions for the trinucleon electrodisintegration are sensitiveMeijgaard89  to final state interactions, precise representations are needed not only for the ground state, but also for the final continuum states.

The calculations of Viviani et al.Viviani00 , made within the pair-correlated hyperspherical harmonics scheme, use the AV18 two-nucleon potential, supplemented by the Urbana-IX three-nucleon interactionPudliner95 . Because the calculations assume a pd final state, they are confined to the narrow excitation region between the two- and three body break-up thresholds at $`E_x=5.5`$ and $`7.7`$ MeV. Both longitudinal and transverse response functions were calculated for 3He at $`q=0.88`$, $`1.64`$, and $`2.47`$ fm-1, three-momentum transfer values that correspond to experimental results obtained by Retzlaff et al.Retzlaff94  at the MIT-Bates accelerator. The effect of MEC is largest at the highest momentum transfer, $`q=2.47`$ fm-1, where the predicted transverse response function $`R_T`$ is doubled by including MEC. Even though the experimental points have $`35`$% uncertainties, the scale of the MEC enhancement is so large that the importance of exchange currents is firmly established.

In this study we present: (1) A new measurement of the threshold transverse response function of 3He, made at the highest momentum transfer probed by Retzlaff et al., but with uncertainties that provide a more rigorous test of the theoretical predictions; (2) Theoretical calculations that include both two- and three-body break-up. These include a self-consistent treatment of final-state interactions and exchange currents; (3) Additional inclusive scattering measurements at $`q=4.4`$ and $`4.7`$ fm-1, a kinematic region sometimes speculated to mark the onset of the transition to quark-gluon dynamics; and (4) A new appraisal of the threshold longitudinal response function, facilitated by the improved information on the transverse.

## II Experiment and Analysis

The new measurements were made at the MIT-Bates Linear Accelerator Center at electron scattering angle $`\theta =160^{}`$, an angle where the cross section is dominantly transverse. Useful continuum data were obtained at effective incident beam energies of $`E_0=263`$, $`506`$, and $`549`$ MeV, corresponding to three-momentum transfers at the two-body break-up threshold of $`q=2.4`$, $`4.4`$, and $`4.7`$ fm-1. Many of the details of this experiment have been previously published in a report on the 3He elastic magnetic form factorNakagawa01 . To recapitulate, the target system contained 4000 STP liters of 3He cooled to $`23`$ K and pressurized to $`50`$ atm. In order to mitigate variation in the 3He density due to beam heating, the gas flow was highly turbulent, an enlarged beam spot was used, and the beam current was held constant at $`19\pm 1\mu `$A. Scattered electrons were detected in a magnetic spectrometer system that included drift chambers for trajectory information, a gas Cerenkov detector and lead-glass shower counter for particle identification, and three layers of plastic scintillators for triggering and timing.

Figure 1 shows the threshold cross section measured with $`263`$ MeV incident electrons, plotted as a function of excitation energy. Experimental backgrounds have been removed, and corrections applied for dead-time losses and detector inefficiencies. The 3He elastic peak and its calculated radiative tail have also been subtracted. As may be seen, in this kinematic region the elastic radiative tail is small. The two- and three-body break-up thresholds are indicated on the figure. Due to the experimental energy resolution, which results mainly from the straggling of electrons travelling different paths within the thick gas target, the continuum cross section begins to rise just before the two-body break-up threshold.

To expedite the comparison of the data to theoretical predictions, corrections were also made for energy lost by the electron in radiative processes, occurring before, after, or during the primary electro-nuclear interaction. In our experiment the overall effect of these processes is to decrease the cross sections measured near threshold.

Radiation corrections were applied using the continuum unfolding procedure of Mo and TsaiMo69  as implemented by MillerMiller71 . Even though our target material is low-Z, the radiative corrections turn out to be large. In principle, these are precisely calculable, but this requires data more extensive than our limited measurements. Hence we have, in part, had to rely on approximations and models to evaluate the radiative corrections. Nevertheless, the uncertainties in the calculated corrections are expected to be always smaller than the statistical uncertainties of the data. Figure 1 shows the result of radiative unfolding for the spectrum measured at $`E_0=263`$ MeV; The table lists the radiation-unfolded spectra for all three beam energies.

## III Results

### III.1 Systematics of the transverse response
function $`R_T`$

In this section we examine the systematic dependence of available experimental information on $`R_T`$ as a function of incident beam and excitation energies. The radiatively-unfolded, inclusive electron scattering cross sections depend on the longitudinal and transverse response functions according to

$`{\displaystyle \frac{d^2\sigma (E_0,E_x)}{d\mathrm{\Omega }dE^{}}}={\displaystyle \frac{d^2\sigma _{Mott}}{d\mathrm{\Omega }dE^{}}}`$
$`\times \left\{R_L(q,E_x)+\left[{\displaystyle \frac{1}{2}}+\mathrm{tan}^2{\displaystyle \frac{\theta }{2}}\right]R_T(q,E_x)\right\}.`$
In the threshold region only Retzlaff et al.Retzlaff94  and Köbschall et al.Kobschall83  have published separated longitudinal and transverse functions for 3He. Nevertheless, due to the large size of the transverse kinematic factor at backward scattering angles, additional information on $`R_T`$ is provided by the spectra measured at $`\theta =180^{}`$ by Jones et al.Jones79  and by our $`\theta =160^{}`$ spectra. Chertok et al.Chertok69  published one additional $`180^{}`$ spectrum, however this lacks continuum radiative corrections, and, in any case, is superseded by the later measurements of Jones et al. made at the same laboratory. Additional spectra were obtained by Kan et al.Kan75 , but at smaller scattering angles where the longitudinal response function is more strongly weighted than the transverse.

For $`180^{}`$ electron scattering the longitudinal response function is negligible and the transverse response can be deduced using

$`R_T(E_0,E_x)\left({\displaystyle \frac{2E_0}{\alpha }}\right)^2{\displaystyle \frac{d^2\sigma (E_0,E_x)}{d\mathrm{\Omega }dE^{}}}.`$
Note that the response obtained in this way is given as a function of incident beam energy $`E_0`$, not the three-momentum transfer $`q`$. However, for the electron beam energies considered here, $`q`$ changes slowly in the vicinity of the break-up threshold.

At the $`160^{}`$ angle of our Bates data, longitudinal contributions to the (e,e) cross section are still small. Based on the results of Retzlaff et al. and our new calculations, longitudinal contributions to the spectrum measured at $`E_0=263`$ MeV are $`<3.5`$%. These estimated components have been subtracted from the data. No similar allowance has been made for longitudinal contributions to data obtained at $`506`$ and $`549`$ MeV, but according to our calculations these are even smaller.

Selected results for $`R_T`$ are compiled as a function of $`q`$ and $`E_x`$ in Fig. 2. As noted above, for the results of Jones et al.Jones79  and our $`160^{}`$ Bates experiment, the value of $`q`$ changes slowly with $`E_x`$, decreasing by about $`0.045`$ fm-1 between the break-up threshold at $`E_x=5.48`$ MeV and $`E_x=15`$ MeV. For these results the $`q`$-values indicated on the plot were calculated at the three-body break-up threshold of $`7.7`$ MeV.

In the threshold region $`R_T`$ peaks at $`q1`$ fm-1, and by $`q=4.5`$ fm-1 has decreased by 4 orders of magnitude. For $`q<1`$ fm-1 there is a tendency for $`R_T`$ to be broadly peaked in the range $`E_x=1020`$ MeV. This resembles the distribution of resonant E1 strength seen in photoabsorption measurements, and which is convincingly explained by very recent Faddeev calculationsGolak02  that use the AV18 NN potential. Indeed, E1 strength will also be large in inclusive electron scattering at $`q<1.5`$ fm-1, although, at low-$`q`$ multipoles other than E1 are predicted to make sizeable near-threshold contributions. For example, early calculations of two-body electrodisintegration by Heimbach et al.Heimbach77  indicated considerable M2 strength in the region $`E_x<20`$ MeV at $`q0.5`$ fm-1. Additional smaller contributions were obtained from the M1 and M3 multipoles.

For $`q>2`$ fm-1, the near-threshold $`R_T`$ increases monotonically with increasing $`E_x`$. At still higher momentum transfers quasifree scattering becomes the dominant reaction mechanism and, notwithstanding resonance effects, final-state interactions, and phase-space suppression close to threshold, the monotonic rise seen in the data taken at $`E_0=506`$ and $`549`$ MeV has the appearance of the high-momentum tail of the quasielastic peak. This tail is of considerable interest since it provides information on elusive high-momentum components of the nuclear wave functionSick80 . In order to test the quasielastic hypothesis we examined the data to see if the cross section scales with $`y`$, the initial momentum component parallel to q that would be carried by a quasielastically scattered nucleon. Such scaling is the signature of quasifree scatteringSick80 ; Day90 . As shown in Fig. 3, even though the momentum transfer is sufficiently high, our results lie outside the band that corresponds to asymptotic $`y`$-scaling. We conclude that, even at the relatively large $`q`$ of our measurements at $`506`$ and $`549`$ MeV, excitation energies of $`25`$ MeV are insufficient to assure the dominance of quasifree scattering.

### III.2 Comparison of $`R_T`$ results with Faddeev calculations that include final-state interactions

The non-relativistic calculations to which we compare our new spectra are similar to those presented in a previous paperGolak95 . As described there, bound and continuum pd and ppn wave functions were obtained by solving Faddeev-like integral equations in momentum space. All final-state interactions are rigorously included. The present calculation is improved in two ways. First, rather than the older-generation Bonn potential, we use the updated Argonne AV18 NN interaction, and second, we include MEC contributions, evaluated using Riska’s prescriptionRiska85 . Most importantly, for these calculations the final-state interactions and exchange currents are fully consistent with the NN interaction potential.

The theoretical calculations were performed on the Cray SV1 of the NIC in Jülich, Germany and the NERSC Computational Facility, USA. Despite the computational power of these facilities, the long CPU times required for the calculations limited what could be achieved. The main results are plotted in Fig. 4. For $`263`$ MeV the agreement with the data could scarcely be better. These new theoretical predictions may also be compared to the pd breakup calculations by Viviani et al.Viviani00  that utilize the same NN potential, but which were carried out by means of pair-correlated hyperspherical harmonics, not by solving Faddeev equations. Motivated by the data of Retzlaff et al.Retzlaff94 , Viviani et al,’s calculations were performed for a momentum transfer $`3.4`$% higher than that of our $`263`$ MeV measurements. Nevertheless, when this difference is taken into account using the $`q`$-dependence given by Viviani et al., it is found that the Faddeev and hyperspherical harmonic calculations, both with and without MEC, are almost indistinguishable in the $`E_x=5.57.7`$ MeV range where comparison is valid. It is reassuring to note that our new Faddeev calculation of the three body final state continues to agree well with the data up to the highest excitation energy. According to an earlier calculationGolak95  based on the previous-generation Bonn B interaction, the three-body contribution to $`R_T`$ grows relatively slowly above threshold: at $`E_x=18`$ MeV and $`q=0.88`$ fm-1, it amounts to just half the two-body contribution.

The agreement of the two calculations, in addition to the agreement with the data, does more than simply confirm the importance of MEC in the 3He break-up—it underscores how accurate these modern calculations can be—to a point. As indicated in Fig. 4, the theoretical prediction is less satisfactory for the higher beam energy of $`506`$ MeV. Here, the relatively flat contribution of the non-relativistic one-body current is modified—in the correct sense—by interference with MEC. Above $`E_x=19`$ MeV the interference is constructive; below $`19`$ MeV it is strongly destructive. (This contrasts strongly with the effect of MEC on the $`263`$ MeV spectrum.) Although this leads to the correct shape, the predicted $`R_T`$ is about a factor of $`2`$ too low throughout the threshold region.

Several factors may account for this discrepancy. For example, the use of non-relativistic dynamics is questionable at such energies. To investigate this we performed an exploratory calculation in which the one-body current only was treated relativistically, and even this in a manner formally inconsistent with the NN interaction. Figure 4 shows that this device leads to an even poorer prediction of the data. Moreover, the momentum transfers of our $`506`$ and $`549`$ MeV spectra probe small wave function components well beyond the Fermi momentum. These components are usually negligible, but their effects can be magnified in scattering at large $`q`$, as shown in an analysis of quasifree scattering by Sick, Day, and McCarthySick80 . From this work it was deduced that the “exact” 3He wave functions obtained from realistic NN interactions have high-momentum components that are too small.

A further indication of theoretical difficulties at large $`q`$ is evident in the observationNakagawa01  that the first diffraction minimum in the elastic magnetic form factor of 3He is located near $`q=4.2`$ fm-1, somewhat higher than predicted by current theories.

### III.3 $`q`$-dependence of $`R_T`$

Viviani et al.Viviani00  have calculated the $`q`$-dependence of $`R_T`$ at a fixed $`1.0`$ MeV excitation above the pd threshold. However, due to the small size of the experimental cross section near threshold, it is unpractical to compare to this prediction. A more reliable comparison may be made by integrating $`R_T`$ in the range $`5.5<E_x<7.7`$ MeV, where break-up is confined to the pd channel evaluated by Viviani. Our integration of the theory assumes a linear dependence for $`R_T`$ on $`E_x`$—an assumption supported by the near-threshold calculations.

The results are shown in Fig. 5. With the exception of the two high-$`q`$ points, excellent quantitative agreement is again obtained. As previously noted, the MEC contribution grows at higher $`q`$. For $`q<3.4`$ fm-1 MEC terms interfere constructively with the one-body matrix elements, raising the prediction by up to a factor of $`3`$. At higher $`q`$ the interference is destructive, an interpretation supported by the two high-$`q`$ points, irrespective of the lack of exact quantitative agreement: near $`q=4.5`$ fm-1 the data lie far below the one-body prediction.

The change in the interference from constructive to destructive agrees with our new theoretical predictions, as indicated in Fig. 4. But note that the destructive interference is confined to low excitation energies: at large $`E_x`$ the interference remains constructive, a prediction confirmed by the data.

### III.4 Longitudinal response function

Unlike the slow and monotonic rise of the transverse response function, the longitudinal response rises abruptly in the first $`2`$ MeV above thresholdRetzlaff94 ; Kobschall83 ; Kan75 , a feature attributedKan75  to a $`{}_{}{}^{2}S^2S`$ Coulomb monopole transition. This behaviour is illustrated in Fig. 6, which reproduces the experimental results of Retzlaff et al.Retzlaff94 , obtained at $`q=2.47`$ fm-1. Our calculation for this momentum transfer has the right shape, but exceeds the data by about a factor of $`2`$. In part, this is attributed to our neglect of the Coulomb barrier that would suppress the emission of low-energy protons. Indeed, as shown in the figure, the calculation of the pd electrodisintegration by Viviani et al.Viviani00 , which includes the Coulomb term, lies closer to the data.

Still better agreement is obtained by including two-body charge operators. As noted by Viviani et al., these operators have relativistic origins and should properly be evaluated by including, in a self-consistent way, relativistic effects in both the interaction models and the nuclear wave functions. Lacking such a method, the only recourse is to perform a model-dependent calculation. This contrasts with the evaluation of the two-body current operators that contribute to the transverse response where, according to the classification scheme of RiskaRiska85 , the main parts of these operators are fixed by current conservation.

Nevertheless, as Fig. 6 shows, even with the inclusion of the Coulomb and two-body charge terms, the best available theoretical result still exceeds the data by roughly $`50`$%. This is a small but notable blemish in what otherwise is a remarkably precise and comprehensive theoretical description of the 3He threshold photo- and electro-disintegration. It emphasizes the potential value of a more rigorous treatment for the two-body charge operators. Further evidence for this comes from the efforts of Schiavilla and collaboratorsCarlson98 ; Marcucci98 , to predict the charge form factors of $`A=3`$ and $`A=4`$ nuclei.

Figure 7 shows the threshold $`q`$-dependence of $`R_L`$. As for Fig. 5, we have integrated the experimental and theoretical response functions in the range $`5.5<E_x<7.7`$ MeV, where break-up is restricted to the pd channel. Our integration in this case takes note of the curvature in the dependence of $`R_L`$ on $`E_x`$. According to the calculations of Viviani et al.Viviani00 , this diminishes at large $`q`$. The experimental points in Fig. 7 include the results of Retzlaff et al.Retzlaff94 , Köbschall et al.Kobschall83 , and Kan et al.Kan75 . Kan et al. were unable to extract $`R_L`$ from their limited measurements, however, by virtue of the now-precise knowledge of $`R_T`$ at the $`q`$-values of their measurements, this separation becomes possible.

As the plot shows, for $`q<2`$ fm-1 the agreement between experiment and theory is generally excellent, while the significance of the discrepancy with Retzlaff’s point at $`q=2.47`$ fm-1 has been discussed in some length. The disagreement for the lowest-$`q`$ point of Kan et al. can perhaps be attributed to uncertainties inherent in a $`>50`$% transverse subtraction. For other points, the transverse contribution is typically $`<20`$%.

## IV Summary

For $`q<3`$ fm-1, exact non-relativistic calculations using the best available NN potentials give a very good description of measurements of the 3He threshold electrodisintegration. At low momentum transfer one-body matrix elements predominate, but as $`q`$ increases the two-body contribution grows significant, particularly in the transverse part of the cross section. For example, at $`q=2.4`$ fm-1, corresponding to one of the three new measurements reported here, the inclusion of two-body terms raises the predicted $`R_T`$ by a factor of two, bringing the theory into close agreement with the data—proof that our understanding of MEC is accurate.

According to the near-threshold theoretical predictions for $`R_T`$, at $`q3.5`$ fm-1 the interference between one- and two-body terms switches from constructive to destructive. This is supported by our two other measurements, made at $`q4.5`$ fm-1. These points fall about a factor of $`7`$ below the one-body prediction, but the destructive interference with two-body current terms lowers the prediction. That it still exceeds the data by a factor of two suggests the need for a more complete interference, but other factors cannot be overlooked.

For example, the high-$`q`$ measurements are in a kinematic region sensitive not only to relativistic effects, but also to high-momentum wave function components, which, as indicated by quasifree scattering results, may be too small in “exact” wave functions obtained from realistic NN interactions. A pointed indication of the importance of relativistic effects is found in the analysis of Viviani et al.Viviani00  of the longitudinal response function at $`q=2.47`$ fm-1. Close to threshold there exists a factor-of-two disagreement between the one-body predictions and experimental values of $`R_L`$. This disagreement is reduced, but not entirely resolved, by considering two-body charge matrix elements. These are equivalent to relativistic corrections. At this time estimates of these corrections are model-dependent: the discovery of a rigorous, self-consistent procedure for evaluating relativistic effects poses a considerable challenge. As has been repeatedly noted, a relativistic formulation of effective hadronic theory is essential for a satisfactory understanding of the transition from hadron to quark regimes.

Additional measurements are needed to guide the theoretical development, particularly above $`q=2`$ fm-1, where few data currently exist. An upcoming experiment in Hall A at Jefferson LaboratoryAniol00  is aimed at measuring the elastic form factors of 3He and 4He, starting at about $`4.5`$ fm-1. Unfortunately, broad energy resolution will limit what can be learned about the threshold break-up.

Especially valuable would be a new measurement of the 3H isobar, for which existing data are very sparse. One useful simplification offered by 3H is the absence of a Coulomb interaction—difficult to incorporate into Faddeev calculations—between the break-up products.

###### Acknowledgements.

This work was supported by U.S. DOE Grant DE-FG02-88ER40415, NSF Grants PHY9114958 and PHY0070858, Grant in Aid for International Research (Nos. 04640281 and 08044050) by the Ministry of Education, Science, and Culture in Japan, and by the Polish Committee for Scientific Research (grant No. 2P03B02818). The theoretical calculations were performed on the Cray SV1 of the NIC in Jülich, Germany and the NERSC Computational Facility, USA.


text:
# Retention of Two–Band Superconductivity in Highly Carbon–Doped MgB2

## Abstract

Tunneling data on MgB1.8C0.2 show a reduction in the energy gap of the $`\pi `$–bands by a factor of two from undoped MgB2 that is consistent with the $`T_c`$ reduction, but inconsistent with the expectations of the dirty limit. Dirty–limit theory for undoped MgB2 predicts a single gap about three times larger than measured and a reduced $`T_c`$ comparable to that measured. Our heavily–doped samples exhibit a uniform dispersion of C suggestive of significantly enhanced scattering, and we conclude that the retention of two–band superconductivity in these samples is caused by a selective suppression of *interband* scattering.

The simple binary compound, MgB2,jNagamatsu01  appears to be an elegant example of a two–band superconductor with a high $`T_c`$ of 39 K, in which superconductivity due to strong electron–phonon coupling in the two–dimensional (2D) $`\sigma `$–bands is induced into the 3D $`\pi `$–bands by a combination of interband electron–phonon coupling and weak interband quasiparticle scattering. Several key signatures of this include the closing of two energy gaps hSchmidt02a ; hSchmidt03 ; pSzabo01 ; fGiubileo01 ; mIavarone02 ; rsGonnelli02b  at the same $`T_c`$, anomalous features in the specific heat and its field dependence fBouquet01a  and the direct observation of interband scattering in the tunneling spectrum.hSchmidt02a

However, two distinct superconductive gaps should only be seen in the clean limit, while in Ref. ayLiu01  a single intermediate gap with a lower $`T_c`$ is predicted in the dirty limit. The observed weak dependence of $`T_c`$ on residual resistivity is not obviously consistent with this latter expectation.iiMazin02b  Mazin *et al.* iiMazin02b  have used supercell band–structure calculations to conclude that the two gaps can be preserved in the presence of enhanced scattering provided the *interband* impurity scattering is extremely weak. As will be shown below, introducing additional scattering centers by doping enables us to reduce the mean free path below the coherence length and thereby offers a direct experimental test of the dirty–limit prediction.

In a two–band superconductor the term *dirty limit* needs to be used carefully, though. The mean free path is reduced by both *intra*– and *interband* scattering, but only *interband* scattering is expected to homogenize the two gaps and result in effective one–band behavior.iiMazin02b  Thus there are two different scenarios for doped MgB2: (i) impurities strongly enhance *interband* scattering, in which case the small gap, $`\mathrm{\Delta }_\pi `$, is expected to increase until it merges with the large gap, $`\mathrm{\Delta }_\sigma `$, at an intermediate value, or (ii) *interband* impurity scattering is minimally affected and two distinct gaps are preserved. In scenario (ii), the ratio $`\frac{2\mathrm{\Delta }_\pi }{k_BT_c}`$ for the $`\pi `$–band is expected to remain small ($`1.5`$ in undopedhSchmidt02a  MgB2), whereas it would increase in scenario (i). Thus determination of $`\frac{2\mathrm{\Delta }_\pi }{k_BT_c}`$ alone allows us to decide between these scenarios, and in this paper we present tunneling data showing that both $`\mathrm{\Delta }_\pi `$ (the smaller gap) and $`T_c`$ are simultaneously reduced by about a factor of two upon $`10`$% C–doping. This represents direct experimental evidence for the selective suppression of *interband* impurity scattering relative to *intraband* scattering.

Several groups wMickelson02 ; raRibeiro03 ; tTakenobu01 ; aBharathi02a ; zhCheng02 ; mParanthaman01a ; iMaurin02 ; kPapagelis02  have investigated the doping of MgB2 with C with conflicting results concerning the C solubility limit, the composition dependence of $`T_c`$ and the possible presence of phase separation. Much of the discrepancy is possibly a synthesis problem due to the slow kinetics of C incorporation. Mickelson *et al.*,wMickelson02  and more recently Ribeiro *et al.*,raRibeiro03  have been able to overcome the diffusion problem by using a C containing starting material, B4C, and obtained samples with sharp diffraction peaks and superconducting transitions. The synthesis using B4C yields multiphase samples containing another borocarbide, MgB2C2, since the solubility limit of C in MgB2 is less than the 20 atom% present in the starting B4C.

Our samples were synthesized from high–purity Mg, graphite, and Eagle–Picher 11B enriched to 99.52%. A low carbon content boron carbide sample with a B0.85C0.15 composition was first made by arc melting the B and graphite together. This material was ground to a coarse powder, mixed with excess Mg and fired at 1000C for 5 hr as described in Ref. dgHinks02 . After the excess Mg was removed by heating to 400C in vacuum, powder neutron diffraction revealed only sharp diffraction peaks (Fig. 1) and the ac susceptibility showed a sharp transition with a midpoint $`T_c`$ of 21.5 K (Fig. 2). Note the absence of broadening of the diffraction lines precludes phase separation that was reported in Ref. iMaurin02 ; kPapagelis02  for highly C–doped MgB2. Rietveld analysis of the neutron diffraction data yield lattice constants $`a=3.0528`$ Å and $`c=3.5219`$ Å with a .07 mole fraction of MgB2C2 impurity. The $`a`$–axis is much reduced from 3.0849 Å for the pure material, while the $`c`$–axis is relatively unchanged. MgB2C2 is still present in our samples (although at a lower concentration then in Refs. wMickelson02 ; raRibeiro03 ) since our starting C content is still larger than the solubility limit of C in MgB2. From the initial carbon content of the boron carbide sample and the final MgB2C2 impurity content, the carbon content of the C–doped MgB2 sample is determined to be MgB1.8C0.2, i.e., a 10 atom% substitution of C for B, with an accuracy of $`\pm 2`$ atom%. This defines the solubility limit that agrees with recent results of Avdeev *et al.*mAvdeev03

Our tunneling data are on superconductor–insulator–superconductor (SIS) break junctions hSchmidt02a ; hSchmidt03  and they commonly identify the smaller gap in the 3D $`\pi `$–bands of C–doped MgB2. An example is shown in Fig. 3 where the current–voltage characteristic, $`I(V)`$, is plotted as the line with open squares along with its differential conductance, $`\frac{dI}{dV}`$, that was generated numerically and is plotted as open circles. Our tunneling apparatus lOzyuzer98  uses Au tips that readily attach to a piece of MgB2 and afterwards the single tunnel junction formed between this attached piece of MgB2 and the bulk sample predominates the spectra, i.e., has by far the highest resistance.hSchmidt02a ; hSchmidt03  The data exhibit rather sharp coherence peaks near $`\pm 3`$ mV that are characteristic for SIS junctions,111Note the qualitative disagreement with an SIN fit using a gap of 2.6 meV and virtually zero lifetime broadening, $`\mathrm{\Gamma }`$. At 4.2 K, thermal effects alone would broaden the peaks to about half their observed height, as shown by the dashed line in Fig. 3. and an overall spectral shape that is well–reproduced by an SIS fit (solid line in Fig. 3) using $`\mathrm{\Delta }_\pi =1.45`$ meV and $`\mathrm{\Gamma }=0.25`$ meV. This ratio of $`\mathrm{\Gamma }/\mathrm{\Delta }_\pi 20`$% is consistent with our tunneling data on undoped MgB2.hSchmidt02a ; hSchmidt03

The Josephson pair–tunneling peak at zero bias is the most prominent difference with the SIS calculation that only considers single–electron, quasiparticle tunneling.222Note that the SIS calculation reveals a small hump at zero bias due to tunneling of thermally–activated quasiparticles. This is seen in high–resistance, low–$`\mathrm{\Gamma }`$ tunnel junctions in undoped MgB2,hSchmidt02a ; hSchmidt03  but like the Josephson current it is a feature that is only found in SIS junctions. The observed zero–bias peak seen in Fig. 3 is too large to be explained by thermal activation, but agrees with the Josephson pair tunneling that is seen in such low–resistance ($`210`$ $`\mathrm{\Omega }`$ in this case) junctions in undoped MgB2. Due to the exceptional strength of this feature in lower–resistance tunnel junctions, and its intimate relation to superconductivity, the disappearance of this zero–bias peak at higher temperature is a very accurate test of the junction–$`T_c`$. We traced this feature to high temperatures in a junction that at 4.2 K showed almost 20 times the normal–state background conductance at zero bias, and find it to merge into a flat background at $`T_c2223`$ K (Fig. 2).333Note that the zero–bias conductance signal is *not proportional* to the Josehson current, $`I_J`$, itself and its temperature dependence is not expected to follow $`I_J(T)`$. However, the disappearance of both quantities coincides and allows for an unambiguous determination of $`T_c`$. This is in excellent agreement with the bulk–$`T_c`$ from magnetization that for comparison is also shown in the same figure, and allows us to draw two conclusions. First, the zero–bias peak is indeed due to Josephson currents, thereby providing solid corroboration of SIS tunneling. Second, the junction–$`T_c`$ is equivalent to the bulk value, confirming that the tunneling results are not affected by surface effects. Higher–resistance junctions do not show the Josephson pair tunneling at zero bias, as is expected. A selection of these, taken on two different samples, are shown in Fig. 4 to indicate the reproducibility of the location of the coherence peak (energy gap) and the general shape of the data.

To explain why we see only the smaller energy gap in the 3D $`\pi `$–band, we follow the reasoning of Ref. hSchmidt02a . Because of the exponential decay of the electron wavefunction in a barrier, tunneling currents are generally dominated by electron momenta that are within a small angle (typically $`15`$ degrees) of the junction normal. Getting within 15 of the $`ab`$–plane in a randomly formed junction is highly improbable; thus while all tunneling data on MgB2 observe the smaller gap in the 3D $`\pi `$–bands, the 2D $`\sigma `$–band is seen only occasionally. When present, contributions from the 2D $`\sigma `$–band are easily distinguished in SIS data hSchmidt03  but are not seen in the data of Figs. 3 and 4, consistent with more than 95% of the junctions we have observed in undoped MgB2. We thus have an unambiguous determination of the smaller gap in the 3D $`\pi `$–bands.

There is a disparity between this small gap value, $`\mathrm{\Delta }_\pi 1.5`$ meV, and a large critical temperature, $`T_c22`$ K, that is crucial experimental evidence for two–band superconductivity. The gap ratio, $`\frac{2\mathrm{\Delta }_\pi }{k_BT_c}1.53`$, is severely reduced from the weak–coupling limit of $`3.53`$, a disparity that would only be increased by strong–coupling effects. This can only be explained by the presence of a second, larger energy gap, i.e., two–band behavior. We note, that $`\frac{2\mathrm{\Delta }_\pi }{k_BT_c}`$ for C–doped MgB2, $`1.53`$, is roughly the same as in undoped MgB2, $`1.49`$, i.e. critical temperature and gap scale together with C–doping.

We therefore conclude that two–band superconductivity persists in C–doped MgB2, consistent with specific heat data that show an anomaly reminiscent of that attributed to two–band effects in undoped MgB2.raRibeiro03  We further conclude that this result gives direct experimental evidence for the absence of *interband* impurity scattering, i.e. scenario (ii). Although a comparable reduction to $`T_c=25.4`$ K was predicted for the dirty limit,aBrinkman02  i.e., strong *interband* scattering, this would be connected to a single intermediate gap, $`\mathrm{\Delta }4.1`$ meV,aBrinkman02  which is about three times the observed value. From the qualitatively different behavior of the small gap we conclude that the lowering of $`T_c`$ in MgB1.8C0.2 is caused by changes in the electronic structure rather than by enhanced scattering.

This reduction of $`T_c`$ may be understood by the model of covalently bonded hexagonal B nets which are ionically bonded by Mg2+ ions.jmAn01 ; jKortus01  The attraction of the $`\pi `$–electrons to the Mg2+ ions lowers the energy of the $`\pi `$– relative to the $`\sigma `$–bands, transferring electrons from the $`\sigma `$– to the $`\pi `$–bands. It is these resultant holes in the $`\sigma `$–bands that are responsible for the high $`T_c`$. Substituting C into MgB2 might be expected to change the effective carrier concentration both by direct charge doping and by indirect charge transfer $`\pi \sigma `$ through a modification of the ionic system. To first order the large decrease in the $`a`$–axis on C doping results from direct electron doping of the bonding $`\sigma `$–bands, causing a net decrease in the hole content. The unchanged $`c`$–axis indicates little effect on the ionic bonding between Mg2+ ions and the B nets. Thus it is probably the lowered density of states (DOS) in the $`\sigma `$–bands resulting from direct charge doping by C that causes the large $`T_c`$ decrease. Medvedeva *et al.* niMedvedeva01b  have calculated the DOS for several different dopants in MgB2. For C they find a reduced DOS and lower $`T_c`$, predicting further that $`T_c`$ should fall to zero at a doping level of about 0.085 electrons/B atom as the $`\sigma `$–band falls below the Fermi energy. Yan *et al.* yYan02b  also find a reduced $`T_c`$ and a smaller B–B interatomic distance on C doping. All these predictions are in qualitative agreement with our data, although we find that superconductivity persists at 22 K for $`0.1`$ electrons/B atom.

The retention of two–band superconductivity at such a high doping level is surprising, as it requires the absence of sizable *interband* impurity scattering. Again, there are essentially two scenarios that allow us to understand this behavior: (a) doping does not increase *any* scattering, or (b) doping increases *intraband* scattering only, while *interband* scattering is unaffected. In scenario (a) the lack of any scattering implies the C forms a superlattice. Scenario (b) requires the selective suppression of *interband* scattering relative to *intraband* scattering, and below we will present evidence that supports scenario (b).

For impurities as dilute as the 10% C for 11B in our samples, the possibility of local ordering, e.g., staged replacement of, on average, every 5th B layer with an ordered (BC)0.5 layer cannot be ruled out, but the long–range order needed for scenario (a) would not be expected. Testing for this is difficult since the contrast between C and 11B is weak by any scattering technique. Electron scattering offers the highest sensitivity to ordering because the diffraction occurs from tiny single crystals in the powder sample. Transmission electron diffraction was performed on crystallites aligned along the  and  zone axes, and no evidence of any ordering of dopant ions was found. We therefore conclude, that C enters the lattice as randomly distributed impurities, and that the sharp $`T_c`$ and the narrow powder diffraction peaks (see Fig. 1) are the result of a solubility limit for C in MgB2 (that in the presence of excess C would be expected to result in a homogenous doping just at this limit) rather than a well–ordered C superlattice.

In the absence of C–ordering, the dopant atoms are expected to act as scattering centers, and we may estimate the effective mean free path from the C–concentration. Using a $`2\times 2\times 1`$ MgB2 supercell containing one C atom, equivalent to $`1/8=12.5`$% impurities, the average C–C in–plane distance is $`0.6`$ nm, about one order of magnitude less than the anisotropic coherence lengths in undoped MgB2 (that are 4 and 10 nm lLyard02 ). Such a significantly reduced mean free path is borne out by the increase in $`\frac{dH_{c2}}{dT}(T_c)`$ by about a factor of two upon similar C–doping.raRibeiro03

We therefore conclude that the retention of two–band superconductivity in the presence of heavy C–doping is a result of the selective suppression of *interband* impurity scattering, scenario (b). (Note that the absence of sizable *interband* scattering is a non–trivial prerequisite for the observation of two–band superconductivity even in undoped MgB2.) Mazin *et al.* iiMazin02b  studied scattering in MgB2 theoretically and argued that *interband* impurity scattering would be suppressed due to the disparity between $`\sigma `$– and $`\pi `$–wavefunctions. Applied to the case of C–doped MgB2, this suggests that the two–band nature of superconductivity (that requires the absence of sizable *interband* scattering) can be retained, while at the same time other properties (that depend on *intraband* scattering) may be suggestive of the dirty limit, e.g., show increased resistivity and $`\frac{dH_{c2}}{dT}(T_c)`$.

In conclusion, our tunneling results on MgB1.8C0.2 demonstrate the retention of two–band superconductivity in the presence of significant scattering, consistent with specific heat measurements of Ribeiro *et al.*raRibeiro03  We explain this by the selective suppression of *interband* scattering, as suggested by Mazin *et al.*iiMazin02b  Our experiment shows that this suppression is a surprisingly robust feature of MgB2 and strong enough to retain the two–band character even in the presence of a C–doping content sufficient to reduce $`T_c`$ by a factor of two.

The authors acknowledge discussions with D.J. Miller, and technical support from S. Short and H. Zheng. This research is supported by the U.S. Department of Energy, Basic Energy Sciences—Materials Sciences, under contract # W–31–109–ENG–38.


text:
# Deformation of crosslinked semiflexible polymer networks

## Abstract

Networks of filamentous proteins play a crucial role in cell mechanics. These *cytoskeletal* networks, together with various crosslinking and other associated proteins largely determine the (visco)elastic response of cells. In this letter we study a model system of crosslinked, stiff filaments in order to explore the connection between the microstructure under strain and the macroscopic response of cytoskeletal networks. We find two distinct regimes as a function primarily of crosslink density and filament rigidity: one characterized by affine deformation and one by non-affine deformation. We characterize the crossover between these two.

The study of biopolymer networks and gels lies at the heart of the understanding of the mechanical properties of the cytoplasm since the mechanical rigidity of the intracellular material is largely governed by the *cytoskeleton*, a complex network of filamentous proteins, cross links, and other associated proteins AlbertsEtc . A key player in this cytoskeleton is F-actin, which exhibits significant rigidity on the cellular scale. The material properties of such semi-flexible polymer networks also pose complex and current problems in polymer physics. Many of the most basic questions concerning these common and important networks, such as how they deform under stress, remain unanswered. The bending rigidity of such polymers introduces a new microscopic elastic parameter that can have consequences for the macroscopic elastic coefficients of the bulk, semi-flexible gel. This changes the traditional rubber elasticity model for the bulk properties of gels composed of cross-linked flexible polymers.

In this letter, we examine a simple model for crosslinked rods that not only allows us to quantitatively test the relationship between the microscopic and macroscopic elastic coefficients of a randomly cross-linked network, but also sheds light on the intimately related issue of the spatial distribution of the network strain. Among the most fundamental properties of polymer networks is the way in which they deform under stress. Since the classical theories of rubber elasticity Rubinstein  it has been suggested that this deformation is *affine*, *i.e.*, that the strain is uniform, as for a sheared Newtonian fluid at low Reynolds number. While this leads to relatively good agreement with experiments for rubber-like networks, much evidence and continuing efforts concern systematic deviations that point to non-affine network strains at a microscopic scale.

The assumption that the deformation field is affine down to length scales comparable to the smallest microscopic scales in the material is a great simplification that allows one to construct quantitative theories relating the macroscopic elastic constants of a gel to the microscopic properties of its constituent polymers. The validity of the affine approximation for this class of semiflexible polymer materials has, however, been the subject of some debate fcm1 ; frey ; satcher . Whether the deformation field is affine or not depends, of course, on length scale; clearly at the scale of the entire sample, all deformations are trivially affine when subject to simple shear. We explore whether this self-averaging property of the deformation field extends to shorter length scales for semiflexible networks. We show that the degree of non-affine strain is a function of length scale and degree of crosslinking. Specifically, we find that the range of non-affine strain can extend well beyond the *mesh size*, or *correlation length*, the typical separation between filaments DeGennes . This occurs near the point of rigidity percolation for the network, or for highly flexible filaments. This results in elastic moduli governed primarily by bending of filaments under non-affine strain, consistent with frey ; satcher . In contrast, we find that these networks become increasingly affine, even down to the smallest scales of the network, *e.g.* the mesh size, at high cross link density, high molecular weight, or for rigid filaments. Here, we find that the bulk elastic moduli converge to those predicted from affine theory fcm1 . We also quantify the degree of non-affine strain and show that this is, indeed, dependent on the length scale.

As we focus on the zero-frequency, or static properties of the system we may ignore the complexities of the network-solvent interaction. To isolate the importance of semi-flexibility upon network properties, we also ignore the complex, nonlinear response of the individual F-actin filaments MarkoSiggia ; fcm1 ; frey . We instead study in detail the dependence of the bulk shear modulus and Young’s modulus of the material upon the cross-link density of the polymer gel as well as the bending and extension moduli of the individual filaments. We model the network filaments via the Hamilitonian per unit length ($`\delta s`$) for a filament

$$\frac{\delta }{\delta s}=\frac{\mu }{2}\left(\frac{\delta l}{\delta s}\right)^2+\frac{\kappa }{2}\left(\frac{\delta \theta }{\delta s}\right)^2.$$
(1)
The first term takes into account the extensional deformation of the filament $`\delta l(s)`$ as a function of arc-length, $`s`$ with modulus $`\mu `$. The second term determines the energy stored in the filament due to bending: the local tangent of the filament makes an angle $`\theta (s)`$ with respect to the $`\widehat{x}`$-axis and the bending modulus of the filament is $`\kappa `$. Note that both of these terms are quadratic. We do not explore nonlinearities, such as buckling when compressed beyond the Euler instability, at the scale of individual filaments. While such nonlinear effects are expected at increasing strains and finite temperature fcm1 , we seek to understand here the fundamental properties of semi-flexible networks and thus address only the linearized version of the problem. We note, however, that thermal fluctuations can result in an effective modulus $`\mu \kappa ^2/(kTl^3)`$ for a segment of length $`l`$. Thus, we consider $`\mu `$ and $`\kappa `$ to be independent parameters, even though they are both determined by single filament elastic properties and geometry at $`T=0`$ LL . Clearly the full exploration of these networks at finite temperature presents an interesting challenge; understanding the $`T=0`$ mechanical properties of these networks, however, is the requisite first step towards this more ambitious program. In addition, we will point out specific aspects of our results that are likely to depend critically on our zero-temperature assumption. After all, it is now understood that for some mechanical properties of networks, the zero temperature presents a singular limit as in the case of rigidity percolation.

Real F-actin networks formed from an actin monomer solution have a complex geometry, arising from the dynamic growth and branching of filaments AlbertsEtc ; real\_networks . For computational efficiency, we ignore such complications and consider static, isotropic networks of monodisperse filaments of length $`L`$. Each filament is represented by a line segment deposited with random position and orientation to a two–dimensional rectangular shear cell. Intersections are identified as permanent, freely–rotating cross links, to mimic e.g. the attachments of double–headed myosin molecules in an ATP–deficient solution active\_exps . The mean distance between cross links is $`l_\mathrm{c}`$ as measured along a filament. Deposition continues until the desired cross link density $`L/l_\mathrm{c}`$ has been reached.

The network is represented by the set of mobile nodes $`\{𝐱_i\}`$ consisting of all cross links and midpoints between cross links (the latter so as to include the dominant bending mode). The total system energy $`(\{𝐱_i\})`$ is then expressed in terms of the $`\{𝐱_i\}`$ using a discrete version of (1). Within our linearized scheme, this $`(\{𝐱_i\})`$ is a high–dimensional paraboloid with a unique global minimum, corresponding to the state of mechanical equilibrium at $`T=0`$. For the initially unstressed network, this minimum corresponds to zero deformation. Depending on whether we wish to measure the shear modulus $`G`$ or the Young’s modulus $`Y`$, a shear or uniaxial strain $`\gamma `$ is applied across the periodic boundaries in a Lees–Edwards manner. This moves the global minimum to a new, non–trivial position, which we numerically find using the preconditioned conjugate gradient method. The stored energy per unit area can then be calculated, which is $`\gamma ^2/2`$ times $`G`$ or $`Y`$ LL  within our linear approximation, and hence the network modulus can be extracted. This procedure is repeated for different network realizations until a reliable estimate of $`G`$ or $`Y`$ has been attained. Further simulation details can be found in long .

Apart from the system size, there are three length scales in the problem: two geometric lengths $`L`$ and $`l_\mathrm{c}`$, and a third material length scale deriving from the stretching and bending modulii, $`l_\mathrm{b}=\sqrt{\kappa /\mu }`$. On dimensional grounds, the modulii can be written in the form $`G=\frac{\mu }{L}f(L,l_\mathrm{c},l_\mathrm{b})`$ (with a similar expression for $`Y`$), provided that the system size is sufficiently large that finite size effects can be ignored, as applies to all of the results presented here.

Results are given in Fig. 2 for the shear modulus $`G`$ and the Poisson ration $`\nu =2Y/G1`$. $`G`$ monotonically increases with the cross link density $`L/l_\mathrm{c}`$ and the ratio $`l_\mathrm{b}/L`$. Both $`G`$ and $`Y`$ simultaneously vanish at the rigidity percolation transition $`[L/l_c]_{\mathrm{crit}}5.9`$, irrespective of $`l_b`$. This is consistent with the more precise value $`[L/l_\mathrm{c}]_{\mathrm{crit}}5.932`$ found by Latva–Kokko et al. latva\_kokko  using the combinatorial “pebble game” method. We find that $`G`$ and $`Y`$ scale near the transition as $`(L/l_\mathrm{c}[L/l_\mathrm{c}]_{\mathrm{crit}})^f`$ with $`f=3.0\pm 0.2`$, placing it in a distinct universality class from both central–force rigidity and bond bending without free rotation at cross links rigidity . The variation of $`\nu `$ is more subtle if at all, varying from $`0.5`$ for high densities to $`0.35\pm 0.1`$ near the critical point (for comparison, standard stability considerations require $`1\nu 1`$ in 2D). Further details of the scaling behavior near the transition will be presented elsewhere long . We remark that at finite temperatures $`G`$ will remain non–zero above the conductivity percolation transition at $`L/l_\mathrm{c}5.42`$, but it is not clear if this small–$`G`$ behavior is experimentally observable for the macromolecules under consideration here.

Away from the critical point, $`G`$ decreases as the bending modulus $`\kappa `$ decreases. Indeed, in the limit $`\kappa 0`$, or equivalently $`l_\mathrm{b}0`$, all filaments freely bend and the model reduces to a random network of Hookean springs, which is already known to have a vanishing $`G`$ for all finite cross link densities zero\_G . In the opposite limit $`l_\mathrm{b}\mathrm{}`$, bending along a filament’s length becomes prohibitively expensive and response becomes dominated by stretching modes. The same is true of the limit $`L/l_\mathrm{c}\mathrm{}`$ with $`l_\mathrm{b}/L`$ fixed, since bending modes of wavelength greater than $`l_\mathrm{c}`$ require an area of network to twist rather than a single filament, and are suppressed. This transition from low $`L/l_\mathrm{c}`$, bending–dominated networks to high $`L/l_\mathrm{c}`$, stretching–dominated ones can be clearly seen Fig. 1.

Pure, affine shear, being a combination of rotation and extension, induces only stretching and compression of filaments. Furthermore, networks dominated by stretching modes must be approximately affine. This follows from the observation that, far above the rigidity transition, there can be no orientational disorder without a corresponding cost in bending energy. It is straightforward to derive analytic expressions for the modulii under an affine strain. A rod of length $`L`$ lying at an angle $`\theta `$ to the $`x`$–axis will undergo a relative change in length $`\delta L/L=\gamma _{\mathrm{xy}}\mathrm{sin}\theta \mathrm{cos}\theta `$ in response to an affine strain field $`\gamma _{\mathrm{xy}}`$ . According to (1) with $`\delta \theta =0`$, the energy cost (after uniformly averaging over all angles $`\theta `$) is $`\delta =\mu L\gamma _{\mathrm{xy}}^2/16`$. To calculate $`G`$, we need to express the number of rods per unit area $`N`$ as a function of $`L`$ and $`l_\mathrm{c}`$. The exact expression is easy to derive, but for current purposes it is sufficient to use the approximate relation $`L/l_\mathrm{c}(\alpha 1)/(12/\alpha )`$ with $`\alpha =2L^2N/\pi `$, valid for $`L/l_\mathrm{c}5`$ and greater (note that the full expression is monotonic). A further correction removes the dangling ends of the rods by renormalizing the rod lengths to $`L2l_\mathrm{c}`$. Then using $`G=2/\gamma _{\mathrm{xy}}^2`$,

$$G_{\mathrm{affine}}=\frac{\pi }{16}\frac{\mu }{L}\left(\frac{L}{l_\mathrm{c}}+2\frac{l_\mathrm{c}}{L}3\right)$$
(2)
so that $`G(\pi /16)\mu /l_\mathrm{c}`$ as $`L/l_\mathrm{c}\mathrm{}`$. Expression (2) is plotted in Fig. 2 and gives a reasonable approximation to the data, with the agreement improving as $`L/l_\mathrm{c}`$ increases. The same calculations can be repeated to give $`Y_{\mathrm{affine}}=3G_{\mathrm{affine}}`$ and hence a Poisson ratio $`\nu _{\mathrm{affine}}=0.5`$, which is also plotted in the figure. Perhaps surprisingly, we find remarkably good agreement between the measured $`\nu `$ and the affine value, even close to the rigidity transition, where the affine approximation fails.

We observe that there is an apparent crossover from a bending-dominated, non-affine regime for either high molecular weight ($`L`$) or for high density. The natural measure for this is the ratio of filament length to average distance between crosslinks. This can be understood by the fact that, unlike networks of flexible polymers, where segments along a single polymer between crosslinks appear to behave as effectively independent network strands, the segments of semiflexible filaments between crosslinks act in series. Thus, segments far from free ends are forced to deform nearly affinely by the many constraints on their neighboring segments imposed by crosslinks. This suggests a physical picture in which non-affine deformations are primarily associated with less constrained free ends. Let $`\lambda `$ denote the range of such non-affine regions near the filament ends. We expect this length $`\lambda (l_c,l_b)`$ to be a function of the local density of filaments (measured by $`l_c`$, the distance between crosslinks) and the material length $`l_b`$. Then, what determines the degree to which the network is affine or not is the relative size of the affine to non-affine regions, or the ratio $`L/\lambda `$. We find that, sufficiently far above the transition, $`G/G_{\mathrm{affine}}`$ collapses to a single master curve with the empirical choice $`\lambda =\sqrt[3]{l_c^4/l_b}`$, as shown in Fig. 3 (the origin of this length scale is discussed elsewhere long ). This confirms the existence of a fundamental length scale $`\lambda `$ for non-affine deformations along the filament backbone, as well as two physically distinct regimes: non-affine behavior for $`L\lambda `$ and affine behavior for $`L\lambda `$. In the first of these, we note that the modulus depends only on $`\kappa `$ and not on $`\mu `$, signalling a bending dominated regime, as predicted by Kroy and Frey frey ; frey2 . Similar results have been found independently frey3 . Note that the crossover can be reached at any fixed density by varying $`l_\mathrm{b}`$ alone, and thus represents distinct physics from the percolation transition, which exists only at a specific density. Hence we use absolute density, rather than its value relative to the transition, in our scaling function.

Our findings demonstrate that the degree of affinity of the deformation field does depend on the length scale on which one looks. This is shown in Fig. 4, where the quantity $`\mathrm{\Delta }\theta ^2(r)`$ is plotted against $`r`$ for different values of $`L/\lambda `$. $`\mathrm{\Delta }\theta ^2(r)`$ is the square deviation of the angle of rotation $`\theta `$ between two points separated by a distance $`r`$, relative to the affine equivalent. This monotonically decreases for increasing $`r`$, suggesting that the deformation appears more affine when viewed on larger length scales. Furthermore, the deviation from affinity at the mesh scale $`r=l_\mathrm{c}`$ is small for networks with $`GG_{\mathrm{affine}}`$, and large for those with $`GG_{\mathrm{affine}}`$, as seen by comparing Figs. 4 and 3. Thus we can reiterate the main findings of our work: (1) there are two qualitatively distinct regimes, one affine and the other non-affine; (2) the physics of the crossover between these is distinct from the rigidity transition; and (3) the crossover is governed by a new length scale $`\lambda `$, where, *e.g.*, affine behavior is seen for filament lengths a few times this length. Since this length is expected to be of order the distance between crosslinks, real networks can be in either regime depending on the length distribution.

AJL would like to acknowledge the hospitality of the Vrije Universiteit. DAH was partly funded by a European Community Marie Curie Fellowship. This work is supported in part by the National Science Foundation under Grant Nos. DMR98-70785 and PHY99-07949.