text:
# Generalized Ward identity and gauge invariance of the color-superconducting gap

## Abstract

We derive a generalized Ward identity for color-superconducting quark matter via the functional integral approach. The identity implies the gauge independence of the color-superconducting gap parameter on the quasi-particle mass shell to subleading order in covariant gauge.

Due to asymptotic freedom of quantum chromodynamics (QCD) the strong coupling constant $`g`$ becomes small at large momentum asymp . In cold, dense quark matter, the typical momentum scale is of the order of the quark chemical potential $`\mu `$. Consequently, at asymptotically large densities, where $`\mu \mathrm{\Lambda }_{\mathrm{QCD}}`$, quark matter is a weakly coupled system. In this case, the dominant interaction between two quarks is single-gluon exchange, which is attractive in the color-antitriplet channel. Due to Cooper’s theorem bcs , at sufficiently low temperatures the quark Fermi surface becomes unstable with respect to the formation of Cooper pairs. Since this is analogous to what happens in ordinary superconductors bcs , this phenomenon was termed color superconductivity bailin . Recently, color superconductivity has been extensively studied within QCD at weak coupling ris03 ; son99 ; pis00 ; ris01 ; miransky , in phenomenological Nambu–Jona-Lasinio-type models alford1 ; Neumann , and in effective theories son00 ; hong-eff ; cas99 ; schafer ; ginzburg-landau .

In a superconductor, exciting particle-hole pairs costs at least an energy amount $`2\varphi _0`$, where $`\varphi _0`$ is the value of the superconducting gap parameter at the Fermi surface for zero temperature. The gap parameter can be computed from a gap equation derived in mean-field approximation which involves single-gluon exchange and bare quark-quark-gluon vertices. Schematically, this gap equation can be written in the form pis00

$$\varphi _0=g^2\left[\zeta \mathrm{ln}^2\left(\frac{\mu }{\varphi _0}\right)+\beta \mathrm{ln}\left(\frac{\mu }{\varphi _0}\right)+\alpha \right]\varphi _0.$$
(1)
In weak coupling, $`g1`$, the solution is SchaferWilczek ; pis00 ; son99 ; miransky

$$\varphi _0=2b\mu \mathrm{exp}\left(\frac{c}{g}\right)\left[1+O(g)\right].$$
(2)
The first term in Eq. (1) contains two powers of the logarithm $`\mathrm{ln}(\mu /\varphi _0)`$: one is the same as in BCS theory bcs  (the so-called “BCS logarithm”) and the other arises from the collinear exchange of almost static magnetic gluons, which is a long-range interaction pis00 ; son99 . The weak-coupling solution (2) implies that this term contributes to the gap equation at order $`O(1)`$. We call this term the leading-order term. The value of the coefficient $`\zeta `$ determines the constant $`c`$ in Eq. (2). The second term in Eq. (1) contains subleading contributions of order $`O(g)`$ to the gap equation, characterized by a single power of the logarithm $`\mathrm{ln}(\mu /\varphi _0)1/g`$. A part of them arises from the exchange of non-static magnetic and static electric gluons pis00 . Another part is due to the contribution from the regular quark self-energy wang ; bro00 . Vertex corrections were reported not to contribute to subleading order bro00 ; mishra . The coefficient $`\beta `$ in Eq. (1) determines the constant $`b`$ in Eq. (2). The third term in Eq. (1) summarizes sub-subleading contributions of order $`O(g^2)`$ with neither a collinear nor a BCS logarithm.

In principle, on the quasi-particle mass shell the gap parameter is an observable quantity, and thus must be independent of the choice of gauge. However, the mean-field approximation to the QCD gap equation may violate this requirement. It was argued in Refs. SchaferWilczek ; pis00 ; ShusterRajagopal  that, in mean-field approximation, gauge-dependent terms enter the QCD gap equation at sub-subleading order. The authors of Ref. pis02  confirmed that this is indeed the case for Coulomb gauge and when taking the gap parameter on the quasi-particle mass shell. In contrast to this result, it was pointed out miransky ; hong03  that, in covariant gauge, the gauge-parameter dependence already shows up at subleading order. Denoting the gauge parameter in general covariant gauges by $`\xi `$, there is then an additional factor $`\mathrm{exp}(3\xi /2)`$ to the prefactor $`b`$ in Eq. (2).

It is a priori not clear why it should depend on the choice of gauge, at which order gauge-parameter dependent terms appear in the QCD gap equation in mean-field approximation. A posteriori, this result is maybe not that puzzling at all, since the mean-field approximation corresponds to a resummation of a particular class of diagrams which may contribute at different orders for different choices of gauge. In any case, in order to remove the gauge-parameter dependence at subleading order in covariant gauge, one has to take into account corrections to the mean-field approximation. Gauge-parameter dependent terms will then only occur at sub-subleading order. The obvious correction so far not taken into account in the mean-field approximation is to replace bare $`qqg`$ vertices by full vertices. It is the purpose of this paper to show how the inclusion of vertex corrections guarantees the gauge-parameter independence at subleading order in covariant gauge.

As a matter of fact, it is not even necessary to compute the vertex corrections explicitly. We shall prove the gauge independence in a rather convenient way by making use of a Ward identity which relates the vertex to the inverse propagator. This approach has been frequently applied to show the gauge independence of physical collective excitations in thermal gauge theories, like hot QCD rebhan ; Pb . However, in a superconductor, the existence of a fermion-fermion condensate necessitates the use of the Nambu-Gor’kov (NG) basis to describe the propagation of quasi-particle excitations. It is therefore desirable to derive the Ward identity in the NG basis. In this paper, we derive this generalized Ward identity for color-superconducting quark matter and apply it to the gap equation, thereby showing that the gap parameter is gauge-independent to subleading order on the quasi-particle mass shell.

We note that our approach is similar in spirit to that of Gerhold and Rebhan gerhold , who used generalized Nielsen identities to give a formal proof that the fermionic quasi-particle dispersion relations in a color superconductor are gauge independent, assuming that the 1PI part of the variation of the effective action with respect to the gauge parameter has no singularities coinciding with those of the quark propagator. Other Ward identities for color-superconducting quark matter have been derived by Miransky, Shovkovy, and Wijewardhana igor , but these identities were not suitable to see the gauge independence of the gap parameter to subleading order.

In order to derive the generalized Ward identity, consider QCD with $`N_f`$ quark flavors and quark chemical potential $`\mu `$. The Lagrangian and the generating functional can be written in the NG basis,

$``$ $`=`$ $`{\displaystyle \frac{1}{2}}\overline{\mathrm{\Psi }}𝒮_0^1\mathrm{\Psi }+{\displaystyle \frac{g}{2}}\overline{\mathrm{\Psi }}𝒯^a\gamma _\mu \mathrm{\Psi }A_a^\mu {\displaystyle \frac{1}{4}}F_{\mu \nu }^aF_a^{\mu \nu }{\displaystyle \frac{1}{2\xi }}(_\mu A_a^\mu )^2,`$ (3)
$`Z[J,\overline{H},H]`$ $`=`$ $`{\displaystyle [dA][d\overline{\mathrm{\Psi }}][d\mathrm{\Psi }]\mathrm{exp}\left\{_X\left[+J_\mu ^aA_a^\mu +\overline{H}\mathrm{\Psi }+\overline{\mathrm{\Psi }}H\right]\right\}}.`$ (4)
The quark fields and their sources are given by

$$\mathrm{\Psi }=\left(\begin{array}{c}\psi \\ \psi _C\end{array}\right),\overline{\mathrm{\Psi }}=(\overline{\psi },\overline{\psi }_C),H=\left(\begin{array}{c}\eta \\ \eta _C\end{array}\right),\overline{H}=(\overline{\eta },\overline{\eta }_C),$$
(5)
where the charge conjugate spinors $`\psi _C,\overline{\psi }_C`$ are defined through $`\psi _C=C\overline{\psi }^T`$, $`\overline{\psi }_C=\psi ^TC`$ with the charge conjugation matrix $`C=i\gamma ^2\gamma _0`$. The inverse free quark propagator (including the chemical potential $`\mu `$) is

$$𝒮_0^1\left(\begin{array}{cc}S_{0,11}^1& 0\\ 0& S_{0,22}^1\end{array}\right),$$
(6)
where $`S_{0,11}^1=i\gamma ^\mu _\mu +\mu \gamma _0m`$, $`S_{0,22}^1=i\gamma ^\mu _\mu \mu \gamma _0m`$. The color part of the $`qqg`$ vertex is

$$𝒯^a\left(\begin{array}{cc}T^a& 0\\ 0& T^{aT}\end{array}\right),$$
(7)
where $`T^a`$ are the generators of $`SU(3)_c`$. The factor $`1/2`$ in front of the first two terms in Eq. (3) accounts for the doubling of the quark degrees of freedom in the NG basis. The gluon field and the field strength tensor are denoted by $`A_\mu ^a`$ and $`F_{\mu \nu }^a`$, respectively. The last term in $``$ is the gauge fixing term, and $`J_\mu ^a`$ is the source term for the gluon field. We denote the space-time integration as $`_X_0^{1/T}𝑑\tau _Vd^3𝐱`$. We neglect ghost terms, as ghosts do not contribute to leading and subleading order in cold, dense quark matter ris03 .

The generating functional $`Z`$ is invariant under the following infinitesimal $`SU(3)_c`$ gauge transformation

$`\delta \mathrm{\Psi }`$ $`=`$ $`i\theta ^a𝒯^a\mathrm{\Psi },`$
$`\delta \overline{\mathrm{\Psi }}`$ $`=`$ $`i\theta ^a\overline{\mathrm{\Psi }}𝒯^a,`$ (8)
$`\delta A_\mu ^a`$ $`=`$ $`{\displaystyle \frac{1}{g}}_\mu \theta ^a+f^{abc}A_\mu ^b\theta ^c.`$
The variation of $`Z`$ with respect to the infinitesimal parameter $`\theta ^a`$ of $`SU(3)_c`$ gauge transformations vanishes, which yields the following identity:

$$\frac{1}{\xi }\left\{^2_\mu A_a^\mu +gf^{abc}\left[\frac{}{}_\mu _\nu D_{bc}^{\mu \nu }(X,X)+A_b^\mu _\mu _\nu A_c^\nu \right]\right\}_\mu J_a^\mu +gf^{abc}J_\mu ^bA_c^\mu +ig\left[\overline{H}𝒯^a\mathrm{\Psi }\overline{\mathrm{\Psi }}𝒯^aH\right]=0,$$
(9)
where the average is taken in the functional sense:

$$F(A,\overline{\mathrm{\Psi }},\mathrm{\Psi })=\frac{1}{Z}[dA][d\overline{\mathrm{\Psi }}][d\mathrm{\Psi }]F(A,\overline{\mathrm{\Psi }},\mathrm{\Psi })\mathrm{exp}\left\{_X\left[+J_\mu ^aA_a^\mu +\overline{H}\mathrm{\Psi }+\overline{\mathrm{\Psi }}H\right]\right\},$$
(10)
and $`D_{bc}^{\mu \nu }=\delta ^2W/\delta J_\mu ^b\delta J_\nu ^c`$ is the two-point Green’s function, where $`W=\mathrm{ln}Z`$ is the generating functional for connected Green’s functions. In the hard-dense-loop (HDL) approximation, the two-point function $`D_{bc}^{\mu \nu }`$ is symmetric in the color indices $`b`$ and $`c`$. Consequently, the term proportional to $`D_{bc}^{\mu \nu }`$ in Eq. (9) vanishes, because the structure constants $`f^{abc}`$ multiplying the two-point function are antisymmetric. In a color superconductor this does not need to be true ris01 , but the violation of symmetry occurs for gluon energies and momenta of order $`\varphi _0`$. In the QCD gap equation, this range of gluon energy and momentum contributes only beyond subleading order ris01 .

