text:
# Isospectrality in Chaotic Billiards

\[

## Abstract

We consider a modification of isospectral cavities whereby the classical dynamics changes from pseudointegrable to chaotic. We construct an example where we can prove that isospectrality is retained. We then demonstrate this explicitly in microwave resonators.

\]

Recently it has been shown that it is possible to construct two drums which have different shapes but sound exactly the same . This answers the famous question asked by Kac  in 1966 “Can you hear the shape of a drum?” the answer being “no”. Gordon, et. al. constructed an example of a pair of two-dimensional domains which had different shapes but had identical eigenvalue spectra for the Laplace operator. Since then, a large number of such isospectral pairs have been obtained.

One common feature of all shapes constructed so far is that they are mostly polygonal. Hence the classical dynamics of a particle in billiards of these shapes is pseudo-integrable. A question of interest then is whether isospectrality can be achieved even for cavities with chaotic dynamics, which is typical of domains that have convex pieces, and hence are non-polygonal. We address this question, viz. are there sound-alike chaotic drums, both theoretically and through experiments using microwave resonators.

Isospectrality is fundamentally a consequence of topology. The essential aspects of isospectrality can be proved using the example of the two isospectral domains $`𝒞1`$ and $`𝒞2`$ shown in Fig. 1. The proof consists in showing that given any eigenfunction in one domain we can construct a corresponding one in the other domain, with the same eigenvalue, and vice-versa. Each domain consists of seven distinct sub-domains each in the form of a triangle. We label these sub-domains in an arbitrary fashion, using numbers $`1,2\mathrm{}7`$ for domain $`𝒞1`$ and the alphabets $`A,B,\mathrm{}G`$ for domain $`𝒞2`$. Note that edges of the triangles are marked differently (By dotted, dashed and solid lines) and this allows us to make a unique correspondence between any pairs of triangles. Consider any wavefunction, $`\psi `$, in domain $`𝒞1`$, which satisfies the eigenvalue equation $`^2\psi =k^2\psi `$ with Dirichlet boundary conditions ($`\psi `$ vanishes on the boundary of the domain). Let us denote by $`\psi _i`$ the restriction of the wavefunction $`\psi `$ in sub-domain $`i`$ (i.e. $`\psi _i(\overline{r})=\psi (\overline{r})`$ if $`\overline{r}`$ is a point in the $`ith`$ sub-domain, else $`\psi _i(\overline{r})=0`$). Similarly we can define the restricted wavefunctions $`\{\psi _A,\psi _B\mathrm{}\psi _G\}`$ from any wavefunction in domain $`𝒞2`$. Starting from the wavefunction $`\psi `$ in $`𝒞1`$ let us construct the following restricted wavefunctions in domain $`𝒞2`$:

$`\psi _A=\psi _2\stackrel{~}{\psi }_1+\psi _7`$ (1)
$`\psi _B=\psi _3+\psi _1+\psi _5`$ (2)
$`\psi _C=\stackrel{~}{\psi }_3+\psi _2+\psi _4`$ (3)
$`\psi _D=\psi _4\stackrel{~}{\psi }_1+\psi _6`$ (4)
$`\psi _E=\psi _5\stackrel{~}{\psi }_2\stackrel{~}{\psi }_6`$ (5)
$`\psi _F=\psi _7\stackrel{~}{\psi }_3+\psi _6`$ (6)
$`\psi _G=\stackrel{~}{\psi }_7+\psi _5\stackrel{~}{\psi }_4`$ (7)
The notation used requires some explanation: to construct $`\psi _A=\psi _2\stackrel{~}{\psi }_1+\psi _7`$, we first move the three domains $`1,2`$ and $`7`$ so that they are on top of each other and all similarly marked edges coincide. This may require us to flip domains about one of the bases and in such cases we have denoted the wavefunction with a tilde (e.g. $`\stackrel{~}{\psi }_1`$). The wavefunction $`\psi _A`$ is then obtained by adding (or subtracting) the values of the three functions at each point. It is easy to see that $`\psi ^{}=\psi _A+\psi _B+\psi _C+\psi _D+\psi _E+\psi _F+\psi _G`$ is an eigenfunction for domain $`𝒞2`$ with the same eigenvalue. For this we notice that:

(1) Laplace’s equation is satisfied in every domain and

(2) It can be verified that the wavefunction vanishes on the boundary and matches smoothly across sub-domains. For example consider the sub-domains $`A`$ and $`B`$ which are separated by a dashed line. The wavefunctions are given by $`\psi _2\stackrel{~}{\psi }_1+\psi _7`$ and $`\psi _B=\psi _3+\psi _1+\psi _5`$. The smoothness follows since from the wavefunction in $`𝒞1`$ we see that $`\psi _2`$ matches smoothly with $`\psi _3`$ across the dashed boundary, similarly $`\psi _7`$ matches $`\psi _5`$ and $`\stackrel{~}{\psi }_1`$ matches $`\psi _1`$.

Similarly one can construct an eigenfunction for $`𝒞1`$ starting from any given eigenfunction in $`𝒞2`$. Thus we have demonstrated a one-to-one correspondence between the states in the two cavities and hence proved isospectrality.

We now modify the domain geometry so as to make the dynamics chaotic. It is expected that making a part of the boundary convex (inwards into the domain) should make the dynamics chaotic. This is related to the fact that on such boundaries, any two particle trajectories which are close to each other, diverge rapidly after being reflected . A well known example of a chaotic billiard is the Sinai billiard obtained by placing a circular scatterer inside a square. In our case, to obtain the modified geometry we first place a scatterer of arbitrary shape inside one of the triangular sub-domains of any one domain and then place one in a similar position in every other triangle. An example with disc-shaped scatterers is shown in Fig. 2. The identification of edges on the two domains makes this construction unique. Thus notice that, in every triangle, the scatterer is placed close to a vertex where a solid and dotted edge meet. The wavefunction in each domain now changes since it has to vanish on and inside the boundary of the scatterers. Our construction of the modified geometry with scatterers is such that the proof for isospectrality given above can be repeated, since the wavefunctions still satisfy the relations given by eq.(1).

A direct physical proof of isospectrality can be obtained by experiments utilizing microwave cavities which provide a simple and powerful method of simulating single particle time-independent quantum mechanics in two dimensions. This follows from the fact that under appropriate geometrical constraints, Maxwell’s equations in a cavity reduces to the Schrödinger equation of a free particle inside a two dimensional domain of arbitrary shape and topology . Infact one can show that, for a cavity with thickness (in the $`z`$ direction, say) small compared to the dimensions in the other transverse directions (in the $`xy`$-plane), the $`z`$-component of the microwave electric field $`\mathrm{\Psi }(x,y)=E_z`$ satisfies the time-independent Schrödinger wave equation $`(_x^2+_y^2)\mathrm{\Psi }=k^2\mathrm{\Psi }`$ (with the identification $`k=2\pi f/c`$, $`f`$ being the frequency and $`c`$ the speed of light) and $`\psi `$ vanishes on the boundary of the domain. This correspondence is exact for all frequencies $`f<c/2d`$ where $`d`$ is the thickness of the cavity. We note that this is also the Helmholtz equation which describes, for example, vibrations of a drum. Using this equivalence, various phenomena (such as quantum chaos) have been studied. Isospectrality has earlier been demonstrated by Sridhar and Kudrolli using microwave cavities shaped as in Fig. 1. In the experiments one obtains the resonance modes of the cavities. Thus the microwave transmission spectra directly yields the eigenvalues of the cavity being measured. The advantage of this approach is that it can be easily applied to arbitrary 2-D domains, for which numerical simulations are very hard  and may sometimes be practically impossible.

In our experiments we consider the same set of cavities (Fig. 1) as the ones considered by Sridhar and Kudrolli  and investigate the question of isospectrality in the presence of scatterers placed in the specified way inside the cavities. A schematic of the experimental cavity is shown in Fig. 3. The desired domain is cut out from a brass plate of thickness $`d=6mm`$. Two other brass plates are placed on top and below the hole to form a closed cavity. As shown in the figure microwaves were coupled in and out using loops terminating coaxial lines that enter through the sides of the cavity. The length of bases of the triangular sub-domains was taken to be $`a=8cm`$ and the thickness of the cavity was $`d=6mm`$. The small thickness of the cavity makes it essentially $`2`$-dimensional and the correspondence between Maxwell’s equations and Schrödinger’s equation is good for frequencies $`f<c/2d=25GHz`$. For all metallic objects in the 2-D space between the plates, Dirichlet boundary conditions apply inside the metal.

All measurements were carried out using an HP8510B vector network analyzer which measured the transmission ($`S_{21}`$) parameters. The typical values of quality factor obtained range from a maximum of $`850`$, at the lower end of the spectrum, to a minimum of $`250`$.

Results: We first attempt to reproduce the results in  for the cavities shown in Fig. 4. We show in Fig 4 the traces of the spectrum for the two cavities in the frequency range $`15GHz`$. The first $`30`$ resonances of the two cavities are listed in Table I and one sees that the eigenvalues match to better than $`1`$%. One sees that each resonance present in one is present in the other. A few lines are missing and this is attributed to the fact that the particular coupling positions we used may not excite some modes. The remaining inaccuracies are due to imperfections in the machining, and in the clamping together of various parts of the cavity. Note, that the amplitudes themselves may be different, as that depends on the location of the coupling and hence to the way the modes are excited.   TABLE I.: The table lists the first $`33`$ resonances in the two cavities.     Resonant frequency Resonant frequency percentage  in $`𝒞1`$(in MHz) in $`𝒞2`$(in MHz) discrepancy  1902.500 1903.750 0.0657  2271.250 2274.375 0.1374  2700.625 2719.375 0.6895  3045.625 3062.500 0.5510  3217.500 3215.000 -0.0778  3612.500 3631.250 0.5163  —- 3892.500 —-  4054.375 —- —-  4184.375 4200.625 0.3868  4303.125 4328.125 0.5776  4488.125 4528.125 0.8834  4743.750 4756.250 0.2628  4898.750 —- —-  5026.875 5031.875 0.0994  5171.875 5194.325 0.4322  5426.250 5474.325 0.8782  5488.750 —- —-  5625.625 5627.500 0.0333  5793.750 5808.750 0.2582  —- 5903.125 —-  5928.750 5940.625 0.1999  6085 —- —-  6222.500 6242.500 0.3204  6312.500 6352.500 0.6297  6497.500 6492.500 -0.0770  6680.000 6707.500 0.4100  6750.000 6775.000 0.3690  —- 6790.000 —-  6855.000 6877.500 0.3272  6930.000 6992.500 0.8938    Thus we have obtained the energy spectrum for the given set of isospectral cavities and verified that each eigenvalue in one is present in the other at the same resonance value.

As a check on the quality of the spectral data we compare the cumulative number of resonance levels as a function of frequency, obtained experimentally, with the Weyl formula for the integrated density of states in a two-dimensional domain :

$`N(k)={\displaystyle \frac{Ak^2}{4\pi }}{\displaystyle \frac{Sk}{4\pi }}+K,`$ (8)
where $`A`$ and $`S`$ are the area and perimeter of the domain, and $`K`$ is a correction term associated with its topology. For a polygonal billiard with inner angles $`\alpha _i`$ this is given by $`K=_i\frac{1}{24}(\frac{\pi }{\alpha _i}\frac{\alpha _i}{\pi })`$. In the present case we find $`K=0.42`$. We show in Table II a comparison between the experimental results with the above formula. The agreement is quite good.   TABLE II.: The table gives a comparison of the cumulative frequency in the cavities with the Weyl estimate     frequency cumulative Weyl’s estimate  (in GHz) resonant frequency   1 0 0  2 1 0.8  3 3 3.4  4 7 7.5  5 13 13.1  6 21 20.4  7 30 29.2  8 38 39.5    Results for the chaotic geometry: We place the scatterers inside the cavity following the prescription outlined above. The scatterers are taken to be metallic cylinders of diameter $`1.0cm`$ and height equal to the thickness of the cavities. The modified spectrum from the two cavities is shown in Fig. 5 and we list the first $`22`$ resonances in Table III. We again find that the eigenvalues in the two cavities match to within $`1`$%. Thus there is clear evidence that isospectrality is retained in the modified chaotic geometry.

