text:
# Random manifolds in non-linear resistor networks: Applications to varistors and superconductors

## Abstract

We show that current localization in polycrystalline varistors occurs on paths which are, usually, in the universality class of the directed polymer in a random medium. We also show that in ceramic superconductors, voltage localizes on a surface which maps to an Ising domain wall. The emergence of these manifolds is explained and their structure is illustrated using direct solution of non-linear resistor networks.

Low energy paths and surfaces in random systems are of broad importance in statistical mechanics and in materials’ theory. Examples include the morphology of flux lines in the presence of random pinning, magnetic domain walls in random-bond Ising magnets, dielectric breakdown paths and fracture surfaces. The directed polymer in a random medium has received an enormous amount of interest, both due to its intrinsic importance and due to its relation with growth processes. Domain walls in random-bond Ising models have also received a good deal of attention, in fact the directed polymer in a random medium was first used to model domain walls in two-dimensional random-bond magnets. More recently it has been realised that random path and surface problems map to classic problems in computer science. In particular, the shortest-path problem is related to the directed polymer in a random media and the minimum-cut problem is related to domain walls in random Ising magnets. The shortest path and minimum cut problems have polynomial time algorithms, and this has enabled physicists to solve path problems with overhangs (i.e. not just the directed case), and domain wall problems in three dimensions, both of which were previously considered difficult. Many more connections between statistical physics and solvable combinatorial optimization problems have recently proven profitable. Here we show that the shortest path and minimum-cut problems also emerge in the study of the onset of current or voltage in non-linear resistor networks and in their applications to varistors and superconductors. This enables a detailed characterization of the manifolds on which current or voltage localize near critical thresholds.

Non-linear flow in random networks is germane to problems ranging from non-Newtonian fluid flow in porous media, to the onset of flow in varistors, and to the onset of voltage in grain-boundary limited superconductors. The emergence of special flow paths in the fluid and varistor cases is well-established and network models with realistic non-linear flow laws have been constructed. Recently lattice models have also been used to model the behavior of polycrystalline superconductors. The networks which model these materials have highly non-linear current-voltage (or pressure-flow) behavior on each bond of a graph. In the case of diodes or varistors each bond has a critical voltage, $`v_c^k`$, which varies from bond to bond . The onset of current flow occurs on special paths and the identification and characterisation of these paths is a central issue in the analysis of varistors. In the superconductor case, each bond has a local critical current, $`i_c^k`$, below which the bond carries current but has zero voltage. Above this threshold current, a local voltage must be applied to the bond in order to increase the current further. The onset of voltage occurs on a surface which we show maps to a minimum cut through the network.

Firstly we define the shortest path and minimum-cut problems. These problems (and the algorithms which solve them) are very basic to algorithmic systems in computer science. They are defined on a graph where each edge is assigned a cost (shortest-path problem) or a capacity (minimum-cut problem). Consider then a graph composed of nodes and edges, where each edge has a cost, $`c_k`$. The shortest path between any two sites, $`s`$ and $`t`$, in the graph is simply the path on which the sum of the edge costs,

$$C_{st}=\mathrm{min}_P(\underset{kϵP}{}c_k),$$
(1)
on the path is smallest. In physics language the shortest path is the lowest energy path. It is clear that this is the same as the problem of a directed polymer in a random media, provided the path is forced to be directed. However if overhangs are allowed this problem can also be solved efficiently using, for example, Dijkstra’s method. It turns out that overhangs do not change the universality class of the problem unless disorder is strong.

The minimum-cut problem is defined on a graph consisting of nodes and edges in which each edge has a flow capacity, $`u_k`$. This is the maximum amount of flow that the edge can carry. There is no energy cost in the computer science definition of this problem, there is just the limiting capacity, $`u_k`$ (However in the end some of the physics applications associate the capacity with an energy). Now consider the amount of flow that can be carried between two sites, $`s`$ and $`t`$. A key theorem states that the maximum flow between $`s`$ and $`t`$ is equivalent to the capacity of the minimum cut between $`s`$ and $`t`$. The minimum cut is the lowest capacity surface of separation in which $`s`$ is on one side of the cut and $`t`$ is on the other. The capacity of the minimum cut, $`U_{st}`$, is the sum of the capacities on this surface of separation, $`S`$,

$$U_{st}=\mathrm{min}_S(\underset{kϵS}{}u_k).$$
(2)
If we associate the bond capacities with the exchange constants in a random-bond Ising magnet, then the minimum cut maps to the lowest energy domain wall. Somewhat surprisingly, the minimum cut can be found efficiently (in almost linear time in the number of edges) and this has enabled detailed study of domain walls in three dimensions.

Now we consider in more detail the emergence of the shortest path and the minimum cut in non-linear random resistor networks. To be concrete, we consider square networks where each bond has a monotonic (single-valued) non-linear V-I behavior. The case of current onset at a critical voltage is applicable to varistors(Fig. 1) and diodes and also to Bingham plastics (replace voltage by pressure and current flow by fluid flow so that flow initiates at a critical pressure). The case of a voltage onset at a critical current is applicable to ideal josephson junctions (Fig. 2) and also to ideal flux-flow laws . Ideal varistor behavior is mathematically described by, $`i(v)=(i_0+\frac{vv_c}{r})\theta (vv_c)`$, while ideal Josephson junction or flux flow behavior is described by, $`v(i)=(v_0+\frac{ii_c}{g})\theta (ii_c)`$, where $`\theta (x)`$ is the step function which is zero for $`x<0`$ and one for $`x0`$. Finding the current or fluid flow through these non-linear networks seems difficult due to the non-smooth behavior at the critical current or voltage. However a variety of regularization methods are available to handle these difficulties. The simplest is to replace the step function by, $`\theta (x)=x/(\xi ^2+x^2)^{1/2}`$ and to take the limit $`\xi 0`$. This form has nice analytic properties which assist in the numerical analysis, so we have used this regularization in our codes. A key observation is that these problems have a related cost function defined by, $`cost=_0^iv(i^{})𝑑i^{}`$, which is convex provided the V-I law on each bond is monotonically increasing. With suitable regularization this is true and solving these non-linear resistor networks then reduces to a convex optimization problem with the constraint of flow conservation at each node. We have used this framework to develop efficient codes for this class of problem. Here we discuss applications to varistors and superconductors.

In the varistor case (see Fig. 1) each bond has a random onset voltage, $`v_c^k`$, and the whole network has the onset voltage, $`V_c`$. At $`V_c`$, current only flows on the shortest path through the network, as we confirmed numerically by finding the path on which $`V_c=_{kϵP}v_c^k`$ is minimal (using Dijkstra’s method). Voltage localization in the superconductor case is illustrated in Fig. 2. In this case each bond has a random critical current, $`i_c^k`$, and the whole network has the critical current, $`I_c`$. We solved the full flow equations directly and then confirmed that at $`I_c`$ voltage localizes on the minimum cut, i.e. the surface on which $`I_c=_{kϵS}i_c^k`$ is minimal.

The easiest way to prove that random manifolds emerge at the macroscopic thresholds, $`I_c`$ and $`V_c`$, is by using the cost function, $`_0^iv(i^{})𝑑i^{}`$. From this cost function, it is evident that if either the current or voltage can be kept to zero, then the cost itself is zero. In the case of a varistor, as the external voltage is increased, the network distributes the voltage drops so as to keep the current at zero. It succeeds in doing this until it is no longer possible and the first time this occurs is on the shortest path, occuring at the critical voltage, $`V_c`$. This argument applies to any cost function which has strictly zero current up to threshold. This means that the behavior of the local $`vi`$ characteristic after threshold has no influence on the value of $`V_c`$. In a similar manner, in the superconductor case, a large network attempts to keep the voltage at zero in all of the bonds until it is impossible to do so. As the applied current is increased, the current is distributed so as to keep the currents in all of the bonds below their threshold value, again in a very co-operative manner. The first applied current, $`I_c`$, at which voltage appears is determined by when a surface of bonds all have reached their local critical currents. This surface is the minimum cut.

To make contact with experiment, note that varistors are materials which are insulating below a critical electric field, $`E_c`$, after which they become highly conducting. The physics of their operation is understood to originate in the grain boundaries of polycrystalline ceramics. The interior of the grains of typical varistors are conducting, but the grain boundaries are insulating at low voltage (below about 3V). The onset of current flow occurs when the applied voltage is large enough to induce a path of grain boundaries to become conducting. This physics is naturally encapsulated in lattice models in which a node represents a grain and a bond represents a grain boundary . However the grains have varying sizes and the grain boundaries have varying properties, with some boundaries having higher onset voltages than others. That is, we must study disordered networks. These networks are also highly non-linear as has been demonstrated by studies of flow onset across individual grain boundaries. The linear dimension, $`L`$ of the lattice models is related to the size of real varistors by $`L=l/g`$, where $`l`$ is the sample size and $`g`$ the grain size. Considerable progress has been made using numerical simulations which have provided good agreement with experiment. Here we make the connection with the statistical physics of random paths.

As discussed above, network models of polycrystalline varistors have a current onset given by the minimum over all paths of $`V_c=_{kϵP}v_c^k`$. This is just the shortest-path problem of Eq. (1), where the voltage threshold on a bond replaces the cost of that bond. The shortest-path problem is usually in the universality class of the directed polymer in a random medium(DPRM) (execept for strong disorder) so that $`V_c`$ behaves like the energy of a DPRM, i.e.

$$V_c=a_1L+a_2L^\theta ,$$
(3)
where $`a_1`$ and $`a_2`$ are independent of sample size and the exponent $`\theta `$ is universal. Its value is known to be exactly $`1/3`$ for paths through two dimensional systems and $`0.248\pm 0.004`$ for paths through three dimensional media. A key feature of (3) is that the threshold electric field of the lattice $`E_c=V_c/L`$ is size independent, in contrast to the size effects produced by rare-fluctuation theories of electrical and dielectric failure, where $`E_c`$ approaches zero logarithmically in the large lattice limit. However the key difference is that in dielectric breakdown, local regions irreversibly make the transition to the conducting state, so that no voltage is required to maintain these regions in the conducting state after failure. In contrast, high quality varistors are reversible so that there is a steady state current at a fixed applied voltage. The path on which current flows is usually self-affine and has length, $`L_p`$, given by, $`L_p=b_1L+b_2L^\zeta `$, where $`b_1`$ and $`b_2`$ are dependent on the disorder distribution, but the exponent $`\zeta `$ is universal. Its value is known to be exactly $`2/3`$ for paths through two dimensional systems and $`0.62\pm 0.01`$ for paths through three dimensional media. The exponents $`\theta `$ and $`\zeta `$ are related by, $`\theta =2\zeta 1`$, which is consistent with the numerical results quoted above. The difference in voltage between the lowest threshold path and the onset voltage of the next filamentary current path scales in the same way as the energy gap in the DPRM problem, which decreases logarithmically with increasing sample size. In the strong disorder limit, the paths become highly tortuous and are no longer in the DPRM class.