We can use the effective action $`\mathrm{\Gamma }[A,\overline{\mathrm{\Psi }},\mathrm{\Psi }]W_X[J_\mu ^aA_a^\mu +\overline{H}\mathrm{\Psi }+\overline{\mathrm{\Psi }}H]`$ to rewrite the identity (9), i.e., we replace $`J_a^\mu `$ by $`\delta \mathrm{\Gamma }/\delta A_\mu ^a`$, $`\overline{H}`$ by $`\delta \mathrm{\Gamma }/\delta \mathrm{\Psi }`$, and $`H`$ by $`\delta \mathrm{\Gamma }/\delta \overline{\mathrm{\Psi }}`$. Note that, from now on, $`A_\mu ^a,\overline{\mathrm{\Psi }}`$, and $`\mathrm{\Psi }`$ denote the expectation values of the gluon and quark fields. Taking the functional derivative $`\delta /\delta \mathrm{\Psi }(X_1)`$ from the right and $`\delta /\delta \overline{\mathrm{\Psi }}(X_2)`$ from the left, and using the fact that the expectation values of the Grassmann-valued quark fields always vanish, $`\overline{\mathrm{\Psi }}=\mathrm{\Psi }0`$, we obtain the following identity

$`{\displaystyle \frac{i}{g}}\left[_\mu ^X\delta ^{ab}gf^{abc}A_\mu ^c(X)\right]{\displaystyle \frac{\delta ^3\mathrm{\Gamma }[A,0,0]}{\delta \overline{\mathrm{\Psi }}(X_1)\delta A_\mu ^b(X)\delta \mathrm{\Psi }(X_2)}}`$ $`=`$ $`\delta ^{(4)}(XX_2){\displaystyle \frac{\delta ^2\mathrm{\Gamma }[A,0,0]}{\delta \overline{\mathrm{\Psi }}(X_1)\delta \mathrm{\Psi }(X_2)}}𝒯^a`$ (11)
$`\delta ^{(4)}(XX_1)𝒯^a{\displaystyle \frac{\delta ^2\mathrm{\Gamma }[A,0,0]}{\delta \overline{\mathrm{\Psi }}(X_1)\delta \mathrm{\Psi }(X_2)}}.`$
The expectation value $`A_\mu ^c`$ of the gluon field on the left-hand side is not necessarily zero in a color superconductor gerhold ; kryjevski ; ris . However, in a two-flavor color superconductor, its value is of order $`\varphi _0^2/(g^2\mu )`$ ris , i.e., much smaller than the derivative $`_\mu ^X`$ which is of the order of the gluon energy, $`p_0\varphi _0`$, or even of the gluon momentum, $`p(g^2\mu ^2\varphi _0)^{1/3}`$ pis00 . Therefore, in the following we may safely neglect the term proportional to $`A_\mu ^c`$ on the left-hand side of Eq. (11). We now rewrite Eq. (11) in momentum space,

$$P_\mu \mathrm{\Gamma }_a^\mu (K,P,K+P)=g[𝒯^a𝒮^1(K+P)𝒮^1(K)𝒯^a],$$
(12)
where $`\mathrm{\Gamma }_a^\mu (K,P,K+P)`$ is the Fourier transform of $`\delta ^3\mathrm{\Gamma }/\delta \overline{\mathrm{\Psi }}(X_1)\delta A_\mu ^a(X)\delta \mathrm{\Psi }(X_2)`$, representing the full $`qqg`$ vertex, while $`𝒮^1(K)`$ is the Fourier transform of $`\delta ^2\mathrm{\Gamma }/\delta \overline{\mathrm{\Psi }}(X_1)\delta \mathrm{\Psi }(X_2)`$, i.e., the full inverse quark propagator. We note that the above identity is a $`2\times 2`$ matrix identity in the NG basis, and thus a generalized Ward identity. Equation (12) is rather similar to the Ward identity derived by Nambu to prove the gauge invariance of the Meissner effect in ordinary superconductors nambu , except that $`𝒯^a`$ assumes the role of the Pauli matrix $`\tau _3`$ in NG space.

The inverse of the quark propagator $`𝒮^1`$ in Eq. (12) is defined as ris03

$$𝒮^1𝒮_0^1+\mathrm{\Sigma },$$
(13)
where $`\mathrm{\Sigma }`$ is the quark self-energy in the NG basis. In momentum space, the Dyson-Schwinger equation for the quark self-energy reads ris03

$$\mathrm{\Sigma }(K)=g_P\mathrm{\Gamma }_\mu ^a𝒮(K+P)𝒯^b\gamma _\nu D_{ab}^{\mu \nu }(P).$$
(14)
The right-hand side of this equation is shown diagrammatically in Fig. 1. Here, $`D_{ab}^{\mu \nu }`$ is the HDL-resummed gluon propagator and $`_PT_nd^3𝐩/(2\pi )^3`$, where the sum runs over all fermionic Matsubara frequencies $`\omega _n=(2n+1)\pi T`$. The first vertex in the above equation is the full one and the second vertex is the bare one.

We now show with the help of the generalized Ward identity (12) that, on the quasi-particle mass shell, the quark self-energy (14) does not depend on the gauge parameter $`\xi `$ in covariant gauge. To this end, consider the contribution from the gauge-dependent part of the gluon propagator to the quark self-energy,

$$\mathrm{\Sigma }_\xi (K)g\xi _P\mathrm{\Gamma }_\mu ^a𝒮(K+P)𝒯^a\gamma _\nu \frac{P^\mu P^\nu }{(P^2)^2}.$$
(15)
Using the generalized Ward identity (12), this becomes

$$\mathrm{\Sigma }_\xi (K)=g^2\xi _P\left[𝒯^a𝒮^1(K+P)𝒮^1(K)𝒯^a\right]𝒮(K+P)𝒯^a\gamma _\nu \frac{P^\nu }{(P^2)^2}.$$
(16)
Now we put the external momentum $`K`$ on the quasi-particle mass shell, which is defined as

$$S^1(K)\mathrm{\Psi }_{\mathrm{on}\mathrm{shell}}=0,\overline{\mathrm{\Psi }}_{\mathrm{on}\mathrm{shell}}S^1(K)=0,$$
(17)
where $`\mathrm{\Psi }_{\mathrm{on}\mathrm{shell}}`$ is the on-shell quasi-particle wave function. Then, when sandwiching Eq. (16) between on-shell wave functions, the second term in brackets in Eq. (16) vanishes. We arrive at

$$\mathrm{\Sigma }_\xi (K)=g^2\xi \frac{4}{3}\mathbf{\hspace{0.17em}1}_P\frac{\gamma _\nu P^\nu }{(P^2)^2}=0,$$
(18)
where $`\mathrm{𝟏}`$ is the $`2\times 2`$ unit matrix in the NG basis. The right-hand side vanishes because the integrand of the $`P`$-integral is an odd function of $`P`$. Hence, we see that the contribution from the gauge-dependent part of the self-energy (14) vanishes on the quasi-particle mass shell. Therefore, the gap parameter as well as the quasi-particle dispersion relation are gauge-independent in a color superconductor.

This result holds strictly speaking only to subleading order in the QCD gap equation. The reason is that the generalized Ward identity (12) is not an exact identity, but only holds under the following approximations. The first approximation is that we neglected the contribution of ghosts. For QCD, which is a non-Abelian gauge theory, in general one has to take into account ghosts and use the BRST transformation and the Slavnov-Taylor identity or, equivalently, generalized Nielsen identities gerhold ; rebhan . However, in dense quark matter, diagrams involving ghosts can be neglected at least to subleading order in the QCD gap equation ris03 . The second approximation is that we neglected the contribution from $`D_{bc}^{\mu \nu }`$ in Eq. (9), arguing that the gluon propagator is symmetric in the color indices. This is true for an HDL-dressed propagator, but not necessarily for the gluon propagator in a color superconductor ris01 . However, terms which violate the symmetry of $`D_{bc}^{\mu \nu }`$ should also be negligible up to subleading order in the gap equation. Finally, the third approximation was to neglect the expectation value of the gluon field in the generalized Ward identity (12). In Ref. ris  it is shown that actually $`A_\mu ^a0`$ in order ensure color neutrality. However, at least in a two-flavor color superconductor, $`A_\mu ^a\varphi _0^2/(g^2\mu )`$, i.e., the expectation value of the gluon field is negligible up to subleading order in the gap equation.

In summary, we derived a generalized Ward identity from QCD for dense, color-superconducting quark matter. The identity implies that, on the quasi-particle mass shell, the gap function and the quasi-particle dispersion relation are independent of the gauge parameter in covariant gauge up to subleading order. There is one potential caveat to this statement. What we have shown is that, to subleading order, the gauge dependence of the quark self-energy (14) arising from the gauge-dependent part of the gluon propagator vanishes on the mass shell. In principle, however, other gauge-dependent terms arise from the gauge dependence of the full vertex $`\mathrm{\Gamma }_\mu ^a`$ when combined with the physical part of the gluon propagator in Eq. (14). In order to show that this is of sub-subleading order in the gap equation, one actually has to compute the gauge-dependent part of the vertex explicitly. This calculation can be simplified considerably with the help of the Ward identity (12), for more details see Ref. hwr . The result is that also this contribution is of sub-subleading order in the gap equation.

Our result shows that in order to obtain a gauge-independent gap function up to subleading order, one has to use the full vertex as well as the full fermion propagator in the NG basis. A consequence is that the prefactor $`\mathrm{exp}(3\xi /2)`$ to the gap function found in the mean-field approximation miransky ; hong03  will be removed by contributions from the full $`qqg`$ vertex when taking the gap function on the quasi-particle mass shell. An explicit diagrammatic proof of this statement will be presented elsewhere hwr .

###### Acknowledgements.

The authors thank R. Pisarski, A. Rebhan, H.-c. Ren, and I. Shovkovy for interesting discussions and a critical reading of the manuscript. D. Hou acknowledges financial support from the Alexander von Humboldt-Foundation and the National Natural Science Foundation of China under grants 1000502 and 10135030. He also appreciates help and support from the Institut für Theoretische Physik of the J.W. Goethe-Universität, and especially from Prof. W. Greiner. The work of Q. Wang is supported by GSI Darmstadt and BMBF.


format: text

text:
# Growth by Random Walker Sampling, and Scaling of the Dielectric Breakdown Model

## I Introduction

The steady state distribution of random walking particles is governed by Laplace’s Equation. As a result Witten and Sander’s model of Diffusion Limited Aggregation (DLA) Witten and Sander (1981), in which a cluster grows by irreversible accretion of dilute diffusing material, has been of double interest: it is readily simulated out to huge cluster sizes Ball and Brady (1985); Tolman and Meakin (1989), whilst at the same time the governance of its growth by Laplace’s equation renders it a landmark mathematical challenge to analyse. The connection to a Laplacian field also underpins the breadth of application of the DLA model, to problems such as viscous fingering in porous media Nittmann et al. (1985) and electrodeposition Matsushita et al. (1984); Brady and Ball (1984). This interest has all been abetted by controversy as to whether the distribution of cluster shapes conforms to simple fractal scaling (see Somfai et al. (1999) and references therein), and interest in the multifractal scaling of the growth measure Amitrano (1989); Ball and Blunt (1990); Jensen et al. (2003); Halsey et al. (1986a).