It may be noted that the introduction of the scatterers changes the topology of the domain from being simply-connected to now being multiply-connected. We now have a polygonal box with $`p=7`$ circular holes. In this case the topology term in Weyl’s formula for the integrated density of states is given by: $`K=_i\frac{1}{24}(\frac{\pi }{\alpha _i}\frac{\alpha _i}{\pi })\frac{p}{6}`$. The comparison with Weyl’s formula is given in Table IV. The agreement is not very good at low frequencies and the number of levels seems to be somewhat higher than that given by the Weyl estimate.   TABLE III.: The first 22 resonant frequencies in the cavities with scatterers    Resonant frequency Resonant frequency percentage  in $`𝒞1`$(in MHz) in $`𝒞2`$(in MHz) discrepancy    2181.875 2175.000 -0.3161  2330.000 2338.875 0.3795  2906.250 2933.125 0.9163  3303.750 3311.250 0.2265  3346.875 3361.875 0.4462  3746.875 3766.250 0.5144  4135.625 —- —-  4175.625 4175.625 0.000  4207.500 4222.500 0.3552  4645.000 4666.250 0.4554  4966.875 4956.250 -0.2144  5056.250 —- —-  5106.250 5101.250 -0.0980  5148.750 5124.375 -0.4757  5540.000 5576.875 0.6612  5603.125 5633.125 0.5326  —- 5889.375 —-  6045.000 6002.500 -0.7080  6435.000 6450.000 0.2325  6490.000 6505.000 0.2306  6850.000 6872.500 0.3274  6955.000 6995.000 0.5718      TABLE IV.:  Comparison of Weyl’s estimate and cumulative frequency obtained experimentally after placing the scatterers     frequency cumulative Weyl’s estimate  (in GHz) resonant frequency   1 0 0  2 0 0  3 3 0.9  4 6 4.5  5 11 9.7  6 17 16.3  7 22 24.5    To make the demonstration more convincing and illustrate the non-triviality of the isospectral construction with scatterers, we consider another geometry (Fig. 6) where the scatterers in the second cavity are placed in a somewhat different manner. The arrangement still seems to follow the folding construction and naively one would expect isospectrality. However on closer inspection one finds that the correct correspondence between the edges of the sub-domains has not been satisfied and the wavefunction matching condition infact no longer holds and so we should not get isospectrality. We plot the spectrum for this case in Fig. 7. We see a marked difference from that in Fig. 5, namely we find that there is no correspondence between the spectral lines from the two cavities. This shows clearly that isospectrality is indeed obtained only for the special arrangement of scatterers in Fig. 2.

In conclusion we have demonstrated that isospectrality is unrelated to the underlying classical dynamics of a particle. We have shown a simple way of introducing scatterers of arbitrary shape into polygonal cavities in such a way that isospectrality is retained. This leads us to a new class of isospectral scatterers and also a better understanding of the essential features necessary for isospectrality.

We thank N. Kumar, Yashodhan Hatwalne and Joseph Samuel for discussions. SS thanks RRI for hospitality while this work was completed. AD acknowledges support from the NSF under grant DMR 0086287. SS was partially supported by NSF 0098801.


format: text

text:
# Solution of the Linear Ordering Problem (NP=P)

We consider the following problem

$`\mathrm{max}`$ $`{\displaystyle \underset{i=1,ij}{\overset{n}{}}}{\displaystyle \underset{j=1}{\overset{n}{}}}c_{ij}x_{ij}`$
s. t. $`0x_{ij}1,`$
$`x_{ij}+x_{ji}=1,`$
$`0x_{ij}+x_{jk}x_{ik}1,ij,ik,jk,i,j,k=1,\mathrm{},n.`$
We denote the corresponding polytope by $`B_n`$. The polytope $`B_n`$ has integer vertices corresponding to feasible solutions of the linear ordering problem as well as non-integer vertices. We denote the polytope of integer vertices as $`P_n`$.

Let us give an example of non-integer vertex in $`B_n`$ and describe an exact facet cut. In what follows we will interested only in generating exact facet cuts.

Fig. 1 shows a graph interpretation of a non-integer vertex ,

$$\begin{array}{c}\text{}\\ \\ \text{Fig. 1}\end{array}$$

where $`\{i_1,\mathrm{},i_m\}`$, $`\{j_1,\mathrm{},j_m\}`$ $`\{1,\mathrm{},n\}`$; $`\{i_1,\mathrm{},i_m\}\{j_1,\mathrm{},j_m\}=\mathrm{}`$; $`x_{i_lj_q}=0`$, $`x_{j_qi_l}=1`$, $`lq`$, $`l,q=1,\mathrm{},m;`$ the other variables that are not shown at the Figure 1 are equal to $`\frac{1}{2}.`$ This is the simplest non-integer vertex of the polytope $`B_n`$. For this vertex all adjacent integer vertices can be written as:

$$\begin{array}{ccc}\alpha _ki_kj_k\beta _k,\hfill & \text{where}\hfill & \alpha _k\text{ is any ordering from the set }\{j_1,\mathrm{},j_m\}\backslash \{j_k\},\hfill \\ & & \beta _k\text{ is any ordering from the set }\{i_1,\mathrm{},i_m\}\backslash \{i_k\},\hfill \\ & & \\ \alpha _{kp}i_kj_ki_pj_p\beta _{kp},\hfill & \text{where}\hfill & \alpha _{kp}\text{ is any ordering from the set }\{j_1,\mathrm{},j_m\}\backslash \{j_k,j_p\},\hfill \\ & & \beta _{kp}\text{ is any ordering from the set }\{i_1,\mathrm{},i_m\}\backslash \{i_k,i_p\},\hfill \end{array}$$

$$kp,k,p=1,\mathrm{},m.$$

All adjacent integer vertices, which number is equal to

$`m\left[(m1)!\right]^2+{\displaystyle \frac{m(m1)}{2}}\left[(m2)!\right]^2`$
lie in one hyperplane

$`f(x)=2{\displaystyle \underset{l=1}{\overset{m}{}}}x_{i_lj_l}{\displaystyle \underset{l=1}{\overset{m}{}}}{\displaystyle \underset{q=1}{\overset{m}{}}}x_{i_lj_q}=1.`$
This hyperplane for $`f(x)1`$ is a facet of the polytope $`P_n.`$

Our aim is to determine exact facet cuts for any non-integer vertex of $`B_n`$ (and not only for them) in an analogous fashion.

Figures 2 and 3 show non-integer vertices of the polytope $`B_n`$:

$$\begin{array}{cc}\text{}& \text{}\\ & \\ \text{Fig. 2}& \text{Fig. 3}\end{array}$$

Noninteger vertex at Figure 2 has an oriented chain $`7583`$ of the length 3, and non-integer vertex at Figure 3 has an oriented $`614`$ chain of the length 2. The oriented chain $`758`$ at Figure 2 is independent, that is if we exchange the chain $`758`$ with any other chain the rest of the graph does not change, while the chains $`81`$ at Figure 2 and $`614`$ at Figure 3 are dependent. We define corresponding dependent and independent oriented chains.

The following Theorems take place.

###### Theorem 1

Let $`x^0`$ be a noninteger vertex in $`B_n`$ and assume that in graph interpretation there is a graph vertex $`i`$ which is the begin or the end of all adjacent arcs. Assume that the vertex $`i`$ can be repeated arbitrarily many times such that each of the new vertices has the same location with the other part of the graph as the vertex $`i`$. Then the new noninteger vertex, corresponding to the new graph, is a noninteger vertex of $`B_n`$, and in the new graph the vertices $`i`$ and new inserted vertices may be put in any order.

###### Theorem 2

Let $`x^0`$ be a noninteger vertex in $`B_n`$. Then there does not exist corresponding dependent oriented chains of the length 4.

The polytope $`B_n`$ has noninteger vertices whose fractional components are equal to $`\frac{l}{r}`$, $`r2,`$ $`l<r`$, as well.

For $`r=2`$ after matrix transformation we get in all cases the following non unimodular minimal standard matrix:

$`\left(\begin{array}{ccc}\hfill 1& \hfill 1& \hfill 0\\ \hfill 1& \hfill 0& \hfill 1\\ \hfill 0& \hfill 1& \hfill 1\end{array}\right).`$ (4)
For $`r=3`$ after matrix transformation we get in all cases the known combination of two minimal standard matrices:

$`\left(\begin{array}{ccccc}\hfill 1& \hfill 1& \hfill 0& \hfill 1& \hfill 0\\ \hfill 1& \hfill 0& \hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 1& \hfill 0& \hfill 0\\ \hfill 1& \hfill 0& \hfill 0& \hfill 0& \hfill 1\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1& \hfill 1\end{array}\right).`$ (10)
For any $`r`$ after matrix transformation we get known combination of $`r1`$ minimal standard matrices.

###### Theorem 3

Let $`x^0`$ be a noninteger vertex of $`B_n`$, which has fractional components $`\frac{l_1}{r_1}`$, $`l_1<r_1`$, $`r_13`$. Then we pass to an adjacent noninteger vertex with fractional components $`\frac{l_2}{r_2},`$ $`r_2<r_1`$, $`l_2<r_2,`$ by changing an equality in a basis (thus changing one or more minimal standard matrices).

Let $`I_1`$, $`I_2`$, …, $`I_s`$ $`\{1,\mathrm{},n\}`$, and assume that each set $`I_p`$, $`p=1,\mathrm{},s,`$ corresponds to noninteger components of a vertex. Then for each $`I_p`$, $`p=1,\mathrm{},s,`$ we can construct a facet cut. If $`s=1`$ a noninteger vertex is called simple. A noninteger vertex is called complex if $`s2.`$

Thus we have given a general description of the polytope $`B_n.`$

###### Theorem 4

Let $`x^0`$ be a simple noninteger vertex of the polytope $`B_n`$ with fractional components $`\frac{1}{2}`$, assume further that there does not exist dependent oriented chains with the length 3. Then all adjacent integer vertices lie in one hyperplane, this hyperplane is a facet of the polytope $`P_n`$, and it can be constructed by a polynomial algorithm.

Now we describe the principle for constructing facets.

Consider a noninteger vertex $`x^0.`$ It can be defined as the solution of the following system of the linear basic equalities

$`x_{i_lj_l}=0,l=1,\mathrm{},q;`$ (11)
$`x_{i_lj_l}+x_{j_lk_l}x_{i_lk_l}=0,l=q+1,\mathrm{},{\displaystyle \frac{n^2n}{2}}.`$
We introduce artificial variables $`x_{j_ln+1}=0,`$ $`x_{i_ln+1}=0,`$ into the first $`q`$ equalities of the system (11):

$`x_{i_lj_l}+x_{j_ln+1}x_{i_ln+1}=0,l=1,\mathrm{},q.`$
With the help of the notation

$`x_{i_lj_l}+x_{j_lk_l}x_{i_lk_l}:=x(i_l,j_l,k_l),l=1,\mathrm{},{\displaystyle \frac{n^2n}{2}},`$
we rewrite the system (11) as follows:

$`x(i_l,j_l,k_l)=0,l=1,\mathrm{},{\displaystyle \frac{n^2n}{2}}.`$
We can determine $`\frac{n^2n}{2}`$ linear independent adjacent integer vertices

$`x^q(i_s,j_s,k_s)=\delta ^q(i_s,j_s,k_s),s,q=1,\mathrm{},{\displaystyle \frac{n^2n}{2}},`$
where $`\delta ^q`$ is either 1 or 0. We can prove that all adjacent integer vertices lie in the hyperplane:

$`f(x)=\left|\begin{array}{cccc}\hfill x(i_1,j_1,k_1)& \hfill \mathrm{}& \hfill x(i_m,j_m,k_m)& \hfill 1\\ \hfill \delta ^1(i_1,j_1,k_1)& \hfill \mathrm{}& \hfill \delta ^1(i_m,j_m,k_m)& \hfill 1\\ & \hfill \mathrm{}& & \\ \hfill \delta ^m(i_1,j_1,k_1)& \hfill \mathrm{}& \hfill \delta ^m(i_m,j_m,k_m)& \hfill 1\end{array}\right|=0`$ (16)
where $`m=\frac{n^2n}{2}.`$

###### Theorem 5

Let $`x^0`$ be a simple noninteger vertex of the polytope $`B_n`$ with fractional components $`\frac{1}{2}`$, assume further that there exist $`\tau `$ dependent oriented chains with the length 3. Then all adjacent integer vertices lie in $`2^\tau `$ hyperplanes, each of them is a facet of the polytope $`P_n`$, and can be constructed by a polynomial algorithm.

###### Theorem 6

Let $`x^0`$ be a simple noninteger vertex of the polytope $`B_n`$ with fractional components $`\frac{l}{r},`$ $`r3`$, $`l<r.`$ Then we can construct all minimal standard matrices and corresponding noninteger vertices with fractional components $`\frac{1}{2}.`$ For every such vertex we can construct facet cuts.

###### Theorem 7

Let $`x^0`$ be a complex noninteger vertex of the polytope $`B_n,`$ and $`I_1`$, $`I_2,`$ …, $`I_s`$ correspond to noninteger values. Then we can construct facet cuts for each set $`I_p`$, $`p=1,\mathrm{},s.`$

Assume we have generated facet cuts. Solving the problem again we get the noninteger vertex $`x^1`$ of the polytope

$`0x_{ij}1,`$
$`x_{ij}+x_{ji}=1,`$
$`0x_{ij}+x_{jk}x_{ik}1,ij,ik,jk,i,j,k=1,\mathrm{},n,`$
$`f_s^1f_s(x)f_s^2,s=1,\mathrm{},q.`$
Without loss of generality we may assume that the noninteger vertex $`x^1`$ satisfies the following linear system:

$`x(i_s,j_s,k_s)=1,i=1,\mathrm{},p,`$
$`f_s(x)=f_s^2,s=1,\mathrm{},q.`$
Now we find all adjacent integer vertices. If all of them lie in one hyperplane and this hyperplane is a facet of $`P_n`$ then we generate this facet and re-solve the problem with a new facet. In case we cannot determine the facet we solve the auxiliary problem:

$$\begin{array}{cc}\hfill \mathrm{max}& \left(_{s=1}^px(i_s,j_s,k_s)+_{s=1}^q\frac{f_s(x)f_s^1}{f_s^2f_s^1}\right),\hfill \\ & 0x_{ij}1,\hfill \\ & x_{ij}+x_{ji}=1,\hfill \\ & 0x_{ij}+x_{jk}x_{ik}1,ij,ik,jk,i,j,k=1,\mathrm{},n,\hfill \end{array}$$

With the solution of this problem we can construct the facet of the polytope $`P_n.`$ In the case of theorems 5 and 6 we can construct the necessary facets by means of a polynomial algorithm.


text:
# Non-existence of normal tokamak equilibria with negative central current

## I Introduction

Tokamaks with reversed central magnetic shear (and thus low core current density) are of interest for at least two reasons: 1) internal transport barriers associated with reduced turbulence are often observed in them, leading to improved energy and particle confinement; and 2) they are the natural result of high beta operation and high bootstrap current fraction used to reduce non-inductive current drive requirements for steady state operation. Both of these features could make reversed magnetic shear operation attractive for a tokamak reactor.

Recent experiments on the Joint European Torus (JET)Hawkes et al. (2001); Stratton et al. (2002a); Hawkes et al. (2002); Stratton et al. (2002b) and the Japan Atomic Energy Research Institute Tokamak-60 Upgrade (JT-60U)Fujita et al. (2001); Miura et al. (2002); Chankin et al. (2003) have pushed the core current density to very low values using off-axis, non-inductive current drive. Large central current holes (regions of nearly zero current density) are produced because off-axis, non-inductive current drive in the same direction as the Ohmic current induces a back electromotive force inside the non-inductive current drive radius that decreases the core current density.

An interesting feature of current hole discharges is that the core current density is approximately zero (within Motional Stark Effect diagnostic measurement errors), even though there is often sufficient current drive power that the core current could in principal go significantly negativeStratton et al. (2002a); Hawkes et al. (2002) (negative relative to the direction of the total plasma current). Recent non-linear toroidal resistive MHD (Magnetohydrodynamic) simulationsStratton et al. (2002a); Breslau et al. (2003) predict that current hole discharges undergo rapid $`n=0`$ reconnection events (axisymmetric sawteeth) that clamp the core current near zero. More generally, this reconnection occurs whenever the current density profile is such that the rotational transform, $`\iota `$, goes to zero on any surface in the plasma (this includes the case where the current density on-axis is positive, but the current profile goes sufficiently negative somewhere off-axis that the total current enclosed by some flux surface vanishes). Reduced MHD simulations in cylindrical geometry have also shown that $`n=0`$ resistive kink instabilities can clamp the core current density at zero when it attempts to go negative.Huysmans et al. (2001) Breslau et al.Breslau et al. (2003) and Stratton et al.Stratton et al. (2002a) stated that a second-order, large aspect ratio expansion of the MHD equations indicates that a normal toroidal equilibrium is not possible if $`\iota `$ crosses through $`0`$ at some radius. (They also stated that a more general proof is needed, which we provide here.) A recent paper by Chu and ParksChu and Parks (2002) used a second order aspect ratio expansion to prove that a normal equilibrium with a peaked pressure profile is not possible with negative core current. They extended the analysis to provide matching conditions at the boundary of a central region with no current and no pressure gradient, showing explicitly that current hole equilibria are theoretically possible (with zero, but not negative, current).

This paper extends some of these results to arbitrary aspect ratio, employing a relatively simple constraint based on a version of the virial theorem to show that a “normal” toroidal MHD equilibrium (with axisymmetric nested flux surfaces around a single magnetic axis, non-singular continuous fields, and a monotonic peaked pressure profile) is not possible with negative core current. Or more generally, a normal equilibrium is not possible if the toroidal current enclosed by any flux surface goes negative relative to the direction of the total plasma current, so that there is an $`\iota =0`$ surface somewhere in the plasma where the poloidal field vanishes (the null surface). Though the starting point of this analysis is based on well-known equations, they are often specialized to large aspect ratio or simplified geometry, while the present analysis is more general.

However, the virial constraint does not necessarily eliminate the possibility of more exotic equilibria, such as with non-nested flux surfaces with islands, or with singular fields and/or off-axis peaks in the pressure profiles that may be spontaneously generated by the plasma near the null surface. Some examples are considered here. In this paper we investigate the consequences of only one integral constraint on equilibria, while there can be other constraints that further limit the types of theoretically possible or experimentally realizable equilibriaChu and Parks (2002).

The non-existence of normal equilibria with negative core current, and/or the rapid axisymmetric sawteeth that are predicted to occur if the enclosed current goes negative, may also explain the results of other experiments, such as the low efficiency seen in some electron cyclotron counter current drive experimentsWesterhof et al. (2001); Miura et al. (2002).

## II Derivation

The MHD equilibrium equation $`p=𝒋\times 𝑩/c`$ can be written as

$$0=\left(p+\frac{B^2}{8\pi }\right)+\frac{1}{4\pi }\left(𝑩\right)𝑩$$
(1)
One common use of the virial theorem is to take the inner product of this equation with the position vector $`𝒙`$ and integrate over all space to show that an isolated MHD equilibrium can not exist by itself (unless there are physical coils or gravity to provide overall force balance).Shafranov (1966); Freidberg (1987) Here we use a version of the virial theorem that can be used to derive the Shafranov shiftRutherford (1981), by focusing on radial force balance of axisymmetric equilibria in cylindrical coordinates $`(R,Z,\varphi )`$. Taking the inner product of Eq. (1) with $`𝑹=R\widehat{𝑹}`$, the radial vector in cylindrical geometry, and integrating over space out to some flux surface of volume $`V`$, gives

$$0=𝑑V𝑹\left(p+\frac{B^2}{8\pi }\right)+\frac{1}{4\pi }𝑑V𝑹\left(𝑩\right)𝑩$$
(2)
For the second integral we use the identity $`𝑹(𝑩)𝑩=(𝑩𝑹𝑩)B_R^2B_\varphi ^2`$. The integral of $`(𝑩𝑹𝑩)`$ vanishes because $`𝑩d𝑺=0`$ on a flux surface. The first integral can be integrated by parts using $`𝑹=2`$, so that Eq. (2) becomes

$`0`$ $`=`$ $`p(\rho ){\displaystyle 𝑑𝑺𝑹}{\displaystyle 𝑑𝑺𝑹\frac{B^2}{8\pi }}`$
$`+2{\displaystyle 𝑑V\left(p+\frac{B^2}{8\pi }\right)}{\displaystyle \frac{1}{4\pi }}{\displaystyle 𝑑V\left(B_R^2+B_\varphi ^2\right)}`$
$`=`$ $`p{\displaystyle 𝑑𝑺𝑹}{\displaystyle 𝑑𝑺𝑹\frac{B^2}{8\pi }}+2{\displaystyle 𝑑V\left(p+\frac{B_Z^2}{8\pi }\right)},`$
where $`p=p(\rho )`$ is the pressure at the surface labeled by $`\rho `$ enclosing the volume $`V(\rho )`$, and $`B_Z`$ is the vertical magnetic field. For the first surface integral we can use Gauss’ theorem to write $`𝑑𝑺𝑹=𝑑V𝑹=2V`$. For the second surface integral, we use $`d𝑺=2\pi R\widehat{\mathit{\varphi }}\times d\stackrel{}{\mathrm{}}`$, where $`d\stackrel{}{\mathrm{}}`$ is a poloidal path length element along the surface, to write $`d𝑺𝑹=2\pi R^2\widehat{𝒁}d\stackrel{}{\mathrm{}}`$ so $`𝑑𝑺𝑹=2\pi R^2\widehat{𝒁}𝑑\stackrel{}{\mathrm{}}=2\pi R^2𝑑Z`$. Since the toroidal field $`B_\varphi 1/R`$ in a flux surface, the $`B_\varphi ^2`$ contribution to this surface integral vanishes, and we have

$$0=p(\rho )\frac{1}{V}𝑑\stackrel{}{\mathrm{}}\widehat{𝒁}\pi R^2\frac{B_{pol}^2}{8\pi }+p+\frac{B_Z^2}{8\pi }$$
(4)
where $`\mathrm{}=𝑑V\mathrm{}/𝑑V`$ denotes a volume average, and $`B_{pol}^2=B_R^2+B_Z^2`$ is the poloidal field strength squared. The poloidal field can be written as $`𝑩_{pol}=\varphi \times \psi =(\widehat{\mathit{\varphi }}\times \rho )(\psi /\rho )/R`$, where $`\rho `$ is a flux surface label. \[While $`\psi `$ is also constant on a flux surface, there can be two surfaces with the same value of $`\psi `$ in the presence of negative central current, so it is convenient to choose another flux surface label $`\rho `$, such as based on the enclosed volume or toroidal flux, to maintain monotonic labeling.\] If the toroidal current near the magnetic axis is in the opposite direction as the total plasma current, then the poloidal field must reverse direction somewhere and there must be a null surface on which the poloidal magnetic field is everywhere zero, as shown in Fig. 1. Another way to see this is to note that the poloidal field is related to the enclosed toroidal current by $`4\pi I_\varphi (\rho )/c=𝑑\stackrel{}{\mathrm{}}𝑩_{pol}=(\psi /\rho )𝑑\stackrel{}{\mathrm{}}\widehat{\mathit{\varphi }}\times \rho /R`$, so $`B_{pol}\psi /\rho =0`$ on any flux surface that encloses zero toroidal current. This is also the flux surface on which the rotational transform $`\iota =0`$ (corresponding to the safety factor $`q=\mathrm{}`$). \[These arguments assume that $`B_{pol}`$ is continuous and finite, we consider a singular exception in the next section.\]