In the case of ceramic superconductors the low-angle grain boundaries have much higher critical currents than the high-angle grain boundaries. As in the varistor case, lattice models take the grain centers to be nodes and the links between grains to represent the grain boundaries. Grain boundaries may have a flux flow character or a Josephson Junction character or be resistive. If a grain boundary is resistive, its critical current is zero. As demonstrated above, the onset of voltage occurs on the minimum cut through such a network and this minimum cut is related to domain walls in Ising magnets . The critical current of these networks then behaves in the same way as the energy of Ising domain walls, i.e.,

$$I_c=c_1L^{d1}+c_2L^\theta $$
(4)
where $`d`$ is the spatial dimension, $`c_1`$ and $`c_2`$ are independent of sample size, and $`\theta =1/3`$ in two dimensions and $`\theta =0.82\pm 0.02`$ in three dimensions. There is thus a sample size independent critical current density in these networks, in contrast to the size effect which occurs in fuse networks. Again this difference is due to the fact that these networks must maintain a current of at least $`i_c^k`$ on the $`k^{th}`$ bond in order for a voltage to appear there. The surface on which voltage localizes at $`I_c`$ (i.e. the minimum cut) is in the universality class of domain walls in the random-bond Ising model and so it is self-affine and has asymptotic roughness, $`w`$, given by, $`w=c_3L^\zeta `$, where $`\zeta =2/3`$ in two dimensions and $`\zeta 0.41\pm 0.01`$ in three dimensions. The exponents $`\theta `$ and $`\zeta `$ are related to each other by the scaling relation, $`\theta =2\zeta +d3`$.

We have demonstrated that current or voltage localize on special manifolds in non-linear networks. As particular examples, we showed that flow in a varistor begins on a path which is equivalent to a shortest path through the random medium, and in superconductors, voltage onset occurs on a minimum cut through the network. The macroscopic current or voltage thresholds are then equivalent to the energy of a random manifold. Since these macroscopic thresholds only depend on the local current or voltage thresholds (and not on the local $`vi`$ characteristics after threshold) the manifolds on which current or voltage localize in varistors and superconductors are described by the universal directed polymer or random surface exponents.

This work has been supported by the DOE under contract DE-FG02-90ER45418, and by a subcontract through Sandia National Laboratories.


format: text

text:
# Novel Methods for Determining Effective Interactions for the Nuclear Shell Model

## I Introduction

The shell-model has long been the standard paradigm for treating the nuclear many-body problem. However, quite recently, the form of the shell model most frequently employed in practice has been criticized by Haxton and coworkers. The chief objection in these papers is that shell-model practitioners often use effective interactions which do not follow in any systematic way from realistic nucleon-nucleon (NN) interactions nor from the inevitable truncations that these entail. They point out that theoretical schemes exist that indicate how “bare”interactions should be modified to account for effects of truncations. Moreover, these schemes emphasize the necessity of using consistently modified operators rather than the bare ones used in many shell-model studies. Finally, they stress that tremendous recent advances in raw computational power and in numerical techniques now make it feasible to implement these theoretical schemes.

During the past few years, two of us (JP and JRS) have extensively studied low-dimensional quantum spin systems, especially the so-called spin-ladder materials . We note that interest in the study of the ladder compounds has been stimulated by suggestions that these deceptively simple materials could exhibit some of the critical behavior believed to be responsible for high-temperature superconductivity . In much of this work, we have made use of a novel “plaquette” basis and of a truncation scheme of the full Hilbert space that is suggested by the structure of this basis. Central to this work has been the use of the COntractor REnormalization (CORE) method of Morningstar and Weinstein. We have relied on CORE for the construction of an effective Hamiltonian and effective operators to be employed in concert with the truncated basis. Our success in applying CORE to important condensed-matter problems has inspired us to attempt adopting the approach to the nuclear-shell model as it directly addresses many of the concerns expressed in the preceding paragraph.

Another important development, the use of effective-field theories in nuclear physics, has been stimulated by novel ideas in modern renormalization theory . Effective theories (ET’s) are low-energy approximations to the exact (and, thus, very complicated) theory. One of the central tenets behind ET’s is that a long-wavelength particle (of wavelength $`\lambda `$) should be insensitive to details of the potential at very short ($`d\lambda `$) distances. Hence, it should be possible to replace the very complicated short-distance structure of the exact theory by a potential with a simpler short-range structure, while leaving the low-energy properties of the theory intact. The aim of this paper is to combine both approaches—CORE and effective theories—in a complimentary manner.

The paper has been organized as follows. In Sec. II we outline the basic steps in the implementation of CORE to the nuclear many-body problem. Particular emphasis is placed in the construction of the similarity transformation that generates the effective interaction. In the same section the basic ideas behind modern effective theories are introduced together with a detailed procedure on how to merge these two (CORE and ET’s) powerful approaches. In Sec. III we illustrate how the effective potential is constructed and the subsequent refinements of the basis that CORE entails. In particular, we focus on predictions of the ground-state energy of the three-body systems using a variety of approximations. Finally, our conclusions and plans for future calculations are presented in Sec. IV.

## II Formalism

In this section we describe in detail how CORE is to be used in the many-body problem. Central to the spirit of CORE is a similarity transformation that ultimately generates a renormalized Hamiltonian to be used in a severely truncated Hilbert space. We follow this discussion with a brief introduction to modern versions of effective theories placing particular emphasis on its use in combination with CORE.

### A CORE for the Shell Model

Here we will briefly outline the CORE method as adapted for the nuclear-shell model. The underlying dynamics is described by a non-relativistic Hamiltonian containing only two-body forces

$$H_A=\underset{n=1}{\overset{A}{}}\frac{𝐩_n^2}{2m}+\underset{m<n=1}{\overset{A}{}}V(𝐫_m𝐫_n),$$
(1)
where $`V(𝐫_m𝐫_n)`$ is a realistic NN interaction. Because of potential problems with center-of-mass motion, the above Hamiltonian is modified by the introduction of a harmonic term that depends exclusively on the coordinate of the center of mass $`𝐑`$:

$$V_A^\mathrm{\Omega }(𝐑)=\frac{1}{2}mA\mathrm{\Omega }^2𝐑^2.$$
(2)
A many-body basis is then constructed from products of single particle harmonic oscillator wave functions characterized by the same oscillator frequency $`\mathrm{\Omega }`$. The modified $`A`$-body Hamiltonian, the standard starting point for most shell-model calculations, is now written as

$$H_A^\mathrm{\Omega }=\underset{n=1}{\overset{A}{}}\left[\frac{𝐩_n^2}{2m}+\frac{1}{2}m\mathrm{\Omega }^2𝐫_n^2\right]+\underset{m<n=1}{\overset{A}{}}\left[V(𝐫_m𝐫_n)\frac{m\mathrm{\Omega }^2}{2A}(𝐫_m𝐫_n)^2\right].$$
(3)
The terms depending on $`\mathrm{\Omega }`$ allow separation of the center-of-mass motion and, after trivial subtractions of the center-of-mass energy, the calculations yield the desired physical spectrum which ultimately is independent of $`\mathrm{\Omega }`$.

The first step in the implementation of CORE is the identification of an “elementary” block and a truncation scheme . In our case the elementary block, or one-body term, consists of the single-nucleon Hamiltonian appearing in Eq. (3). This one-body term, henceforth denoted by $`h_1`$, defines the basis via products of one-body harmonic-oscillator eigenstates. Subsequently, the truncation scheme is defined by retaining only those states with energy less than some maximum (cutoff) value $`E_{\mathrm{max}}`$, or equivalently oscillator quanta less than some maximum value $`N`$. In this way it is guaranteed, by construction, that the spectrum of the effective Hamiltonian $`h_1`$ reproduces identically the low-energy properties of the exact one-body Hamiltonian.

The next step involves computation of the effective Hamiltonian for the two-body system. To begin, the two-body problem must be solved exactly. This is accomplished by using a very large basis, i.e., a basis as large as can be afforded computationally and one characterized by a cutoff $`N_{\mathrm{}}`$ which is effectively infinite. For the method to have any utility we evidently require $`N_{\mathrm{}}N`$. CORE now instructs us that the effective two-body Hamiltonian be given in terms of a similarity transformation:

$$H_{\mathrm{eff}}(1,2)SE_{\mathrm{diag}}S^T.$$
(4)
where $`E_{\mathrm{diag}}`$ is the diagonal matrix containing the lowest $`𝒩`$ exact eigenvalues of the two-body system where $`𝒩`$ is the number of two-body states defined by the cutoff $`N`$. Evidently, the structure of $`H_{\mathrm{eff}}(1,2)`$ guarantees, independent of the details of the similarity transformation, the low-energy spectrum of the full theory will be preserved. We now discuss the particular form of $`S`$ prescribed by CORE.

By elementary linear algebra, the columns of the similarity transformation are the eigenvectors of $`H_{\mathrm{eff}}(1,2)`$. CORE demands that the first column of $`H_{\mathrm{eff}}(1,2)`$, i.e., the eigenvector associated with the lowest eigenvalue, has a non-zero overlap with the lowest exact eigenvector of the full Hamiltonian. A simple way to implement this demand is by choosing the first eigenvector in the model space to be proportional to the lowest exact eigenvector after it has been truncated and then properly normalized. To construct the second column of the similarity transformation one starts with the second exact eigenvector, truncates it and then, via a Gram-Schmidt procedure, subtracts any component parallel to the first column. This procedure guarantees the fulfillment of another important CORE requirement: the $`n^{\mathrm{th}}`$ eigenvector in the model space must have a non-zero overlap with the $`n^{\mathrm{th}}`$-lowest exact eigenvector yet no such overlap is allowed with the previous $`n1`$ eigenstates. The third column is determined in a similar way and must have no component parallel to the first two. This process is repeated until the $`𝒩`$ columns of $`S`$ are completely determined. Note that, in practice, Gram-Schmidt is numerically unstable and more sophisticated—but, in principle, equivalent—techniques must be used. Finally the effective two-body interaction is obtained by subtracting the one-body terms from $`H_{\mathrm{eff}}(1,2)`$:

$$V_{\mathrm{eff}}(1,2)=H_{\mathrm{eff}}(1,2)h_1(1)h_1(2).$$
(5)
Hence, the effective 2-body Hamiltonian to be used in the truncated model space is given by

$$H_{\mathrm{eff}}^{(2)}=\underset{n=1}{\overset{A}{}}h_n+\underset{m<n=1}{\overset{A}{}}V_{\mathrm{eff}}(m,n).$$
(6)
It is possible to construct three-body, four-body, etc, interactions by similar methods. Indeed, this is essential if many-body properties are to be reproduced exactly. Note that the induced many-body interactions appear as a consequence of truncations even when the original nuclear Hamiltonian contains at most two-body terms as in Eq. (1). The usefulness of the approach depends, of course, on being able to stop at a low order in the cluster expansion while retaining acceptable accuracy.

### B CORE and Effective Theories

Effective Theories (ET’s) are now widely employed in many fields of physics including nuclear theory. The philosophical similarities between CORE and ET’s are obvious: in both cases one attempts to incorporate in a systematic fashion the contributions from short-ranged (high momentum/energy) physics into effective interactions and operators that are defined in a Hilbert space which excludes the high momentum/energy states. Indeed, Bogner and Kuo have recently examined the connection between CORE-like approaches, such as Lee-Suzuki  and Krenciglowa-Kuo  schemes on one hand and ET’s on the other. We are presently undertaking similar studies using concepts from CORE and renormalization-group (RG) methods outlined by Birse and collaborators . Our results will appear in a future publication.