Physical analogies and the mathematical connections have led to interest in other models where the growth is governed *non-linearly* by a Laplacian field. In particular Niemeyer, Pietronero and Wiesmann introduced the family of Dielectric Breakdown Models (DBM) Niemeyer et al. (1984) where

$$v_n\left|_n\varphi \right|^\eta ,^2\varphi =0,\text{ }\varphi _{\text{interface}}0$$
(1)
and $`\eta `$ is (for interest) a positive parameter. Interest in this is further prompted by Ball and Somfai’s proposal Ball and Somfai (2002, 2003) that growth proportional to field with non-trivial spatial cutoff maps onto simple DBM but with $`\eta 1`$. It is important that diffusion controlled growth which is locally limited by the capillary energy associated with surface ramification is a case in point. Computationally these non-linear models have been a challenge, as random walkers only directly sample the harmonic measure linearly, realising only the $`\eta =1`$ case.

First in this paper we show that random walkers can be exploited to sample non-linear moments of the harmonic measure. In this way we obtain results for the active portion of the multifractal spectrum of DLA far beyond existing results. The key strength of these methods is that no explicit solving of the Laplace equation is involved.

We then exploit this to establish a new method of growing DBM clusters by random walker accretion. This also entails adopting the noise reduction strategies lately introduced in Ref. Ball et al. (2002), and enables us to explore the DBM class out to unprecedentedly large clusters. We show that in two dimensions this largely resolves how the exponents of the DBM model depend on $`\eta `$: in particular superuniversality of the tip scaling exponent alpha is strongly refuted, in favour of a continuous variation of exponents which also confirms the hypothesised upper critical value Hastings (2001a, b); Halsey (2002); Sanchez et al. (1993) $`\eta _c=4`$.

Our walker-DBM results are supported by extensive computations using established iterative conformal mapping methods due to Hastings and Levitov (HL) Hastings and Levitov (1998) and also by direct integration of the Shraiman-Bensimon Equations Shraiman and Bensimon (1984) exploiting the mappings of Ball and Somfai Ball and Somfai (2002, 2003). For exponents, all three agree within statistical errors. Below the level of the errors there is a systematic difference between walker-DBM and HL, and separate results for the relative penetration depth suggest that it is the walker-DBM clusters which are more converged to asymptotic behaviour. Unlike HL and our Shraiman-Bensimon integrations, the new walker-DBM technique is not limited to two dimensions and so the way forward appears open to a full exploration of the DBM class in three dimensions.

## II Random Walker Sampling

We consider first the problem of sampling the harmonic measure of an equipotential surface. The harmonic measure is given on the surface by $`\frac{d\mu }{ds}=\frac{\varphi }{n}`$ where outside the surface the scalar field $`\varphi `$ obeys a Poisson equation $`^2\varphi =S(x)`$ with $`S(x)`$ the source density; typically we will be interested in cases where the source is concentrated at points, particularly $`\mathrm{}`$. We can sample this measure by introducing random walkers at the source points, and tracing their Brownian trajectories to the point of first encounter with the surface, whereupon any given walker is discarded. The points of first encounter uniformly sample the measure $`\mu `$.

By firing a large enough sample of random walkers and collecting frequency counts of their hit distribution, one can build up an approximation to the entire harmonic measure. This has been successfully used by Somfai *et al.* Somfai et al. (1999). However such methods are expensive and give differing quality estimates across the support.

Here we focus on recurrence times defined as follows. We first divide the support into (many) small partitions (hereafter termed “sites”) for each of which we aim to estimate the corresponding hit probability $`\mu `$. We then fire (independent) walkers sequentially at the surface and when each walker hits, the number of walkers fired since the previous hit on that site we will call the “age” $`a`$, and this provides a simple estimate of the hit probability of that site $`\mu _1=1/a`$. This is a standard way to estimate frequency of uncorrelated events, by their recurrence time. The probability distribution of the estimator $`\mu _1`$, given the true underlying value $`\mu `$ for that site, is given by simple Poisson statistics as

$$p_1(a|\mu )=\mu e^{\mu a},$$

assuming that $`\mu 1`$ so that we can approximate age as a continuous variable. We can generate more reliable estimates by using the age $`a_k`$ since the $`k`$’th previous hit, which is distributed according to

$$p_k(a_k|\mu )=\mu e^{\mu a_k}(\mu a_k)^{k1}/(k1)!,$$
(2)
in terms of which our estimate is $`\mu _k=k/a_k`$.

In applications considered below we will be interested in calculating moments of the measure, $`M_q=\mu ^q`$, where the summation is over the sites over which we partitioned the support. This can be interpreted as $`M_q=\mu ^{q1}`$ where the average is weighted by the measure $`\mu `$ itself which in turn is sampled by where random walkers hit. At each impact we can use the site age in our estimate for the factors $`\mu ^{q1}`$, leading to the estimate

$$M_{qk}=A(k,q)(a_k/k)^{1q},$$

where the average is now over all random walkers hitting the surface, we use $`k`$’th ages, and $`A(k,q)`$ is a trivial numerical front factor.

Random walkers being cheap, our main concern is when this estimate converges. Using the distribution (2) one can readily check that $`M_{qk}`$ converges to $`M_q`$ provided $`k>q1`$, and we use the appropriate prefactor $`A(k,q)=\frac{k^{1q}(k1)!}{(kq)!}`$. Thus we can compute moments of finite degree simply by using high enough order ages. In practice we will see that for problems of interest the most significant limitation comes not at high $`q`$ but rather at highly negative $`q`$, for which random walkers give an inefficient sample of the relevant parts of the measure.

## III Moments under Growth

The simple ideas above become rather useful when extended to compute moments of the measure as a surface grows. The harmonic measure of Diffusion Limited Aggregation provides a well studied (but not entirely resolved) test case. For a cluster of $`N`$ added particles the conventional multifractal scaling would lead to $`M_qN^{\tau (q)/D}`$, where $`D`$ is the fractal dimension relating radius $`R`$ to $`N`$ through $`NR^D`$. Summing these moments over $`N`$ with weight $`N^{t1}`$gives us a partition function $`Z(q,t)`$ which we can estimate (ignoring numerical prefactors) as

$$Z(q,t,N)=\underset{n=1}{\overset{N}{}}a^{1q}n^{t1},$$

where the sum is now over the particles used to grow the cluster and the corresponding ages of the sites where they hit. Following the spirit of how Halsey *et al.* Halsey et al. (1986b) generalised the identification of multifractal spectrum, we can now identify that $`\tau (q)/D`$ separates the values $`t<\tau (q)/D`$ for which $`Z(q,t,N)\mathrm{}`$ as $`N\mathrm{}`$ from the values $`t>\tau (q)/D`$ for which $`Z(q,t,N)0`$.

The above definition does not restrict the behaviour of $`Z(q,t,N)`$ on the locus $`t=\tau (q)/D`$ but the simple expectation is of a logarithmic divergence with $`N`$. Then a numerical strategy is to choose $`t`$ such that $`Z(q,t,N)Z(q,t,ϵN)=\underset{n=ϵN+1}{\overset{N}{}}a^{1q}n^{t1}`$ becomes independent of $`N`$ as $`N`$ becomes large with fixed $`ϵ`$.

To obtain results at $`q2`$ we have to use higher order ages. The first order age is naturally thought of as the age of the parent to a given new site and a corresponding estimate of $`a_2`$ is given by the age of its grandparent (the parent of its parent) and so on. It is of course a concern that the more such generations one looks back for the age, the more the cluster geometry will have evolved in the meantime: we will see that overall cluster screening cuts in below a threshold value of $`q`$ in a well understood way. To alleviate the worry above threshold we have studied noise reduced clusters where each particle deposited advances the local geometry by only a small fraction $`A`$ of one notional particle size, so that we can be confident that the local geometry does not evolve significantly on looking back of order $`1/A`$ generations.

Figure 1 shows the moment spectrum $`\tau (q)`$ resulting from imposing the condition $`Z(q,t,N)Z(q,t,ϵN)=Z(q,t,ϵN)Z(q,t,ϵ^2N)`$ on a single off-lattice DLA cluster. As expected from our discussion in section II, the results using ages of order $`k`$ break away at $`q=k+1`$. An important check is given by the measured value $`\tau (3)/D=1.01`$ compared to the prediction of unity from the Electrostatic Scaling Law.

For more definitive results we have studied a sample of $`10`$ clusters out to $`N=10^7`$ using $`A=0.1`$. Even allowing for the noise reduction which means that of order $`1/A=10`$ particles were deposited per local unit of advance in the growth, it is unprecedented to probe the harmonic measure out to such large scales. Because of the volume of age data, it was cumbersome to seek stationarity of $`Z(q,t,N)Z(q,t,ϵN)`$ and instead we extracted $`t(q)`$ from the simple scaling expectation that $`Z(q,t^{},N)Z(q,t^{},ϵN)N^{t^{}t(q)}`$ where we simply chose $`t^{}`$ to assure reasonably uniform weighting of the data. We then took $`\tau (q)=Dt(q)`$ using $`D=1.71`$, where any error in this value is not significant to the accuracy of $`t(q)`$. The resulting $`f(\alpha )`$ spectrum is shown in Fig. 2 for ages of order $`k=1,3,9,27`$ and compared to one earlier report of the $`f(\alpha )`$ spectrum which is distinguished by carrying full error bars Ball and Spivack (1990) and one more recent spectrum Jensen et al. (2002) which claimed better convergence than all previous.

The success of our data lies in the region $`q1`$ corresponding to $`\alpha 1`$, where it passes two important tests by meeting the (dashed) tangent lines corresponding to Makarov’s TheoremMakarov (1985) at $`q=1`$ and (for the higher $`k`$ as appropriate) Halsey’s electrostatic scaling lawHalsey (1987) at $`q=3,`$ and in these respects it clearly improves on the earlier published data.

A significant limitation of the $`f(\alpha )`$ spectrum obtained is that it substantially undershoots at large $`\alpha `$, but this is qualitatively consistent with expectations. It is obvious that we cannot expect to probe probabilities to hit a site which are smaller than $`1/N`$ as these are unlikely to be sampled by the $`N`$ walkers used to grow the cluster. Thus we certainly cannot probe the spectrum for $`\alpha >D`$.

A more careful argument suggests that in principle we should (in the limit of large enough $`N`$) be able to probe all the way up to $`\alpha =D`$, as follows. Consider a typical site $`X`$ “born” with hit probability $`R^\alpha `$, $`\alpha <D`$. It is expected that the probability for growth within distance $`r`$ of this point varies as $`p(r)(r/R)^\alpha `$, and when this vicinity has received $`r^D`$ walkers the local structure up to length scale $`r`$ will have been completely regrown. This requires a number of walkers added to the cluster given by $`\delta N(r/R)^\alpha r^D`$ leading to $`\delta NR^\alpha r^{D\alpha }`$. The argument is consistent for $`\alpha <D`$ because it shows that for length scales larger than $`r`$ the regrowth threshold will not have been met, validating our retention of the hit probability $`p(r)(r/R)^\alpha `$. From the point of view of site $`X`$ significant reorganisation of the cluster happens first on the smallest scales, and requires $`\delta NR^\alpha `$ by which point we can expect to have hit site $`X`$ itself with non-zero probability. For $`\alpha >D`$ the scenario is of screening being dominated by distant growth and hitting the site is unlikely. A more detailed discussion of the screening dynamics is given by Ball and Blunt Ball and Blunt (1990).

Figure 3 shows how the $`f(\alpha )`$ spectrum from first order ages develops as a function of size up to $`N=10^7`$. As expected the spectrum builds up on the RHS as we increase $`N`$ consistent with our expectation that ultimately we can probe as far as $`\alpha =D`$. Quantitatively the convergence is poor beyond about $`\alpha 1.5`$.