On such a flux surface where the poloidal field vanishes, the second term of Eq. (4) vanishes and we are left with constraint

$$0=p+p+\frac{B_Z^2}{8\pi }$$
(5)
Since the last two terms are positive definite, the only way this equation can possibly be satisfied is if the pressure at this flux surface, $`p`$, is larger than the volume-averaged pressure inside that flux surface, $`p`$. I.e., the pressure profile must be hollow at least in some region, and can not be a monotonically decreasing function of $`\rho `$ at all radii as usual pressure profiles do.

This is in agreement with the result of Chu and ParksChu and Parks (2002), who also found that a normal equilibrium with a negative central current is not possible if the pressure profile is monotonically decreasing. These earlier results used a second-order large aspect ratio ordering while our derivation is valid for arbitrary aspect ratio. In other ways, their calculation goes beyond ours, as they have investigated additional constraints that can further limit the class of accessible equilibria.

## III Possible alternate solutions

Here we consider several possible alternate solutions for satisfying force balance in equilibria. Each differs from “normal” equilibria in a different way.

The argument in the previous section applies rigorously only for nested flux-surfaces where the fields are continuous and finite. Within the framework of ideal MHD, in principle there could be a singular poloidal field on a flux surface such that $`𝑑\stackrel{}{\mathrm{}}𝑩_{pol}=0`$ so that this surface encloses zero toroidal current, but $`𝑑\stackrel{}{\mathrm{}}\widehat{𝒁}\pi R^2B_{pol}^2/(8\pi )`$ is still finite and can contribute to Eq. (4) so that integrated force balance can be satisfied. An example of such a field might be the limiting case $`B_{pol}|\psi |C\mathrm{exp}(\mathrm{}^2/w^2)/\sqrt{w}`$. Then in the limit as $`w0`$ we have zero toroidal current enclosed while still giving a finite contribution to the second term of Eq. (4) (due to a singularity in the field chosen to be at $`\mathrm{}=0`$ in this example). Although $`B_{pol}`$ is becoming infinite at some point on the flux surface, it is an integrable singularity containing a finite amount of energy, and so could formally be considered as an admissible solution of ideal MHD. Since the spacing between two nearby flux surfaces labeled by poloidal flux $`\psi _1`$ and $`\psi _2`$ is given by $`\mathrm{\Delta }(\psi _2\psi _1)/|\psi |`$ (except where second derivatives have to be considered), there will be flux surfaces with finite separation at some places which will approach one another at the singular point where $`|\psi |B_{pol}\mathrm{}`$. The topology of this configuration is illustrated in Fig. 2. \[This sketch is intended only to illustrate the topology of a possible solution which satisfies the integral force balance constraint, Eq. (4). An actual detailed solution that would satisfy force balance locally at all points is left for future work, and may involve highly distorted flux-surface shapes with boundary layersPark et al. (1980); Rosenbluth et al. (1973); Cowley et al. (1991).\] Note that not only is the poloidal field infinite at the singular point, it must also flip signs from $`+\mathrm{}`$ to $`\mathrm{}`$ in the limit as the singular point is approached radially from opposite directions. Of course this singular configuration is strongly susceptible to magnetic reconnection when resistivity is included, consistent with the interpretation of rapid reconnection observed in simulationsStratton et al. (2002a); Breslau et al. (2003). If finite resistivity is included, then the singularity in the field is smoothed out and loss of equilibrium can help drive reconnecting flows.

Another possible way of satisfying integrated force balance with negative central current, while keeping the field finite and continuous, is if the flux surfaces are non-nested. An example is shown in Fig. 3, which is similar to the intermediate configurations observed during some toroidal simulations of axisymmetric reconnection in negative central current plasmas (for example, Fig. 11 of Ref. Stratton et al., 2002a, though in other cases they see islands with higher poloidal mode numbers). This case has an axisymmetric islandand the poloidal field is non-zero almost everywhere (except at the two magnetic axes and at the X-point) and so in principle can be arranged to give a negative contribution to the second term of Eq. (4) to satisfy integrated force balance. This is related to the role in a normal equilibrium of the Shafranov shift, which provides a larger value of $`R^2B_{pol}^2`$ on the outer part of a flux surface than on the inner part, so that the second term in Eq. (4) is negative. The X-point of a non-nested equilibrium might not be on the low-field side, and another possible equilibrium might be obtained by rotating Fig. 3 by 180 and shifting the spacing between flux surfaces so that the integral of $`R^2B_{pol}^2`$ on the outer part of the flux surface is again larger than on the inner part. One or the other of these configurations may be an unstable equilibrium and prefer to flip to the other orientation. \[TakizukaTakizuka (2002) earlier proposed another possible non-nested equilibria with negative central current, involving $`(m=2,n=0)`$ islands, while the example we discussed here has an $`(m=1,n=0)`$ island.\]

Another way of thinking about ideal MHD equilibria is to modify the time-dependent ideal MHD equations to include viscosity and parallel thermal conduction while retaining the ideal Ohm’s lawFreidberg (1987). Since viscosity should eventually damp the velocity $`𝒖`$ to zero, and parallel thermal conduction will lead to $`𝑩p=0`$, the dissipative terms vanish in a stationary steady state and the solution is also an ideal MHD equilibrium. One can then start with any arbitrary initial configuration of the magnetic field (which can be assumed to be nested flux surfaces) with arbitrary initial profiles, as functions of toroidal flux $`\mathrm{\Phi }`$, for the rotational transform $`\iota (\mathrm{\Phi })`$ and the adiabatic parameter $`\mu (\mathrm{\Phi })=p/n^\mathrm{\Gamma }`$ (where $`p(\mathrm{\Phi })`$ is the pressure profile, $`n(\mathrm{\Phi })`$ is the density profile, and $`\mathrm{\Gamma }`$ is the ratio of specific heats). Since this initial condition is not necessarily an equilibrium, flows will be driven and the plasma will move about, perhaps oscillating for a while. But it seems reasonable to assume that the viscosity will eventually damp out the oscillations and the plasma will settle into an equilibrium configuration while conserving $`\iota (\mathrm{\Phi })`$ and $`\mu (\mathrm{\Phi })`$. (The motions are assumed to be constrained to be axisymmetric to find such an equilibria. This approach to equilibria of course does not address the issue of stability, and these equilibria might then be unstable to symmetry-breaking perturbations.)

This was the logic that motivated the flux-conserving tokamak equilibria conceptsFreidberg (1987) that showed that there is formally no equilibrium limit on the pressure in a tokamak, since as the pressure increases, the Shafranov shift and the current can also increase to provide sufficient poloidal magnetic field on the outboard side to provide force balance.

Presumably this procedure would also find an equilibrium even if the rotational transform changed sign so that there was a null flux surface where $`\iota =0`$. In some cases with certain initial conditions, it might be possible for the plasma to spontaneously adjust flux surfaces near the null $`\iota `$ surface to produce a local peak in the pressure profile that can satisfy Eq. (5). \[We have focussed on the consequences of only one integral constraint that rules out “normal” equilibria with negative central current, and there can be other constraints that would further limit the practical accessibility of such non-monotonic pressure equilibriaChu and Parks (2002).\] The more typical case is probably that the equilibrium that is approached will have a singular structure, such as in Fig. 2, in order to satisfy Eq. (4). \[This is similar to studies showing that the nonlinear saturation of an internal kink mode approaches a neighbouring equilibrium state with singular currentsPark et al. (1980); Rosenbluth et al. (1973).\] These singular or near-singular states will be subject to strong reconnection if a small amount of resistivity is introduced, and the change in topology may dominate what happens. Realistic evaluations of what happens may depend on fully including various dissipation mechanisms (thermal, momentum, and particle anisotropic transport driven by small scale turbulence and collisional effects, as well as resistivity, current drive, heating and loss mechanisms). We leave detailed investigation of these issues to other work.

## IV Relation to the Shafranov shift

For completeness, we show the relation of Eq. (4) to usual expressions for the Shafranov shiftShafranov (1966); Freidberg (1987). The second term of Eq. (4) can be written as

$$T_2=\frac{\pi }{V}𝑑\stackrel{}{\mathrm{}}\widehat{𝒁}R^2\frac{B_{pol}^2}{8\pi }=\frac{\pi }{V}\frac{(\psi /\rho )^2}{8\pi }𝑑\stackrel{}{\mathrm{}}\widehat{𝒁}|\rho |^2$$
(6)
At this point, many previous calculations specialize to a large aspect ratio expansion and/or to specified shapes for the flux surfaces. For example, assume shifted circular flux surfaces with $`R(\rho ,\theta )=R_0\mathrm{\Delta }(\rho )+\rho \mathrm{cos}\theta `$ and $`Z(\rho ,\theta )=\rho \mathrm{sin}\theta `$, where $`\rho `$ has now been chosen to be the minor radius of the flux surface, and $`\mathrm{\Delta }`$ is the Shafranov shift. It can be shown that $`|\rho |^2=1/(1\mathrm{\Delta }^{}\mathrm{cos}\theta )^2`$, where $`\mathrm{\Delta }^{}=d\mathrm{\Delta }/d\rho `$. Defining $`\psi /\rho =B_p(\rho )(R_0\mathrm{\Delta })`$, we have $`B_{pol}=B_p[R_0\mathrm{\Delta }]/[R(1\mathrm{\Delta }^{}\mathrm{cos}\theta )]`$ (this would be exact if the flux surfaces really were shifted circles), and $`T_2`$ becomes

$$T_2=\frac{B_p^2}{8\pi }\frac{R_0\mathrm{\Delta }}{r}\frac{\mathrm{\Delta }^{}}{(1\mathrm{\Delta }^{})^{3/2}}$$
(7)
where we have used $`𝑑\theta \mathrm{cos}\theta /(1\mathrm{\Delta }^{}\mathrm{cos}\theta )^2=2\pi \mathrm{\Delta }^{}/(1\mathrm{\Delta }^{})^{3/2}`$. Inserting this into Eq. (4) yields

$$\frac{\mathrm{\Delta }^{}}{(1\mathrm{\Delta }^{})^{3/2}}=\frac{r}{R_0\mathrm{\Delta }}\frac{8\pi }{B_p^2}\left(pp+\frac{B_Z^2}{8\pi }\right)$$
(8)
Evaluating this in the large aspect ratio limit at the plasma edge where $`p=0`$ gives the familiar form $`\mathrm{\Delta }^{}=(r/R_0)(\beta _{pol}+\mathrm{}_i/2)`$, where $`\beta _{pol}`$ is the poloidal beta and $`\mathrm{}_i`$ is the internal inductance per unit length. \[Note that $`\mathrm{}_i/2=B_Z^2/B_p^2(a)=B_{pol}^2cos^2(\theta )/B_p^2(a)=B_{pol}^2/2B_p^2(a)`$, so that the factor of $`1/2`$ comes from the fact that only the vertical magnetic field contributes to the numerator in the shift.\] The nonlinear form of the left-hand side of Eq. (8) has the nice property of insuring that the flux surfaces are well behaved and don’t cross ($`|\mathrm{\Delta }^{}|<1`$) for arbitrarily high $`\beta _{pol}`$, though this equation only rigorously applies if the flux surfaces remained shifted circles, which breaks down at high beta.

## V Summary