In the present work we examine how CORE and the methods of ET can be applied to the shell-model problem in a complementary manner. The basic idea is as follows: the bare $`NN`$ interaction (see, e.g., Ref. ) is pathological in the sense that it is characterized by a very strong short-range repulsion. Obtaining accurate many-body energies using the bare interaction therefore requires huge bases. CORE and other similar methods are used to remedy this problem. We propose to solve this problem by first constructing an effective $`NN`$ potential of a convenient form for use in shell-model calculations. The main theoretical underpinning behind effective theories is that a long-wavelength particle should be insensitive to details of the short-range interaction. Thus, one adds new (smoother) interactions to “mimic” the effects of the unknown short-distance physics. Since no explicit short-range contributions will appear in the effective potentials, one expects smaller bases to be sufficient. Moreover, there is no reason that CORE cannot be used together with an effective potential to further refine the basis, a strategy that may result in substantially increased computational efficiency.

Our determination of the effective potential follows from the procedures outlined by Lepage  and those later adapted by Steele and Furnstahl  to treat the $`NN`$ interaction. For simplicity we work in one dimension and treat all particles as as spinless and distinguishable; they are assumed to have identical masses of $`m=939`$ MeV. We assume a bare $`NN`$ interaction with the same pathologies of a realistic interaction. That is, it is given by the sum of a strong short-range repulsive exponential and a medium-range attractive exponential:

$$V(x)=V_\mathrm{s}e^{m_\mathrm{s}|x|}+V_\mathrm{v}e^{m_\mathrm{v}|x|}.$$
(7)
The two masses were chosen to be equal to $`m_\mathrm{s}=400`$ MeV and $`m_\mathrm{v}=783`$ MeV, respectively. The strengths of the potentials, $`V_\mathrm{s}=506`$ MeV and $`V_\mathrm{v}=+1142.49`$ MeV, were chosen to give a binding energy and point root-mean-square (rms) radius for the symmetric (“deuteron”) state of $`E_\mathrm{b}=2.2245`$ MeV and $`r_{\mathrm{rms}}=1.875`$ fm and $`\sqrt{r^2}=1.875`$ fm, respectively. Note that for the $`NN`$ system an odd-parity state bound by $`0.180358`$ MeV is also found and that the three body ground state for this potential is bound by $`6.32`$ MeV.

Employing an option suggested by Lepage, we choose a gaussian form for our effective potential:

$$V_{\mathrm{eff}}(x)=\frac{1}{a}\left(c+d\frac{^2}{\xi ^2}+e\frac{^4}{\xi ^4}+\mathrm{}\right)\mathrm{exp}(\xi ^2);\xi x/a.$$
(8)
The range parameter has been fixed at $`a`$=1.16 fm as it reproduces the deuteron binding energy exactly. Yet, as indicated in Fig. 1, variations of this quantity—with the corresponding adjustment of the effective parameters to reproduce the low-energy phase shifts—have very little effect on the deuteron binding energy and, indeed, on any of our calculations of low-energy properties. Note that the parameters of the effective potential ($`c`$, $`d`$, and $`e`$) are fixed separately for the even and odd channels. They are constrained to reproduce the low-energy scattering phase shifts for the appropriate channel, as we explained below.

## III Results

Following the discussions in Ref. , we fit the parameters of $`V_{\mathrm{eff}}`$ to quantities which yield the effective range expansion when expanded in even powers of the wavenumber, $`k`$. In one dimension we find for the even channel

$$k\mathrm{tan}\delta _\mathrm{e}(k)=+\frac{1}{a_\mathrm{e}}+\frac{r_\mathrm{e}}{2}k^2+𝒪(k^4).$$
(9)
where $`\delta _\mathrm{e}`$ is the even-channel phase shift and $`a_\mathrm{e}`$ and $`r_\mathrm{e}`$ are the scattering-length and effective-range parameters, respectively. This apparently unorthodox form of the effective-range expansion emerges from the fact that, in one dimension, the symmetric wave-function does not vanish at the origin. For the odd channel we find the familiar expression

$$k\mathrm{cot}\delta _\mathrm{o}(k)=\frac{1}{a_\mathrm{o}}+\frac{r_\mathrm{o}}{2}k^2+𝒪(k^4).$$
(10)
We have found that including only the the first two terms of Eq. (8) for $`V_{\mathrm{eff}}`$ gives, as expected, an excellent reproduction of both scattering-lengths and effective-range parameters. The parameters of the effective potential are shown in Table I along with the binding energies obtained using the bare and effective potentials. The resulting $`V_{\mathrm{eff}}`$’s for various orders of Eq. (8) are compared with the bare potential in Fig. 2. Insofar as the low-energy properties of the theory are concern, the smooth $`V_{\mathrm{eff}}`$ is practically indistinguishable from the exact, and very pathological, $`V_{\mathrm{bare}}`$. As Table I indicates, the effective potentials reproduce the exact binding energies very accurately. That the effective potentials do so without possessing the “pathologies” of the bare $`NN`$ potential alluded to earlier is demonstrated in Fig. 3. The exact binding energy for the even-parity state is determined accurately by solving the Schrödinger equation using the Numerov algorithm. This energy is also determined by evaluating the matrix elements of the potential in a harmonic-oscillator basis and then diagonalizing the matrix. Fig. 3 shows that, as suggested above, such calculations using $`V_{\mathrm{eff}}`$ converge much faster than with the bare interaction where, eventually, round-off errors cause discrepancies to grow as the basis size increases.

Finally, we turn to the three-body system which represents the testing ground for our model. Using the bare potential with very large basis sizes ($`N100`$) and then extrapolating to $`N\mathrm{}`$, we find an exact ground-state energy for the three-body system of $`6.32`$ MeV. Table II summarizes approximate calculations of this quantity employing severely truncated bases; i.e., $`N30`$. (Note that the number of three-body states scales roughly like $`N^3`$.) The same information also appears in Fig. 4. Four types of calculations are reported: calculations using either the bare or effective $`NN`$ interaction both with and without CORE. The results show clearly that CORE corrections are important, especially for small values of $`N_{\mathrm{basis}}`$. Yet these corrections are significantly more important for the bare potential than for the effective one. This is to be expected as, in constructing the effective interaction, we have already “softened” the potential to such a degree that most of its large momentum components are absent. Put another way, the smoothed effective potential does not explicitly contain the “pathological” short-range contributions which CORE must treat. At the same time, the convergence to the exact ground-state energy is much faster overall for the effective interaction. Note that this is true even when one considers that the exact binding energy of the three-body system for the effective potential ($`6.89`$ MeV) differs slightly from the exact three-body energy for the bare interaction ($`6.32`$ MeV). Presumably, this difference is due to the omission of induced three-body interactions. This relatively small discrepancy suggests that such three-body interactions are tractably small.

## IV Conclusions

A contractor renormalization (CORE) approach has been applied in concert with modern versions of effective theories (ET’s) to the nuclear many-body problem. CORE is a “model-space” approach similar to existing techniques, such as the Lee-Suzuki method, that have been applied already with some success to the nuclear many-body problem. The central concept behind CORE is a similarity transformation responsible for generating a renormalized Hamiltonian with identical low-energy properties to the exact Hamiltonian. The initial study presented here has been limited, for simplicity, to one spatial dimension. Yet the bare nucleon-nucleon potential employed throughout the paper possesses the same “pathologies”—the short-range repulsive core—as any realistic interaction. To remove these pathologies we proceeded as follows. First, the high-momentum components of the potential were removed through the introduction of an effective potential that was constrained – in effect – by low-energy scattering phase shifts. Such an effective potential predicts a two-body (“deuteron”) bound-state energy identical to that obtained with the bare interaction. Second, in calculating the ground-state energy of the three-body system we depended on CORE to compensate for severely truncating the size of the basis employed in the calculation. In computing the ground-state energy of the three-body system we observed that, as in most implementations of model-space techniques, CORE corrections improve significantly the convergence towards the exact ground-state energy. However, use of the effective two-body potential significantly accelerated the rate of convergence while simultaneously reducing the magnitude of the CORE corrections.

Overall these results are gratifying as they suggest that ET’s can be profitably combined with more traditional model-space methods for use with the nuclear shell model. For us, much remains to be done. Apart from the obvious extension to more realistic systems, we must determine how to optimally mesh ET and CORE as well as how to correct operators in a manner consistent with the construction of the effective potential. Also, a preliminary comparison between CORE and a commonly used  approach for evaluating effective interactions—the Lee-Suzuki method—has been made by us. We find that, at least for the one-dimensional three-body system treated here, the two methods give similar results, with ground-state energies computed using CORE being slightly—but consistently—closer to the exact value. (See Fig.5 for relevant comparison with relatively small basis sizes.) In a future study we intend to elucidate the relationships between CORE, Lee-Suzuki, and other related approaches such as the Krenciglowa-Kuo scheme . The formalism developed by Andreozzi  may prove useful in this regard. Finally, as mentioned earlier, it would be very useful and informative to find a set of RG-like equations which would lead directly to a sufficiently accurate version of $`V_{\mathrm{eff}}`$. Illuminating early steps in this direction may be found in Refs. .

###### Acknowledgements.

We are grateful to the Institute for Nuclear Theory at the University of Washington for their hospitality and support during the workshop on Effective Field Theories and Effective Interactions. It is also a pleasure to acknowledge very useful conversations with B. Barrett, W. Haxton, W. Leidemann, P. Navratil, and G. Orlandini. This work was supported by the DOE under Contracts Nos. DE-FC05-85ER250000, DE-FG05-92ER40750 and DE-FG03-93ER40774.


text:
# Thermal emission from low-field neutron stars

## 1 Introduction

Modern X-ray observatories detect the thermal surface emission of a number of neutron stars. In analogy to the classic model atmosphere analysis of normal stars, such observations permit the direct measurement of fundamental neutron star properties, such as their effective temperatures, atmospheric abundances and surface gravities. The latter point is especially of great importance, as an accurate measurement of the surface gravity of a neutron star is (with the distance known) equivalent to a measurement of its mass/radius ratio. The observational confirmation/rejection of the predicted neutron star mass-radius relations is a fundamental test of our understanding of the physics of matter above nuclear densities.

The first neutron star model atmospheres involving realistic opacities were computed by Romani (1987, henceforth R87), using atomic data from the Los Alamos Opacity Library and employing a simple angle-averaged radiation transfer. As a major result, R87 could show that the thermal emission of neutron stars differs substantially from a Planck spectrum. For low-metallicity (helium) atmospheres, the emitted spectrum is harder than the corresponding blackbody spectrum. The spectra emitted from high-metallicity atmospheres (carbon, oxygen or iron) are closer to a blackbody distribution, but show strong absorption structures in the energy range observable with X-ray telescopes.