## IV Dielectric Breakdown Model

### IV.1 The random walker model

Niemeyer, Pietronero and Wiesmann Niemeyer et al. (1984) introduced the generalisation of diffusion controlled growth in which the local advance velocity is set by $`\mu ^\eta `$. In this model $`\eta =1`$ corresponds to DLA, while growth with higher values of $`\eta `$ was proposed to model the evolution of dielectric breakdown patterns. More recently it has been argued by two of us Ball and Somfai (2002, 2003) that non-trivial values of $`\eta `$ can also be induced when mapping between different types of ultraviolet cutoff, even when the underlying growth is strictly proportional to flux. The computational challenge of the DBM is that random walkers sample the harmonic measure proportional to $`\mu `$, so for $`\eta 1`$ an explicit calculation of the measure appears to be required. Relaxation methods Amitrano (1989); Sanchez et al. (1993) and more recently (for two dimensions) conformal mapping techniques Hastings (2001b) have been used, but the resulting computations are dramatically more expensive than for simple DLA simulation.

We propose that the following random walker growth model is equivalent to DBM. We fire random walkers at the cluster, whose first encounter linearly samples $`\mu `$. At each encounter we use the $`k`$’th order age of the site hit to give the estimate $`\mu _k=k/a_k`$ and accordingly advance the growth locally by an amount proportional to $`\delta m=Aa_k^{1\eta }`$, which is interpreted as the mass added to the cluster. There are two clear constraints, the first is that the local advance must converge in the mean: $`a_k^{1\eta }_\mu \mu ^{\eta 1}`$, requiring that $`k>\eta 1`$. The second is that we desire that the $`k`$’th order ages should look back no further than a particle size, so individual advance steps should be bounded by $`1/k`$. Naïvely this requires $`A<1/k`$ for $`\eta >1`$ and $`A<N^{\eta 1}`$ for $`\eta <1`$, where $`N`$ is the total number of walkers, but we will discuss below how we can greatly improve on these constraints by adjusting $`A`$ dynamically.

### IV.2 DBM with optimal growth steps

We now consider allowing the growth step prefactors $`A`$ to vary as the growth proceeds. It is convenient to interpret $`A=\delta Q`$ as the charge borne by each walker adding to the growth, and the growth is governed by the local advance rate $`\left(r/Q\right)_n=\mu ^\eta `$ where $`Q`$ is the cumulative charge added. It was noted by two of us Ball and Somfai (2003) that if we characterise the extremal tips of the growth as having $`\mu _{\text{tip}}R^{\alpha _{\text{tip}}}`$ that the growth law at the tips trivially leads to $`QR^{1+\eta \alpha _{\text{tip}}}`$ as a generalisation of a standard relation for the fractal dimension of DLA, and less trivially it is expected that $`1+\eta \alpha _{\text{tip}}=\tau (2+\eta )`$ which we will use below to render some expressions less cumbersome. Within all this framework we are now free to choose the charge increments per walker $`\delta Q`$ to vary systematically with growth of a cluster (but not biased by where each particular walker hits), with the corresponding mass increments given by $`\delta m=\delta Qa_k^{1\eta }`$.

We focus on the case $`\eta >1`$ for which the concern lies with anomalously low ages which would give large growth steps. At given $`\mu =R^\alpha `$, the probability that a single sampling of the $`k`$’th age is less than $`a_k`$ can be found from Eq. (2) and varies as $`(R^\alpha a_k)^k`$ for $`\mu a_k1`$. If we grow the cluster by a significant fraction of its size then the number of walkers sampling this distribution is given by $`(Q/\delta Q)R^{f(\alpha )\alpha }`$ where the first factor is the total number of walkers added and the second is the probability of one walker landing on a site with $`\mu R^\alpha `$. Therefore the number of times ages below $`a_k`$ are expected to be sampled is given by

$$\underset{\alpha }{\mathrm{max}}\left[Q/\delta QR^{f(\alpha )\alpha }(R^\alpha a_k)^k\right],$$
(3)
where we have taken the maximally contributing value of $`\alpha `$. To find the smallest likely value of $`a_k`$ we then set this expression to unity. Maximising the $`\alpha `$ dependence gives $`\mathrm{max}_\alpha \left[R^{f(\alpha )\alpha }R^{k\alpha }\right]R^{\tau (k+1)}`$ and it is convenient to similarly substitute $`QR^{\tau (2+\eta )}`$, leading to expected minimum age given by

$$\mathrm{min}\{a_k\}\left(R^{\tau (k+1)\tau (2+\eta )}\delta Q\right)^{1/k}.$$
(4a)
The above calculation is valid when it yields an increasing function of $`R`$; otherwise the scaling assumption in Eq. (3) does not hold, and we have
$$\mathrm{min}\{a_k\}=1.$$
(4b)
We have tested the predicted minimum age by looking at the extremal ages sampled during DLA cluster growth, where we have $`\delta Q=1`$, $`\eta =1`$ and $`\tau (3)=D`$. Figure 4 shows scatter plots of all the ages as a function of cluster growth in terms of $`NR^D`$, using ages of different orders and comparing with the result above.

We can now exploit the expected minimum age to set an optimal choice of $`\delta Q`$ in DBM growth: we substitute $`\mathrm{min}\{a_k\}`$ into the restriction that the growth steps be suitably limited $`\delta m1`$, leading to

$$\delta QR^{(\eta 1)\frac{\tau (k+1)\tau (2+\eta )}{k+1\eta }}.$$
(5)
Given that we are restricted to $`k+1>\eta `$, the exponent of $`R`$ in the expression for $`\delta Q`$ increases monotonically with $`k`$ towards asymptote $`(\eta 1)\alpha _{\text{min}}`$ as $`k\mathrm{}`$.

In practice we are inhibited from setting $`k`$ too large because $`k`$-dependent prefactors compete with optimising the power law for finite $`R`$, but taking $`k=10`$ puts us quite close to the limit. Also, we cannot use Eq. (5) explicitly because the values of $`\tau (k+1)`$ and $`\tau (2+\eta )`$ are not known a priori for a generic $`\eta `$. Instead, armed with the knowledge that such an asymptotic power law exists, we use the empirical formula $`\delta Q=A_0k^{\eta 2}N^\beta `$ for which the parameters are chosen such that the constraints are met. For the parameter values of Table 1, used in the simulations below, the condition $`\delta m<1`$ always holds whilst the distance looked back by the $`k`$’th order age only exceeds one particle diameter a few times per million walkers and never exceeds two particle diameters.

Our algorithm is most competitive for $`\eta `$ near 1, including the important region $`1<\eta <2`$ corresponding to surface tension regularization Ball and Somfai (2002, 2003), where it can generate cluster masses inaccessible by other methods. For too small or large $`\eta `$ (particularly $`\eta 3`$) it becomes less advantageous.

In principle for $`k\mathrm{}`$ the number of walkers required to grow a cluster is given by

$$N=\frac{dQ}{\delta Q}R^{1+\alpha _{\text{min}}+\eta (\alpha _{\text{tip}}\alpha _{\text{min}})}.$$

The result is familiar in that for $`\eta =1`$ it recovers the DLA result (where $`N`$ is also the cluster mass). There has been recent controversy about the precise value of $`\alpha _{\text{tip}}\alpha _{\text{min}}`$ through which our method is explicitly sensitive to $`\eta `$ Jensen et al. (2003), although certainly this difference is very small at $`\eta =1`$, where $`0\alpha _{\text{tip}}\alpha _{\text{min}}0.03\pm 0.005`$, see Fig. 5.

### IV.3 Numerical results

As a test of the new walker-DBM method, we compared the relative penetration depth (the growth zone width, normalized by the average deposition radius) with measurements on clusters grown by the more established HL method, for $`\eta =2`$. The asymptotic behavior of the penetration depth can be considered a sensitive measurement, as for DLA it caused considerable controversy in the past. As seen on Fig. 6, both methods converge to the same asymptotic value of the relative penetration depth, and of the two the walker-DBM converges slightly faster. This can be understood in terms of the higher effective noise reduction levels.

To probe the scaling, we measured the tip scaling exponent $`\alpha _{\text{tip}}`$ and the fractal dimension $`D`$ for several values of $`\eta `$. On Fig. 7 we compare these results with measurements obtained with different methods: HL and by integrating the Shraiman-Bensimon equation; the latter described in more detail in the next section. The data agree within the stated uncertainties, showing a smooth monotonic transition between the extremes $`\eta =0`$ and $`\eta \eta _c`$, while below the level of significance point-by-point one can see some differences. One difference is that for HL, $`\alpha _{\text{tip}}`$ shows very little variation between $`1\eta 2`$ (seemingly agreeing with an earlier prediction Ball and Somfai (2002, 2003)), while the walker-DBM does not show this feature. Another difference is that for high $`\eta `$ the walker-DBM yields lower values for both $`\alpha _{\text{tip}}`$ and $`D`$ than HL. Although for walker-DBM it was not practical to obtain values for $`\eta >3`$, its extrapolation is more consistent with $`\eta _c=4`$.

## V Shraiman-Bensimon Equation

Shraiman and Bensimon observed Shraiman and Bensimon (1984) that for DLA in two dimensions the evolution of the complex potential reduces to a non-local problem in one space dimension (plus time). This is based on the cluster boundary $`z(\theta )`$ expressed as a complex function of the imaginary part $`\theta `$ of the complex potential, corresponding to cumulative harmonic measure or more loosely, charge. In Refs. Ball and Somfai (2002, 2003) we adapted this to the DBM class, arguing that changing from a fixed cutoff scale in physical space to a fixed cutoff in charge space could be offset by adjustment of the DBM parameter from $`\eta _0`$ to $`\eta _1=\alpha \eta _0.`$ Then in terms of

$$(iz/\theta )^{(1+\eta _1)/2}=\psi (\theta ,t)=1+\underset{k=1}{\overset{K}{}}\psi _k(t)e^{ik\theta }$$

the dynamics for DBM growth along a channel of width $`2\pi `$ with periodic boundary conditions reduces to

$`{\displaystyle \frac{d}{dt}}\psi _k=k\psi _k{\displaystyle \underset{m=0}{\overset{K}{}}}\psi _m\overline{\psi _m}`$ (6)
$`+{\displaystyle \underset{j=1}{\overset{k}{}}}{\displaystyle \underset{m=0}{\overset{Kj}{}}}((3+\eta _1)j2k)\psi _{kj}\psi _{j+m}\overline{\psi _m}`$
with $`\psi _0=1`$. The fixed bandwidth limit $`kK`$ corresponds qualitatively to imposing a cutoff of fixed minimal charge on the growing tips, and inevitably means that the most screened regions of the physical growth are not tracked. Scaling arguments lead to the expectation that $`\psi _m\overline{\psi _m}m^{1+\eta _1(11/\alpha )}`$ for $`1mK`$, but give no indication of how wide a range of $`m`$ is required to see the scaling behaviour and hence determine the tip scaling exponent $`\alpha `$, and we will see that this is a major issue below.

The trilinear form of the RHS is efficiently evaluated by Fourier methods costing of order $`K\mathrm{ln}K`$ per timestep for the whole system, and in the results presented here we used $`K=1023`$ and a simple predictor-corrector timestepping scheme. We set the timestep adaptively such that for every $`\psi _k`$ at each timestep: either the predicted and corrected updates $`\delta \psi _k`$ agreed within 10% or else $`\psi _k`$ was being updated by less than 20%. As predictor-corrector is a second order method, the resulting worst case error is of order 5%: even this sounds generous but as it was applied to the worst case of $`1023`$ variables the typical precision achieved was very much higher, and we chose 20% on the basis of obtaining results statistically indistinguishable under refinement.