We have presented a relatively simple integral constraint on toroidally axisymmetric MHD equilibrium that shows that a normal equilibrium (with nested magnetic flux surfaces, non-singular fields, and a peaked pressure profile that falls monotonically with radius) can not exist if the toroidal current inside any flux surface is negative relative to the total plasma current. This generalizes previous resultsBreslau et al. (2003); Stratton et al. (2002a); Chu and Parks (2002) to arbitrary aspect ratio.

However, the integral constraint does not necessarily prevent negative central current equilibria with non-nested or singular magnetic flux surfaces. Possible examples of this are shown in Figs. 2 and 3. A plasma with nested flux surfaces and negative central current that is initially out of equilibrium could presumably move towards an equilibrium, though it seems most likely that this new equilibrium would have singular or near-singular fields and thus would be subject to strong reconnectionBreslau et al. (2003); Stratton et al. (2002a) if finite resistivity is introduced, changing the topology. There can also be other constraints that limit the accessible class of alternate equilibriaChu and Parks (2002). One might be able to understand the structure of some of these possible solutions in the vicinity of the null surface as a boundary layer analysis of a shock-like solution. But a realistic evaluation of such scenarios would require including finite cross-field transport, viscosity, resistivity, and FLR effects. We leave these issues to future work. Other interesting questions to consider are whether such “non-normal” MHD equilibria are stable to ideal and/or resistive MHD modes and/or are experimentally accessible.

###### Acknowledgements.

We thank Drs. Joshua Breslau, Ming-Sheng Chu, Nathaniel Fisch, Nikolai Gorelenkov, Don Monticello, Paul Parks, Emilia Solano, and Leonid Zakharov for helpful discussions on these topics. GWH also thanks Dr. Paul Rutherford for teaching a course that covered a virial theorem approach to the Shafranov shift. Supported by DOE Contract # DE-AC02-76CH03073. NOTE: After submitting this paper for publication, we became aware of recent related work by Martynov et al.Martynov et al. (2003), who study island equilibria with negative central current in much more detail.


text:
# The expected photospheric emission of GRBs in the internal shock model

## 1. The photosphere of a relativistic wind

We do not discuss here the nature of the source responsible for the initial energy release leading to the GRB. We assume that a relativistic wind has emerged at radius $`r_{\mathrm{acc}}`$ with a variable distribution of the Lorentz factor $`\mathrm{\Gamma }(t)`$ (with $`\overline{\mathrm{\Gamma }}>100`$). The energy and mass injection rates are $`\dot{E}(t)`$ and $`\dot{M}(t)=\dot{E}(t)/\mathrm{\Gamma }(t)c^2`$. The wind production lasts for a total duration $`t_\mathrm{w}`$ (total ejected mass : $`M_\mathrm{w}`$). At $`r=r_{\mathrm{acc}}`$, the wind is still optically thick, whereas at the radius $`r_{\mathrm{IS}}`$ where the internal shock phase starts, it has to be already transparent, as indicated by the dominant non-thermal nature of the observed GRB spectra. The radius at which each layer of the relativistic wind becomes transparent can be estimated (see Daigne & Mochkovitch 2002) : photons emitted by a layer at $`r`$ will successively cross all the layers emitted earlier by the source before they can escape from the relativistic wind at radius $`r_{\mathrm{esc}}(r)`$ : as $`\mathrm{\Gamma }1`$, we usually have $`r_{\mathrm{esc}}(r)r`$. The optical depth of a layer at radius $`r`$ is approximatively given by $`\tau (r)\kappa \dot{E}/8\pi c^3\mathrm{\Gamma }^2r`$, where the opacity $`\kappa `$ is mainly due to the ambient electrons (Thomson opacity). Therefore the photospheric radius ($`\tau 1`$) of the layer is given by : $`r_{\mathrm{ph}}\kappa \dot{E}/8\pi c^3\mathrm{\Gamma }^2`$. Fig. 1 shows this radius for an example where $`\dot{E}=10^{52}\mathrm{erg}.\mathrm{s}^1`$, $`t_\mathrm{w}=10\mathrm{s}`$ and the Lorentz factor is given in the same figure.

## 2. Time profile and spectrum of the photospheric emission

The photospheric luminosity and temperature of each layer can be estimated in the context of the fireball model (see e.g. Piran 1999) :

$$L_{\mathrm{ph}}\dot{E}\left(\frac{r_{\mathrm{ph}}}{r_{\mathrm{acc}}}\right)^{2/3}\mathrm{and}kT_{\mathrm{ph}}kT^0\left(\frac{r_{\mathrm{ph}}}{r_{\mathrm{acc}}}\right)^{2/3},$$
(1)
where the initial fireball temperature is given by $`kT^0=1.3\dot{E}_{52}^{1/4}\mu _1^{1/2}\mathrm{MeV}`$ ($`\mu _1`$ is the mass of the central black hole in units of $`10\mathrm{M}_{}`$). The emitted photons will closely follow a Planck distribution, so that the phostopheric emission leads to a spectrum and a count rate in an energy band $`[E_1,E_2]`$ :

$$\frac{dn^{\mathrm{ph}}(E)}{dEdt}=\frac{\left(1+z\right)L_{\mathrm{ph}}}{_{\mathrm{Planck}}\left(kT_{\mathrm{ph}}\right)^2}\times \frac{x^2}{\mathrm{exp}x1}\mathrm{and}C_{12}^{\mathrm{ph}}=\frac{L_{\mathrm{ph}}}{kT_{\mathrm{ph}}}\frac{1}{_{\mathrm{Planck}}}_{x_1}^{x_2}\frac{x^2dx}{\mathrm{exp}x1},$$
(2)
where $`x=(1+z)E/kT_{\mathrm{ph}}`$, $`x_{1,2}=(1+z)E_{1,2}/kT_{\mathrm{ph}}`$ and $`_{\mathrm{Planck}}=\pi ^4/15`$.

## 3. Comparison with the emission from internal shocks

The internal shock luminosity is assumed to be $`L_{\mathrm{IS}}=f_\gamma \dot{E}`$ where $`f_\gamma `$ is the efficiency for the conversion of kinetic energy into $`\gamma `$-rays by internal shocks. This leads to a spectrum and a count rate in energy band $`[E_1,E_2]`$ given by

$$\frac{dn^{\mathrm{IS}}(E)}{dEdt}=\frac{L_{\mathrm{IS}}}{\left(1+z\right)E_\mathrm{p}^2_{\mathrm{Band}}}\left(\frac{E}{E_\mathrm{p}}\right)\mathrm{and}C_{12}^{\mathrm{IS}}=\frac{L_{\mathrm{IS}}}{\left(1+z\right)E_\mathrm{p}}\frac{1}{_{\mathrm{Band}}}_{x_1}^{x_2}(x)𝑑x$$
(3)
where $`E_\mathrm{p}`$ is the observed peak energy of the internal shock spectrum ($`x_{1,2}=E_{1,2}/E_\mathrm{p}`$) and $``$ is the usual Band spectral shape (with $`_0^+\mathrm{}x(x)𝑑x=_{\mathrm{Band}}`$), characterized by a low and high-energy slope $`\alpha `$ and $`\beta `$ (Band et al. 1993). Fig. 2 shows the spectrum and the time profile obtained in the example of fig. 1 including both the internal shock and the photospheric emission (with $`\alpha =1.0`$, $`\beta =2.25`$ and $`z=1`$). The internal shock emission is computed with the model developped by Daigne & Mochkovitch (1998). The global observed peak energy of the non-thermal photons is 200 keV.

## 4. Discussion

The ratio of the photospheric over internal shock count rates is given by :

$$R_{12}=\frac{C_{12}^{\mathrm{ph}}}{C_{12}^{\mathrm{IS}}}=1.6f_{\gamma }^{}{}_{0.1}{}^{1}\dot{E}_{52}^{1/4}\mu _1^{1/2}\frac{(1+z)E_\mathrm{p}}{200\mathrm{keV}}\frac{_{\mathrm{Band}}}{_{\mathrm{Planck}}}\frac{_{x_1}^{x_2}\frac{x^2dx}{\mathrm{exp}x1}}{_{x_1}^{x_2}(x)𝑑x}.$$
(4)
This expression indicates that the photospheric thermal emission is usually dominant in the X- and $`\gamma `$-ray bands (see fig. 2) ! This is clearly in disagreement with the observations showing that the GRB spectrum is non-thermal in this spectral range. From eq. 4, there are two possible solutions to recover GRBs dominated by the non-thermal internal shock component : (i) increase of $`f_\gamma `$, i.e. more efficient internal shocks. However, our calculations show that even for $`f_\gamma `$ close to unity, the thermal component would still be easily detectable ; (ii) decrease of $`L_{\mathrm{ph}}`$ and $`kT_{\mathrm{ph}}`$, i.e. a less hot and bright photosphere, compared to the prediction of the standard fireball model. If we define $`\lambda `$ as the fraction of the energy released by the source which is initially under thermal form, the photospheric luminosity and temperature are given by :

$$L_{\mathrm{ph}}\lambda \dot{E}\left(\frac{r_{\mathrm{ph}}}{r_{\mathrm{acc}}}\right)^{2/3}\mathrm{and}kT_{\mathrm{ph}}\lambda ^{1/4}kT^0\left(\frac{r_{\mathrm{ph}}}{r_{\mathrm{acc}}}\right)^{2/3},$$
(5)
so that the expression of $`R_{12}`$ in eq. 4 has now to be multiplied by a factor $`\lambda ^{3/4}`$. In the standard fireball model $`\lambda =1`$ but much smaller values can be expected if a large fraction of the energy realeased by the source is initially under magnetic form. Fig. 3 shows exactly the same example as fig. 1, the only difference being that the photospheric emission has been computed with $`\lambda =0.01`$. A non-thermal burst is recovered in agreement with the observations. A small thermal precursor is still present in X-rays. It can be more intense for intermediate values of $`\lambda `$, so that the photospheric emission predicted in this case could explain the X-ray precursors observed in a few GRBs.

## 5. Conclusion

We have computed the photospheric thermal emission of a relativistic wind with a highly variable initial distribution of Lorentz factor. We find that : (1) The photosphere predicted in the standard fireball model for GRBs is usually too hot and luminous. The thermal photospheric emission dominates the non-thermal internal shock emission, even for large internal shock efficiencies. This is in disagreement with the observations. (2) A less luminous photosphere is expected if the energy released by the source is initially injected mainly under magnetic rather than thermal form. In this case, we recover GRB spectra and temporal profiles dominated by the non-thermal internal shock emission. (3) We then obtain X-ray precursors which may correspond to the precursor activity observed in a few GRBs by GINGA and WATCH/GRANAT. These resultats are presented in more details in Daigne & Mochkovitch (2002).

## References

Band, D. et al. 1993, ApJ, 413, 281.

Daigne, F., & Mochkovitch, R. 1998, MNRAS, 296, 275.

Daigne, F., & Mochkovitch, R. 2002, MNRAS, 336, 1271.

Piran, T. 1999, Physics Reports, 314, 575.


text:
# Designability and Thermal Stability of Protein Structures

## I Introduction

Two remarkable features of natural proteins are the simple fact that they fold and the limited number of distinct folds they adopt. Random amino-acid sequences typically do not fold to a unique structure. Rather they have many competing configurations of similar minimum free energy. Nature has evolved sequences that do fold stably, but it is estimated that the total number of qualitatively different folds is only about 1,000 .

To attempt to explain these remarkable features of natural proteins, we have proposed a principle of designability . The designability of a structure is the number of sequences which have that structure as their unique lowest-energy configuration . In a wide range of models, structures are found to differ dramatically in designability, and sequences associated with highly designable structures have unusually high thermal stability . We refer to this connection between the designability of a structure and the stability of its associated sequences as the designability principle. In model studies, highly designable structures are rare. As a result, thermally stable sequences are also rare. We conjecture that the designability principle also applies to real proteins, and that natural protein folds are exceptional, highly designable structures.