Rajagopal & Romani (1996, henceforth RR96) computed hydrogen, solar abundance, and iron model atmospheres for low-field neutron stars, using improved opacity and equation of state data from the OPAL project (Iglesias & Rogers 1996), but employing the same radiation transfer as R87. Contemporaneously, a similar set of low-field atmospheres, partially based on the same atomic OPAL data, but employing a more sophisticated radiation transfer, was presented by Zavlin et al. (1996, henceforth ZPS96). The RR96 and ZPS96 models, broadly confirming the results of R87, were applied to X-ray observations of presumably low-field neutron stars, i.e. millisecond pulsars (RR96, Zavlin & Pavlov 1998), isolated neutron stars (Pavlov et al. 1996), and transiently accreting neutron stars in quiescent LMXBs (e.g., Rutledge et al. 1999, 2000, 2001).

In the case of strongly magnetic neutron stars ($`B10^{12}10^{14}`$ G), both the observation and the computation of the thermal surface emission is significantly more difficult, as non-thermal magnetospheric radiation dilutes the surface emission and as little atomic data is readily available (for a review of the properties of matter in strong magnetic fields, see Lai 2001). Magnetic models for hydrogen and iron atmospheres were computed by Pavlov et al. (1995) and Rajagopal et al. (1997), respectively.

The wealth of high-quality neutron star observations expected from Chandra and XMM justifies the computation of modern model atmospheres. In this paper, we present new model atmosphere grids for low-field neutron stars that overcome a number of shortcomings and errors in the RR96 and ZPS96 calculations and that are made available to the community.

## 2 The model atmospheres

The new grid of neutron star atmospheres was computed combining a fully frequency- and angle-dependent radiation transfer code (Gänsicke et al. 1995) and the OPAL opacities from RR96. The employed radiation transfer is part of a model atmosphere code that has previously been used for the computation of white dwarf atmospheres which have been widely applied to observations of cataclysmic variables and white dwarf/main-sequence binary stars.

The model atmospheres were created under the following assumptions. (1) Plane-parallel geometry. The atmosphere of a neutron star has a vertical extension of a few centimeters at most, compared to a radius of $`10`$ km. Its curvature is hence completely negligible. (2) Hydrostatic equilibrium. For the extreme surface gravity of neutron stars and the relatively low temperatures that we consider here, the atmosphere is static. The situation is different in accreting neutron stars during X-ray bursts, where the nuclear burning atmosphere considerably expands on short time scales. (3) Radiative equilibrium. The atmosphere contains no source of energy, but acts only as a blanket, through which the thermal energy of the underlying core leaks out. In Sect. 3, we will discuss how abandoning this assumption will impact the spectrum of the neutron star. (4) Local thermal equilibrium (LTE). For the high densities and the low temperatures in the considered atmospheres collisional ion-ion interactions dominate over interactions between matter and the radiation field throughout most of the atmosphere. Deviations from LTE in the outermost tenuous layers of the atmosphere may somewhat affect the depth of the absorption lines, but considering the overall level of uncertainty in the atomic physics involved, and the relatively poor quality of the available spectroscopy, the assumption of LTE appears to be a reasonable approximation. Below, we discuss a quantitative comparison to a NLTE model.

With these approximations, the computation of a model atmosphere can be separated into two independent parts: calculating the structure of the atmosphere from the integration of the hydrostatic equation (Sect. 2.1) and solving the radiation transfer (Sect. 2.2). The temperature structure of the atmosphere, $`T(z)`$, is a free parameter in the overall process and is adjusted iteratively until radiative and hydrostatic equilibrium are satisfied (Sect. 2.3).

### 2.1 Atmosphere structure

The vertical scale of the atmosphere is transformed from the geometrical depth $`z`$ to the optical depth $`\tau `$ by $`dz=d\tau /\rho \chi `$, with $`\chi `$ the total mean opacity (see below) and $`\rho `$ the density. The structure of the atmosphere is then obtained by integrating the equation of hydrostatic equilibrium

$$\frac{dP_g}{d\tau }=\frac{g}{\chi }$$
(1)
where $`P_g`$ is the gas pressure, from an outer boundary $`\tau _{\mathrm{min}}`$ to an inner boundary $`\tau _{\mathrm{max}}`$. Considering the extremely high surface gravity encountered in a neutron star atmosphere, we neglect the effects of radiation pressure. We compute a starting value for the integration of eq. (1), $`P_g(\tau =\tau _{\mathrm{min}})`$, assuming that the outer layers of the atmosphere are isothermal and that the degree of ionization is constant. With $`\chi =\mathrm{const}`$,

$$P_g(\tau _{\mathrm{min}})=\frac{g}{\chi (\tau _{\mathrm{min}})}\tau _{\mathrm{min}}.$$
(2)
$`\chi `$ being itself a function of $`(T,P_g)`$, we use $`\chi (\tau _{\mathrm{min}})=1`$ as an estimate and solve eq. (2) with a Newton-Raphson method. From $`P_g(\tau =\tau _{\mathrm{min}})`$ we proceed with the integration of eq. (1) to $`\tau _{\mathrm{max}}`$.

This integration involves at each depth the evaluation of the equation of state (EOS) as well as the computation of the absorption and scattering coefficients. For the EOS and the radiative opacities, we use the OPAL data (Iglesias & Rogers 1996) from RR96, that cover three different chemical compositions: hydrogen, solar abundances, and iron. The opacities are tabulated as a function of the temperature $`T`$, $`R=\rho /T_6^3`$ with $`T_6`$ the temperature in $`10^6`$ K, and $`u=E_\gamma /k_BT`$, with $`E_\gamma `$ the photon energy (for additional details on the OPAL opacity tables, see RR96). For a given depth in the atmosphere, defining $`T`$ and $`\rho `$ and at a given energy, the radiative absorption coefficient $`\kappa _{\mathrm{rad}}(T,\rho ,E)`$ is calculated from the opacity tables using a bilinear interpolation in $`\mathrm{log}T`$ and $`\mathrm{log}R`$ and a linear interpolation in $`E`$. In the construction of the model atmospheres, the radiative absorption coefficient $`\kappa _{\mathrm{rad}}(T,\rho ,\delta E)`$ for a finite energy interval $`\delta E`$ is given by the harmonic mean of $`n`$ individual evaluations of $`\kappa _{\mathrm{rad}}(T,\rho ,E_n)`$, with $`E_n`$ within $`\delta E`$, and $`n`$ sufficiently large to sample well the resolution of the OPAL tables. Thomson scattering is an important source of opacity only in low-Z atmospheres at high temperatures and high energies. The OPAL tables implicitly include Thomson scattering as radiative opacity. $`\kappa _{\mathrm{rad}}`$ is, hence, the sum of true absorption plus Thomson scattering. We will come back to this issue in Sect. 2.2. Finally, we compute a radiative Rosseland mean opacity $`\overline{\kappa }_{\mathrm{Ross}}`$ from the energy dependent $`\kappa _{\mathrm{rad}}(E)`$.

Following RR96, we include a conductive opacity $`\kappa _{\mathrm{con}}=(2.5\times 10^4/n_e)\times (Z^2/A)(T_7^{1/2}/\rho )\mathrm{cm}^2\mathrm{g}^1`$ (Cox & Giuli 1968) to account for energy transport by electron conduction, with $`n_e`$ the number of electrons per a.m.u., $`A`$ the average atomic mass, $`Z=An_e`$ the average ionic charge, and $`T_7`$ the temperature in $`10^7`$ K. The radiative (Rosseland) opacity and the conductive opacity are harmonically added. The total opacity is, hence, given by $`\chi =(\overline{\kappa }_{\mathrm{Ross}}^1+\kappa _{\mathrm{con}}^1)^1`$.

Once eq. (1) is integrated to $`\tau _{\mathrm{max}}`$, a complete macroscopic description of the atmosphere structure is at hand.

### 2.2 Radiation transfer

For the assumption given above, the radiation transfer equation takes the form

$$\mu \frac{dI_\nu (\mu )}{d\tau _\nu }=I_\nu (\mu )S_\nu ,$$
(3)
with $`\mu =\mathrm{cos}\theta `$ the cosine of the angle between the normal to the atmosphere and the direction of the considered radiation beam, $`\tau _\nu `$ the frequency-dependent optical depth, $`I_\nu (\mu )`$ the frequency- and angle-dependent specific intensity, and $`S_\nu `$ the source function.

As mentioned in Sect. 2.1, Thomson scattering is implicitely included in the OPAL opacity tables. Considering that the contribution of Thomson scattering to the total opacity is negligible in high-Z atmospheres, we restrict an explicit treatment of Thomson scattering to our hydrogen models (Sect. 2). We compute Thomson scattering, $`\sigma _e=N_e\mathrm{\hspace{0.17em}6.65}\times 10^{25}/\rho `$, but taking into account the cross-section reduction by collective effects (Boercker 1987), with $`N_e`$ the electron density in $`\mathrm{cm}^3`$, calculated under the assumption of full ionisazion, and the true frequency-dependent absorption $`\kappa _\nu =\kappa _{\mathrm{rad}}(\nu )\sigma _e`$, where $`\kappa _{\mathrm{rad}}(\nu )`$ is the radiative frequency-dependent OPAL opacity from Sect. 2.2.

With the isotropic and, to first order approximation, coherent Thomson scattering term, the source function is given by

$$S_\nu =\frac{\kappa _\nu }{\kappa _\nu +\sigma _\mathrm{e}}B_\nu +\frac{\sigma _\mathrm{e}}{\kappa _\nu +\sigma _\mathrm{e}}J_\nu $$
(4)
with $`B_\nu `$ the Planck function and $`J_\nu `$ the mean intensity. We use as boundary conditions for eq. (3) $`I(\mu <0)=0`$ at $`\tau _{\mathrm{min}}`$ and the diffusion approximation at $`\tau _{\mathrm{max}}`$. Equation (3) is solved using the well-documented Rybicki method (e.g. Mihalas 1978).

For the high-Z composition atmospheres (solar abundances and iron), we treat Thomson scattering as true absorption, directly using the OPAL $`\kappa _{\mathrm{rad}}`$.

From the resulting angle-dependent specific intensities, the flux at the surface of the atmosphere is computed,

$$F_\nu =\underset{0}{\overset{1}{}}\mu I_\nu (\mu )𝑑\mu |_{\tau =\tau _{\mathrm{min}}}.$$
(5)
Radiative and conductive equilibrium requires

$$\underset{0}{\overset{\mathrm{}}{}}F_\nu 𝑑\nu +F_{\mathrm{con}}=\frac{\sigma T_{\mathrm{ns}}^{}{}_{}{}^{4}}{\pi }$$
(6)
at each depth in the atmosphere, where $`\sigma `$ is the Stefan-Boltzmann constant, $`F_{\mathrm{cond}}=(16\sigma T^3)/(3\rho \kappa _{\mathrm{con}})dT/dz`$ is the conductive energy flux, as detailed in RR96, and $`T_{\mathrm{ns}}`$ is the effective temperature of the neutron star. If the deviations from the equilibrium (6) are larger than a given error, correction terms $`\mathrm{\Delta }T(\tau )`$ are computed from linearization of the angle-integrated form of eq. (3) (e.g. Feautrier 1964; Auer & Mihalas 1969, 1970; Gustafsson 1971).