Figure 8(a) shows the observed behaviour of $`\left(_{j<k}\left|\psi _j\right|^2\right)^{1/\eta _1}`$ vs $`k`$ with increasing time for a representative value of $`\eta _1`$. On logarithmic scales the expectation is a straight line with slope $`11/\alpha `$ and a difficulty is immediately apparent, in that two slope behaviour is clearly observed. It is clear that it is the slope emerging from lower $`k`$ values which predominates at large times, although different slope is preserved over a limited upper range of $`k`$. Figure 8(b) shows results at long times for various $`\eta _1`$ from which it is clear that whilst the high $`k`$ slopes are insensitive to $`\eta `$, the (we claim) asymptotic slopes at lower $`k`$ exhibit a systematic variation.

From the present data the conclusion has to be that the slopes in the lower range of $`k`$ give our best estimate of $`11/\alpha `$ and it is these results which are compared (and agree) with the results by other methods in Fig. 7(a). The relatively constant slopes at high $`k`$ are what dominated our earlier numerical conclusions in Ball and Somfai (2002, 2003) which had a factor of 10 less range in $`k`$, and appeared to lend support to the claim that the slope, and hence $`\alpha `$, might be independent of $`\eta `$. That earlier suggestion is clearly refuted, although we cannot rule out more surprises from simulations at decades larger $`K`$.

## VI Conclusions

We have shown that recurrence times between random walkers provide estimates of the local harmonic measure which are highly effective for estimating moments and for simulating non-linear growth models.

For the scaling of multifractal moments this technique is limited to the regions of the growth which are sufficiently active that nearby hits dominate the screening of a given site, but in that “active range” our results for DLA in two dimensions clearly surpass the quality of previous data and offer more conclusive support to the theoretical exponent relations of Makarov and Halsey. Given that our method is also applicable in higher dimensions, we believe we have established it as the technique of choice to determine harmonic measure in the active range.

Our walker-based simulations of the Dielectric Breakdown Model have greatly consolidated knowledge of how DBM exponents vary with the nonlinearity parameter $`\eta `$ in two dimensions. In particular we can now be much more confident of the continuous variation of fractal dimension with $`\eta `$, looking at the trend in terms of $`\alpha _{\text{tip}}`$ sharpens the issue, and we can be fairly confident that DBM becomes trivial for $`\eta >\eta _c=4`$ as conjectured theoretically Hastings (2001a, b); Halsey (2002); Sanchez et al. (1993).

The clearest results for the scaling of the DBM model in two dimensions over the full range of $`\eta `$ now come, ironically, from our integration of the Shraiman-Bensimon equation. However unravelling the confusing two slope behaviour which this approach exhibits for the extraction of the exponent $`\alpha _{\text{tip}}`$, would not have been prompted without the walker-DBM results. Overall it is the agreement of both methods with earlier results which establish the definitive picture.

Finally we note that the most crucial role of walker-based DBM will be in three dimensions, where the techniques based on complex analysis have nothing to offer and relatively little is known about the exponent behaviour. We look forward to exploring this in a subsequent paper.

###### Acknowledgements.

This research has been supported by the EC under Contract No. HPMF-CT-2000-00800. The computing facilities were provided by the Centre for Scientific Computing of the University of Warwick, with support from the JREI. We thank Joachim Mathiesen and Anders Levermann for providing us with the data of Ref. Jensen et al. (2002).


text:
# Does the detection of X-ray emission from SN1998bw support its association with GRB980425?

## 1. Introduction

Recent Chandra observations of SN1998bw (Kouveliotou et al., 2004) have allowed to identify one of the X-ray sources detected earlier by BeppoSAX (Pian et al., 2000) as associated with the supernova. We discuss here the implications of the detected X-ray emission, in particular to the suggested association of SN1998bw with GRB980425.

The association of GRBs with type Ib/c supernovae is motivated by the coincidence of GRB980425 and SN1998bw (Galama et al., 1998), and by the identification of a SN1998bw-like spectrum in the optical afterglow of GRB030329111See, however, Katz (1994), who suggests that a SN like emission may result from the impact of the relativistic GRB debris on a nearby dense gas cloud. (Stanek et al., 2003; Hjorth et al., 2003). The $`\gamma `$-ray luminosity $`L_{\gamma ,\mathrm{Iso}.}10^{51}\mathrm{erg}/\mathrm{s}`$ of GRB030329, inferred from the redshift z=0.1685 of its host galaxy, lies within the range of typical cosmological GRBs, $`L_{\gamma ,\mathrm{Iso}.}10^{52\pm 1}\mathrm{erg}/\mathrm{s}`$ (e.g. Schmidt, 2001). The subscript ”Iso.” indicates luminosity derived assuming isotropic emission. The association of GRB980425 with SN1998bw sets the distance to this burst to 38 Mpc (for $`H_0=65\mathrm{k}\mathrm{m}/\mathrm{s}\mathrm{Mpc}`$), implying that its luminosity is nearly 5 orders of magnitude lower than that typical for cosmological GRBs (Pian et al., 2000).

Two hypotheses are commonly discussed, that may account for the orders of magnitude difference in luminosity. First, it may be that SNe Ib/c produce two different classes of GRBs, with two different characteristic luminosities. It is now commonly believed that long duration, $`T>2`$ s, cosmological, $`L_{\gamma ,\mathrm{Iso}.}10^{52}\mathrm{erg}/\mathrm{s}`$, GRBs are produced by the collapse of SN Ib/c progenitor stars. It is assumed that the stellar core collapses to a black-hole, which accretes mass over a long period, $`T`$, driving a relativistic jet that penetrates the mantle/envelope and then produces the observed GRB (Woosley, 1993; Paczyński, 1998; MacFadyen & Woosley, 1999). This scenario is supported by the association of GRB030329 with SN2003dh, and by additional evidence for optical supernovae emission in several GRB afterglows (Bloom, 2003). The origin of a second, low-luminosity, class is unknown. It may be, e.g., due to supernova shock break-out (Colgate, 1968; Woosley & Weaver, 1986; Matzner & McKee, 1999; Tan et al., 2001).

A second possibility is that GRB980425 was a typical, cosmological GRB jet viewed off-axis (Nakamura, 1998; Eichler & Levinson, 1999; Woosley, Eastman & Schmidt, 1999; Granot et al., 2002; Yamazaki, Yonetoku & Nakamura, 2003). Due to the relativistic expansion of jet plasma, with Lorentz factor $`\mathrm{\Gamma }100`$ during $`\gamma `$-ray emission (Krolik & Pier, 1991; Baring, 1993), $`\gamma `$-rays are concentrated into a cone of opening angle comparable to the jet opening angle $`\theta _j`$ (assuming $`\theta _j>1/\mathrm{\Gamma }`$). Thus, if the jet is viewed from a direction making an angle larger than $`\theta _j`$ with the jet axis, the $`\gamma `$-ray flux may be strongly suppressed. In this scenario, strong radio emission, $`L_\nu 10^{30}\nu _{10\mathrm{G}\mathrm{H}\mathrm{z}}^{1/2}\mathrm{erg}/\mathrm{s}\mathrm{Hz}`$, is expected at $`1`$ yr delay (Frail, Waxman & Kulkarni, 2000; Livio & Waxman, 2000) as the jet decelerates to sub-relativistic speed and its emission becomes nearly isotropic (for a more detailed discussion see Waxman, 2004, and references therein).

We have recently shown (Waxman 2004) that the low radio luminosity of SN1998bw, compared to that expected from a decelerated GRB jet at $`1`$ yr delay, may be consistent with the off-axis jet interpretation of GRB980425 provided that either the magnetic field energy fraction behind the shock driven by the jet into the wind is atypically low, $`ϵ_B10^4`$, or that the density of the wind surrounding the progenitor is lower than typically expected, $`\dot{m}(\dot{M}/10^5M_{}\mathrm{yr}^1)/(v_\mathrm{w}/10^3\mathrm{km}\mathrm{s}^1)0.1`$. Lower values of $`\dot{m}`$ and of $`ϵ_B`$ reduce the specific luminosity at the transition to sub-relativistic expansion. A low value of $`\dot{m}1`$ further delays the time at which the flow approaches sub-relativistic expansion. We consider the low $`ϵ_B`$ scenario less likely, since we expect $`ϵ_B`$, which is determined by the shock micro-physics, to be similar for different bursts, and $`ϵ_B`$ close to equipartition is inferred from radio observations of other bursts. The latter scenario, $`\dot{m}0.1`$, is consistent with the constraints imposed on $`\dot{m}`$ by the observed radio emission from the supernova (Kulkarni et al., 1998), which is interpreted as due to synchrotron emission from a shock driven into the wind by a mildly relativistic shell ejected by the supernova explosion (Kulkarni et al., 1998; Waxman & Loeb, 1999; Chevalier & Li, 1999). In this scenario, transition to sub-relativistic expansion is expected over $`10`$ yr time scale (Waxman, 2004).

The recent identification of the X-ray emission of SN1998bw was argued to support the interpretation that GRB980425 was an off-axis ”standard” GRB jet associated with SN1998bw (Kouveliotou et al., 2004). We show here that the observed X-ray emission may be naturally explained as synchrotron emission from a shock driven into the wind surrounding SN1998bw progenitor by the mildly relativistic shell that was inferred to exist from earlier radio observations. In §2 we construct a simple model that describes the emission from a shock driven into the wind by a sub-relativistic shell, based on the analysis of Waxman (2004). In §3 we show that the model can account for both X-ray and radio observations. While X-ray observations, or radio observations, alone do not allow to determine all model parameters, combined X-ray and radio data over constrain the model, and therefore enhance our confidence in its validity. Our conclusions are summarized in §4.

## 2. Model

### 2.1. Dynamics

Let us consider a sub-relativistic shell ejected by the supernova explosion, with mass $`M`$, total energy $`E`$ and initial velocity $`v=\beta c`$. We assume the shell’s thermal energy is dominated by its kinetic energy, since the thermal energy deposited in the shell by the passage of the supernova shock leads to expansion and post-shock acceleration which converts the thermal energy to kinetic energy (As pointed out by Matzner & McKee, 1999, this post-shock acceleration is essential for achieving mildly relativistic shell velocity). The shell is assumed to propagate into a $`r^2`$ density profile created by stellar mass loss. We assume a density profile

$$\rho =Kr^2,K\frac{\dot{M}}{4\pi v_\mathrm{w}},$$
(1)
where $`\dot{M}`$ is the mass loss rate and $`v_\mathrm{w}`$ is the wind velocity.

As the shell expands, it drives a shock wave into the surrounding wind, and a reverse shock is driven backward (in the shell frame) into the shell. It is straightforward to show that as long as the mass of shocked wind plasma, $`M_w(r)=4\pi Kr`$, is small compared to the shell’s mass, the shock driven backward into the shell does not lead to significant deceleration. Denoting by $`u<v`$ the velocity of shocked wind (and shell) plasma, pressure balance between the forward and reverse shocks implies $`\rho _wu^2=\rho _s(vu)^2`$, where $`\rho _w=M_w/4\pi r^3`$ and $`\rho _s=\eta _sM/4\pi r^3`$ denote wind and shell density respectively. Here $`r/\eta _s`$ is the thickness of the ejected shell. For $`M_wM`$ we find $`(vu)/v=(M_w/\eta _sM)^{1/2}`$. Thus, at small radii, where $`M_wM`$, the shell expands with a time independent velocity, $`v=\beta c`$. At this stage, the velocity of the shock propagating into the wind is also time independent and given by $`(\widehat{\gamma }+1)v/2`$, where $`\widehat{\gamma }`$ is the adiabatic index of the wind plasma.