In this article, we review the designability principle. We start from a minimal model of protein structure in which the designability of a structure can be understood geometrically as the size of its basin of attraction in sequence space. More detailed models, including all 20 amino-acid types and off-lattice backbone configurations, reinforce the basic principle and provide a framework for the design of qualitatively new protein folds.

## II Purely Hydrophobic (PH) Model

Generally, the folding of proteins relies on the formation of a hydrophobic core of amino acids. Consideration of hydrophobicity alone leads to a very simple description of proteins–the “purely hydrophobic” (PH) model . In this model, sequences consist of only two types of amino acids, hydrophobic and polar . Structures are compact walks on a cubic or square lattice. An example of a $`6\times 6`$ square structure is shown in Fig. 1(a). As indicated, each structure consists of core sites surrounded by surface sites. The energy of a particular sequence folded into a particular structure is the number of hydrophobic amino acids occupying core sites, multiplied by $`1`$,

$$E=\underset{i=1}{\overset{N}{}}s_ih_i.$$
(1)
A binary string $`\{s_i\}`$ represents each folded structure: $`s_i=1`$ if site $`i`$ along the chain is in the core, and $`s_i=0`$ if the site is on the surface. Similarly, a binary string $`\{h_i\}`$ represents each sequence: $`h_i=1`$ if the $`i`$th amino acid in the sequence is hydrophobic, and $`h_i=0`$ if the amino acid is polar.

Within the PH model, structures differ dramatically in their designability. In practice, the designability is obtained by sampling a large number of binary sequences, and, for each sequence, recording the unique lowest-energy structure if there is one. Finally, the number of sequences which map to, i.e. “design”, each structure is summed to give the designability $`N_S`$ of the structure. Fig. 1(b) shows a histogram of designability $`N_S`$ for compact $`6\times 6`$ structures. There are 57,337 structures, with 30,408 distinct binary strings. Most structures have a designability around 50, but a small number of structures have designabilities more than 10 times this high. If sequences were randomly assigned to structures, the result would the Poisson distribution which is shown for comparison, and there would be no structures with such high designability.

Importantly, the PH model has a simple geometrical representation that explains both the wide range of designabilities and the close connection between designability and thermal stability. To find the relative energies of different structures for a given sequence, the energy in Eq. (1) can be replaced by

$$E=\underset{i=1}{\overset{N}{}}(h_is_i)^2.$$
(2)
This replacement is allowed because the extra term $`h_i^2`$ is a constant for a given sequence, and the other extra term $`s_i^2`$ is also a constant, equal to the number of core sites, for all compact structures. Eq. (2) indicates that the energy of a sequence folded into a particular structure is simply the Hamming distance  between their respective binary strings. So, for a given sequence, the lowest-energy structure is simply the closest structure. The designability of a structure is therefore the exclusive volume of binary strings (sequences) that lie closer to it than to any other structure, as shown schematically in Fig. 2.

The wide range of structure designabilities can be traced to the nonuniform density of structures in the space of binary strings. Most structures are found in dense regions, i.e. in clusters of structures with similar patterns of surface and core sites. Structures found in these crowded regions have small exclusive volumes, and so, by definition, have small designabilities. In fact, many groups of distinct structures share an identical surface-core pattern (binary string) and therefore have zero designability. In contrast, a few structures fall in low-density regions, that is they have unusual surface-core patterns, and so have large exclusive volumes. These are the highly designable structures. In Fig. 3(a), we plot the number of structures $`n(d)`$ at a Hamming distance $`d`$ from a structure with low, intermediate, and high designability, respectively. It shows that both low- and high-density neighborhoods typically have a large spatial extent, reaching nearly halfway across the space of binary strings.

The geometrical representation of the PH model makes clear the connection between thermal stability and designability. A sequence is considered to be thermally unstable if it has a small or vanishing energy gap between its lowest-energy structure and all other structures, and if there are many such competing structures. In the PH model, the energy of a sequence folded into a particular structure is the distance between their binary strings. A sequence which folds to a structure in a dense region (cf. Fig. 2) will necessarily lie close to many other structures, and will therefore have many competing low-energy conformations. Even if the sequence perfectly matches the structure, with hydrophobic amino acids at all core sites and polar amino acids at all surface sites, the high surrounding density of structures with similar surface-core patterns implies a large number of competing folds. This is the hallmark of thermal instability. Therefore, the low-designability structures, which are found in high-density regions, will have associated sequences which are thermally unstable.

In contrast, if a sequence folds to a structure in a low-density region, there will be relatively few nearby structures, and so relatively few competing folds. These sequences will be thermally stable. Therefore, the highly designable structures, from low-density regions, will have associated sequences of high thermal stability. This is the designability principle in a nutshell–high designability and thermal stability are connected because both arise from low-density regions in the space of binary strings which represent folded structures.

A measure of the “neighborhood density” of structures around a particular structure is the variance $`\gamma `$ of the quantity $`n(d)`$ shown in Fig. 3(a). The variance $`\gamma `$ is directly related to the thermal stability–smaller $`\gamma `$ implies lower neighhood density and hence higher thermal stability. In Fig. 3(b) we plot this variance as a function of designability. It shows a strong correlation between the designability and the thermal stability.

Since low-energy competing structures also represent kinetic traps, one expects the thermally stable sequences associated with highly designable structures to be fast folders as well. This has been tested for a lattice model closely related to the PH model .

## III Miyazawa-Jernigan (MJ) Matrix Model

Natural proteins contain 20 amino acids, not two, and their interactions are more complicated than simple hydrophobic solvation. Some of these real-world features are captured in the Miyazawa-Jernigan (MJ) matrix model. The MJ matrix is a set of amino-acid interaction energies inferred from the propensities of different types of amino acids to be neighbors in natural folded structures . The model assigns the appropriate energy from the MJ matrix to every pair of amino acids that are on neighboring lattice sites, but are not adjacent (covalently bonded) on the chain, as indicated in Fig. 4(a) . In studies using the MJ-matrix model, there are generally too many possible sequences (20N) to enumerate, but the relative designabilities of structures can be obtained accurately by random sampling of sequences.

Fig. 5(a) shows a histogram of designability for compact $`6\times 6`$ structures obtained using the MJ matrix of interaction energies. The form of the distribution is very similar to the PH-model histogram, including the tail of highly designable structures (Fig. 1(b)). There is also a strong correlation between thermal stability and designability in the MJ-matrix model . For thermal stability, one can use some measure of “neighborhood” density of states. We find that in the MJ-matrix model, the thermal stability of a sequence folded into a structure is well correlated with the local energy gap  between the lowest-energy structure and the next lowest. Fig. 5(b) shows the energy gap averaged over sequences which fold to structures of a given designability $`N_S`$. With increasing designability, there is a clear increase in the average gap, and hence in the thermal stability of associated sequences. The results of the MJ-matrix model are very similar to those obtained with the PH model. Indeed, the same structures are found to be highly designable in both models, including the same top structure shown in Fig. 4(a). The most designable $`3\times 3\times 3`$ structure is shown in Fig. 4(b). Qualitatively, the results of the MJ-matrix model are the same for three-dimensional structures (Fig. 6) as for two-dimensional ones.

Why are the results of the purely hydrophobic model and the Miyazawa-Jernigan-matrix model so similar? In fact, both models are dominated by hydrophobic solvation energies. The interaction energy between any two amino acids $`i`$ and $`j`$ in the MJ matrix can be well approximated by $`(h_i+h_j)`$, where $`h_i`$ is an effective hydrophobicity for each amino acid . This implies that the energy of formation of a non-covalent nearest-neighbor pair is simply the desolvation energy of shielding one face of each amino acid from the surrounding water. As a result, the MJ-matrix model can be viewed as a variant of the PH model in which there are 20 possible values of hydrophobicity instead of just two. An additional distinction is that the MJ-matrix model has a range of different site types (core, edge, and corner in two dimensions; core, edge, face, and corner in three dimensions) rather than just surface and core as in the PH model. Overall, these differences are not enough to alter the basic designability principle, or even to change the set of highly designable $`6\times 6`$ structures.

For the MJ-matrix model, one can still construct a space of structure strings, now including several levels of solvent exposure between 0 (most exposed) and 1 (most buried). As in the PH model, some regions of this space are dense with structures and some have few structures. Structures with similar surface-exposure strings compete for sequences. As a result, structures in high-density regions have small basins of attraction for sequences, and structures in low-density regions have large basins. Moreover, sequences associated with structures in low-density regions have few competing conformations and are the most thermally stable. Therefore, the designability principle holds in the MJ-matrix model for the same reason it does in the PH model: high designability and high thermal stability are connected because both arise in low-density regions in the space of strings, i.e. the space of surface-exposure patterns of structures.

Lattice-protein models in which hydrophobic solvation does not dominate may show different behavior. For example, Buchler and Goldstein reported results of a variant of the MJ-matrix model in which the dominant hydrophobic term $`(h_i+h_j)`$ had been subtracted out . They found a set of highly designable structures different from that obtained with the full MJ matrix, and similar to the set obtained for a random pairing potential between amino acids.

## IV Off-lattice models

Natural proteins fold in three dimensions, and their main degrees of freedom are bond rotations. Does the designability principle extend to off-lattice models with more realistic degrees of freedom? One model for which designability has been studied off-lattice is a 3-state discrete-angle model, of the type introduced by Park and Levitt . The results strongly confirm the designability principle, and suggest the possibility of creating new, highly designable folds in the laboratory .

The main degrees of freedom of a protein backbone are the dihedral angles $`\varphi `$ and $`\psi `$. Certain pairs of $`\varphi `$-$`\psi `$ angles are preferred in natural structures, since they lead to conserved secondary structures such as $`\alpha `$-helices and $`\beta `$-strands . Discrete-angle models for protein structure take advantage of these preferences by allowing only certain combinations of angles. For an $`m`$-pair model, the total number of backbone structures grows as $`m^N`$. With $`m=3`$, it is possible to computationally enumerate all structures up to roughly $`N=30`$ amino acids.

Figure 7(a) shows an example of a protein backbone of length $`N=23`$ generated using a 3-state model with ($`\varphi `$,$`\psi `$) = (-95,135), (-75,-25), and (-55,-55), where the first pair of angles corresponds to a $`\beta `$-strand and the second two correspond to variants of $`\alpha `$-helices. Structures are decorated with spheres representing sidechains, as shown in Fig. 7(b). Only compact self-avoiding structures are considered as possible protein folds.

To assess designability of these off-lattice structures, the solvent-exposed area of each sidechain sphere is evaluated. An energy of hydrophobic solvation is defined as in Eq. (1) by $`E=_{i=1}^Ns_ih_i`$, where now $`s_i`$ is the fractional exposure to solvent of the $`i`$th sidechain, and the $`h_i`$ are amino-acid hydrophobicities. Figure 8(a) shows a histogram of designability for the 3-state model. There is a wide range of designabilities, with a tail of very highly designable structures. A strong correlation exists between the designability of a structure and the thermal stability of its associated sequences, as shown in Fig. 8(b).

The designability principle evidently applies to the 3-state model as well as to the lattice models discussed earlier. This is not surprising, because folding in the 3-state model is also driven by hydrophobic solvation. Each structure in the model is represented by a string of sidechain solvent exposures, represented by real numbers between 0 and 1. Again, the space of these strings has high- and low-density regions, with the, by now familiar, relation between low density and high designability and thermal stability leading to the designability principle.

A major advantage of the 3-state model is that it addresses structures that a real polypeptide chain can adopt. Among the highly designable folds, one recovers several recognizable natural structures, including an $`\alpha `$-turn-$`\alpha `$ fold and a zincless zinc-finger. In addition, some of the highly designable folds, including a $`\beta `$-$`\alpha `$-$`\beta `$ structure, have not been observed in nature as independent domains. Results of our effort to create this fold in the laboratory are encouraging .

## V Discussion and Conclusion