### 2.3 Model atmosphere computation

The actual computation of a model atmosphere and the emergent spectrum is done in two separate stages.

(i) Starting from the temperature structure of a grey atmosphere (Chandrasekhar 1944), we iteratively construct (Sect. 2.1 and Sect. 2.2) an atmosphere structure that satisfies radiative and hydrostatic equilibrium. At this stage, we use an energy grid of 1000 logarithmic equidistant bins covering the range $`10^{3.4}10`$ keV. For each energy bin, 30 evaluations of the OPAL opacity tables are harmonically added for the computation of the radiative absorption coefficient. The atmosphere structure is constructed on an optical depth grid with 100 points covering $`\tau _R=10^65\times 10^3`$ (except for the coldest models, which were computed on a grid covering $`\tau _R=10^610^3`$ due to limitations in the available OPAL data). The angular dependency of the radiation field is sampled by three Gaussian points in $`\mu `$. Radiative equilibrium (eq. 6 is satisfied to better then $`10^5`$ at each optical depth) is reached after $`515`$ iterations. Generally, the hydrogen atmospheres converged the slowest, as the steep drop of the opacity $`\nu ^3`$ makes the atmosphere transparent for high-energy photons, thus radiatively coupling the deep hot layers with the surface layers. Electron conduction carries a few percent of the total flux at the largest optical depths (Rosseland optical depths $`\tau _{\mathrm{Ross}}>100`$) in the coldest iron and solar atmospheres. For $`\mathrm{log}(T_{\mathrm{ns}})>5.75`$, the conductive flux is less than 1%, which is in good agreement with the estimates of ZPS96. Including the conductive opacity in our atmosphere models has, hence, no noticeable effect on the emergent spectra.

(ii) From a given atmosphere structure, we recompute the emergent spectrum solving the radiation transfer (Sect 2.2) on a much finer energy grid (10 000 bins, 10 OPAL evaluations per bin), and using 6 Gaussian points in $`\mu `$ to resolve the angular dependence of the specific intensity. The higher number of angular points has practically no influence on the angle-averaged flux from the neutron star surface, but may be used to account for limb darkening in the computation of spectra from neutron stars with non-homogeneous temperature distributions.

### 2.4 The model spectra

We computed grids of model atmospheres for the three different compositions covered by the OPAL opacity/EOS tables of RR96: hydrogen, solar abundances, and iron. We restricted ourselves to the canonical neutron star configuration, $`M_{\mathrm{ns}}=1.4`$$`M_{}`$ and $`R_{\mathrm{ns}}=10`$ km, corresponding to a surface gravitational acceleration of $`\mathrm{log}g=14.386`$. The three grids cover the temperature range $`\mathrm{log}T_{\mathrm{ns}}=5.16.5`$ in steps of 0.05, i.e. a total of 29 spectra per abundance grid.111The angle-averaged fluxes are available from CDS, for the angle-dependent specific intensities, or models covering additional temperatures/surface gravities, please contact the authors.

Figure 1 shows the emergent spectra for the three different compositions. As already evident in the previous neutron star atmosphere calculations (R87, RR96, ZPS96), the model spectra differ significantly from blackbody distributions. The pure hydrogen models show a strong flux excess over the blackbody distributions at energies above the peak flux. The free-free and bound-free opacity in these atmospheres drops off rapidly at high energies, leaving Thomson scattering as the dominant interaction between the atmosphere matter and the radiation field. As a consequence, the atmosphere is highly transparent to the hard X-ray photons from deep hot layers. In contrast to the hydrogen models, the iron and solar abundance models are overall closer to the blackbody distributions because of the milder energy dependence of their opacities. They show, however, a substantial amount of absorption structures, lines and edges, which are especially pronounced in the $`0.110`$ keV range well observable with most of the present and past X-ray satellites.

Noticeable absorption structures with equivalent widths of $``$ Å are present in the iron spectra of moderately cool neutron stars also in the ultraviolet. The (non-magnetic) hydrogen models contain only the extremely pressure broadened $`\mathrm{Ly}\alpha `$ line. In a magnetic field $`\mathrm{Ly}\alpha `$ splits in three Zeeman-components, with the $`\sigma ^{+/}`$ components asymmetrically shifted by $`100`$ Å for $`B10^9`$ G (Ruder et al. 1994). Because the hydrogen in these hot atmospheres is ionised to a large extent, the $`\mathrm{Ly}\alpha `$ Zeeman components are expected to be rather weak, but they may possibly be detected with future large aperture ultraviolet telecopes, allowing a direct measurement of the magnetic field strength. Figure 2 compares the redshifted222 $`z=(12GM_{\mathrm{ns}}/R_{\mathrm{ns}}c^2)^{1/2}1=0.306`$ for an assumed $`R_{\mathrm{ns}}=10`$ km and $`M_{\mathrm{ns}}=1.4M_{}`$. iron and hydrogen spectra in the far ultraviolet.

Figure 3 shows the angle-dependent intensity $`I_\nu (\mu )`$ for a $`T_{\mathrm{ns}}=10^6`$ K iron atmosphere. It is apparent that the emission from the neutron star surface is highly anisotropic, as previously discussed by ZPS96. Calculations of the emission of a neutron star with a non-homogeneous surface temperature distribution must take this anisotropy into consideration. An application of the angle-dependent intensities will be discussed elsewhere.

### 2.5 Comparison to earlier neutron star atmosphere models

#### 2.5.1 Rajagopal & Romani (1996)

The new model spectra agree quite well with the models of RR96. The only systematic differences are found for the hydrogen models, where the new models have a significantly lower flux at energies $`1`$ keV. This difference is due to the inclusion of Thomson scattering in the new models, which increases the opacities at large optical depths in the hydrogen atmospheres. For the heavy element atmospheres, this effect is negligible, as the opacity is dominated by bound-bound and bound-free absorption.

Small (5-15 %) differences in the emergent flux of the heavy element atmospheres are found in the low-temperature ($`\mathrm{log}T_{\mathrm{ns}}5.75`$) models in the energy range 1–10 keV. The specific intensities at these high energies are $`1525`$ orders of magnitude below the peak intensities, and we believe that the differences between the new models and those of RR96 are due to the different numerical treatment of the radiation transfer. At higher energies, the agreement between the two different model generations is better than 1%, except in strong lines where different energy sampling naturally leads to somewhat larger discrepancies.

#### 2.5.2 Zavlin et al. (1996)

ZPS96 computed grids of non-magnetic neutron star LTE model atmospheres for three different compositions, hydrogen, helium, and iron. The hydrogen and helium opacities and non-ideal EOS employed in these models were calculated by ZPS96 using the occupation probability formalism described by Hummer & Mihalas (1988). The iron models were computed using the same OPAL data used also by RR96 and in the present paper. The iron models show a significant difference both in the overall shape of the continuum as well as in the details of the absorption features. We traced this to an erroneous use of the opacity grid in ZPS96. This has been confirmed (Pavlov, priv. comm.)333 After we alerted G. Pavlov to their errors in employing the opacities, their group has computed revised models in reasonable agreement with other work. An example revised spectrum appears in (Pavlov & Zavlin 2000).. The ZPS96 hydrogen model is somewhat harder than both our new spectrum and that of RR96. The flux excess of the ZPS96 model (or the flux deficiency of our model) is not fully understood, but likely stems from small differences in the absorption/scattering opacities.

#### 2.5.3 Werner & Deetjen (2000)

Werner & Deetjen (2000) presented the first NLTE calculation for a neutron star atmosphere using bound-free opacities from the Opacity Project (Seaton et al. 1994) and explicitly treated line blanketing by millions of lines from the Kurucz & Bell (1995) line list. Figure 4 compares the $`T_{\mathrm{ns}}=\mathrm{200\hspace{0.17em}000}`$ K pure iron NLTE model from Werner & Deetjen (2000) with a corresponding model from our LTE calculations. While the NLTE model clearly shows more detail because of the much larger number of considered transitions, the overall agreement between the two models is very good. Werner & Deetjen (2000) quote a flux difference between their LTE and NLTE models of $`10`$ % in the line cores, and much less in the continuum. For higher temperatures, the NLTE effects might be stronger, but so far no detailed atomic data are available.

## 3 Beyond “classic” stellar atmospheres

Part of the ongoing effort in computing neutron star atmospheres has been carried out for the analysis of X-ray observations of isolated neutron stars (e.g. Pons et al. 2002). However, recent high-resolution X-ray spectroscopy of the brightest isolated neutron star candidates RX J1856.5–3754 (Burwitz et al. 2001) apparently excludes all the neutron star atmosphere models presented so far: while the overall spectral energy distribution (SED) rules out low-$`Z`$ atmospheres, the featureless XMM and Chandra spectra are inconsistent with the strong structures expected from a heavy element atmosphere.

As in all previous calculations (R87, RR96, ZPS96, Werner & Deetjen 2000), the model atmospheres presented here assume radiative equilibrium. It is worth remembering that when the model fits suggest that the flux is dominated by hot polar caps with a small fraction of the full neutron star surface area, some sort of local heating is likely implicated. For instance, anisotropic interior conductivities can produce smooth variations in the surface effective temperture; e.g. Heyl & Hernquist (1998). For slow, long $`P`$ neutron stars, local heating can be caused by low $`\dot{M}`$ accretion. Zampieri et al. (1995) and Zane et al. (2000) have produced emergent spectra for ionized H stars accreting at low rates. With proton stopping depths of $`20`$$`\mathrm{g}\mathrm{cm}^2`$, much of the energy is deposited fairly deep in the atmosphere, but shocks in the accretion flow can apparently heat the outer layers, providing excess flux on the Rayleigh-Jeans tail. For an active magnetosphere, local heating without accretion arises from precipitating $`e^\pm `$ (Arons 1981) or from illumination by energetic (polar cap or outer gap) photons.

An analogous situation – pole caps heated by accretion – is well-documented in polars, a subtype of cataclysmic variables in which a strongly magnetic white dwarf accretes from a late type donor star. In theses systems, the white dwarf atmosphere near the magnetic poles is heated by strong irradiation with cyclotron radiation and thermal bremsstrahlung from a stand-off shock. As a consequence of this accretion heating, the spectra of these heated pole caps are almost devoid of spectral structures. The Lyman lines, typically very strong absorption features in the ultraviolet spectra of white dwarfs, are almost completely flattened out (Gänsicke et al. 1998), and also the X-ray spectra of polars are very close to blackbody distributions, unlike the expected emission from a high-gravity photosphere strongly enriched with heavy elements (Mauche 1999). While the physical processes involved in heating the pole caps of accreting white dwarfs and accreting neutron stars differ markedly, the observed phenomenology is very similar: RX J0720.4–3125, another bright isolated neutron star shows a quasi-sinusoidal X-ray light curve suggesting the presence of large heated pole caps, whereas the X-ray spectrum contains no significant absorption structure (Haberl et al. 1997; Paerels et al. 2001) reminiscent of the ultraviolet observations of the heated white dwarf in the polar AM Her (Gänsicke et al. 1998). It appears likely, therefore, that heating effects may affect the line spectra of heavy element neutron star atmospheres.