Significant deceleration of the shell begins at a radius where the shocked wind mass, $`M_w(r)`$, is similar to the shell’s mass $`M`$. At larger radii the flow approaches the Sedov-von Neumann-Taylor solutions describing expansion of a spherical strong shock wave into a $`r^2`$ density profile (e.g. Chapter XII of Zel’dovich & Raizer 2002). In these solutions, the shock radius is given by

$$R=\xi (\widehat{\gamma })\left(\frac{E}{K}t^2\right)^{1/3}.$$
(2)
$`\xi (\widehat{\gamma })`$ is a dimensionless parameter of order unity, which may be approximated by $`\xi (\widehat{\gamma })=(3/2)[(\widehat{\gamma }+1)^2(\widehat{\gamma }1)/2\pi (7\widehat{\gamma }5)]^{1/3}`$ (Waxman, 2004). For $`\widehat{\gamma }=5/3`$, appropriate for sub-relativistic flow, we have $`\xi =0.73`$. At the radius where $`M_w(r)M`$, the velocity given by the Sedov-von Neumann-Taylor for the appropriate energy $`E`$ is similar to the initial shell velocity $`v`$. We therefore define the deceleration radius, $`R_{\mathrm{dec}}`$, as the radius at which the Sedov-von Neumann-Taylor post-shock fluid velocity equals $`v`$, $`2\dot{R}/(\widehat{\gamma }+1)=v`$,

$$R_{\mathrm{dec}}=\frac{16\pi }{9}\xi ^3\frac{E/[(\widehat{\gamma }+1)v/2]^2}{\dot{M}/v_w}=2.1\times 10^{15}\frac{E_{49}}{\beta ^2\dot{m}}\mathrm{cm},$$
(3)
where $`E=10^{49}E_{49}`$ erg. The deceleration time, $`t_{\mathrm{dec}}`$, is defined as

$$t_{\mathrm{dec}}\frac{R_{\mathrm{dec}}}{(\widehat{\gamma }+1)v/2}=0.62\frac{E_{49}}{\beta ^3\dot{m}}\mathrm{day}.$$
(4)
For $`tt_{\mathrm{dec}}`$ we therefore expect shell expansion at constant velocity, $`v`$, and a shock driven into the wind at a constant speed, $`(\widehat{\gamma }+1)v/2`$. For $`tt_{\mathrm{dec}}`$, we expect the expansion to follow the Sedov-von Neumann-Taylor solution, Eq. (2), with shock velocity $`\dot{R}t^{1/3}`$.

### 2.2. Synchrotron emission

Let us consider synchrotron emission from electrons accelerated to relativistic energies by the collisionless shock driven into the wind. We assume that a constant fraction $`ϵ_e`$ ($`ϵ_B`$) of the post-shock thermal energy is carried by relativistic electrons (magnetic field), and that the electron distribution function follows a power-law, $`\mathrm{d}\mathrm{ln}n_e/\mathrm{d}\mathrm{ln}\gamma _e=p`$ for $`\gamma _e\gamma _m`$. In what follows, we assume $`p=2`$, as observed for non-relativistic (Blandford & Eichler, 1987) as well as for relativistic shocks (Waxman, 1997). This power law is produced by Fermi acceleration in collisionless shocks (Blandford & Eichler, 1987; Bednarz & Ostrowski, 1998; Achterberg et al., 2001), although a first principles understanding of the process is not yet available.

At times $`tt_{\mathrm{dec}}`$ the flow is well described by the spherical non-relativistic self-similar solution, where the shock radius is given by Eq. (2). For simplicity, we assume that at this stage the shocked plasma is concentrated into a thin shell behind the shock, $`R/\eta R,`$ within which the plasma conditions are uniform (the Chernyi-Kompaneets approximation). This implies, in particular, $`\eta =(\widehat{\gamma }+1)/(\widehat{\gamma }1)`$. The synchrotron emission obtained under the above assumptions was derived in (Waxman, 2004). The scaling of magnetic field amplitude $`B`$ and of $`\gamma _m`$ with time, at $`tt_{\mathrm{dec}}`$, is (Waxman, 2004)

$$Bt^1,\gamma _mt^{2/3},$$
(5)
and for $`p=2`$ the specific luminosity $`L_\nu `$ is

$`\nu L_\nu `$ $``$ $`10^{42}(3ϵ_e)(3ϵ_B)^{3/4}\dot{m}^{3/4}E_{49}`$ (6)
$`\times `$ $`\left({\displaystyle \frac{h\nu }{1\mathrm{k}\mathrm{e}\mathrm{V}}}\right)^{1/2}\left({\displaystyle \frac{t}{100\mathrm{d}\mathrm{a}\mathrm{y}}}\right)^{3/2}\mathrm{erg}/\mathrm{s}.`$
Eq. (6) holds at frequencies well above the characteristic synchrotron frequency of the lowest energy electrons ($`\gamma _e=\gamma _m`$), and below the cooling frequency, $`\nu _c`$. $`\nu _c`$, the characteristic synchrotron frequency of electrons for which the synchrotron cooling time is comparable to the shock expansion time, $`t`$, is

$$\nu _c2\times 10^{13}(3ϵ_B)^{3/2}\dot{m}^{3/2}\left(\frac{t}{100\mathrm{d}\mathrm{a}\mathrm{y}}\right)\mathrm{Hz}.$$
(7)
Let us now consider synchrotron emission from the shock driven into the wind at times $`tt_{\mathrm{dec}}`$. At frequencies $`\nu >\nu _c`$, the specific luminosity (for $`p=2`$) is given by $`L_\nu =L_m(\nu _c/\nu _m)^{1/2}(\nu /\nu _c)^1=L_m(\nu _c\nu _m)^{1/2}/\nu `$. Here, $`L_m`$ is the specific intensity at $`\nu =\nu _m`$, the characteristic synchrotron frequency of the lowest energy electrons ($`\gamma _e=\gamma _m`$). In order to determine the time dependence of $`L_\nu `$, we first note that $`L_m`$ is proportional to the product of the number of radiating electrons, $`N_e`$, and the magnetic field strength, $`B`$, $`L_mN_eB`$, and that $`\nu _m\gamma _m^2B`$. The Lorentz factor $`\gamma _c`$ of electrons that cool on a time $`t`$ is given by comparing the synchrotron cooling time, $`t_{\mathrm{syn}}1/\gamma _cB^2`$ to $`t`$, yielding $`\gamma _c1/tB^2`$ and $`\nu _c\gamma _c^2Bt^2B^3`$. The dependence of $`\nu _m`$, $`\nu _c`$, and $`L_m`$ on time is now obtained by noting that the time independent shock velocity implies time independent thermal energy per shocked particle, i.e. $`\gamma _mt^0`$, thermal energy density $`r^2t^2`$ (since the density drops as $`1/r^2`$), i.e. $`Bt^1`$, and $`N_et`$. Thus, $`L_mt^0`$, $`\nu _m\gamma _m^2Bt^1`$, $`\nu _ct`$ and $`L_\nu t^0`$ for $`\nu >\nu _c`$. Synchrotron emission at times $`tt_{\mathrm{dec}}`$ is therefore described by