The designability principle has been explored in a number of models for proteins, including all 20 amino acids and realistic backbone conformations. In these models, the strong link between designability and thermal stability can be traced to the dominance of the hydrophobic solvation energy. Whenever hydrophobicity is dominant, each structure can be reduced to its pattern of solvent exposure. In the same vein, each sequence can be reduced to its pattern of hydrophobicity. Sequences will fold so as to best match their hydrophobic residues to the buried sites of structures. Both designability (number of sequences per structure) and thermal stability depend on a competition among structures with similar patterns of solvent exposure. Highly designable structures are those with unusual patterns of surface exposure, and therefore with few competitors. This lack of competitors also implies that the sequences folding to highly designable structures are thermally stable.

Since hydrophobicity is generally accepted to be the dominant force for folding of real proteins, the designability principle may provide a guide to understanding the selection of natural protein structures. Of course, real proteins are held together by forces other than hydrophobicity. Next to hydrophobicity, the formation of hydrogen bonds is the most important factor in determining how a typical protein folds. The backbone hydrogen bonds of $`\alpha `$-helices and $`\beta `$-sheets help stabilize particular folds. These secondary structures can be incorporated within the framework of designability as a favorable energy bias for formation of $`\alpha `$-helices and $`\beta `$-sheets.

One way to incorporate hydrogen bonding in the design of new protein folds is to specify in advance the secondary structure of the protein. This approach has the added advantage of greatly reducing the number of degrees of freedom. The desired secondary structures can be designed into the sequence via the propensities of particular amino acids to form $`\alpha `$-helices and $`\beta `$-strands.

This approach to design was recently carried out for four-helix bundles . Compact, self-avoiding structures consisting of four tethered 15-residue $`\alpha `$-helices were generated and assessed for designability. Figure 9 shows the four most designable distinct folds, which closely correspond to natural four-helix bundles. As shown in Fig. 10, the histogram of designability for the four-helix model has the characteristic long tail of highly designable structures.

The principle of designability has been motivated here in terms of hydrophobic solvation. More generally, the dependence of both designability and thermal stability on a competition among structures broadens the application of the principle. For example, designability and thermal stability have been found to correlate in non-solvation models including random-interaction models  and folding of two-letter RNA . In the future, we hope that designability will provide a guide to the design of new structures both for polymers other than proteins and for solvents other than water.

We gratefully acknowledge the contributions of many coworkers in developing the notion of designability, in particular Eldon Emberly, Robert Helling, Régis Mélin, Jonathan Miller, Tairan Wang, and Chen Zeng.


text:
# Hadron formation and attenuation in deep inelastic lepton scattering off nuclei

## Abstract

We investigate hadron formation in deep inelastic lepton scattering on $`N`$, $`Kr`$ and $`Xe`$ nuclei in the kinematic regime of the HERMES experiment. The elementary electron-nucleon interaction is described within the event generator PYTHIA while a full coupled-channel treatment of the final state interactions is included by means of a BUU transport model. We find a good agreement with the measured charged hadron multiplicity ratio $`R_M^h`$ for $`N`$ and $`Kr`$ targets by accounting for the deceleration and absorption of the primarily produced particles as well as for the creation of secondary hadrons in the final state interactions.

journal: Physics Letters B

, , ,

Hadron production in deep inelastic lepton-nucleus scattering (DIS) offers a promising tool to study the physics of hadronization Kop . The reaction of the exchanged virtual photon (energy $`\nu `$, virtuality $`Q^2`$) with a bound nucleon leads to the production of several hadrons. While the primary production is determined by the fragmentation function – in medium possibly different from that in vacuo – the number of ultimately observed hadrons and their energy distribution depends also on their rescattering in the surrounding nuclear medium. Consequently, the particle spectrum of a lepton-nucleus interaction will differ from that of a reaction on a free nucleon. In order to explore such attenuation effects the HERMES collaboration has investigated the energy $`\nu `$ and fractional energy $`z_h=E_h/\nu `$ dependence of the charged hadron multiplicity ratio

$$R_M^h(z,\nu )=\left(\frac{N_h(z,\nu )}{N_e(\nu )}\right)_A/\left(\frac{N_h(z,\nu )}{N_e(\nu )}\right)_D$$
(1)
in DIS off $`N`$ HERMES1  and $`Kr`$ HERMES2  nuclei. Here $`N_h(z,\nu )`$ represents the number of semi-inclusive hadrons in a given ($`z,\nu `$)-bin and $`N_e(\nu )`$ the number of deep inelastically scattered leptons in the same $`\nu `$-bin. It was suggested in Ref. HERMES1 , that a phenomenological description of the $`R_M^h`$ data can be achieved if the formation time, i.e. the time that elapses from the moment when the photon strikes the nucleon until the reaction products have evolved to physical hadrons, is assumed to be proportional to $`(1z_h)\nu `$ in the target rest frame. This $`(1z_h)\nu `$ dependence of the formation time $`\tau _f`$ is compatible with the gluon-bremsstrahlung model of Ref. Kop . In the investigations of Ref. HERMES1  any interaction of the reaction products with the remaining nucleus during this formation time has been neglected. After the formation time the hadrons could get absorbed according to their full hadronic cross section. Another interpretation of the observed $`R_M^h`$ spectra – as being due to a combined effect of a rescaling of the quark fragmentation function in nuclei due to partial deconfinement as well as the absorption of the produced hadrons – has recently been given by the authors of Ref. Acc02 . Furthermore, calculations based on a pQCD parton model Wang ; Arl03  explain the attenuation observed in the multiplicity ratio solely by partonic multiple scattering and induced gluon radiation completely neglecting any hadronic final state interactions (FSI). It has already been pointed out by the authors of Ref. Acc02  that a shortcoming of the existing models is the purely absorptive treatment of the FSI. We will avoid this problem by using a semi-classical coupled-channel transport model Eff99  which already has been employed to describe high energy photo- Photo  and electroproduction Fal02  off nuclei.

We point out that the formation time also plays an important role in studies of ultra-relativistic heavy-ion reactions. For example, the observed quenching of high transverse momentum hadrons in $`Au+Au`$ reactions relative to $`p+p`$ collisions is often thought to be due to jet quenching in a quark gluon plasma. However, the attenuation of high $`p_T`$ hadrons might at least partly be due to hadronic rescattering processes Gal03 ; Cas03 .

Unfortunately, the average hadron formation time is not well known and the number of $`\tau _f1`$ fm/c as used commonly in the Bjorken estimate for the energy density bjorken  is nothing but an educated guess. The nonperturbative nature of this number – due to time scales of $``$ 1 fm/c and hadronic size scales of 0.5–1 fm – excludes perturbative evaluation schemes; it is hard to calculate $`\tau _f`$ from first principles and formation times cannot be addressed in present lattice QCD simulations. One might expect that the rather successful string models LUND  shed some further light on this number, since the intrinsic time scale $`\tau _0`$ for the $`q\overline{q}`$ formation vertex can be related to the fragmentation function and string tension, respectively Anders . However, the actual parameters employed in current transport codes are not unique, with hadron formation times ranging from 0.3–2 fm/c Cas99 ; UrQMD1 , depending on the flavor, momentum and energy of the created hadrons. In fact, the rapidity and transverse mass spectra from relativistic nucleus-nucleus collisions are not very sensitive to the formation time $`\tau _f`$ Cas03 . It is therefore essential to check whether these times are compatible with constraints extracted from reactions, where the collision geometry is much better under control.

To this aim the attenuation of antiprotons produced in $`p+A`$ reactions at the AGS energies of 12.3 GeV and 17.5 GeV has been investigated on various nuclear targets in Ref. Cas02  and a range of values for $`\tau _f=0.40.8`$ fm/c has been extracted in comparison to the data from the E910 Collaboration E910 .

In our present approach the lepton-nucleus interaction is split into two parts: 1) In the first step the virtual photon is absorbed on a bound nucleon of the target; this interaction produces a bunch of particles that in step 2) are propagated within the transport model. Coherence length effects in the entrance channel, that give rise to nuclear shadowing, are taken into account as described in Ref. Fal02 . The virtual photon-nucleon interaction itself is simulated by the Monte Carlo generator PYTHIA v6.2 pythia  which well reproduces the experimental data of Ref. HERMES3  on a hydrogen target. Depending on whether the photon interacts directly with the nucleon or via one of its hadronic fluctuations ($`\rho ^0,\omega ,\mathrm{\Phi },J/\mathrm{\Psi }`$ or a perturbative $`q\overline{q}`$ fluctuation) the reaction leads to the excitation of one or more hadronic strings which in our approach are assumed to fragment very rapidly into colorless prehadrons.

The time that the reaction products need to evolve to physical hadrons, i.e., the production time $`t_p`$ of the prehadrons plus the time needed to build up the hadronic wave function, we denote as formation time $`t_f`$ in line with the convention in transport models. For simplicity we assume that the formation time is a constant $`\tau _f`$ in the rest frame of each hadron and that it does not depend on the particle species. Due to time dilatation the formation time $`t_f`$ in the laboratory frame is then proportional to the particle’s energy

$$t_f=\gamma \tau _f=\frac{z_h\nu }{m_h}\tau _f.$$
(2)
The size of $`\tau _f`$ can be estimated by the time that the constituents of the hadrons need to travel a distance of a typical hadronic radius (0.5–0.8 fm).

To illustrate the situation shortly after the photon nucleon reaction Fig. 1 shows the excitation and fragmentation of a hadronic string in a deep inelastic scattering process. For simplicity we do not show any gluon bremsstrahlung of the struck quark in this figure. Note, however, that the possibility of such final state gluon radiation plus subsequent $`q\overline{q}`$ splitting is included in the PYTHIA part of our model and leads to the creation of additional strings pythia . It has been emphasized in Ref. Cio02  that the string propagating through the nucleus is a rather short (white) object of length $`1`$ fm since the slow end of the string is accelerated very fast. When the primary string fragments – due to the creation of $`q\overline{q}`$ pairs from the vacuum – new colorless prehadrons are produced, which we propagate in space-time.

As discussed in Ref. Kop03  the production time $`t_p`$ of these prehadrons has to be distinguished from the total formation time $`t_f`$ of the final hadrons which is dilated according to Eq. (2). It has been confirmed Bia87  within the Lund model that the production time $`t_p`$ vanishes for $`z_h0`$ and $`z_h1`$. In the present numerical realization of our model we first approximate this behavior by setting the production time $`t_p`$ to zero for all prehadrons, but will also discuss the effect of a finite production time at the end of this study.

Right after the photon nucleon interaction the primary string should interact with a hadronic cross section because its transverse size is essentially that of the original nucleon. Motivated by the constituent quark model we assume that this hadronic cross section is shared by the quark/diquark at the string ends and after the fragmentation by the so called leading prehadrons that contain this quark or diquark. Our PYTHIA simulations show that in most cases the prehadrons with $`z_h1`$ are such leading hadrons since they contain constituents (valence- or sea-quarks) from the target nucleon or the hadronic component of the photon. They can therefore interact directly after the photon-nucleon interaction with a constant effective cross section which we denote as $`\sigma _{\mathrm{lead}}`$. The cross sections of the other prehadrons, that solely contain quarks and antiquarks created from the vacuum, emerge at intermediate $`z_h`$. They are assumed to be non-interacting until $`t_f`$. This assures that the summed cross section of the complete final state right after the photon-nucleon interaction is approximately that of the original nucleon. Each time when a new hadron has formed, the summed cross section rises just like in the approach of Ref. Cio02 . After the hadron formation time $`t_f`$ by definition all hadrons interact with their full hadronic cross section $`\sigma _h`$. Note, that our concept of leading hadrons is in accordance with those of other transport models for high energy reactions Cas03 ; Cas99 ; Geiss ; URQMD2 .