As a simple illustration, we ignored the details of particular physical heating effects and computed a number of spectra under the assumption that an energy $`\alpha \sigma T_{\mathrm{ns}}^4`$ is deposited at a characteristic depth in the atmosphere, $`\tau _{\mathrm{dep}}`$ (Fig. 5). In practice, the emergent spectrum of a “heated” neutron star is computed by (1) modifying the thermal structure of an undisturbed atmosphere $`T(\tau )`$ to account for the deposited energy, (2) recomputing the atmosphere structure based on the modified temperature profile (Sect. 2.1), and finally solving the radiation tansfer within the “heated” atmosphere (Sect. 2.2). Figure 6 shows the emergent fluxes from a $`T_{\mathrm{ns}}=10^6`$ K iron atmosphere for two different sets of parameters, $`\alpha =0.82`$ and $`\tau _{\mathrm{dep}}=1`$, equivalent to $`\mathrm{\Sigma }=3.2\times 10^3`$$`\mathrm{g}\mathrm{cm}^2`$, and $`\alpha =0.035`$ and $`\tau _{\mathrm{dep}}=0.01`$, equivalent to $`\mathrm{\Sigma }=2.2\times 10^4`$$`\mathrm{g}\mathrm{cm}^2`$.

Two effects are evident. The flatter temperature in the outer layers of the heated atmospheres supresses the strong absorption edge around 1 keV, returning the Wien tail colour quite close to that of a simple blackbody. Further, the equivalent widths of the lines are, of course, strongly affected. Heating at $`\tau 1`$ strongly supresses the $``$ keV absorption features; surface heating in fact drives the line features on the Rayleigh-Jeans tail into emission.

We stress that this illustrative approach is by no means a self-consistent model of a heated neutron star atmosphere, as we include neither any physical assumption on the actual heating mechanims, e.g., there are currently no plausbile physical processes known which are capable of heating such shallow layers, nor a proper energy balance of the heating/cooling processes. It is merely intended as a possible explanation for the apparent paradox that the SEDs of isolated neutron stars (in particular RX J1856.5–3754) are best-fitted with heavy-element model spectra (Pavlov et al. 1996; this paper; Burwitz et al. 2001; Pons et al. 2002), whereas their X-ray spectra show at best very weak spectral structures.

More detailed future work will need to carefully examine the efficiency of the various possible heating mechanisms, in particular to establish the depth-dependent energy deposition, and to account properly for the energy balance in the neutron star atmosphere.

## 4 Conclusions

We have presented a set of angle-dependant low field neutron star atmosphere spectra. These radiative and hydrostatic equilibrium atmospheres should be useful to researchers pursuing exploratory fits of soft X-ray/UV/optical data. However, it is important to inject a note of caution into the discussion. At a minimum, non-uniform surface temperature, combined with limb-darkening and gravitational focussing will have subtle, but important effects which can be modeled using the present grid.

Many other physical effects can, of course, strongly perturb neutron star spectra. The effects on heavy element absorption lines and edges can be particularly strong. Magnetic fields characteristic of young radio pulsars will have a dramatic effect on the opacities and the emergent spectra. When considering the strength of the absorption lines in spectra from magnetic atmospheres with heavy elements (Rajagopal et al. 1997) it is important to remember that the line energies depend sensitively on the B-field. Even simple dipole $`B`$ variation across a polar cap can shift line energies by $`10`$%, strongly decreasing the equivalent width of spectral features in phase-averaged spectra. Convection has also been considered (eg. RR96), although even small magnetic fields may supress convective transport and keep atmospheres close to the hydrostatic solution.

On the whole, while we feel that our atmosphere grid should be of use for fitting of observed data; complications are certainly expected. In particular strong $`B`$-field variations and surface re-heating can decrease the equivalent width of heavy element line features in neutron stars with active magnetosphers; this may explain in part the difficulty in finding such features in early Chandra/XMM data (e.g. Paerels et al. 2001; Burwitz et al. 2001).

###### Acknowledgements.

BTG was supported in part by a travel grant of the Deutscher Akademischer Auslandsdienst (PKZ: D/99/08935) and by the DLR under grant 50 OR 99 03 6. We thank K. Werner and G. Pavlov for providing their model spectra for a quantitative comparison, and B. Rutledge for comments on an earlier draft. We thank the referee, Dr. Zavlin, for useful comments.


text:
# The Quest for Light Sea Quarks: Algorithms for the Future

## 1 Introduction

The past six years have witnessed a number of large simulations of full QCD with two degenerate sea quarks and Wilson like actions. All of them were based on the hybrid Monte Carlo algorithm (HMC) which has been systematically improved over the years. These simulations, however, albeit being milestones for lattice QCD in the pre-Teracomputing age , suffer from severe restrictions that prevent them to be truly realistic: we need to go to lighter quark masses and operate with three flavors with masses eventually approaching the physical mass.

This situation has been motivating research and development of alternative algorithms of the multibosonic variety (MBA) that promise in principle to overcome some of the restrictions of HMC mentioned. It has been shown that MBA can indeed be geared to reach the efficiency of HMC  and encouraging results have been achieved with MBA for the treatment of finite density QCD , supersymmetric field theories  as well as the case of three dynamical fermion flavors .

In this contribution we wish to describe our benchmarking of MBA, performed in the quest to achieve lighter quark masses than previously possible in the SESAM and other full QCD simulations. For details regarding the HMC algorithm used see ref. .

## 2 Gauging performance

The efficiency of an algorithm is most suitably quoted as the computational cost, $`E_{\text{ind}}`$, for producing a statistically independent vacuum gauge field configuration within the Markov process, expressed in units of the number of necessary Dirac-matrix vector $`M\varphi `$ multiplies:

$$E_{\text{ind}}=2\tau _{\text{int}}E_{\text{traj}},$$
(1)
with $`\tau _{\text{int}}`$ being the integrated autocorrelation time measured in units of update sequences called “trajectories” each having a cost of $`E_{\text{traj}}`$.

Being particularly interested in the quark mass dependence of the algorithms, we chose for our analysis two different quark mass settings, as given in table 1. Note that both quark masses are from the regime of current full QCD simulations.

## 3 Defining the MBA algorithm

The parameter space of MBA algorithms is large w.r.t. the HMC situation: MBA optimization proceeds (a) by finding a not too high order polynomial approximating the inverse of the fermion matrix at the given working point, (b) by devising an efficient “trajectory” in form of an updating sequence for the degrees of freedom and (c) by exploiting a method to correct for the remaining systematic error. Today, all MBAs use a noisy correction step for (c). For a detailed discussion see e.g. ref. .

For the construction of appropriate polynomials we followed the prescription proposed in ref.  which employs the GMRES algorithm to find a non-hermitian approximation of the inverse Wilson matrix with ultraviolet-filtering on a small sample of thermalized gauge field configurations. This was done to minimize the polynomial order and hence the number of boson fields. Quadratically optimized polynomials  need about $`23`$ times more boson fields for similar acceptance rates.

The other key element for the optimization of an MBA is the choice of “trajectory”, $`T`$, i.e., the basic updating pattern. For that matter, one has to devise an appropriate blend of gauge and boson field updates. We have chosen the “trajectory” to be

$$T=CO_g\left(O_bO_g\right)^5H_b\left(O_gO_b\right)^5O_g,$$
(2)
with $`C`$ being the final accept/reject correction step, $`O_g`$ ($`O_b`$) a gauge (boson) field overrelaxation sweep and $`H_b`$ a boson field global quasi-heat bath .

The number of boson fields is tuned such as to ensure a good acceptance rate in the ball-park of $`6070\%`$. As a result we obtained the values quoted in table 2.

## 4 Results

The determination of accurate autocorrelation times, $`\tau _{\text{int}}`$, needs very long Monte Carlo time histories. The number of trajectories underlying the present analysis together with our results for the efficiencies, $`E_{\text{ind}}`$, are collected in table 4. At $`\kappa =0.160`$, we have also used the topological charge, $`Q`$, as a monitor for autocorrelations.

As a check for our simulation code we have computed the average plaquette and (non-singlet) pseudo-scalar and vector meson mass values, as shown in table 3. They agree within errors.

$`E_{\text{ind}}(\text{plaq})`$ has been obtained from direct integration of the plaquette autocorrelation functions while for the masses the jackknife blocking technique has been applied. The latter method only allows for a rather crude estimate of autocorrelation times, so that we concentrate mainly on the plaquette in the following, where we have safe control over statistical errors:

The statistics at the point $`\kappa =0.159`$ covers about $`320`$ autocorrelation times for the HMC and $`140`$ for the MBA. At the working point $`\kappa =0.160`$, the HMC series has a length of $`140`$ and the MBA of $`160`$ times $`\tau _{\text{int}}`$. Hence, the statistics is sufficient for the results to be conclusive.

From table 4 one finds that $`E_{\text{ind}}(\text{plaq})`$ is about equal for MBA and HMC at the larger quark mass. The total effort in the case of the HMC increases by a factor of $`8`$ when stepping down in quark mass, while at the same time the cost for MBA ‘only’ rises by a factor of $`2.7`$. Consequently, the MBA has become almost a factor of three more efficient than HMC.

## 5 Conclusions

We have demonstrated in a realistic setting that MBA-type sampling techniques are not at all inferior to the state-of-the-art HMC techniques and appear to show superior scaling behavior in quark masses. We have established here a better scaling behavior of MBA, leading to a gain factor of almost three in favor of the particular MBA variant at currently reachable quark masses.

It remains to be seen to what extent this encouraging result will pertain in the regime of smaller quark masses. The ultimate aim is to push full QCD simulation with Wilson like fermions below the point $`m_\pi /m_\rho 0.5`$.

Another point to check is whether the MBA variant considered in this paper with its non-hermitian approximation will remain operational for practical simulations in the deeper chiral regime. Because the GMRES method suffers from numerical instabilities for increasing orders $`n`$, one might have to take recourse — even on $`32\times 16^3`$ lattices — to $`128`$-bit precision. A more extended MBA study in this direction is in preparation .

## Acknowledgments

N.E. is supported under DFG grant Li701/3-1, M.P. by the European Community’s Human potential program under HPRN-CT-2000-00145 Hadrons/Lattice QCD, and W.S. by the DFG Graduiertenkolleg “Feldtheoretische Methoden in der Elementarteilchentheorie und Statistischen Physik”. The numerical productions were run on the APE-100 systems installed at NIC Zeuthen.


text:
# Enhanced Gilbert Damping in Thin Ferromagnetic Films

## Abstract

The precession of the magnetization of a ferromagnet is shown to transfer spins into adjacent normal metal layers. This “pumping”of spins slows down the precession corresponding to an enhanced Gilbert damping constant in the Landau-Lifshitz equation. The damping is expressed in terms of the scattering matrix of the ferromagnetic layer, which is accessible to model and first-principles calculations. Our estimates for permalloy thin films explain the trends observed in recent experiments.

The magnetization dynamics of a bulk ferromagnet is well described by the phenomenological Landau-Lifshitz-Gilbert (LLG) equation Gilbert (1955)