$$\nu L_\nu \{\begin{array}{cc}\nu ^{1/2}t^{1/2},\hfill & \nu <\nu _c(t)\text{;}\hfill \\ \nu ^0t^0,\hfill & \nu >\nu _c(t)\text{,}\hfill \end{array}$$
(8)
with $`\nu _c`$ given by Eq. (7). Combined with Eqs. (4) and (6), this determines the synchrotron emission predicted by the model at $`t<t_{\mathrm{dec}}`$.

The following line of arguments shows that synchrotron emission from the shock driven backward into the expanding shell, which exists at $`t<t_{\mathrm{dec}}`$, is dominated by the emission from the shock driven into the wind, under the assumption that $`ϵ_e`$ and $`ϵ_B`$ are similar for the post-shock plasma behind the two shocks. Since the thermal energy density is the same in the shocked wind and shell plasmas (pressure equilibrium), the magnetic field is the same in the two regions (assuming similar $`ϵ_B`$ for both shocks) and the product of the electron number density and their minimum Lorentz factor, $`n_e\gamma _m`$, which is proportional to the thermal energy density, is the same in the shocked wind and shell plasmas. Similar $`B`$ implies similar $`\nu _c`$ for both shocks, and hence the ratio of the specific synchrotron luminosity of the shocked shell plasma and shocked wind plasma, $`L_{\nu ,s}/L_\nu `$, is $`L_{\nu ,s}/L_\nu =L_{m,s}\nu _{m,s}^{1/2}/L_m\nu _m^{1/2}=N_{e,s}\gamma _{m,s}/N_e\gamma _m=\mathrm{\Delta }_s/\mathrm{\Delta }`$, where the subscript $`s`$ denotes quantities related to the shocked shell plasma, and $`\mathrm{\Delta }`$ ($`\mathrm{\Delta }_s`$) is the width of shock wind (shell) plasma. The thickness of the region of shocked wind/shell plasma is proportional to the velocity of the shock propagating into the wind/shell plasma, measured in the shocked plasma rest frame, $`\mathrm{\Delta }_s/\mathrm{\Delta }=(vu)/u`$ \[see discussion following Eq. (1)\]. Thus, for $`M_wM`$ we find $`L_{\nu ,s}/L_\nu =(M_w/\eta _sM)^{1/2}1`$.

The following point should be emphasized here. The model described above is highly simplified. For example, the wind density is assumed to follow a pure $`1/r^2`$ law and to be free of inhomogeneities, and the expanding shell is assumed uniform with no internal (density, velocity) gradients. In comparing model predictions with observations one should therefore expect an approximate, rather than exact, agreement. Construction of a more detailed, and more accurate, model requires a more detailed knowledge of, for example, the wind density distribution and its homogeneities and of the internal structure of the shell. We find that such a more detailed analysis is not warranted by the present data, since the number of free parameters in a more detailed model would be larger than the observational constraints.

## 3. Implications to SN1998bw

The luminosity of the X-ray source associated with SN1998bw, given in Fig. 3 of Kouveliotou et al. (2004), may be approximately described as

$$L_X6\times 10^{40}\mathrm{erg}/\mathrm{s}\times \{\begin{array}{cc}1,\hfill & t<\text{110 day;}\hfill \\ (t/110\mathrm{d}\mathrm{a}\mathrm{y})^{3/2},\hfill & \text{otherwise.}\hfill \end{array}$$
(9)
The $`t^{3/2}`$ decline at late times is consistent with the prediction of Eq. (6), describing synchrotron emission from the shock driven into the wind at the self-similar stage of expansion, $`t>t_{\mathrm{dec}}`$ \[see Eq. (4)\], at frequencies lower than the cooling frequency, given by Eq. (7). The time independent flux at earlier time is consistent with synchrotron emission from the shock driven into the wind at the stage of expansion at constant speed, $`t<t_{\mathrm{dec}}`$, at frequencies higher than the cooling frequency \[see Eq. (8)\]. The observed X-ray luminosity may therefore be interpreted as due to the shock driven by a sub-relativistic shell into the wind, for which deceleration starts at $`t_{\mathrm{dec}}110`$ day and the cooling frequency passes through the X-ray band at $`tt_{\mathrm{dec}}`$, i.e. $`h\nu _c(t=110\mathrm{day})1`$ keV.

Using Eq. (7), the requirement $`h\nu _c(t=110\mathrm{day})1`$ keV implies

$$(3ϵ_B\dot{m})^{3/4}10^2.$$
(10)
Using this result in Eq. (6), the observed luminosity at $`t=t_{\mathrm{dec}}110`$ implies

$$3ϵ_eE_{49}6.$$
(11)
Using this result in Eq. (4), the requirement $`t_{\mathrm{dec}}110`$ day implies

$$\beta ^3\dot{m}0.03/(3ϵ_e).$$
(12)
Our model for the X-ray emission of SN1998bw is consistent with its radio emission, which is also interpreted as synchrotron emission from a shock wave driven into the wind by a mildly relativistic shell ejected by the supernova explosion. Within the framework of this model, radio observations imply electron energy fraction close to equipartition, $`3ϵ_e1`$ (see Eq. (17) of Waxman & Loeb, 1999). Near equipartition value of $`ϵ_e`$ is also inferred from radio observations of GRB970508 for the shock driven by the fireball into surrounding plasma at late stages, where the shock becomes mildly relativistic (Frail, Waxman & Kulkarni, 2000). We therefore adopt $`3ϵ_e1`$, which implies, using Eq. (11), $`E_{49}6`$.

Radio observations further allow to determine the values of $`ϵ_B`$, $`\beta `$ and $`E`$ as functions of $`\dot{m}`$. The lowest allowed value of $`\dot{m}`$ (and of $`E`$) is that corresponding to magnetic field equipartition, with larger values of $`\dot{m}`$ implying lower values of $`ϵ_B`$ (and larger values of $`E`$). Model predictions were shown (Kulkarni et al., 1998; Waxman & Loeb, 1999; Chevalier & Li, 1999) to be consistent with radio observations for, e.g., {$`\dot{m}=0.04`$, $`ϵ_B0.1`$, $`E_{49}3`$, $`\beta 0.8`$}, and for {$`\dot{m}=6`$, $`E_{49}50`$, $`ϵ_B10^6`$, $`\beta 0.6`$}. The energy inferred from X-ray observations, Eq. (11), implies that the former choice of parameters, {$`\dot{m}=0.04`$, $`ϵ_B0.1`$, $`E_{49}3`$} with $`\beta 0.8`$, should be chosen for consistency with observations.

X-ray and radio observations therefore over constrain the model: All model parameters are determined from radio observations with the additional X-ray constraint given by Eq. (11). The inferred parameter values are consistent with the two additional independent X-ray constraints, Eqs. (10) and (12).

## 4. Discussion

Simple analytic expressions are given in § 2 for the specific luminosity emitted by a shock wave driven into a wind surrounding a supernova progenitor by a sub-relativistic shell ejected by the supernova explosion. We have shown in § 3 that the observed X-ray and radio emission from SN1998bw are consistent with such a model. In particular, the approximately time independent X-ray flux at early time, $`t<100`$ day, and the $`t^{3/2}`$ decline of the flux at late time \[see Eq. (9)\] are characteristic for the emission from a shock driven into a wind by a shell which suffers significant deceleration at $`t100`$ day \[see Eqs. (6), (8)\]. Combined X-ray and radio data allow to determine all model parameters: The shell kinetic energy is $`E10^{49.7}`$ erg and its initial velocity is $`\beta 0.8`$; The normalized mass loss rate of the wind, $`\dot{m}(\dot{M}/10^5M_{}\mathrm{yr}^1)/(v_\mathrm{w}/10^3\mathrm{km}\mathrm{s}^1)`$, is $`\dot{m}0.04`$; The post-shock magnetic field is not far below equipartition, $`ϵ_B0.1`$. Combined X-ray and radio data over constrain the model, and therefore enhance our confidence in its validity. We note, that the inferred energy in mildly relativistic ejecta, $`\mathrm{\Gamma }\beta >1`$ where $`\mathrm{\Gamma }(1\beta ^2)^{1/2}`$, is an order of magnitude higher than the value obtained from analyses of shock emergence from the assumed progenitor of SN1998bw, $`10^{48.5}`$ erg (Tan et al., 2001).

X-ray observations do not provide therefore evidence for existence of an off-axis ”standard” GRB jet associated with SN1998bw, that may have produced GRB980425. The lack of observational signatures typically expected to be produced by such an off-axis jet on a 1 yr time scale may be due, however, to the low value of $`\dot{m}`$, $`\dot{m}0.1`$, which implies that an off-axis jet will become observable only on $`10`$ yr time scale (Waxman, 2004).

This research was supported in part by AEC, Minerva, and ISF grants.


text:
# Enhancement of solid-state proton NMR via SPINOE with laser-polarized xenon

## I Introduction

Many difficulties appear in scaling the liquid-state Nuclear Magnetic Resonance (NMR) quantum computers to multi-qubit, since the signal intensity is low sensitive and decreases exponentially with increasing the number of qubits. Therefore, the maximal number of qubits of quantum computer using the current liquid-state NMR technology might be limited to 10 Warren . Although there are some methods developed for enhancing polarization, such as lower temperature and higher magnetic fields, Dynamic Nuclear Polarization (DNP) Abragam ; Bates , Nuclear Overhauser Effect (NOE) Overhauser , cross polarization Slichter , etc., unfortunately these approaches can provide only limited relief or need rigid conditions for polarization enhancement. A promising route to enhance the nuclear polarization is the utilization of optical pumping and spin exchange Happer84 ; Zeng ; Happer , which can increase the nuclear polarization by four or five orders of magnitude over thermal equilibrium. The hyperpolarized 129Xe and 3He gases are greatly effective in the development of magnetic resonance imaging Albert , surface science Pietrass , probing of biological system Rubin ; Cherubini , polarized targets Xu , neutron polarization Jones , precision measurement Bear , polymer scienceNagasaka , quantum computation Chuang , etc..

Hitherto, proton sensitivity enhancements are not very large in all experiments Song ; Navon ; Fitzgerald ; Room  except for surface enhancements Happer93 , and all of these experiments have almost been implemented by using gaseous or liquid hyperpolarized xenon. The 1H polarization has been enhanced by a factor of 0.1 to 2 on a 4.2T NMR spectrometer at room temperature via cross-relaxation between dissolved hyperpolarized gaseous 129Xe and 1H of the liquid benzene solvent, which has also been firstly called the SPINOE by Pines’ groupNavon . Through dissolving compounds in hyperpolarized liquid xenon, the enhanced signal of over 45 for 1H has also been observed at 1.4T and 200K by Happer’s group Fitzgerald . In Xe ice, larger nuclear polarization of hyperpolarized 129Xe has been transferred to 13CO2 by low-field thermal mixing Bowers . However, as far as we know, solid-state proton polarization enhancement has not been reported in the literature.

Solid-state quantum computers, which could execute the 106 qubits operation brought out by DiVincenzo et al., were widely noticed DiVincenzo . Considering that hydrogen chloride (HCl) and xenon have almost the same melting points between 160 and 170K, and the longitudinal relaxation time (T1) of hyperpolarized 129Xe in solid-state (4K) is about 500 hours Gatzke , we put forward the idea of laser-enhanced solid-state quantum computer (LESSQC) by mixing hyperpolarized 129Xe and 1HCl.

In this letter, we firstly demonstrate that the polarization of hyperpolarized solid-state 129Xe produced by spin exchange is transferred to the proton of solid-state 1HCl via SPINOE without using low-field thermal mixing or Hartmann-Hahn matching condition. This method yields the proton enhancement of 6 by comparison of that without optical pumping on a Bruker SY-80M NMR spectrometer (1.87T) at 142 K, leading us to take a step towards realizing LESSQC.

## II Experiment and Results

The experimental apparatus, similar to that in Ref. 26, is shown schematically in Figure 1. Briefly, the whole system consists of an optical pumping system (the right part of the valve K4), which is used to produce the hyperpolarzied gaseous 129Xe, and a cross-relaxation and detection system (the left part of the valve K4). The two parts are connected by a cylindrical Pyrex tube and separated by stopcocks. The whole system was evacuated with K2, K5 valves closed and K1, K3 and K4 valves opened. When the vacuum reached to 1.5$`\times `$10-5 Torr for several hours, the valves K1, K4 and K5 were closed and K2 and K3 opened. The cylindrical pump cell containing a few drops of metal Cs, with a volume of 600 cm3, was loaded 760 Torr natural isotopic xenon gas at room temperature and an atmosphere, then all valves were closed. The pump cell, which was placed in a 25 Guass magnetic field generated by Helmholtz coils, was maintained at approximately 333$`\pm `$4 K by a resistance heater (not shown) during the optical pumping. The inner surfaces of the cylindrical Pyrex tube and the pump cell were coated with silane in order to slow down the relaxation of the 129Xe upon collision with the tube wall. As soon as all valves were closed, laser light from a 15 W cw tunable semiconductor-diode laser array (Opto Power Co. Model OPC-D015-850-FCPS) at wavelength $`\lambda `$=852.1 nm was introduced to the system. After passing through beam expander, Glan prism, $`\lambda `$/4 plate and convex lens, the laser light became circularly polarized, and illuminated almost 4/5 of the pump cell volume. The propagation of the laser paralleled to the orientation of the magnetic field B. The circularly polarized laser resonates with the Cs D2 absorption line and induces an electron spin polarization in the Cs atoms via a standard optical pumping process Zeng ; Happer . After about 25 minutes, the hyperpolarized 129Xe gas was produced by spin-exchange collision with the Cs atoms.

At the part of cross-relaxation and detection system, there were a liter of 1HCl gas at room temperature and an atmosphere in the tank. K4 and K5 valves were opened to allow the hyperpolarized 129Xe and 1HCl to be mixed, then transferred to the NMR sample probe (10 mm diameter, 8 mm inner diameter), pre-cooled to 172$`\pm `$2 K, of the Bruker SY-80M NMR spectrometer. The temperature of the probe was subsequently reduced to 142 K in an evacuated glass dewar by flowing cold nitrogen gas from liquid nitrogen tank, which controlled by the Bruker variable temperature unit. This temperature was kept, and the variation was less than 2 K at all times so that the mixture could be in the solid-state. Since there exist imbalance of polarization, the hyperpolarized 129Xe and the thermally polarized proton, in the mixed system, the interaction between the 129Xe and the proton can make spin exchange occur via SPINOE.

The time dependence of the 129Xe NMR signal intensity (solid circle) observed after blending hyperpolarized 129Xe and 1HCl is shown in Figure 2(a). The observed spin-lattice relaxation time of 129Xe in the presence of 1HCl is 29.6$`\pm `$0.6 min. The initial rise of signal manifests the accumulation of solid-state hyperpolarized 129Xe in the probe, similar to that in Ref. 19. At the peak signal, the increase of accumulated 129Xe magnetization reaches a balance with the decay of 129Xe magnetization via SPINOE with 1HCl. After this peak point, 129Xe magnetization decays as a result of SPINOE. So the solid line, a fit to data after the peak, represents the time dependence of 129Xe residual magnetization after a SPINOE experiment.

Figure 2(b) displays the time dependence of the solid-state proton NMR signal (solid circle) from the same run. Both the accumulation of the condensed state mixture and the cross-relaxation between the proton and the hyperpolarized 129Xe contribute to the initial rise of the proton signal. The proton NMR signal reaches its peak vaule at a time t=280 s after the mixing, and decays towards its thermal equilibrium value (dashed line) at the same rate as 129Xe signal. The solid line represents a fit to data to guide eyes. We used pulse flip angles of 40 and 90 for 129Xe and 1H respectively, and 129Xe was performed with a home-bulit probe.

Typcial spectrum for the enhanced proton is shown in Figure 3. In order to obtain the large 1H NMR signal, we used 900 pulse angle for a single acquisition at a time t=280 s. Figure 3(b) presents a typical enhanced 1H NMR spectrum comparing with the one at thermal equilibrium (Figure 3(a)). Due to the relation of the magnetization versus the nuclear spin polarization given by Abragam Abragam , the enhancement factor of the solid-state 1H was about 6, which corresponds to the proton polarization of 8.55$`\times `$10$`^5,`$ on the basis of compariosn of the integrated intensity of the laser-enhanced signal with that at equilibrium.

## III Theoretical Analysis and Discussion

The interactions are complicated in solid-state, but in this paper, we will limit our discussion to only spins-$`\frac{1}{2}`$ nuclei coupled via the direct nuclear dipolar interaction:

$$H_D=\frac{\gamma __S\gamma __I}{r^3}[SI3(rS)(rI)],$$
(1)
where $`I`$ and $`S`$ are the spin angular momentum operators, $`\gamma _I`$ and $`\gamma _S`$ are the gyromagnetic ratios of I and S spins, respectively, and $`r`$ is the distance between two spins. Due to the perturbation of Hamiltonian $`H_D,`$ the transition probabilities $`W`$ between the eigenstates can be given Solomon :

$$\begin{array}{c}W_0^{IS}=\frac{2\delta }{20}J(\omega _I\omega _S),\\ W_{1I}^{IS}=\frac{3\delta }{20}J(\omega _I),\\ W_{1S}^{IS}=\frac{3\delta }{20}J(\omega _S),\\ W_2^{IS}=\frac{12\delta }{20}J(\omega _I+\omega _S),\end{array}$$
(2)
with

$$\delta =\frac{\mathrm{}^2\gamma _S^2\gamma _I^2}{r^6},$$
(3)
$$J(\omega )=\frac{\tau _c}{1+\omega ^2\tau _c^2},$$
(4)
here, $`\mathrm{}`$ is the Planck constant divided by 2$`\pi `$, and $`\tau _c`$ is the correlation time of spin systems.

Accroding to the Solomon equations, the evolution of this two-spin systems, 1H (I=$`\frac{1}{2}`$) and 129Xe (S=$`\frac{1}{2}`$), could be described by Solomon ; Abragam :

$$\frac{d}{dt}\left(\begin{array}{c}I_z\\ S_z\end{array}\right)=\left[\begin{array}{cc}\rho _I& \sigma _{IS}\\ \sigma _{SI}& \rho _S\end{array}\right]\left(\begin{array}{c}I_zI_0\\ S_zS_0\end{array}\right),$$
(5)
where Iz and Sz are z components of the spins I and S respectively, I0 and S0 are their equilibrium values, $`\rho _I`$ and $`\rho _S`$ are the autorelaxation rates of the 1H and the 129Xe spins, and $`\sigma _{IS}`$ and $`\sigma _{SI}`$ are the corresponding cross-relaxation rates. The elements of cross-relaxation matrix can be expressed by the transition probabilities W resulting from II, SS and IS interactions Ernst :

$$\begin{array}{c}\rho _I=2(n_I1)(W_1^{II}+W_2^{II})+n_S(W_0^{IS}+2W_{1I}^{IS}+W_2^{IS}),\\ \rho _S=2(n_S1)(W_1^{SS}+W_2^{SS})+n_I(W_0^{IS}+2W_{1S}^{IS}+W_2^{IS}),\\ \sigma _{IS}=n_S(W_2^{IS}W_0^{IS}),\\ \sigma _{SI}=n_I(W_2^{IS}W_0^{IS}),\end{array}$$
(6)
here, $`n_I`$ and $`n_s`$ are magnetically equivalent I and S spins. A solution of equation (5) of particular interest in this system is the one corresponding to the initial conditions:

$$\begin{array}{c}(I_zI_0)_{t=0}=0,\\ (S_zS_0)_{t=0}=S_i.\end{array}$$
(7)
So the solution can be given by:

$$\begin{array}{c}I_z(t)=I_0+C(\mathrm{exp}(\lambda _1t)\mathrm{exp}(\lambda _2t)),\\ S_z(t)=S_0+C[r_1(\mathrm{exp}(\lambda _1t)r_2\mathrm{exp}(\lambda _2t))],\end{array}$$
(8)
where $`\lambda _1`$ and $`\lambda _2`$ are given by:

$$\begin{array}{c}\lambda _1=\frac{(\rho _I+\rho _S)\sqrt{(\rho _I\rho _S)^2+4\sigma _{IS}\sigma _{SI}}}{2},\\ \lambda _2=\frac{(\rho _I+\rho _S)+\sqrt{(\rho _I\rho _S)^2+4\sigma _{IS}\sigma _{SI}}}{2},\end{array}$$
(9)
and

$$\begin{array}{c}r_1=\frac{\rho _S\sqrt{(\rho _I\rho _S)^2+4\sigma _{IS}\sigma _{SI}}}{2\sigma _{IS}},\\ r_2=\frac{\rho _S+\sqrt{(\rho _I\rho _S)^2+4\sigma _{IS}\sigma _{SI}}}{2\sigma _{IS}},\end{array}$$
(10)
$$C=\frac{S_i}{r_1r_2}.$$
(11)
Therefore, the 1H enhancemet via cross-relaxation comparing with that at thermal equilibrium is:

$$\frac{I_z(t)I_0}{I_0}=\frac{\sigma _{IS}}{\rho _I}\frac{\gamma _S}{\gamma _I}\frac{S_z(t)S_0}{S_0},$$
(12)
where $`[S_z(t)S_0]/S_0`$, the enhancement of hyperpolarized solid-state 129Xe, is about 6000 in our experiment, which corresponds to the 129Xe polarization of 2.16% Zhou . Inserting the values (2) and (6) into (12) we get:

$$\frac{I_z(t)I_0}{I_0}=\frac{\gamma _S(S_z(t)S_0)n_S[6\delta J(\omega _I+\omega _S)\delta J(\omega _I\omega _S)]}{\gamma _IS_0\{(n_I1)[3\delta ^{}J(\omega _I)+12\delta ^{}J(2\omega _I)]+n_S[\delta J(\omega _I\omega _S)+3\delta J(\omega _I)+6\delta J(\omega _I+\omega _S)]\}},$$
(13)
with

$$\delta ^{}=\frac{\mathrm{}^2\gamma _I^4}{r^6},$$
(14)
here $`r^{}`$ is the distance between two neighboring I spins.

In our experiment, $`\omega _I`$ (1H) and $`\omega _S`$ (129Xe) are 80.13MHz and 22.16MHz (on Bruker AC-80 spectrometer). The concentrations of 1H and 129Xe are 34.6mmol/cm3 and 5.46mmol/cm3 respectively. We can image a model that one 129Xe nucleus is surrounded by six 1H nuclei on average. The autorelaxation rate of the proton $`\rho _I`$ of (66 s)-1 is measured in the experiment so that one can obtain the cross relaxation $`\sigma _{IS}`$ by the equation (12). Therefore, the correlation time $`\tau _c`$ is estimated to be 6.18$`\times `$10-5 s. Using the equation (6), we can calculate that autorelaxation rates $`\rho _I`$ and $`\rho _S`$ are (81 s)-1 and (62 min)-1, and cross relaxation rates $`\sigma _{IS}`$ and $`\sigma _{SI}`$ are 5.53$`\times `$10-5 s-1 and 3.32$`\times `$10-4 s-1, respectively. Thus from equation (8), the time evolution of the proton polarization can be theoretically writen as:

$$I_z(t)=I_0+C(\mathrm{exp}(t/3744)\mathrm{exp}(t/81)).$$
(15)
The theoretical simulation (dot line) is visualized in Figure 2(b) for clear comparison to experimental results. Although only the direct nuclear dipolar interaction is considered in the above discussion, in fact, there are so many factors, such as spin rotation, paramagnetisim, relaxation with the wall, inhomogenuous magnetic field, etc., which can induce the relaxation. Consequently the theoretical relaxation rates are smaller than the experimental ones, and predictions of the time dependence of the proton polariztion are larger than experimental results.

Assuming $`r`$ is equal to $`r^{}`$ for simplification in our theoretical model, one can calculate the maximum solid-state proton enhancement of 7.1 using equation (13), which is in general agreement with the measured value. Although it is substantially smaller than the 129Xe enhancement, this can not indicate the possibility that only a fraction of the total number of 1HCl molecules interact with the hyperpolarized 129Xe by diploar-diploar. Because only the natural xenon (26.4% enriched 129Xe, 21.2% enriched 131Xe) has been used, at the temperature of 142K, the cross-relaxation between 129Xe and the isotope 131Xe can decrease the efficiency of polarization transfer from 129Xe to 1H Gatzke . On the other hand, dipolar relaxation from diffusing vacancies dominates at this temperature for 129Xe spin-lattice relaxation Cates . If we take into account of all above factors, together with the inhomogenuous magnetic field due to the fluctuation of eletromagnet and the polarization loss during the phase transition, etc., all of these losses maybe counterbalance the difference between experimental values and theoretical ones.

According to equation (13), for further increasing the proton enhancement, one should first obtain the greatest 129Xe enhancement via employing the high power and narrow bandwidth laser, and/or increase the gas pressure of the pump cell in order to enhance the optical-pumped absorbed power when using the wide bandwidth laser. A small proton number density by using partially deuterated sample and a large 129Xe number density by using Xe isotopically enriched with 129Xe also can boost transfer efficiency.

## IV Conclusion

In conclusion, we firstly obtained the solid-state laser-enhanced NMR proton signal of 1HCl via SPINOE by cooling hyperpolarzied 129Xe to solid-state at 1.87T and 142K in the flow system. The amplification factor of the proton nuclear polarization is about 6. The theoretical calculated enhancement is in agreement with the measured value. Furthermore, the hyperpolarized 129Xe has long lifetime in the solid-state, which means that 1H could keep the longer coherence time with hyperpolarized 129Xe. Therefore this method may be useful to overcome one of the difficult problems in liquid-state NMR quantum computer and to establish solid-state quantum computers.

Solid-state signal enhancement via SPINOE with hyperpolarized 129Xe is not limited to the proton, i.e., it can be expanded to other nuclei. Although we have focused on quantum computers, this method should be readily extended to material science and determination of the 3D structure of large biomolecules, since NOE can provide unique information on molcular structure which can not be obtained with any other known technique. Traditional NOE needs to irradiate one nucleus in order to observe anther. However, it is very convenient and available for probing interactions between nuclei without any additional conditions by using this method.

## V Acknowledgment

This work is supported by the National Natural Science Foundation of China under Grant No. 10234070, National Science Fund for Distinguished Young Scholars under Grant No. 29915515 and National Fundamental Research Program under Grant No. 2001CB309306.

Figure captions

Fig. 1. Schematic diagram of the experimental setup.

Fig. 2. (a) Time dependence of the solid-state 129Xe NMR signal (solid circle) observed after the blend of hyperpolarized 129Xe and 1HCl at 142 K and 1.87 T. The solid line, a fit to data after the peak, represents how much 129Xe magnetization is left after a SPINOE experiment. (b) Time dependence of the solid-state proton NMR signal (solid circle) from the same sample, and it relaxes towards its thermal equilibrium value (dashed line). The dot line represents the theoretical predictions for comparison to experimental results (solid line).

Fig. 3. (a) The solid-state proton NMR spectra of 1HCl under the conditions of thermal equilibrium and (b) laser-enhanced signal via SPINOE with hyperpolarized 129Xe.