Since the lighter (intermediate $`z_h`$) hadrons have large formation times in the target frame (see Eq. (2)) they may escape the nucleus without being attenuated if they are non-leading. However, many ($`2/3`$) of the observed hadrons with intermediate $`z_h`$ are not directly produced in the string fragmentation but stem from decays of the much heavier vector mesons with correspondingly shorter formation times. These vector mesons may therefore form inside the nuclear volume and thus be subjected to FSI. The effect of the FSI, finally, will depend dominantly on the nuclear geometry, i.e. the size of the target nucleus.

In our present study the FSI are described by a coupled-channel transport model based on the Boltzmann-Uehling-Uhlenbeck (BUU) equation which describes the time evolution of the phase space density $`f_i(\stackrel{}{r},\stackrel{}{p},t)`$ of particles of type $`i`$ that can interact via binary reactions. These particles involve nucleons, baryonic resonances and mesons ($`\pi `$, $`\eta `$, $`\rho `$, $`K`$, …) that are produced either in the primary reaction or during the FSI. In this work we also have to account for the prehadrons emerging from the string fragmentation. For a particle species $`i`$ the BUU equation takes the form

$$\left(\frac{}{t}+\frac{H}{\stackrel{}{p}}\frac{}{\stackrel{}{r}}\frac{H}{\stackrel{}{r}}\frac{}{\stackrel{}{p}}\right)f_i(\stackrel{}{r},\stackrel{}{p},t)=I_{coll}[f_1,\mathrm{}f_i,\mathrm{},f_M],$$
(3)
where the Hamilton function $`H`$ includes a position and momentum dependent mean-field potential for baryons. The collision integral on the right hand side accounts for the creation and annihilation of particles of type $`i`$ in a collision as well as elastic scattering. The transition rates are determined from the particular (vacuum) cross sections. For fermions Pauli blocking is taken into account in $`I_{coll}`$ via blocking factors. The prehadrons are treated like ordinary hadrons except for their modified interaction cross section and the fact that they are not allowed to decay during the formation time. The BUU equations of each particle species $`i`$ are coupled via the collision integral and the mean field in case of baryons. The resulting system of coupled differential-integral equations is solved via a test particle ansatz for the phase space density. For further details of the transport model we refer the reader to Ref. Eff99 .

Most of the hadronic FSI happen with invariant energies $`\sqrt{s}2.2`$ GeV and are described within the Lund string formation and decay scheme LUND  as also implemented in the transport approaches Eff99 ; Photo ; Fal02  as well as Cas99 ; Cas02 ; Geiss . The important difference between a purely absorptive treatment of the FSI and the coupled-channel description provided by the BUU model is that in an interaction with a nucleon a hadron might not only be absorbed or recreated but also be decelerated in an elastic or inelastic collision. Furthermore, it may in addition produce several low energy particles. In the case of electroproduction of hadrons these interactions leads to a redistribution of strength from the high $`z_h`$ part of the hadron energy spectrum to lower values of the energy fraction $`z_h`$.

In our calculation we employ all kinematic cuts of the HERMES experiment as well as the geometrical cuts of the detector. In actual numbers: we require for the Bjorken scaling variable $`x=\frac{Q^2}{2m_N\nu }>0.06`$ (with $`m_N`$ denoting the nucleon mass), for the photon virtuality $`Q^2>1`$ GeV2 and for the energy fraction of the virtual photon $`y=\nu /E_{beam}<0.85`$. In addition, the PYTHIA model introduces a lower cut in the invariant mass of the photon-nucleon system at $`W=4`$ GeV that is above the experimental constraint $`W>2`$ GeV. This limits our calculations to minimal photon energies of $`\nu _{min}=8.6`$ GeV as compared to $`\nu _{min}=7`$ GeV in the HERMES experiment and leads to a suppression of high $`Q^2`$ events at energies below $`\nu 15`$ GeV.

In Fig. 2 we present the average values of $`Q^2`$ and $`\nu `$ ($`z_h`$) for our simulated event samples on $`N`$ and $`Kr`$ as a function of $`z_h`$ ($`\nu `$) in comparison to the experimental quantities private . In order to compare with the $`\nu `$ dependence of $`Q^2`$ and $`z_h`$ we only account for hadrons with $`z_h>0.2`$ as in the HERMES experiment. For both the $`N`$ and the $`Kr`$ target the average kinematical variables are well reproduced by our model. The underestimation of $`Q^2`$ at low photon energies $`\nu <12`$ GeV is due to the PYTHIA cut in $`W4`$ GeV which suppresses higher values of $`Q^2`$ at low photon energies $`\nu `$. This is also the reason why the average value of $`Q^2`$ in the $`z_h`$ spectrum comes out slightly too low within our model.

We proceed with a discussion of the actual results of our PYTHIA + BUU simulations. In Fig. 3 we show the calculated multiplicity ratio $`R_M^h`$ for $`N`$ and $`Kr`$ targets using an ’estimated’ formation time of 0.5 fm/c and different values for the leading prehadron cross section $`\sigma _{\mathrm{lead}}`$. The data have been taken from Refs. HERMES1 ; HERMES2 . Since the particles with $`z_h`$ close to 1 are predominantly leading hadrons we can use the high $`z_h`$ part in the fractional energy spectrum to obtain information on $`\sigma _{\mathrm{lead}}`$. The data for both nuclei indicate that $`\sigma _{\mathrm{lead}}`$ has to be in the range 0–0.5 $`\sigma _h`$ with $`\sigma _h`$ ($`h=\pi ^\pm ,K^\pm ,p,\mathrm{}`$) taken from PDG . For the heavier nucleus, $`Kr`$, we clearly underestimate the hadron attenuation with $`\sigma _{\mathrm{lead}}`$=0 since most of the particles, especially those with large energies, escape the nucleus due to time dilatation. If one wants to describe the strong attenuation of hadrons at large fractional energy $`z_h`$ without any prehadronic interactions one would need an unphysical short formation time $`\tau _f<0.1`$ fm/c. This, however, is ruled out by the measured $`\nu `$ dependence of $`R_M^h`$ since a vanishing formation time leads to a multiplicity ratio $`R_M^h`$ which is considerably too low (see dash-dotted and dash-dot-dotted curves in Fig. 4).

The dotted line in Fig. 3 shows the result of a calculation where the leading prehadrons interact with the full hadronic cross section $`\sigma _h`$; this leads to a too strong attenuation of charged hadrons. A good agreement with the data is achieved for $`\sigma _{\mathrm{lead}}=0.33\sigma _h`$ during the formation time $`\tau _f`$ as can be seen in Fig. 4. We note, that this value for $`\sigma _{\mathrm{lead}}`$ represents an average value over time from the virtual photon-nucleon interaction to the actual hadron formation time. For a detailed investigation we refer the reader to a forthcoming study F2 .

In Fig. 4 we investigate the influence of different formation times $`\tau _f`$ on $`R_M^h`$ using the effective cross section $`\sigma _{\mathrm{lead}}=0.33\sigma _h`$. We find, that formation times $`\tau _f0.3`$ fm/c are needed to describe the experimental data with little sensitivity to higher values. This is compatible to the range of values extracted from the antiproton attenuation studies in Ref. Cas02 . The steep rise of $`R_M^h`$ for $`z_h0.2`$ is caused by the energy loss and the production of low energy secondary particles in elastic and inelastic FSI. We mention that in models dealing with purely absorptive FSI a ratio much larger than one can only be explained by a drastic change of the fragmentation function.

We note that the photon energy dependence of the ratio $`R_M^h`$ for the $`Kr`$ target is less well reproduced for energies below 14 GeV. The reason for this can be traced back to the $`W4`$ GeV cut of the PYTHIA model. Our simulations show that the number of leading hadrons decreases with $`Q^2`$ for fixed energy, but the W$`4`$ GeV cut discards larger values of $`Q^2`$ at energies $`\nu <15`$ GeV. For higher values of $`Q^2`$ the importance of DIS events rises as compared to events, where the photon interacts via a vector meson fluctuation (VMD events). In the latter case one has initially 5 constituent quarks and antiquarks and thus gets more leading prehadrons than in a DIS event. This leads to a mismatch of leading and nonleading prehadrons in the region $`z_h0.4`$. Since the number of prehadrons created in an electron-nucleon collision decreases exponentially with increasing $`z_h`$, this region contributes dominantly to the $`z_h`$ integrated $`\nu `$ spectrum. This deficiency of the present model will be cured in a more detailed upcoming work F2 .

In the pQCD parton model of Ref. Wang  the multiple parton scattering leads to a modification of the fragmentation function and predicts a hadron attenuation $`A^{2/3}`$. In our approach, however, the fragmentation is assumed to be decided on time scales of the nucleon dimension itself such that only the ’free’ fragmentation function enters. All attenuation effects then are attributed to FSI of the leading and secondary (pre-)hadrons. In order to distinguish experimentally between the different concepts, it is thus important to get the scaling with target mass $`A`$. To this aim Fig. 4 also shows predictions for a $`Xe`$ target. In accordance with the authors of Ref. Acc02  we predict only a small change in the multiplicity spectra compared to the $`Kr`$ target such that the scaling exponent is lower than $`2/3`$.

Finally, we discuss the effect of a finite production time $`t_p`$ (in the lab frame) for prehadrons which we adopt, for consistency, from the Lund model Bia87 :

$$t_p=\left(\frac{\mathrm{ln}(1/z^2)1+z^2}{1z^2}\right)\frac{z\nu }{\kappa }.$$
(4)
Here $`\kappa 1`$ GeVfm-1 denotes the string tension. In the following we assume no interaction before $`t_p`$ and the full hadronic cross section for prehadrons for $`tt_p`$, i.e., there is no dependence on the formation time $`t_f`$ anymore. The dashed line in Fig. 5 represents the result of a Glauber-like treatment of the FSI where every time a prehadron interacts with another particle it is removed from the outgoing channel. As the authors of Ref. Acc02  we get a good description of the $`z`$-dependence of the multiplicity ratio. However, the $`\nu `$ spectrum is not attenuated strongly enough at the higher $`\nu `$ (r.h.s. of Fig. 5). The solid curves show the effect of the coupled channels, i.e., a particle is not only absorbed in a collision but produces a bunch of low energy particles, thereby shifting strength to the low $`z`$ part of the spectrum and thereby underestimating the attenuation at low $`z`$. Similarly, the attenuation is also too weak in the $`\nu `$-spectrum. Since an additional formation time with reduced cross sections would further enhance these discrepancies, we conclude that the data cannot be described with the production time (4). For a detailed investigation of finite production times and alternative models we refer the reader to a forthcoming work F2 .

In summary, we have shown that one can describe the experimental data of the HERMES collaboration for hadron attenuation on nuclei without invoking any changes in the fragmentation function due to gluon radiation. In our dynamical studies, that include the most relevant FSI, we employ only the ’free’ fragmentation function on a nucleon and attribute the hadron attenuation to the deceleration of the produced (pre-)hadrons due to FSI in the surrounding medium. We find that in particular the $`z`$-dependence of $`R_M^h`$ is very sensitive to the interaction cross section of leading prehadrons and can be used to determine $`\sigma _{\mathrm{lead}}`$. The interaction of the leading prehadrons during the formation time could be interpreted as an in-medium change of the fragmentation function, which however could not be given in a closed form. The extracted average hadron formation times of $`\tau _f0.3`$ fm/c are compatible with the analysis of antiproton attenuation in $`p+A`$ reactions at AGS energies Cas02 . In an upcoming work we will investigate in detail the spectra for different particle species ($`\pi ^\pm ,K^\pm ,p,\overline{p}`$) to examine, if the formation times of mesons and antibaryons are about equal. In addition we will improve our model to describe the primary photon-nucleon reaction below the PYTHIA threshold of $`W`$ 4 GeV.

The authors acknowledge valuable discussions with N. Bianchi, E. Garutti, C. Greiner, V. Muccifora and G. Van Steenhoven. This work was supported by DFG.