$$\frac{d𝐦}{dt}=\gamma 𝐦\times 𝐇_{\text{eff}}+\alpha 𝐦\times \frac{d𝐦}{dt},$$
(1)
where $`𝐦`$ is the magnetization direction, $`\gamma `$ is the gyromagnetic ratio, and $`𝐇_{\text{eff}}`$ is the effective magnetic field including the external, demagnetization, and crystal anisotropy fields. The second term on the right-hand side of Eq. (1) was first introduced by Gilbert Gilbert (1955) and the dimensionless coefficient $`\alpha `$ is called the Gilbert damping constant. For a constant $`𝐇_{\text{eff}}`$ and $`\alpha =0`$, $`𝐦`$ precesses around the field vector with frequency $`\omega =\gamma H_{\text{eff}}`$. When damping is switched on $`\alpha >0`$, the precession spirals down to a time independent magnetization along the field direction on a time scale of $`1/\alpha \omega `$. Study of $`\alpha `$ in bulk metallic ferromagnets has drawn a significant interest over several decades. Notwithstanding the large body of both experimental Bhagat and Lubitz (1974) and theoretical Korenman and Prange (1972) work, the damping mechanism in bulk ferromagnets is not yet fully understood.

The magnetization dynamics in thin magnetic films and microstructures is technologically relevant for, e.g., magnetic recording applications at high bit densities. Recent interest of the basic physics community in this topic is motivated by the spin-current induced magnetization switching in layered structures Slonczewski (1996); Myers et al. (1999); Katine et al. (2000). The Gilbert damping constant was found to be $`0.04<\alpha <0.22`$ for Cu-Co and Pt-Co Back et al. (1999); Myers et al. (1999), which is considerably larger than the bulk value $`\alpha _00.005`$ in Co Schreiber et al. (1995); Katine et al. (2000). Previous attempts to explain the additional damping in magnetic multilayer systems involved an enhanced electron-magnon scattering near the interface Berger (1996) and other mechanisms Wegrowe (2000), both in equilibrium and in the presence of a spin-polarized current.

In this Letter we propose a novel mechanism for the Gilbert damping in normal-metal–ferromagnet (N-F) hybrids. According to Eq. (1), the precession of the magnetization direction $`𝐦`$ is caused by the torque $`𝐦\times 𝐇_{\text{eff}}`$. This is physically equivalent to a volume injection of what we call a “spin current”. The damping occurs when the spin current is allowed to leak into a normal metal in contact with the ferromagnet. Our mechanism is thus the inverse of the spin-current induced magnetization switching: A spin current can exert a finite torque on the ferromagnetic order parameter, and, vice versa, a moving magnetization vector loses torque by emitting a spin current. In other words, the magnetization precession acts as a spin pump which transfers angular momentum from the ferromagnet into the normal metal. This effect can be mathematically formulated in terms of the dependence of the scattering matrix of a ferromagnetic layer attached to normal metal leads on the precession of $`𝐦`$, analogous to the parametric charge pumping in nonmagnetic systems Brouwer (1998). The damping contribution is found to obey the LLG phenomenology. Enhancement of the damping constant $`\alpha ^{}=\alpha \alpha _0`$ can be expressed in terms of the scattering matrix at the Fermi energy of a ferromagnetic film in contact with normal metal reservoirs, which can be readily obtained by model or first-principles calculations. Our numerical estimates of $`\alpha ^{}`$ compare well with recent experimental results Mizukami et al. (2001a). Earlier experiments reported in Ref. Heinrich et al. (1987) can also be understood by our model Tserkovnyak et al. (2002).

We consider a ferromagnetic film sandwiched between two paramagnetic layers as shown in Fig. 1. Spin pumping is governed by the ferromagnetic film and the vicinity of the N-F interfaces. The normal metal layers are, therefore, interpreted as reservoirs attached to non-magnetic leads. The quantity of interest is the $`2\times 2`$ current matrix in spin space $`\widehat{I}=\widehat{1}I_c/2\widehat{𝝈}𝐈_se/\mathrm{}`$ for the charge ($`I_c`$) and spin flow ($`𝐈_𝐬`$) from the magnetic film into adjacent normal metal leads, where $`\widehat{1}`$ is the unit matrix and $`\widehat{𝝈}`$ the vector of Pauli spin matrices.

When no voltages are applied and the external field is constant, the charge current vanishes. Two contributions to the spin current $`𝐈_s`$ on either side of the ferromagnet may be distinguished, viz. $`𝐈_s^{\text{pump}}`$ and $`𝐈_s^{(0)}`$. $`𝐈_s^{\text{pump}}`$ is the spin current pumped into the normal metal to be discussed below, whereas $`𝐈_s^{(0)}`$ is the current which flows back into the ferromagnet. The latter is driven by the accumulated spins in the normal metal and gives, e.g., rise to the spin-current induced magnetization switching Slonczewski (1996); Myers et al. (1999); Katine et al. (2000). Here, we model the normal metal as an ideal sink for the spin current, such that a spin accumulation does not build up. This approximation is valid when the spins injected by $`𝐈_s^{\text{pump}}`$ decay and/or leave the interface sufficiently fast, e.g., when the dimensionless conductance of the N-F interface is smaller than $`h/\tau _{\text{sf}}\delta `$ Tserkovnyak et al. (2002). Here, $`\tau _{\text{sf}}`$ is the spin-flip relaxation time and $`\delta `$ is the energy level spacing at the Fermi surface of a normal metal film with a thickness which is the smaller one of the geometrical film thickness and the spin-flip diffusion length.

The current $`\widehat{I}(t)`$ pumped by the precession of the magnetization into the right and left paramagnetic reservoirs, connected to the ferromagnet by normal metal leads ($`R`$) and ($`L`$), may be calculated in an adiabatic approximation since the period of precession $`2\pi /\omega `$ is typically much larger than the relaxation times of the electronic degrees of freedom of the system. The adiabatic charge-current response in nonmagnetic systems by a scattering matrix which evolves under a time-dependent system parameter $`X(t)`$ has been derived in Büttiker et al. (1994); Brouwer (1998). Adopting Brouwer’s notation Brouwer (1998), the generalization to the $`2\times 2`$ matrix current (directed into the normal metal lead $`l=R`$ or $`L`$) reads

$$\widehat{I}(t)^{\text{pump}}=e\frac{\widehat{n}(l)}{X}\frac{dX(t)}{dt},$$
(2)
where the matrix emissivity into the lead $`l`$ is

$$\frac{\widehat{n}(l)}{X}=\frac{1}{4\pi i}\underset{mnl^{}}{}\frac{\widehat{s}_{mn,ll^{}}}{X}\widehat{s}_{mn,ll^{}}^{}+\text{H.c.}$$
(3)
and $`\widehat{s}`$ is the $`2\times 2`$ scattering matrix of the ferromagnetic insertion. $`m`$ and $`n`$ label the transverse modes at the Fermi energy in the normal metal leads and $`l^{}=R,L`$. Spin-flip scattering in the contact is disregarded. $`\widehat{s}`$ depends on the magnetization direction $`𝐦`$ of the ferromagnet through the projection matrices $`\widehat{u}^{}=\left(\widehat{1}+\widehat{𝝈}𝐦\right)/2`$ and $`\widehat{u}^{}=\left(\widehat{1}\widehat{𝝈}𝐦\right)/2`$ (Ref. Brataas et al. (2000)): $`\widehat{s}_{mn,ll^{}}=s_{mn,ll^{}}^{}\widehat{u}^{}+s_{mn,ll^{}}^{}\widehat{u}^{}`$. The spin current pumped by the magnetization precession is obtained by identifying $`X(t)=\phi (t)`$, where $`\phi `$ is the azimuthal angle of the magnetization direction in the plane perpendicular to the precession axis. The resulting current is traceless, $`\widehat{I}^{\text{ pump}}=(e/\mathrm{})\widehat{𝝈}𝐈_s^{\text{pump}}`$, i.e., charge current indeed vanishes, and

$$𝐈_s^{\text{pump}}=\frac{\mathrm{}}{4\pi }\left(A_r𝐦\times \frac{d𝐦}{dt}A_i\frac{d𝐦}{dt}\right),$$
(4)
where the interface parameters are

$`A_r`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{mn}{}}\{\left|r_{mn}^{}r_{mn}^{}\right|^2+\left|t_{mn}^{}t_{mn}^{}\right|^2\},`$ (5)
$`A_i`$ $`=`$ $`\text{Im}{\displaystyle \underset{mn}{}}\{r_{mn}^{}(r_{mn}^{})^{}+t_{mn}^{}(t_{mn}^{})^{}\}.`$ (6)
Here, $`r_{mn}^{}`$ \[$`r_{mn}^{}`$\] is the reflection coefficient for spin-up \[spin-down\] electrons in the $`l`$th lead and $`t_{mn}^{}`$ \[$`t_{mn}^{}`$\] is the transmission coefficient for spin-up \[spin-down\] electrons into the $`l`$th lead. (See Fig. 1 for $`l=R`$.) Using unitarity of the scattering matrix for each spin direction, we can summarize Eqs. (5) and (6) by $`A_r+iA_i=g^{}t^{}`$, where $`g^{\sigma \sigma ^{}}=_{mn}\{\delta _{mn}r_{mn}^\sigma (r_{mn}^\sigma ^{})^{}\}`$ is the (DC) conductance matrix Brataas et al. (2000); Waintal et al. (2000) and $`t^{}=_{nm}t_{mn}^{}(t_{mn}^{})^{}`$. The spin current (4) trivially vanishes for the steady state,  i.e., when $`d𝐦/dt=0`$, and for unpolarized contacts $`s_{mn,ll^{}}^{}=s_{mn,ll^{}}^{}`$.

Per revolution, the precession pumps an angular momentum into an adjacent normal metal layer which is proportional to $`A_r`$, in the direction of the (averaged) applied magnetic field, and decaying in time. At first sight, it is astonishing that a pump can be operated by a single parameter varying in time, whereas the “peristaltic” pumping of a charge current requires at least two periodic parameters Brouwer (1998). However, there are actually two periodic parameters (out of phase by $`\pi /2`$) hidden behind $`\phi (t)`$, viz. the projections of the unit vector defined by $`\phi `$ in the plane perpendicular to the axis of precession.

By conservation of angular momentum, the spin torque on the ferromagnet resulting from the spin pumping into the nonmagnetic leads gives an additional term to the LLG equation (1). After including this term, Eq. (1) remains valid, but the gyromagnetic ratio and the damping constant are renormalized:

$`{\displaystyle \frac{1}{\gamma }}`$ $`=`$ $`{\displaystyle \frac{1}{\gamma _0}}\{1+g_L[A_i^{(L)}+A_i^{(R)}]/4\pi M\},`$ (7)
$`\alpha `$ $`=`$ $`{\displaystyle \frac{\gamma }{\gamma _0}}\{\alpha _0+g_L[A_r^{(L)}+A_r^{(R)}]/4\pi M\}.`$ (8)
Here, $`g_L`$ is the Landé factor and $`M`$ is the total magnetic moment (in units of $`\mu _B`$) of the ferromagnetic film; subscript $`0`$ denotes the bulk values of $`\gamma `$ and $`\alpha `$; superscripts $`(L)`$ and $`(R)`$ denote parameters evaluated on the left and right side of the F layer, respectively. Eqs. (7) and (8) are the central result of this paper. $`A_r`$ and $`A_i`$ affect, e.g., ferromagnetic resonance experiments as a shift of the resonance magnetic field via $`A_i^{(L)}+A_i^{(R)}`$, whereas $`A_r^{(L)}+A_r^{(R)}`$ increases the relative resonance linewidth.

From now on we focus on ferromagnetic films which are thicker than the coherence length $`\lambda _{\text{fc}}=\pi /(k_{}k_{})`$, where $`k_{}`$ are the spin-up and spin-down Fermi wavevectors, i.e., thicker than a few monolayers in the case of transition metals. In this regime, spin-up and spin-down electrons transmitted or scattered from one N-F interface interfere incoherently at the other interface, $`t^{}`$ vanishes and the mixing conductance $`g^{}`$ is governed by the reflection coefficients of the isolated N-F interfaces.

$`A_i=\text{Im}g^{}`$ vanishes for ballistic and diffusive contacts as well as nonmagnetic tunnel barriers Brataas et al. (2000). First-principles calculations find very small $`A_i`$ for Cu-Co and Fe-Cr Xia et al. (2001a). It is, therefore, likely that $`A_i`$ may be disregarded in many systems. If $`A_i`$ does vanish on both sides of the ferromagnetic film, it follows from Eqs. (7) and (8) that the resonance frequency is not modified $`\gamma =\gamma _0`$ and the enhancement of the Gilbert damping is given by $`\alpha ^{}=g_L[A_r^{(L)}+A_r^{(R)}]/4\pi M`$.

The coefficient $`A_r`$ can be estimated by simple model calculations Brataas et al. (2000). For ballistic (point) contacts, $`A_r^B=(1+p)g`$ with the polarization $`p=(g^{}g^{})/(g^{}+g^{})`$ and the average conductance $`g=(g^{}+g^{})/2`$. For diffusive N-F hybrids, $`A_r^D=g_N`$, the conductance of the normal metal part. A nonmagnetic tunneling barrier between F and N suppresses the spin current exponentially. The magnetization precession of a magnetic insulator can also emit a spin current into a normal metal, since $`g^{}`$ does not necessarily vanish because the phase shifts of reflected spin-up and spin-down electrons at the interface may differ Xia et al. (2001a).

Let us now estimate the damping coefficient $`\alpha ^{}`$ for thin films of permalloy (Ni80Fe$`_{20},`$ Py), a magnetically very soft material of great technological importance. Mizukami et al. Mizukami et al. (2001a) measured the ferromagnetic resonance linewidth of N-Py-N sandwiches and discovered systematic trends in the damping parameter as a function of Py layer thickness $`d`$ for different normal metals. The spin polarization of electrons emitted by Py has been measured to be $`p0.4`$ in point contacts Soulen et al. (1998), the magnetization per atom is $`f1.2`$, and Landé factor–$`g_L2.1`$ Mizukami et al. (2001a). The interface conductance of metallic interfaces with Fe or Co is of the order of 1015 $`\mathrm{\Omega }^1`$m-2, with significant but not drastic dependences on interface morphology or material combination Xia et al. (2001b). This corresponds to roughly one conducting channel per interface atom. Assuming the Fermi surface of the normal metal is isotropic, we arrive at the estimate $`\alpha ^{}1.1/d(\text{Å})`$. The factor $`1/d`$ does not reflect an intrinsic effect; a reduced total magnetization is simply more sensitive to a given spin-current loss at the interface. Comparing with the intrinsic $`\alpha _00.006`$ of permalloy Patton et al. (1975); Mizukami et al. (2001a) the spin-current induced damping becomes important for ferromagnetic layers with thickness $`d<100`$ Å. We can refine the estimate by including the significant film-thickness dependence of the magnetization measured by the same group Mizukami et al. (2001a). We, therefore, improve our above estimate as

$$\alpha ^{}(d)\kappa \times \frac{1.1}{d(\text{Å})}\times \frac{f_0}{f(d)},$$
(9)
where $`f_0`$ and $`f(d)`$ are the atomic magnetization of the permalloy bulk and films. $`\kappa `$ is an adjustable parameter representing the number of scattering channels in units of one channel per interface atom, which should be of the order of unity.

The experimental results for the damping factor $`\alpha `$ and the relative magnetization $`f/f_0`$ for N-Py-N sandwiches with  N=Pt, Pd, Ta, and Cu are shown in the insets of Fig. 2. Our estimate (9) appears to well explain the dependence of $`\alpha `$ on the permalloy film thickness $`d`$ (see Fig. 2) for reasonable values of $`\kappa `$. First-principles calculations are called for to test these values.

The lack of a significant thickness dependence of damping parameter of the Cu-Py system requires additional attention. An opaque interface might be an explanation, but it appears more likely that due to long spin-flip relaxation times in Cu, the 5 nm thick buffer layers in Mizukami et al. (2001a) do not provide the ideal sink for the injected spins as assumed above. This means that a nonequilibrium spin accumulation on Cu opposes the pumped spin current and nullifies the additional damping when $`h/\tau _{\text{sf}}\delta `$ is comparable or smaller than the conductance $`g`$. For 5 nm Cu buffers, $`g\delta /h10^{13}`$ s-1, whereas $`1/\tau _{\text{sf}}10^{12}`$ s-1 Meservey and Tedrow (1978). It follows that Cu is indeed a poor sink for the injected spins and the Gilbert damping constant is not enhanced. On the other hand, Pt, Ta, and Pd are considerably heavier than Cu and, since $`1/\tau _{\text{sf}}`$ scales as $`Z^4`$ Abrikosov and Gor’kov (1962), where $`Z`$ is the atomic number, have much larger spin-relaxation rates and our arguments hold.

A physical picture of the effect of magnetization precession in layered systems has been proposed earlier by Hurdequint et al. Monod et al. (1972) in order to explain ferromagnetic and conduction electron spin resonance experiments. These authors realized that the precessing magnetization is a source of a nonequilibrium spin accumulation which diffuses out of the N-F interfaces into the adjacent normal metal layers where it can dissipate by spin-flip processes. Enhanced Gilbert damping in thin ferromagnetic films in contact with normal metal has also been discussed by Berger Berger (1996) for a ballistic N-F interface in a spin-valve configuration. His expression for the damping coefficient \[Eq. (20) in Ref. Berger (1996)\] scales like ours as a function of layer thickness, but differs as a function of material parameters. E.g., in contrast to our result, Berger’s expression does not vanish with vanishing exchange splitting.

In conclusion, we demonstrated that the Gilbert damping constant is enhanced in thin magnetic films with normal metal buffer layers by a spin-pump effect through the N-F contact. The damping is significant for transition metal films thinner than about 10 nm. Recent experiments on permalloy films Mizukami et al. (2001a) are well explained.

We are grateful to L. Berger, B. I. Halperin, B. Heinrich, D. Huertas-Hernando, and Yu. V. Nazarov for stimulating discussions. This work was supported in part by the DARPA Award No. MDA 972-01-1-0024, the Norwegian Research Council, the NEDO International Joint Research Grant Program “Nano-magnetoelectronics”, NSF Grant No. DMR 99-81283, and the Schlumberger Foundation.

Note added.—After submission of this paper, we became aware of new exciting experimental results Urban et al. (2001) on enhanced Gilbert damping in ultrathin iron films which support the spin-pumping mechanism proposed in this Letter Tserkovnyak et al. (2002).


text:
# Tracing Cosmic Evolution with the Las Campanas Distant Cluster Survey

## 1. Summary

The past few years have been a period of tremendous progress in detection of distant clusters. As recently as five years ago no single survey contained more than ten cluster candidates at $`z`$$`>`$0.5, but a series of efforts have since led to a several order of magnitude increase in the number of known high-redshift systems (Postman et al. 1996, Scodeggio et al. 1999, Gladders & Yee 2001, Gonzalez et al. 2001). The recently completed Las Campanas Distant Cluster Survey (LCDCS) is our contribution to this field. Clusters in the LCDCS are detected as regions of excess optical surface brightness relative to the background sky, a technique that permits detection of clusters to $`z`$$``$0.9 with shallow, drift-scan imaging. The resulting statistical catalog consists of 1073 cluster candidates at $`z`$$`>`$0.3. To enable efficient use of these catalogs we use follow-up imaging and spectroscopy to characterize the properties of the sample, including the contamination rate as well as estimates of the redshifts and relative velocity dispersions of the candidates (see Gonzalez et al. 2001$`a`$ for details).

One application of this sample is determination of the cluster correlation function (Gonzalez et al. 2001$`b`$). For roughly dispersion-limited subsamples at $`z`$$``$0.5 the observed angular correlation functions are well fit by the power law $`\omega (\theta )=A_\omega \theta ^{1\gamma }`$ with $`\gamma =2.1`$. Using the angular correlation data, we derive the spatial correlation length ($`r_0`$) as a function of mean intracluster separation ($`d_c`$) via the cosmological Limber inversion. We find correlation lengths for the LCDCS that are similar to those observed for local catalogs, and observe a dependence of $`r_0`$ upon $`d_c`$ that is comparable to the results of Croft et al. (1997) for the APM catalog. These results are consistent with theoretical expectations, which predict only modest evolution in the clustering amplitude since $`z`$=0.5.

Several programs are also being conducted with LCDCS clusters to study the evolution of cluster galaxies. Nelson et al. (2001$`a,b`$) are studying the structural and luminosity evolution of BCG’s; Nelson et al. (2001$`c`$) are studying evolution in the luminosity function and color of the red envelope. Meanwhile, the EDisCS collaboration is conducting a deep multicolor imaging and spectroscopic survey of 20 of the most massive LCDCS candidates at $`z`$=0.5 and $`z`$=0.8, with the aim of constraining cluster galaxy morphologies, scaling relations, cluster dynamics, and the evolution of cluster galaxies at these epochs.

## References

Croft, R. A. C., Dalton, G. B., Efstathiou, G., Sutherland, W. J., & Maddox, S. J. 1997, 291, 305

Gladders. M., & Yee H. K. C. 2001 (this volume)

Gonzalez, A. H., Zaritsky, D., Dalcanton, J. J., & Nelson, A. E. 2001, astro-ph/0106055

Gonzalez, A. H., Zaritsky, D., & Wechsler, R. H. 2001, submitted to ApJ

Nelson, A. E., Dalcanton, J. J., Simard, L., Zaritsky, D., & Gonzalez, A. H. 2001$`a`$, submitted to ApJ

Nelson, A. E., Gonzalez, A. H., Zaritsky, D., & Dalcanton, J. J. 2001$`b`$, submitted to ApJ

Nelson, A. E., Gonzalez, A. H., Zaritsky, D., & Dalcanton, J. J. 2001$`c`$, astro-ph/0108279

Postman, M., Lubin, L. M., Gunn, J. E., Oke, J. B., Hoessel, J. G., Schneider, D. P., & Christensen, J. A. 1996, AJ, 111, 615

Scodeggio, M. et al. 1999 A&AS, 137, 83