text:
# “TOY MODELS” OF TURBULENT CONVECTION AND THE HYPOTHESIS OF THE LOCAL ISOTROPY RESTORATION

## Abstract

A brief review is given of recent results devoted to the effects of large-scale anisotropy on the inertial-range statistics of the passive scalar quantity $`\theta (t,𝐱)`$, advected by the synthetic turbulent velocity field with the covariance $`\delta (tt^{})|𝐱𝐱^{}|^\epsilon `$. Inertial-range anomalous scaling behavior is established, and explicit asymptotic expressions for the structure functions $`S_n(𝐫)[\theta (t,𝐱+𝐫)\theta (t,𝐱)]^n`$ are obtained; they are represented by superpositions of power laws with universal (independent of the anisotropy parameters) anomalous exponents, calculated to the first order in $`\epsilon `$ in any space dimension. The exponents are associated with tensor composite operators built of the scalar gradients, and exhibit a kind of hierarchy related to the degree of anisotropy: the less is the rank, the less is the dimension and, consequently, the more important is the contribution to the inertial-range behavior. The leading terms of the even (odd) structure functions are given by the scalar (vector) operators. The small-scale anisotropy reveals itself in odd correlation functions: for the incompressible velocity field, $`S_3/S_2^{3/2}`$ decreases going down towards to the depth of the inertial range, while the higher-order odd ratios increase; if the compressibility is strong enough, the skewness factor also becomes increasing.

This paper is intended for the memorial issue of “Zapiski Nauchnykh Seminarov POMI” dedicated to the memory of Anatolij Georgievich Izergin.

Dedicated to the memory of Anatolij Georgievich Izergin

The investigation of intermittency and anomalous scaling in fully developed turbulence remains essentially an open theoretical problem. Both the natural and numerical experiments suggest that the deviation from the predictions of the classical Kolmogorov–Obukhov theory  is even more strongly pronounced for a passively advected scalar field than for the velocity field itself; see, e.g.,  and literature cited therein. At the same time, the problem of passive advection appears to be easier tractable theoretically: even simplified models describing the advection by “synthetic” velocity fields with given Gaussian statistics reproduce many of the anomalous features of genuine turbulent heat or mass transport observed in experiments. Hence the term “toy models.” Therefore, the problem of passive scalar advection, being of practical importance in itself, may also be viewed as a starting point in studying anomalous scaling in the turbulence on the whole.

The most remarkable progress has been achieved for the so-called rapid-change model, which dates back to Obukhov and Batchelor and was explicitly introduced by Kraichnan : for the first time, the anomalous exponents have been calculated on the basis of a microscopic model and within regular perturbation expansions; see, e.g., \[6–10\] and references therein.

Within the “zero-mode approach,” developed in \[6–9\], nontrivial anomalous exponents are related to the zero modes (unforced solutions) of the closed exact equations satisfied by the equal-time correlations. Within the approach based on the field theoretic renormalization group (RG) and operator product expansion (OPE), the anomalous scaling emerges as a consequence of the existence in the model of composite operators with negative critical dimensions, which determine the anomalous exponents \[10–13\].

An important question recently addressed is the effects of large-scale anisotropy on the inertial-range statistics of passively advected fields \[8,9,12–18\] and the velocity field itself .

According to the classical Kolmogorov–Obukhov theory , the anisotropy introduced at large scales by the forcing (boundary conditions, geometry of an obstacle etc) dies out when the energy is transferred down to smaller scales owing to the cascade mechanism. A number of recent works confirms this picture for the even correlation functions, thus giving some quantitative support to the aforementioned hypothesis on the restored local isotropy restoration of the inertial-range turbulence for the velocity and passive fields \[12–19\]. More precisely, the exponents describing the inertial-range scaling exhibit universality and hierarchy related to the degree of anisotropy, and the leading contribution to an even function is given by the exponent from the isotropic shell \[12,15–19\]. Nevertheless, the anisotropy survives in the inertial range and reveals itself in odd correlation functions, in disagreement with what was expected on the basis of the cascade ideas. The so-called skewness factor decreases down the scales much slower than expected , while the higher-order odd dimensionless ratios (hyperskewness etc) increase, thus signalling of persistent small-scale anisotropy . The effect seems rather universal, being observed for the scalar  and vector  fields, advected by the Gaussian rapid-change velocity, and for the scalar advected by the two-dimensional Navier-Stokes (NS) velocity field .

Below we briefly review the most recent results obtained on the example of the rapid-change model using the RG and OPE techniques; more detailed discussion can be found in the original papers \[10–13\]. The detailed description of the RG approach to the models of fully developed turbulence is given in the book  and review paper . For simplicity, we restrict ourselves with the rapid-change models; the case of a finite-correlated velocity field is considered .

The advection of a passive scalar field in the presence of an imposed linear gradient is described by the equation

$$_t\theta +(𝐯\mathbf{})\theta =\nu _0^2\theta (\mathrm{𝐡𝐯}).$$
(1)
Here $`\theta (x)\theta (t,𝐱)`$ is the random (fluctuation) part of the total scalar field $`\mathrm{\Theta }(x)=\theta (x)+(\mathrm{𝐡𝐱})`$, $`𝐡`$ is a constant vector that determines distinguished direction, $`\nu _0`$ is the molecular diffusivity coefficient, $`_t/t`$, $`_i/x_i`$, and $`^2_i_i`$ is the Laplace operator.

The velocity field $`𝐯(x)=\{v_i(x)\}`$ obeys a Gaussian distribution with zero mean and correlator

$$v_i(x)v_j(x^{})=\delta (tt^{})K_{ij}(𝐫)$$
(2)
with

$$K_{ij}(𝐫)=D_0\frac{d𝐤}{(2\pi )^d}\frac{P_{ij}(𝐤)}{k^{d+\epsilon }}\mathrm{exp}[\mathrm{i}𝐤𝐫],𝐫𝐱𝐱^{},$$
(3)
where $`P_{ij}(𝐤)=\delta _{ij}k_ik_ȷ/k^2`$ is the transverse projector, $`𝐤`$ is the momentum, $`k|𝐤|`$, $`D_0>0`$ is an amplitude factor, and $`0<\epsilon <2`$ is a free parameter. The IR regularization is provided by the cut-off in the integral (3) from below at $`km`$, where $`m1/L`$ is the reciprocal of the integral turbulence scale; the precise form of the cut-off is not essential. For $`0<\epsilon <2`$, the difference

$$S_{ij}(𝐫)K_{ij}(\mathrm{𝟎})K_{ij}(𝐫)$$
(4)
has a finite limit for $`m0`$:

$$S_{ij}(𝐫)=Dr^\epsilon \left[\left(d+\epsilon 1\right)\delta _{ij}\epsilon \frac{r_ir_j}{r^2}\right],$$
(5)
with

$`D={\displaystyle \frac{D_0\mathrm{\Gamma }(\epsilon /2)}{(4\pi )^{d/2}2^\epsilon (d+\epsilon )\mathrm{\Gamma }(d/2+\epsilon /2)}},`$

where $`\mathrm{\Gamma }(\mathrm{})`$ is the Euler gamma function (note that $`D>0`$). It follows from Eq. (5) that $`\epsilon `$ can be viewed as a kind of Hölder exponent, which measures the roughness of the velocity field. In the RG approach, the exponent $`\epsilon `$ plays the same role as the parameter $`\epsilon =4d`$ does in the RG theory of critical phenomena . The relations

$$g_0D_0/\nu _0\mathrm{\Lambda }^\epsilon $$
(6)
define the coupling constant $`g_0`$ (i.e., the expansion parameter in the ordinary perturbation theory) and the characteristic ultraviolet (UV) momentum scale $`\mathrm{\Lambda }`$.

The term $`(\mathrm{𝐡𝐯})`$ in Eq. (1) acts as a driving force; it maintains the steady state and, simultaneously, introduces the large-scale anisotropy. In the original (isotropic) formulation of the model this term is replaced with an artificial Gaussial random force $`f(x)`$ with the covariance

$$f(x)f(x^{})=\delta (tt^{})C(mr),r|𝐱𝐱^{}|,$$
(7)
where $`C(mr)`$ is a function finite as $`mr=0`$.

The issue of interest is, in particular, the behavior of the equal-time structure functions

$$S_n(𝐫)[\theta (t,𝐱)\theta (t,𝐱^{})]^n$$
(8)
in the inertial range, specified by the inequalities $`\mathrm{\Lambda }>>1/r>>m`$. The dimensionality considerations give

$$S_n(𝐫)=\nu _0^{n/2}r^nR_n(\mathrm{\Lambda }r,mr,z),$$
(9)
where $`z`$ is the cosine between the directions of $`𝐫`$ and $`𝐡`$ and $`R_n`$ are some functions of dimensionless parameters. In principle, they can be calculated within the ordinary perturbation theory (i.e., as series in $`g_0`$), but this is not useful for studying inertial-range behavior: the coefficients are singular in the limits $`\mathrm{\Lambda }r\mathrm{}`$ and/or $`mr0`$, which compensate the smallness of $`g_0`$ (assumed in ordinary perturbation theory), and in order to find correct IR behavior we have to sum the entire series. The desired summation can be accomplished using the RG and OPE; see Refs. .

The RG analysis consists of two main stages. On the first stage, the multiplicative renormalizability of the model is demonstrated and the differential RG equations for its correlation functions are obtained. The asymptotic behavior of the functions like (8) for $`\mathrm{\Lambda }r>>1`$ and any fixed $`mr`$ is given by IR stable fixed points of the RG equations and has the form

$$S_n(r)=D_0^{n/2}r^{n(1\epsilon /2)}R_n(mr,z),\mathrm{\Lambda }r>>1,$$
(10)
with certain, as yet unknown, scaling functions $`R_n(mr,z)`$.

On the second stage, the small $`mr`$ behavior of the functions $`R_n`$ is studied within the general representation (10) using the OPE. It shows that, in the limit $`mr0`$, the functions $`R_n`$ have the asymptotic forms

$$R_n(mr,z)=\underset{F}{}C_F(mr,z)(mr)^{\mathrm{\Delta }_F},$$
(11)
where $`C_F`$ are coefficients regular in $`mr`$. In general, the summation is implied over all possible renormalized composite operators $`F`$ allowed by the symmetry; $`\mathrm{\Delta }_F`$ being their critical dimensions.

The peculiarity of the models describing turbulence is the existence of the so-called “dangerous” composite operators with negative critical dimensions . Their contributions into the OPE give rise to a singular behavior of the scaling functions for $`mr0`$, and the leading term is given by the operator with minimal $`\mathrm{\Delta }_F`$. In the isotropic case only scalar operators have nonvanishing mean values and contribute to the right hand side of Eq. (11). The odd functions $`S_{2k+1}`$ vanish, while the leading contributions to $`S_{2k}`$ are determined by scalar gradients $`F_k=(_i\theta _i\theta )^k`$ and have the form

$$S_{2k}(r)D_0^kr^{k(2\epsilon )}(mr)^{\mathrm{\Delta }_k},$$
(12)
where the critical dimensions $`\mathrm{\Delta }_k`$ of the operators $`F_k`$ are given by

$$\mathrm{\Delta }_k=2k(k1)\epsilon /(d+2)+O(\epsilon ^2)=2k(k1)\epsilon /d+O(1/d^2).$$
(13)
The expression (12) agrees with the results obtained earlier in  using the zero-mode techniques; the $`O(\epsilon ^2)`$ contribution to $`\mathrm{\Delta }_k`$ is reported in . Note that the operators $`F_l`$ with $`lk`$ (whose contributions would be more important) do not appear in the right hand side of the OPE for $`S_{2k}`$; this is a consequence of the linearity of the original equation (1) in $`\theta (x)`$.

In the anisotropic case, an important role is played by the irreducible tensor operators of the form

$$F[n,p]_{i_1}\theta \mathrm{}_{i_p}\theta (_i\theta _i\theta )^l+\mathrm{},$$
(14)
where $`p`$ is the number of the free vector indices and $`n=p+2l`$ is the total number of the fields $`\theta `$ entering into the operator; the vector indices of the symbol $`F[n,p]`$ are omitted. The dots stand for the subtractions (needed for $`p2`$) which ensure that the resulting expression is traceless with respect to any contraction, for example, $`F[2,2]=_i\theta _j\theta \delta _{ij}F[2,0]/d`$.

For $`𝐡0`$, the operators (14) acquire nonvanishing mean values, built of the vectors $`𝐡`$ and delta symbols, and give nontrivial contributions to the OPE for $`S_n`$; in particular, this gives rise to nonvanishing odd functions.

The critical dimensions associated with the operators (14) have the form  (see also Refs. ):

$$\mathrm{\Delta }[n,p]=\frac{2n(n1)(d+1)(np)(d+n+p2)}{2(d1)(d+2)}\epsilon +O(\epsilon ^2).$$
(15)
The straightforward analysis of the explicit one-loop expression (15) shows that for any fixed $`n`$, the dimension $`\mathrm{\Delta }[n,p]`$ decreases monotonically with $`p`$ and reaches its minimum for the minimal possible value of $`p=p_n`$, i.e., $`p=0`$ if $`n`$ is even and $`p=1`$ if $`n`$ is odd:

$$\mathrm{\Delta }[n,p]>\mathrm{\Delta }[n,p^{}]\mathrm{if}p>p^{}.$$
(16)
Furthermore, this minimal value $`\mathrm{\Delta }[n,p_n]`$ is negative and decreases monotonically as $`n`$ increases:

$$0>\mathrm{\Delta }[2k,0]>\mathrm{\Delta }[2k+1,1]>\mathrm{\Delta }[2k+2,0].$$
(17)
Finally, we note that for any fixed $`p`$, the dimension (15) decreases monotonically as $`n`$ increases:

$$\mathrm{\Delta }[n,p]>\mathrm{\Delta }[n^{},p]\mathrm{if}n<n^{}.$$
(18)
The mean value of the operator with the dimension $`\mathrm{\Delta }[k,p]`$ is a traceless (for $`p2`$) tensor built of the vector $`𝐡`$ and Kronecker delta symbols; the contraction with the vector indices of the corresponding coefficient $`C_F`$ in Eq. (11) gives rise to $`p`$-th order Legendre polynomial $`P_p(z)`$, so that the expansion (11) is analogous to the decomposition in irreducible representations of the rotation group, used in Refs. \[15–19\].

The inequalities (16), (17) show that the contributions of the tensor operators (14) into the asymptotic expression (11) exhibit a kind of hierarchy: the less is the rank, the more important is the contribution.

The leading term of the expression (11) for the even (odd) function $`S_n`$ is determined by the scalar (vector) composite operator consisting of $`n`$ factors $`\theta `$ and has the form

$$S_n(r)(hr)^n(mr)^{\mathrm{\Delta }[n,p_n]}.$$
(19)
Let us now explore the consequences of our results for the issue of the local isotropy restoration in this problem.

As a measure of small-scale anisotropy one usually uses the dimensionless ratios of the structure functions:

$$R_n(𝐫)S_n(𝐫)/S_2^{n/2}(𝐫).$$
(20)
From Eq. (19) it then follows that in inertial range we have:

$$R_{2k+1}(mr)^{\mathrm{\Delta }[2k+1,1](2k+1)\mathrm{\Delta }[2,0]/2},$$
(21)
$$R_{2k}(mr)^{\mathrm{\Delta }[2k,0]k\mathrm{\Delta }[2,0]}.$$
(22)
The dependence on the Péclet number, $`Pe(\mathrm{\Lambda }/m)^\epsilon `$ can be estimated by replacing $`r`$ with $`1/\mathrm{\Lambda }`$; see Ref. . Using explicit $`O(\epsilon )`$ expressions for $`\mathrm{\Delta }[n,p]`$ we then obtain:

$$R_{2k+1}Pe^{(d+24k^2)/[2(d+2)]},$$
(23)
$$R_{2k}Pe^{2k(k1)/(d+2)}.$$
(24)
Since the leading terms of the even functions (19) are determined by the exponents of the isotropic shell (i.e., those related to the scalar composite operators), the inertial-range behavior of the even ratios (22), (24) is the same as in the isotropic model. This gives a quantitative support to the universality of leading anomalous exponents for the even functions with respect to different classes of forcing. On the other hand, the odd quantities (21), (23) appear sensitive to the anisotropy: $`R_3`$ in (23) slowly decreases for $`Pe\mathrm{}`$, while ratios $`R_{2k+1}`$ with $`k2`$ increase with $`Pe`$, thus signalling of the persistent small-scale anisotropy.

The picture outlined above seems rather general. Indeed, it is compatible with that recently established for the NS turbulence , the scalar field passively advected by a NS velocity in the two-dimensional inverse energy cascade , and for the vector (magnetic) field advected by the velocity of the type (3), see Refs. .

To conclude with, let us briefly discuss the effects of compressibility onto the small-scale anisotropy persistence. The transversal projector in the correlator (3) is then replaced by the combination

$$P_{ij}(𝐤)+\alpha Q_{ij}(𝐤),$$

where $`Q_{ij}(𝐤)=k_ik_j/k^2`$ is the longitudinal projector and $`\alpha >0`$ is a free parameter. The model (1) then corresponds to the advection of a tracer (entropy, temperature, concentration of a pollutant), while the nonlinearity $`\theta ^{}\mathbf{}(𝐯\theta )`$ describes advection of a density.

Equation (15) for a tracer becomes

$$\mathrm{\Delta }[n,p]=\frac{2n(n1)(1\alpha )(np)(n+p+d2)(d+1+\alpha )}{2(d+2)(d1+\alpha )}\epsilon ,$$
(25)
the $`O(\epsilon ^2)`$ correction is calculated in  for all $`n`$, $`p`$, $`d`$ and $`\alpha `$.

Although the hierarchy relations (16)–(18) remain valid for all values of $`\alpha >0`$ (in particular, $`\mathrm{\Delta }[n,p]/p>0`$), the corrections become closer to leading terms as $`\alpha `$ increases:

$$^2\mathrm{\Delta }[n,p]/p\alpha <0.$$

Furthermore, the expressions (23) take on the form

$$R_{2k+1}(mr)^{\epsilon [(d1+\alpha )(d+24k^2)8\alpha k^2]/2(d+2)(d1+\alpha )}.$$
(26)
It follows from (26) that $`R_3`$ diverges for $`mr0`$ provided $`\alpha `$ is large enough \[namely, $`\alpha >(d1)(d+2)/(10d)+O(\epsilon )`$\], while the divergence of the higher ratios $`R_{2k+1}`$ enhances as $`\alpha `$ increases.

This means that compressibility enhances the penetration of the large-scale anisotropy towards the depth of the inertial range. This fact also seems universal, being observed in the model of the passively advected magnetic field .

Some remarks are now in order.

The critical dimensions of all composite operators (15) are independent of the forcing: they remain unchanged, when the isotropic stirring force (7) is replaced in Eq. (1) with the anisotropic term $`(\mathrm{𝐡𝐯})`$. The difference is that for the former case, only scalar operators contribute to the representations like (11), while for the latter case the irreducible tensor operators acquire nonzero mean values and also contribute to Eq. (11). This is easily understood in the language of the zero-mode approach: the noise $`f`$ and the term $`(\mathrm{𝐡𝐯})`$ do not affect the differential operators in the equations satisfied by the equal-time correlations functions; the zero modes (homogeneous solutions) coincide in the two cases, but in the latter case the modes with nontrivial angular dependence should be taken into account. In the language of the RG (which is also applicable to the case of a finite-correlated velocity, see Refs. ) this is explained as follows: the stirring force or field $`𝐡`$ do not enter into the diagrams that determine the renormalization of the operators (15). However, the anomalous exponents of the passive scalar become nonuniversal and acquire the dependence on the anisotropy parameters if the velocity field is taken to be strongly anisotropic also at small scales .

The picture outlined above and in Refs.  for passively advected fields (a superposition of power laws with universal exponents and nonuniversal amplitudes) seems rather general, being compatible with that established recently in the field of NS turbulence, on the basis of numerical simulations of channel flows and experiments in the atmospheric surface layer, see Refs.  and references therein. In those papers, the velocity structure functions were decomposed into the irreducible representations of the rotation group. It was shown that in each sector of the decomposition, scaling behaviour can be found with apparently universal exponents. The amplitudes of the various contributions are nonuniversal, through the dependence on the position in the flow, the local degree of anisotropy and inhomogeneity, and so on .

This is rather surprising because the equations for the correlation functions in such cases are neither closed nor isotropic and homogeneous. Although the hierarchy similar to Eq. (16) is demonstrated by the critical dimensions of certain tensor operators in the stirred NS turbulence, see Ref.  and Sec. 2.3 of , the relationship between them and the anomalous exponents is not obvious there. One can thus speculate that the anomalous scaling for the genuine turbulence can also appear a linear phenomenon. It is worth recalling here that the so-called “additive fusion rules,” hypothesized for the NS turbulence in a number of papers, Refs. , and characteristic of the models with multifractal behavior (see Ref. ), arise naturally in the context of the rapid-change models owing to their linearity. The existing results for the Burgers turbulence can also be interpreted naturally as a consequence of similar fusion rules, where only finite number of dangerous operators contributes to each structure function, see Ref. .

Of course, one should not insist too much on this bold assumption.

To conclude, let us compare briefly the situation for the passively advected fields with the case of weak acoustic turbulence, where the spectra can be obtained as solutions of the linear kinetic equations (see Refs. ). For weakly dispersive waves (e.g., with the dispersion law $`\omega (k)k^{1+\delta }`$ with $`\delta <<1`$), the anisotropy introduced by the large-scale forcing enhances going down towards to the depth of the inertial range . The hierarchy of the exponents related to the Legendre decomposition is opposite to that established below and in Refs. : anisotropic corrections decrease slower for larger $`p`$’s . On the contrary, for the nondispersive waves ($`\delta =0`$) the hierarchy of the exponents is similar to that in our case, the anisotropic corrections decay faster and faster with $`p`$ and the spectrum tends to become isotropic at small scales . To the best of our knowledge, no information is available for the higher-order correlation functions for such models. One can thus conclude that turbulent systems can exhibit essentially different types of behavior with respect to the small-scale isotropy restoration.

Acknowledgments. I am thankful to L. Ts. Adzhemyan, M. Hnatič, J. Honkonen, M. Lässig, A. Mazzino, P. Muratore-Ginanneschi, M. Yu. Nalimov and A. N. Vasil’ev for discussions. The work was supported in part by the Russian Foundation for Fundamental Research (Grant No. 99-02-16783) and the Grant Center for Natural Sciences of the Russian State Committee for Higher Education (Grant No. 97-0-14.1-30). A part of this work was performed during my stay at the Department of Physics of the University of Helsinki, financed by the Academy of Finland.


format: text

text:
# Pressure–induced residual resistivity anomaly in CeCu5Au

## Abstract

The electrical resistivity of the magnetically ordered CeCu5Au has been investigated under pressure up to 8.5 GPa. In the magnetically ordered region ($`p<3.4`$ GPa) the residual resistivity $`\rho _0`$ shows a pronounced maximum as a function of pressure. Even in the nonmagnetic region $`\rho _0`$ decreases monotonically by more than a factor of three. These two effects can be qualitatively explained in terms of the interplay of pressure, magnetism and disorder in a strongly correlated electron system with weak disorder.

The behavior of the electrical resistivity $`\rho `$ of metals when approaching zero temperature $`T`$ has been of interest ever since the beginning of solid state physics. While it is clear that the residual resistivity $`\rho _0`$ of conventional metals is largely governed by lattice defects and impurities, the situation is much less clear for metals with strong electronic correlations and/or magnetic order at low $`T`$. In the independent-electron approximation, $`\rho _0=m^{}/e^2n\tau =\frac{3}{2}\pi (h/e^2)/(kF^2l)`$, where $`m^{}`$ is the effective mass, $`\tau `$ the scattering time arising from electron scattering by defects and impurities, $`l`$ the corresponding mean free path, $`n`$ the electron density and $`k_\mathrm{F}`$ the Fermi wave-number. In heavy-fermion (HF) systems, $`m^{}`$ is enhanced by up to a factor of 1000 due to the Kondo effect arising from the exchange coupling of conduction electrons and nearly localized $`f`$-electrons. Nevertheless, $`\rho _0`$ is rather low in stoichiometric compounds: the quasi particles obey Bloch’s theorem and scatter at $`T=0`$ only from rare defects. The simple independent-electron expression above suggests that hydrostatic pressure $`p`$ should affect $`\rho _0`$ only to a minor extent since both $`k_\mathrm{F}`$ and $`l`$ should depend on $`p`$ only through the (small) change of electron density and interatomic distances.

In this letter, we report $`\rho _0`$ measurements under hydrostatic pressure on the stoichiometric HF compound CeCu5Au where heavy quasiparticles, leading to a Sommerfeld coefficient $`\gamma =0.64`$ J/molK2 for $`T0`$, coexist with incommensurate antiferromagnetic (AF) order . The magnetic ordering temperature $`T_\mathrm{N}`$ in this system can be suppressed by pressure . Surprisingly, a strong dependence of $`\rho _0`$ on $`p`$ is found: starting from $`\rho _0=28\mu \mathrm{\Omega }\mathrm{cm}`$ for $`p`$ = 0, $`\rho _0`$ passes over a pronounced maximum of 58 $`\mu \mathrm{\Omega }\mathrm{cm}`$ at 1.8 GPa and decreases to 20 $`\mu \mathrm{\Omega }\mathrm{cm}`$ at $``$ 3 GPa. Even in the nonmagnetic Fermi-liquid state ($`p3.4`$ GPa), $`\rho _0`$ decreases substantially upon further pressure increase, reaching 6 $`\mu \mathrm{\Omega }\mathrm{cm}`$ at 8.5 GPa.

The single crystal of CeCu5Au (space group Pnma) used in this study was grown with the Czochralski method from the starting constituents Ce (5N), Cu (5N), and Au (4N). A part of the crystal was used for neutron experiments which revealed a very good quality of the single-phase sample . Two small rectangular pieces ($`35\times 73\times 861\mu `$m3 and $`28\times 85\times 561\mu `$m3) cut from the same crystal were used in two pressure experiments. A clamped-anvil high-pressure device , capable of reaching 10 GPa, was used to measure $`\rho (T)`$ by the four-point method for a current $`I`$ along the crystallographic $`b`$ direction, i.e. perpendicular to the magnetic ordering vector $`𝐐`$ . The sample and the pressure manometer (a thin Pb foil) were embedded in a soft pressure-transmitting medium (steatite) to ensure homogeneous pressure conditions. Both samples showed the same $`\rho _0(p)`$ behavior, ruling out an error due to a possible rearrangement of the contact leads under pressure.

The $`\rho (T)`$ data of the larger sample are shown in Fig. 1. The AF phase transition at ambient $`p`$ occurs at $`T_\mathrm{N}=2.35`$ K, in perfect agreement with specific-heat and magnetization measurements . Reducing the unit–cell volume by external pressure tunes $`T_\mathrm{N}`$ to zero . This behavior is well known for a number of HF systems and is attributed to an increase of hybridization between conduction electrons and local moments, thus favoring the Kondo effect over the RKKY interaction. The detailed analysis of these data  reveals a critical pressure $`p_c`$ 3.4 GPa for the quantum-critical point (QCP) where $`T_\mathrm{N}0`$.

Figure 2 shows the pressure dependence of $`\rho _0`$, determined in the $`T`$ range 30 mK $`<T<50`$ mK. The particularly striking observations are (i) the strong nonmonotonic variation of $`\rho _0`$ within a few GPa where $`T_\mathrm{N}`$ is monotonically suppressed to zero, and (ii) the strong decrease of $`\rho _0`$ even in the paramagnetic state (for $`p>3.4`$ GPa).

The maximum in $`\rho _0(p)`$ is centered at $`p1.8`$ GPa, i. e., well inside the magnetic phase. The $`\rho _0(p)`$ in the magnetic region is in a superficial regard qualitatively similar to the $`\rho _0(x)`$ variation in CeCu6-xAux . However, $`\rho _0(x)`$ with a maximum for $`x=0.5`$ arises from $`x`$-dependent disorder and negative lattice pressure while external pressure is not expected to change the number of impurities in CeCu5Au.

Pressure is not the only parameter to change $`\rho _0`$ drastically. Likewise, an external magnetic field $`B`$, applied along the (easy) crystallographic $`c`$ direction, yields a reduction of $`\rho `$ at low $`T`$. Figure 3 shows $`\rho (T)`$ curves at $`p`$ = 2.2 and 3.5 GPa for $`B`$ = 0 and 7.5 T. Well inside the magnetic region (for $`p`$ = 2.2 GPa, $`T_\mathrm{N}=1.35`$ K), $`B`$ = 7.5 T destroys the AF order and causes a decrease of $`\rho _0`$ by more than 60%. In the nonmagnetic phase ($`p`$ = 3.5 GPa), $`B`$ affects $`\rho (T)`$ as well, but here the influence becomes weaker as $`T0`$. Hence, the “magnetic” contribution to $`\rho _0`$ decreases as the system is tuned away from the magnetic-nonmagnetic instability.

The impressive $`\rho _0(p)`$ dependence of CeCu5Au should be regarded in the context of the $`p`$-induced $`\rho _0`$ changes in other HF systems. For CeCu2Ge2 a maximum is found well inside the nonmagnetic phase as for CeCu2Si2 . The magnitude of the $`\rho _0`$ peak in superconducting CeCu2Ge2 at a pressure close to the $`T_\mathrm{c}(p)`$ maximum, can be associated with charge fluctuations  or even a valence transition. For YbCu2Si2 only a fraction of the Yb ions orders magnetically at $`p_\mathrm{c}`$  and the broad and shallow $`\rho _0(p)`$ maximum may correspond to full ordering at higher $`p`$. In the case of CeRu2Ge2 the $`\rho _0(p)`$ variation at low $`p`$ might be related to a change of the magnetic ordering vector between low and intermediate $`p`$ , but

the effect is rather small.

The following discussion treats our major observations for CeCu5Au, i.e., the (i) nonmonotonic $`\rho _0(p)`$ dependence in the AF phase and (ii) the strong $`\rho _0`$ decrease for $`p>p_c`$, in terms of models that are generic to strongly correlated systems with some disorder.

We first focus on the AF phase and argue that the strong $`\rho _0(p)`$ dependence is a consequence of the quasi-particle scattering from local variations in the staggered magnetization. The qualitative picture is as follows: impurities change the chemical environment of the Ce atoms and lead to a large spatial variation of the local susceptibilities $`\chi _\mathrm{L}(𝐫)1/T_\mathrm{K}(𝐫)`$ with a “Kondo temperature” $`T_\mathrm{K}(𝐫)`$ which depends exponentially on the local environment. In the AF phase this leads to a variation of the local magnetization $`\mathrm{\Delta }M_\mathrm{L}(𝐫)f[H_{\mathrm{AF}}/T_\mathrm{K}^{\text{imp}}]f[H_{\mathrm{AF}}/T_\mathrm{K}^{\text{bulk}}]`$, where $`H_{\mathrm{AF}}(𝐫)`$ is a staggered (exchange) field proportional to the AF order parameter and $`f[x]`$ describes the crossover from weak magnetization $`f[x0]x`$ to a saturation of the magnetic moment for large $`H_{\mathrm{AF}}`$ (upper inset of Fig. 4). The electrons scatter off the variations of the (local) magnetization $`\mathrm{\Delta }M_\mathrm{L}`$. Close to the QCP, $`|\mathrm{\Delta }M_\mathrm{L}|`$ increases with increasing staggered bulk magnetization $`M_{\mathrm{AF}}`$. Deep in the AF phase, the local magnetization saturates, $`H_{\mathrm{AF}}\mathrm{max}[T_\mathrm{K}^{\text{imp}},T_\mathrm{K}^{\text{bulk}}]`$, and $`|\mathrm{\Delta }M_\mathrm{L}|`$ decreases in size. This explains qualitatively the pronounced maximum of $`\rho _0(p)`$. It also corresponds to the observation that close to the $`\rho _0(p)`$ maximum, $`T_\mathrm{N}`$ is of the order of the coherence temperature, measured by the maximum in $`\rho (T)`$ . This signals a smooth crossover from a spin-density wave with small magnetic moments just below $`p_c`$ to local-moment AF order with large mo-

ments as $`p`$ approaches zero.

To substantiate this qualitative picture we have performed a microscopic calculation  for an Anderson-lattice model in the limit $`U\mathrm{}`$: $`H=_{ij,\sigma }t_{ij}c_{i\sigma }^{}c_{j\sigma }+Vc_{i\sigma }^{}f_{i\sigma }+ϵ_f^if_i^{}f_i+Uf_i^{}f_if_i^{}f_i`$. We mimic changes of $`p`$ by a variation of the hybridization $`V`$. Increasing the values of $`V`$ corresponds to a $`p`$ increase. The local physics of this strongly correlated system was treated in the dynamical mean-field theory (MFT), the resulting self-consistent Anderson model was solved within slave-boson MFT . While this approximation fails to describe the inelastic processes correctly, it is expected to give a qualitatively correct description for low $`T`$. The disorder is modeled by a small density of $`f`$-electron sites with a different local energy $`ϵ_f^i=ϵ_f^{\text{imp}}ϵ_f^{\text{bulk}}`$. For simplicity, the long-range order is described by a staggered magnetic field with ordering wave vector $`𝐐`$ = (1/2 0 0). In Fig. 4, $`\mathrm{\Delta }\rho _0`$ in a direction perpendicular to $`𝐐`$ (like in the experiment) as a function of $`M_{\mathrm{AF}}`$ shows indeed the maximum which we expected from the arguments given above. The fact that we obtain a sharp cusp instead of a smooth maximum in $`\mathrm{\Delta }\rho _0(M_{\mathrm{AF}})`$ is a well known artefact of the MFT treatment of the single-impurity Anderson model where a spurious phase transition is induced when the magnetization is close to saturation.

The approximation used above neglects non-local effects which should be important in the AF phase close to the QCP, when the correlation length $`\xi `$ is large. We estimate this effect in perturbation theory by considering the scattering from non-local variations of the staggered magnetization $`\mathrm{\Delta }M(𝐫)\chi (𝐫𝐫^{})\mathrm{\Delta }M_\mathrm{L}(𝐫^{})/\chi _\mathrm{L}𝑑𝐫^{}`$, where $`\chi `$ is the susceptibility of the AF phase, characterized by $`\xi `$.

For independent impurities, small $`\mathrm{\Delta }M`$ and $`\mathrm{\Delta }\rho _0\rho _0`$ , we find an enhancement of $`\mathrm{\Delta }\rho _0`$ for large $`\xi `$, $`\mathrm{\Delta }\rho _0\xi |\mathrm{\Delta }M_\mathrm{L}|^2`$. Unfortunately, there are no data on the $`p`$ dependence of the magnetization $`M_{\mathrm{AF}}`$ or $`\xi `$. From MFT, we expect $`\mathrm{\Delta }M_\mathrm{L}\sqrt{p_cp}`$ and $`\xi 1/\sqrt{p_cp}`$ and therefore $`\mathrm{\Delta }\rho _0\sqrt{p_cp}`$ close to $`p_c`$. This is not observed as $`\rho _0(p)`$ appears to rise linearly with $`p_cp`$ and may be attributed to the fact that MFT strongly overestimates the increase of the size of $`M_{\mathrm{AF}}`$ as suggested by the $`x`$-dependence of $`M_{\mathrm{AF}}`$ in CeCu6-xAux . Also the assumption of independent impurities breaks down close to the QCP. Furthermore, $`\xi `$ may be too short in the relevant $`p`$ regime (inelastic neutron scattering at the QCP of CeCu6-xAux reveals relatively broad structures ).

Other potentially important effects like the opening of a gap in some parts of the Fermi surface and an associated reduction of scattering rates  have been included in our model calculation. We do not expect them to dominate the $`p`$ dependence of $`\rho _0`$ measured perpendicular to $`𝐐`$. Likewise, our qualitative picture should hold also for the incommensurate ordering vector $`𝐐`$ = (0.56 0 0) for CeCu5Au . A semi-quantitative check for the local picture is a comparison with the magnetoresistivity at the QCP, where according to our local approximation a similar rise of $`\rho _0`$ is expected for not too strong fields. Indeed, in CeCu5.9Au0.1, $`\rho _0`$ (measured for $`I`$ along the $`b`$ direction) rises by approximately $`20\%`$ for a local magnetization of $`0.3\mu _\mathrm{B}`$ . This is of the expected order of magnitude. A more precise comparison to our pressure data is not possible due to the different nature of the impurities in CeCu5.9Au0.1. Finally, we note that the strong field dependence in the range where $`\rho _0(p)`$ is large, is in qualitative agreement with our model.

We now turn to the pronounced $`\rho _0(p)`$ dependence in the paramagnetic region of the phase diagram which is indeed surprising from a theoretical point of view since in a Kondo lattice with nonmagnetic impurities, $`\rho _0`$ is only weakly renormalized by interactions . It is tempting to interpret the increase of $`\rho _0`$ upon approaching the QCP from the paramagnetic side as the signature of the spontaneous formation of magnetic domains of size $`\xi `$ around certain disorder configurations. Such an interpretation seems not viable as the same strong variation of $`\rho _0`$ is observed at a considerable distance from the QCP, where the number of those domains should be exponentially small. We have therefore theoretically investigated the question whether a strong $`p`$ dependence may arise from purely local effects.

Each impurity, i.e., a Kondo ion with a local environment differing from the bulk due to disorder, contributes according to Friedel’s sum rule a factor $`\mathrm{\Delta }\rho \mathrm{sin}^2\pi \mathrm{\Delta }n`$ to $`\rho _0`$ if $`s`$-wave scattering dominates and the density of impurities is small. $`\mathrm{\Delta }n`$ is the total charge accumulated (in the $`s`$-wave channel) around the impurity. For the following argument we consider an impurity site where the $`f`$-electron is more strongly coupled to the environment than a typical bulk site ($`ϵ_f^{\text{bulk}}<ϵ_f^{\text{imp}}<0`$). For high $`p`$ or large $`V`$, the occupation of the $`f`$ orbital at the impurity site is therefore reduced and $`\mathrm{\Delta }n`$ is negative. For small $`V`$ both the impurity site and bulk sites are deep in the Kondo regime and the number of $`f`$-electrons per site and spin is fixed at $`1/2`$. Nevertheless, extra charge is accumulated around sites with the larger local $`T_\mathrm{K}`$. For a single-impurity Kondo model this effect is of order $`T_\mathrm{K}/E_\mathrm{F}`$, with $`E_\mathrm{F}`$ the Fermi-energy, and therefore negligible. For a lattice $`\mathrm{\Delta }n`$ is of order $`1`$ and positive for a particle-like Fermi surface. In a generic situation one therefore expects a sign change in $`\mathrm{\Delta }n`$ (for a particle-like Fermi surface) and correspondingly a strong $`p`$ dependence of $`\mathrm{\Delta }\rho _0\mathrm{sin}^2\pi \mathrm{\Delta }n`$ as is shown in the lower inset of Fig. 4. A similar argument holds in the case $`ϵ_f^{\text{bulk}}>ϵ_f^{\text{imp}}`$. Remarkably, in the limit $`T_\mathrm{K}^{\text{imp}}T_\mathrm{K}^{\text{bulk}}`$ (or $`T_\mathrm{K}^{\text{imp}}T_\mathrm{K}^{\text{bulk}}`$) $`\mathrm{\Delta }n`$ takes on a universal value which depends only on bulk properties .

With this mechanism, it is possible to explain the unexpected large variation of $`\rho _0(p)`$ in the paramagnetic phase as an intrinsic property of a weakly disordered Kondo lattice. We obtain changes of $`\rho _0`$ by a factor of $`3`$ accompanied by changes of $`V`$ by a few percent and variations of $`T_\mathrm{K}`$ by a factor of 2-3, in reasonably good agreement with experiment. Note that our analysis does not contradict the results of Kotliar and Varma  who find that in a Kondo lattice $`\rho _0`$ is only weakly renormalized by interactions due to nonmagnetic impurities: the strong $`\rho _0(p)`$ variations are indeed “weak” as they are not of order $`1/T_\mathrm{K}`$. A prediction of this scenario is that $`\rho _0`$ rises again for even higher $`p`$ upon entering the mixed-valence regime. From the fact that $`\rho _0(p)`$ in CeCu6 decreases for $`p`$ up to 8.5 GPa , we deduce that in CeCu5Au such an increase, if present, might occur at pressures well above $`10`$ GPa.

In conclusion, we have observed a surprisingly strong variation of the residual resistivity with pressure in both the antiferromagnetic and paramagnetic phase of CeCu5Au. We argued that the pronounced $`\rho _0(p)`$ maximum in the AF phase can be attributed to scattering from inhomogeneities in the magnetization. The explanation for the strong decrease of $`\rho _0`$ in the paramagnetic phase by a factor of three is less obvious: a possible mechanism can nevertheless be found even in a purely local picture where the sign of the charge, which is accumulated in a $`s`$-wave channel, changes under typical conditions

This work was supported by the Swiss National Science Foundation and by the Deutsche Forschungsgemeinschaft through grant Lo 250/20-1 and the Emmy Noether program.

$``$ Present Address: H. H. Wills Physical Laboratory, University of Bristol, Tyndall Avenue, Bristol, BS8 1TL, Great Britain.


text:
# Numerical Results for the Hubbard Model: Implications for the High 𝐓_𝑐 Pairing Mechanism

## Abstract

Numerical studies of the Hubbard model and its strong-coupling form, the $`t`$-$`J`$ model, show evidence for antiferromagnetic, $`d_{x^2y^2}`$-pairing and stripe correlations which remind one of phenomena seen in the layered cuprate materials. Here, we ask what these numerical results imply about various scenarios for the pairing mechanism.

Dedicated to Martin C. Gutzwiller on the occasion of his 75th birthday.

In his 1963 article on the Effect of Correlation on the Ferromagnetism of Transition Metals , Gutzwiller suggested that it was “instructive” to consider a model Hamiltonian of the form

$$H=\underset{ks}{}ϵ_kn_{ks}+U\underset{i}{}n_in_i.$$
(1)
He noted that the first term arose from solving Schrödinger’s equation for one electron in the effective periodic potential of the lattice while the second term described the Coulomb repulsion between two electrons of opposite spin which happen to be in the same orbit at site $`i`$. In this same year, in a paper entitled Electron Correlations in Narrow Bands , Hubbard also studied this model. He wrote $`H`$ in the form

$$H=\underset{ijs}{}t_{ij}\left(c_{i,s}^+c_{j,s}+c_{j,s}^+c_{i,s}\right)+U\underset{i}{}n_in_i$$
(2)
which is now called the Hubbard model. Here $`c_{i,s}^+`$ created an electron of spin $`s`$ in an orbital at site $`i`$ and $`n_{is}=c_{i,s}^+c_{i,s}`$ was the orbital occupation number for electrons of spin $`s`$ on site $`i`$. For a 2D square lattice with $`t_{ij}=t`$ for near-neighbor sites only, $`ϵ_k`$ of Eq. (1) is $`2t(\mathrm{cos}k_x+\mathrm{cos}k_y)`$. In this case, two parameters, $`U/t`$ and the average site occupation $`n`$ determine the state of the system. In his 1963 paper, Gutzwiller also proposed a new trial wave function for large $`U/t`$. This wave function was obtained by applying a projection operator

$$P_G(\eta )=\underset{i}{\Pi }\left(1(1\eta )n_in_i\right)$$
(3)
to the usual Slater determinant of Bloch functions. In particular, for $`\eta =0`$, $`P_G(0)`$ projects out all configurations with doubly occupied sites. The $`t`$-$`J`$ model is the large $`U/t`$ limit of the Hubbard model obtained by expanding to order $`t^2/U`$ and dropping three-site terms which are proportional to the doping. Using the Gutzwiller projection operator, the $`t`$-$`J`$ model has the form

$$H=P_G(0)\left[\underset{ij,s}{}t\left(c_{is}^+c_{js}+c_{js}^+c_{is}\right)+J\underset{ij}{}\left(\stackrel{}{S}_i\stackrel{}{S}_j\frac{1}{4}n_in_j\right)\right]P_G(0)$$
(4)
with the exchange interaction $`J=4t^2/U`$. Here the state of the system is determined by $`J/t`$ and the site occupation $`n`$.

Both Gutzwiller and Hubbard were interested in the strongly-correlated case in which $`U/t`$ was large and in understanding what conditions favored ferromagnetism. Now, almost forty years later, we are still seeking to understand the behavior of the deceptively simple looking Hubbard and $`t`$-$`J`$ models. In particular, we would like to know if these models contain the essential physics of the high $`T_c`$ cuprate superconductors, and if so what they tell us about the basic pairing mechanism. The premise of this brief report is that numerical calculations have shown that the Hubbard and $`t`$-$`J`$ models exhibit a number of the phenomena which have been observed in various cuprate materials. Nevertheless, as we will discuss, the question regarding the nature of the underlying pairing mechanism remains open with a diverse range of proposals set forth.

Turning first to our premise, we know from Monte Carlo calculations that the groundstate of the half-filled 2D Hubbard model has long-range antiferromagnetic order  as does the 2D Heisenberg model , which is just the undoped $`t`$-$`J`$ model. It also has been shown from finite temperature Monte Carlo calculations  that the low energy spin-fluctuations of the insulating cuprates are well described by the Heisenberg model. In particular, the instantaneous spin-spin correlation length calculated for the 2D $`S=1/2`$ Heisenberg model with a near-neighbor exchange is in excellent agreement with neutron scattering measurements of $`La_2CuO_4`$ above its Neel Temperature . Thus, the 2D Hubbard and $`t`$-$`J`$ models provide an excellent description of the low energy magnetic behavior of the undoped layered cuprates. Likewise, DMRG calculations on half-filled $`n`$-leg Hubbard and $`t`$-$`J`$ ladders have shown that the even-leg ladders are spin-gapped while the odd-leg ladders are not . This is in agreement with experimental data on $`SrCu_2O_3`$ and $`Sr_2Cu_3O_5`$ which contain weakly-coupled arrays of 2- and 3-leg ladders respectively .

When a 2-leg Hubbard or $`t`$-$`J`$ ladder is doped away from half-filling, the added holes form $`d_{x^2y^2}`$-like pairs with power law pair-field and CDW correlations and there are physically reasonable parameter regimes in which the pair-field correlations are dominant . Under pressure, the 2-leg ladder material $`(SrCa)_{14}Cu_{24}O_{41}`$ has a superconducting transition  with $`T_c12`$K, however transport evidence suggests that the pressure has made this system into an anisotropic 2D material rather than a weakly Josephson coupled array of 2-leg ladders. On doped 3-leg ladders, calculations indicate that the holes first go into the odd band and for doping $`x.05`$ the system exhibits properties similar to a 1-leg, $`t`$-$`J`$ ladder, which is known to have charge-spin separation . At higher doping in the 3-leg ladder, pairs and domain walls form. When two holes are added to 4- and 6-leg ladders, they form $`d_{x^2y^2}`$-like pairs and at finite doping, domain walls form. On doped 8-leg, $`t`$-$`J`$ ladders, DMRG calculations  find domain walls with a linear hole density of 0.5 per $`Cu`$ separating $`\pi `$-phase shifted antiferromagnetic regions, similar to the stripes observed by neutron scattering  from $`La_{1.6x}Nd_{0.4}Sr_xCuO_4`$. When a small next-near-neighbor hopping $`t^{}`$ is added, $`d_{x^2y^2}`$-pairing correlations are favored over the striped domain wall configurations in the DMRG calculations .

Nevertheless, considerable difference of opinion remains about the nature of the correlations in the 2D system. While early RPA calculations for the doped 2D Hubbard model found that the exchange of antiferromagnetic spin fluctuations could lead to $`d_{x^2y^2}`$-superconductivity , competing channels which arise from nesting and van Hove singularities lead to a breakdown of conventional perturbation theory. A number of authors  have carried out renormalization group analyses which treat the various particle-hole and particle-particle channels on an equal footing. For the weak-coupling doped Hubbard model with a near-neighbor hopping $`t`$, a $`d_{x^2y^2}`$-superconducting state is inferred from the divergence of the one-loop $`d_{x^2y^2}`$-pairing susceptibility. More generally, the one-loop calculation can lead to divergencies of vertices associated with various channels. In particular, with a next-near-neighbor hopping $`t^{}`$, different instabilities arising from the position of Van Hove points on the fermi surface can compete. For example, from an analysis of the renormalization group flow, Furukawa et. al.  have recently suggested that an Umklap-gapped, spin-liquid state might exist for a 2D Hubbard model with appropriate values of $`t^{}`$ and doping.

The original finite temperature Monte Carlo calculations were unable to reach low temperatures for the doped Hubbard model because of the fermion determinantal sign problem. They have, however, been used to study the momentum and Matsubara frequency dependence of the effective pairing interaction at temperatures of order one-half to one-third of the exchange energy scale $`4t^2/U`$. The leading eigenvalue in the singlet channel associated with this interaction was shown to have $`d_{x^2y^2}`$ symmetry . However, zero temperature constrained path Monte Carlo calculations  find only short-range pairing correlations in the groundstate of the 2D Hubbard model.

For the 2D $`t`$-$`J`$ model, questions of phase separation, stripe formation, and pairing arise. Our study of long $`n`$-leg $`t`$-$`J`$ ladders with $`n`$ up to 6 support the notion that the $`t`$-$`J`$ model does not phase-separate in the physically relevant $`J/t`$ region . This conclusion disagrees with Green’s function Monte Carlo calculations of Hellberg and Manousakis  but is in agreement with other Green’s function Monte Carlo calculations with stochastic reconfiguration by Calandra and Sorella . These latter authors also find evidence for long-range superconducting order for dopings $`x>0.1`$ rather than stripes. Thus, while questions remain regarding the existence of long-range superconducting order, stripes, and phase separation, it is clear from a variety of numerical work that the Hubbard and $`t`$-$`J`$ model have low-lying states that exhibit a number of features which are seen in the cuprate materials.

So what can we conclude from this about the basic pairing mechanism? Certainly, in a general way, the fact that these models exhibit antiferromagnetism, $`d_{x^2y^2}`$-pairing correlations and stripes imply that these phenomenon can arise naturally in systems with strong short-range Coulomb interactions in which there is a competition between kinetic and exchange energies: $`d_{x^2y^2}`$ pairs as well as domain walls form as the doped holes locally arrange themselves so as to satisfy the competing requirements of minimizing both their kinetic energy and the disturbance of the background exchange interactions. However, beyond this general picture, forced by the minimalist nature of these models which after all have only $`U/t`$ or $`J/t`$ along with the band filling $`n`$ to parameterize the physics, there remain a remarkably diverse set of views regarding the underlying pairing mechanism.

In part, this diversity reflects the fact that for the parameter regime appropriate to the cuprates ($`U/t1`$ or $`J/t<0.5`$) these models describe a strongly-coupled system which is delicately balanced. Thus, changes in the doping, the lattice or the next-near-neighbor hopping can alter the nature of the groundstate. While we believe that this delicacy in fact provides further evidence that these models contain much of the essential physics, it poses special problems for the theorist. One difficulty with a strongly-interacting system is that when a particular type of mean-field order is assumed, it has a good chance of being self-consistently “found”. That is, strongly-interacting systems can support a self-fulfilling ansatz, keeping the true nature of the groundstate hidden. One goal of the numerical work is to avoid this. However, it is easy to slip back into this framework in interpreting the numerical results. That is, one’s description of “the mechanism” may depend upon one’s starting ansatz. The delicate balance of the system suggests that one must simultaneously treat a variety of channels. In addition, one must ask if the weak-coupling regime of the Hubbard model is continuously connected to the strong-coupling regime. That is, what is the nature of the state out of which the superconducting state is formed as the temperature is lowered or at zero temperature as the doping is changed. Depending upon how one interprets the limited numerical data or what additional interactions or degrees of freedom one imagines including, different pairing scenarios arise.

Here, with this in mind, we’ll examine some of the numerical results. Consider the schematic phase diagram for the 2-leg Hubbard ladder shown in Fig. 1. Here the average site filling $`n`$ is plotted along the horizontal axis and the ratio of the rung hopping $`t_{}`$ to the leg hopping $`t`$ is plotted along the vertical axis. For the non-interacting system, the antibonding band is lifted above the bonding band for $`t_{}/t>2`$. Thus, at half-filling, $`n=1`$, for large values of $`t_{}/t`$ the ladder is simply a typical band insulator with an even number of electrons per unit cell (2 per rung) filling the bonding band. In the absence of the interaction $`U`$, the bonding and antibonding bands of the half-filled system would overlap and it would become a metal for $`t_{}/t<2`$. However, with $`U`$ present it remains a spin-gapped insulator at half-filling. The important point is that for $`n=1`$, the spin-gapped small $`t_{}/t`$ insulating state is adiabatically connected to the large $`t_{}/t`$ Bloch band-insulating state.

The behavior of this system under doping is, of course, quite different depending upon the size of $`t_{}/t`$. In the shaded 2-band region, doped holes enter as quasi-particles with charge $`e`$ and spin $`s=1/2`$. Above a value of $`t_{}/t`$ which, for the interacting system depends upon $`U`$, the doped holes enter a single band Luttinger liquid and exhibit spin-charge separation. In the shaded region of the phase diagram, the doped holes form pairs consisting of superpositions of rung and leg singlets with the $`d_{x^2y^2}`$-like phasing shown in Fig. 2. Here, power-law pairing correlations compete with CDW correlations . How should these results be interpreted? One might have thought that the 2-leg ladder geometry is ideally suited for an RVB description. In the limit in which the rung exchange $`J_{}`$ is large, one can imagine that the groundstate of the undoped ladder contains a set of rung valence-bond singlets with a spin gap $`\frac{3}{4}J_{}`$. Upon hole doping, these latent rung singlet pairs are free to move and power law pairing correlations develop. In the isotropic case in which the leg and rung exchange couplings are equal, the undoped groundstate still has a robust spin gap of order $`J/2`$ and one can again picture that the pairing correlations arise when the system is hole-doped and the singlet pairs (shown in of Fig. 2) are free to move. Thus, there is a continuous transition from the spin-gapped Mott insulating phase to the $`d_{x^2y^2}`$-superconducting phase as the chemical potential exceeds a critical value. However, contrary to some of the original RVB tenents , there are no spinon excitations in the undoped, 2-leg ladder and the type of spin-charge separation found in the one-leg system is absent. In fact, for the 2-leg Hubbard ladder the weak and strong-coupling regions appear to be adiabatically connected.

Another view argues that the exchange of short-range antiferromagnetic spin fluctuations mediate the pairing . Indeed, Monte Carlo calculations  of the effective pairing interaction $`V(q,\omega _m)`$ show that it resembles the spin susceptibility $`\chi (q,\omega _m)`$ as suggested by weak-coupling RPA calculations. Fig. 3 shows Monte Carlo results for $`V(q,0)`$ and $`\chi (q,0)`$ versus $`q_x`$ for $`q_y=\pi `$. As the temperature is lowered and the short-range, spin-spin correlations develop, the pairing interaction $`V`$ shown in Fig. 3(a) exhibits a $`q`$ and $`T`$ dependence similar to $`\chi `$ shown in Fig. 3(b). This same behavior is seen on $`8\times 8`$, 2D clusters  However, because of the fermion sign problem these Monte Carlo calculations are limited to temperatures greater than of order $`J/3`$. Moreover, the simple RPA framework clearly fails to describe the spin-gapped nature of the undoped ladder and an approach such as the RNG-bozonization  method is needed. In this approach, a weak-coupling renormalization group calculation is used which treats all channels on an equal footing until the dominant couplings emerge and can be treated by abelian bosonization. Lin, Balents, and Fisher , and Arrigoni and Hanke  have shown that in weak coupling a generic ladder model, including both longer-range Coulomb interactions and hoppings , flows to a manifold with $`SO(5)`$ symmetry. This is the symmetry group originally proposed by Zhang  in which the antiferromagnetic and superconducting order parameters are combined into a five-dimensional superspin. For an $`SO(5)`$ symmetric ladder, the magnon dispersion about $`(q_x=\pi ,q_y=\pi )`$ should be identical to the hole pair dispersion about $`(q_x=\pi ,q_y=\pi )`$. This was shown to be the case for an explicitly constructed $`SO(5)`$ ladder . However, this need not be the case for the Hubbard and $`t`$-$`J`$ models in the physically relevant region of parameter space. The point is that the RNG results were obtained in the weak-coupling regime. Thus, for example, for the Hubbard ladder, the question is whether there is room for a sufficient renormalization group flow to approach the region in which $`SO(5)`$ provides a useful description of the low energy properties when $`U`$ is of order the bandwidth. A similar question, of course, also arises for the $`t`$-$`J`$ ladder. One test of $`SO(5)`$ is to compare the magnon and pair dispersions which would be identical for an $`SO(5)`$ ladder. The question is whether they are close so that $`SO(5)`$ provides a good starting point for describing the collective modes of the $`t`$-$`J`$ ladder, or whether they are significantly different. With the open boundary conditions on a $`2\times L`$ ladder used in the DMRG calculations, the coefficient of the $`q^2`$ dispersion measured relative to $`(\pi ,\pi )`$ for the magnon and $`(0,0)`$ for the pair, is given by the slope of the $`L^2`$-dependence of the energy $`E_m=E_0(S=1)E_0(S=0)`$ to add a magnon and the energy to add a pair $`E_p=E_0(n_n=2)E_0(n_n=0)`$, respectively. As shown in Fig. 4, for a 2-leg $`t`$-$`J`$ ladder with $`J/t=0.5`$ these two slopes differ by a factor of 2. Thus, in the strong-coupling limit the 2-leg Hubbard model, or here the $`t`$-$`J`$ model, does not exhibit an exact $`SO(5)`$ symmetry. However, recent work on the Hubbard model and $`SO(5)`$ symmetry  has shown that in intermediate to strong coupling, the presence of the on-site $`U`$ requires one to implement a projected $`SO(5)`$ symmetry in which the high-energy double occupancies, responsible for the Mott-Hubbard gap, are Gutzwiller-projected out. This restores $`SO(5)`$ symmetry for static correlation functions but introduces changes in dynamic ones , such as the differences in the magnon and pair dispersions discussed above.

Turning next to the 8-leg, $`t`$-$`J`$ ladder results, the DMRG calculations  for the doped system find evidence for the formation of charged stripes separating $`\pi `$-phase shifted antiferromagnetic regions as illustrated in Fig. 5. Both site-centered “one-leg” and bond-centered, “two-leg” stripes with a linear charge density $`\rho _L=0.5`$ have been found in these calculations, depending upon the length of the open end $`8\times L`$ system. While the existence of striped domain walls was found in early Hartree-Fock calculations , Emery and Kivelson  have argued that stripes arise when phase separation in the $`t`$-$`J`$ model is frustrated by long-range Coulomb forces. In addition, they have suggested that the holes on the stripes exhibit charge-spin separation, at least over some suitable length scale, so that spinons can fluctuate into the “insulating spin-gapped” regions between the stripes. In this spin-gap proximity effect scenario, these spinons then transfer the spin-gap from the “insulating spin-gapped” region back onto the charged stripe leading to enhanced pairing correlations . Ultimately the pairing correlations on different stripes couple through a Josephson coupling to form the superconducting state.

The DMRG calculations certainly do find striped domain walls in a variety of $`n`$-leg ladder calculations. Furthermore, there is a significant spin gap of order $`J/2`$ on an undoped 2-leg ladder which could make such regions attractive candidates for inducing a “spin-gap proximity effect”. Nevertheless, our numerical calculations on the $`t`$-$`J`$ model imply that it does not phase-separate in the physically relevant region of parameter space . Rather, we find that stripe formation in the $`t`$-$`J`$ model is driven by the competition between the kinetic and exchange energies. Furthermore, when we turn on a next-near-neighbor hopping $`t^{}`$, the stripes and stripe-stripe correlations decrease as the pairing correlations increase in strength. The pairing correlations are present when the static stripes have disappeared as well as any evidence for the fluctuation stripe correlations in the density-density correlation function. It appears from these calculations that in this region the stripes have evaporated .

Still other theories propose additional order parameters such as the $`d`$-wave density state which has antiferromagnetic orbital currents. Monte Carlo calculations on a Gutzwiller projected BCS $`d`$-wave variational state are reported to have a power law decay of the orbital currents . It has also been suggested that the pseudogap phase of the cuprates is associated with a $`d`$-wave density phase which is prevented from fully developing by disorder . However, we find that the antiferromagnetic orbital currents decay exponentially on the 2-leg $`t`$-$`J`$ ladder  and also decay rapidly on $`6`$-leg $`t`$-$`J`$ ladders. In addition, were a phase transition to a $`d`$-wave density phase to occur for the lightly-doped 2D Hubbard model at a temperature $`T^{}`$ well above the superconducting transition temperature, one would have expected to see some indication of this in the finite temperature Monte Carlo results.

There has also been a recently proposed $`Z_2`$ gauge theory in which electrons can fractionalize into separate charge and spin degrees of freedom in dimensions greater than one . Here the spin of the electron is carried by a neutral fermionic excitation, the spinon, and the charge is carried by a bosonic excitation, the chargon. Central to this fractionalization is an underlying topological order. If this $`Z_2`$ topological order is present, the groundstate of such a system on a 2D periodic lattice should exhibit a 4-fold degenerate groundstate. Here, we comment on the relationship of such a theory to what is known numerically about the 2D Hubbard and $`t`$-$`J`$ models. For the 2D Heisenberg model with periodic boundary conditions, exact diagonalization calculations on finite sized clusters find a “Neel tower of states” with the groundstate energy $`E_0`$ varying with the total spin as $`S(S+1)`$. This implies that the groundstate of the half-filled $`t`$-$`J`$ model is not a “topological antiferromagnetic” state which would have had an $`S=0`$, 4-fold degenerate groundstate in the limit of a large 2D lattice with periodic boundary conditions. Similarly, in the weak-coupling limit of the half-filled 2D Hubbard model, renormalization group calculations imply that the system has a Neel antiferromagnetic groundstate and electron-like quasi-particles above a gap. Since, in strong coupling, the Hubbard model goes to the $`t`$-$`J`$ model, which at half-filing is Heisenberg-like, it seems likely that the groundstate of the half-filled Hubbard model is the traditional antiferromagnetic state for all $`U/t`$. It may be that additional 4-site ring exchange interactions, or other interactions laying outside the simple near-neighbor exchange-coupled Heisenberg model, can lead to such a topological state. This remains to be studied.

For the doped 2D Hubbard and $`t`$-$`J`$ models, it is possible that a topological phase exists. This would, of course, be a key finding and would lay the groundwork for a spin-charge fractionalization mechanism in which the electron separates into a spinon which carries the electron’s spin and fermi statistics and a chargon which can bose condense giving superconductivity . While it seems unlikely that this can happen at infinitesimal doping, it could be that it occurs at some finite doping. At present, our numerical work on the Hubbard and $`t`$-$`J`$ models has failed to find evidence of such a topological state. However, these studies are limited and there are of course additional interactions present. For example, we know that the cuprates are actually charge transfer insulators rather than Mott-Hubbard insulators. While we believe that this enhances the pairing tendencies that one sees in the Hubbard and $`t`$-$`J`$ models due to the extra exchange paths and increased pair mobility , it is possible that these additional interactions produce more profound effects.

To summarize, the Hubbard model and its strong-coupling limit, the $`t`$-$`J`$ model, have certainly proven “instructive” as orginally suggested by Gutzwiller. The numerical results exhibit a variety of antiferromagnetic, $`d_{x^2y^2}`$-pairing and stripe correlations reminding us of the phenomena seen in the cuprates and providing tests for various pairing scenarios. The numerical results suggest that (1) the competing requirements of minimizing the kinetic energy and the disturbance of the exchange energy background underlie both pairing and stripe domain wall formation, (2) the momentum and temperature dependence of the effective pairing is similar to the spin susceptibility, and (3) stripes compete with superconductivity.

###### Acknowledgements.

We are grateful to Leon Balents and Matthew Fisher for useful discussions and insights. D.J. Scalapino acknowledges support from the NSF under grant # DMR98-17242. S.R. White acknowledges support from the NSF under grant # DMR98-70930.


text:
# Light Curve and Spectral Models for the Hypernova SN 1998bw associated with GRB980425

## 1. INTRODUCTION

On 1998 April 25, GRB980425 triggered the Narrow Field Instrument (NFI) and Wide Field Camera (WFC) detectors on board BeppoSAX (Soffita et al. 1998). SN 1998bw was discovered within the WFC error box in the optical (Galama et al. 1998a) and radio wavelength bands (Kulkarni et al. 1998a) only 0.9 and 3 days after the date of the GRB, respectively. The X-ray afterglow detected by the NFI within the error box of GRB980425 is consistent with SN 1998bw (Pian et al. 1999). The small likelihood of finding a supernova and a GRB in such a small field over such a short interval of time suggests that SN 1998bw and GRB980425 are related (Galama et al. 1998b; Kulkarni et al. 1998b).

The early optical spectra of SN 1998bw apparently lack dominant line features, displaying only some broad emissions at 500, 620 and 800 nm (Galama et al. 1998b; Stathakis et al. 2000; Patat et al. 2000). These were shown to be the emission components of P-Cygni profiles of very broad line blends. They are caused mostly by FeII lines in the blue, by SiII near 600 nm and by OI + CaII near 720 nm (Iwamoto et al. 1998, hereafter IMN98). The absence of any hydrogen line after subtraction of the galaxy background and the fact that the Si II line at 615 nm has a very large velocity indicate that the supernova is neither a Type II nor a normal Type Ia (Galama et al. 1998b). Following the conventional classification scheme (Filippenko 1997), the lack of strong He I features led Patat & Piemonte (1998) to conclude that SN 1998bw is likely to be a Type Ic supernova (SN Ic) rather than a Type Ib.

The light curve of SN 1998bw (Galama et al. 1998b) showed a very early rise and reached a peak at $``$ 17 days (in the V band) after the explosion, and then has been declining exponentially with time (McKenzie & Schaefer 1999; Patat et al. 2000). This clearly indicates that the light curve is not a typical optical afterglow of a gamma-ray burst, but it is powered by the radioactive decay of 56Ni and 56Co as in usual supernovae (IMN98). The distance modulus to SN 1998bw is estimated as $`\mu =32.89`$ mag, so that the peak luminosity of SN 1998bw is $`1\times 10^{43}`$ ergs sec-1, which is about ten times brighter than typical SNe Ib/Ic (Clocchiatti & Wheeler 1997). To achieve such a high luminosity, a large amount of 56Ni must have been synthesized in SN1998bw (IMN98; Woosley, Eastman & Schmidt 1999), again about ten times as much as that produced in typical core-collapse-induced supernovae. The very broad spectral features and the light curve shape led various groups to the conclusion that SN 1998bw had a very large kinetic energy of explosion $`E_\mathrm{K}`$ (IMN98; Woosley et al. 1999; Branch 2000).

IMN98 constructed models of the core-collapse-induced explosion of C+O cores of initially massive stars that had lost their hydrogen and helium-rich layers before the explosion. Among those models, the energetic explosion ($`E_\mathrm{K}(25)\times 10^{52}`$ ergs) of a C+O star of 13.8 $`M_{}`$ successfully fit the first 60 days of the light curve and the early spectra, although the synthetic spectral lines were still narrower than the very broad observed features (IMN98). Since the kinetic energy is more than one order of magnitude larger than the energy of typical supernovae, SN 1998bw was called a “hypernova” (IMN98), a term we use to describe events with $`E_\mathrm{K}\text{ }>10^{52}`$ erg without specifying whether the central engine is a collapser (MacFadyen & Woosley 1999), a magnetar (Nakamura 1998; Wheeler et al. 2000) or a pair-instability supernova.

Interestingly, photometry after $`60`$ days showed that SN 1998bw declined significantly more slowly than the rate predicted by the model of IMN98 (McKenzie & Schaefer 1999; Patat et al. 2000). Also the bolometric light curve has been constructed up to day 500 (Patat et al. 2000; Sollerman et al. 2000). Therefore, we have recomputed the light curve of SN 1998bw using progenitors of different masses and explosions of different energies, and provide better estimates for these parameters.

The hydrodynamical models for SN 1998bw are described in §2. Sections 3 and 4 are devoted to light curve models and synthetic spectra models, respectively. Explosive nucleosynthesis in SN 1998bw is discussed in §5. Finally, the nature of this peculiar supernova is summarized in §6, with the emphasis being placed on possible evidence that the explosion was aspherical. Preliminary results have already been reported by Nakamura et al. (1999a, 2000) and Nomoto et al. (2000ab).

## 2. Hydrodynamical models for SN 1998bw

Hydrodynamical models are constructed as follows. C+O stars are chosen as progenitors as in IMN98. Light curves and spectra are computed for various C+O star models with different values of the kinetic energy $`E_\mathrm{K}`$ and the ejecta mass $`M_{\mathrm{ej}}`$. These parameters can be constrained by comparing calculated light curves, synthetic spectra, and photospheric velocities with the observational data of SN 1998bw.

1. In the ordinary, low energy SN Ic model (model CO138E1), a C+O star with a mass $`M_{\mathrm{CO}}=13.8M_{}`$ (which is the core of a 40 $`M_{}`$ main-sequence star; Nomoto & Hashimoto 1988; Nomoto et al. 1997) explodes with $`E_\mathrm{K}=1.0\times 10^{51}`$ ergs and $`M_{\mathrm{ej}}=M_{\mathrm{CO}}M_{\mathrm{cut}}12M_{}`$. $`M_{\mathrm{cut}}`$ (= 2 $`M_{}`$ in this case) denotes the position of the mass cut, which corresponds to the mass of the compact star remnant. This is either a neutron star or a black hole, depending on $`M_{\mathrm{cut}}`$.
2. For the hypernova models CO138E50, CO138E30, and CO138E7, the progenitor C+O star of $`M_{\mathrm{CO}}=13.8M_{}`$ is the same as in CO138E1. These models have different explosion energies: $`E_\mathrm{K}=5\times 10^{52}`$ erg (CO138E50), $`3\times 10^{52}`$ erg (CO138E30), and $`7\times 10^{51}`$ erg (CO138E7). The ejecta mass is $`M_{\mathrm{ej}}1011.5M_{}`$, i.e., $`M_{\mathrm{cut}}`$ 2.5 - 4 $`M_{}`$. The parameters of the models are summarized in Table 1. The position of the mass cut is chosen so that the ejected mass of 56Ni is the value required to explain the observed peak brightness of SN 1998bw by radioactive decay heating. The compact remnant in these hypernova models may well be a black hole, because $`M_{\mathrm{cut}}`$ can exceed, sometimes significantly, the maximum mass of a stable neutron star.

The hydrodynamics at early phases was calculated using a Lagrangian PPM code (Colella & Woodward 1984). All models are spherically symmetric. The explosion is triggered by depositing thermal energy in a couple of zones just below the mass cut so that the final kinetic energy has the required value. A strong shock wave forms and propagates toward the surface. Explosive nucleosynthesis takes place behind the shock wave. Radioactive 56Ni is produced in the deep, low velocity layers of the ejecta.

In SN II and SN Ib models, it has been demonstrated that Rayleigh-Taylor instabilities develop at the H/He and He/C+O interfaces and induce mixing of elements in the ejecta (e.g. Arnett et al. 1989; Hachisu et al. 1991, 1994; Iwamoto et al. 1993). Although bare C+O stars, the progenitors of SNe Ic, lack composition interfaces with such pronounced density jumps, Rayleigh-Taylor instabilities can be driven by neutrino heating and develop at the Ni+Si/O interface (Kifonidis et al. 2000). In addition, polarization measurements suggested that SNe Ic, and hypernovae in particular, are asymmetric explosions (Wang et al. 1999; Danziger et al. 1999; Patat et al. 2000). The asymmetry of the explosion may reinforce the instability and bring heavy elements up to high velocity layers. This is particularly true of jet-like explosions, where it is likely that extensive mixing takes place in velocity space. In view of the large uncertainties in our knowledge of the mixing process, we assume that the material ejected is uniformly mixed out to a velocity $`v=v_{\mathrm{mix}}`$.

The hydrodynamical models become homologous ($`vr`$) at $`t250`$ sec, and are then used as input for a radiation transfer code. The lines in Figure 1 show the density–velocity distribution of the homologously expanding ejecta (top) and the enclosed mass $`M_r`$ as a function of velocity (bottom) for models CO138E50 (solid), CO138E30 (long-dashed), CO138E7 (short-dashed), and CO138E1 (dash-dotted).

## 3. Light curve models

Synthetic light curves are computed with a radiative transfer code (Iwamoto et al. 2000) which takes into account the balance between photo-ionizations and recombinations and includes a simplified treatment of line opacity. The width of the light curve peak, $`\tau _{\mathrm{LC}}`$, depends on $`E_\mathrm{K}`$ and $`M_{\mathrm{ej}}`$ approximately as $`\tau _{\mathrm{LC}}(\kappa /c)^{1/2}M_{\mathrm{ej}}^{3/4}E_\mathrm{K}^{1/4}`$, where $`\kappa `$ and $`c`$ are the optical opacity and the speed of light, respectively (Arnett 1982). This means that the light curve can be reproduced with different explosion models that have the same values of $`M_{\mathrm{ej}}^3E_\mathrm{K}^1`$. However, these parameters can be further constrained from both photospheric velocities and spectra because the velocity scales roughly as $`M_{\mathrm{ej}}^{1/2}E_\mathrm{K}^{1/2}`$. The light curve shape depends also on the distribution of the radioactive heating source 56Ni, for which we examine the dependence on $`v_{\mathrm{mix}}`$2).

Our calculations are compared with the bolometric light curve of SN 1998bw constructed by Patat et al. (2000), which uses a redshift distance $`\mu =32.89`$ mag ($`d=37.8`$ Mpc, $`H_0=65`$ km s-1 Mpc-1) and $`A_V=0.05`$ mag. The time of core collapse is set at the detection of GRB 980425. Figure 2 shows the bolometric light curves of the energetic models CO138E50 (solid) and CO138E30 (long-dashed), and Figure 3 the less energetic models CO138E7 (short-dashed) and CO138E1 (dash-dotted).

Figure 4 shows the evolution of the calculated photospheric velocities of CO138E50 (solid line), CO138E30 (long-dashed line), CO138E7 (short-dashed line), and CO138E1 (dash-dotted line) compared with those obtained from spectral models (filled circles), and with the observed velocity of the Si II 635.5 nm doublet measured in the spectra at the absorption core (open circles; Patat et al. 2000), and that of the Ca II H + K doublet measured in the spectrum of May 23 (square; Patat & Piemonte (1998)).

### 3.1. Early Phase

The early part of the light curve ($`t\text{ }<25`$ days) is well reproduced by the two energetic models, CO138E50 and CO138E30. The difference between the two curves is too small to see on the scale of Figure 2 because the light curve width scales as $`E_\mathrm{K}^{1/4}`$. The early part of the light curve also depends on the 56Ni distribution in the ejecta. To reproduce the early sharp rise of the light curve, extensive mixing of 56Ni, out to $`v_{\mathrm{mix}}=22,000`$ km s-1 for CO138E50 and $`v_{\mathrm{mix}}=30,000`$ km s-1 for CO138E30, is required. A better fit can be obtained if we adopt the following ad hoc 56Ni distribution: X(56Ni) = 0.79 at $`v11,000`$ km s-1, 0.011 at $`v=`$ 11 - 17,000 km s-1, and 0.032 at $`v=`$ 17 - 40,000 km s-1, where X(56Ni) denotes the mass fraction of 56Ni. The result is shown in Figure 5. Such a non-uniform distribution of 56Ni might be due to the mixing caused by Rayleigh-Taylor instability (§2), or reflect a complicated structure of a possible jet-like ejecta (§6). The 56Ni mass is determined to be 0.4 $`M_{}`$ from the fitting of the maximum luminosity and the date of maximum. (Note that IMN98 slightly overestimated the mass of 56Ni because they assumed that the bolometric correction is negligible and adopted a different absorption and distance.)

In contrast, the early light curves of the less energetic models (CO138E7 and CO138E1) evolve too slowly, reaching the maximum too late as compared with the observations (Figure 3), even if 56Ni is distributed uniformly throughout the ejecta.

The photospheric velocities $`v_{\mathrm{ph}}`$ provide clearer diagnostics to distinguish between CO138E50 and CO138E30. In particular, in the early phase ($`t<`$ 20 days), $`v_{\mathrm{ph}}`$ of CO138E50 is in good agreement with the observed $`v_{\mathrm{ph}}`$, while $`v_{\mathrm{ph}}`$ of CO138E30 is clearly too low to be consistent with the observations (Fig. 4).

In the early phase, therefore, model CO138E50 shows the best agreement with both the light curve and the photospheric velocities of SN 1998bw. The synthetic spectra also require the most energetic model CO138E50 as will be shown in §4. CO138E50 has the same mass as the best model in IMN98 but a larger $`E_\mathrm{K}`$, which is necessary especially to improve the fit to the spectra (§4).

### 3.2. Intermediate Phase

At intermediate phases ($`t25`$ \- $`200`$ days), the decline rate of the light curve is determined mainly by the fraction of the $`\gamma `$-rays emitted by 56Co decay which are trapped in the ejecta. The optical depth of the ejecta to the $`\gamma `$-rays scales as $`\kappa _\gamma \rho RMR^2M^2E_\mathrm{K}^1t^2`$. Thus the behaviour of the various models can be seen more easily than at earlier phases as follows:

1. The light curve of CO138E50 is consistent with SN1998bw until day 50 but declines faster than the observation afterwards (Figure 2; McKenzie & Schaefer 1999; Patat et al. 2000). The light curve of CO138E30 declines more slowly than that of CO138E50, and is in better agreement with SN1998bw, but it still declines faster than the observations.
2. The photospheric velocity $`v_{\mathrm{ph}}`$ shows a similar tendency. Figure 4 shows that $`v_{\mathrm{ph}}`$ of CO138E50 is the best fit to the data for the first 20 days, but afterwards $`v_{\mathrm{ph}}`$ of CO138E30 gives as good a fit. The other two models do not fit $`v_{\mathrm{ph}}`$ at all.
3. Figure 3 shows that the apparently exponential decline of SN1998bw after day 60 (McKenzie & Schaefer 1999; Patat et al. 2000) is well reproduced by the lower energy model CO138E7. In this model, $`\gamma `$-ray trapping is more efficient than CO138E50 and CO138E30, but the observed low flux level at this phase requires a reduced 56Ni mass of $`0.28M_{}`$ (dotted line). Note, however, that $`0.28M_{}`$ of 56Ni is too small to reproduce the observed light curve maximum. Note also that because of the interplay of $`E_K`$ and $`M`$(Ni), it is difficult to establish those parameters from the intermediate phase light curve alone.
4. The normal SN Ic model CO138E1 in Figure 3 is slow enough to trap most of the $`\gamma `$-rays emitted from the 56Co decay. The decline of the light curve is therefore too slow compared with the observed rate.

These comparisons between SN 1998bw and the model light curves and $`v_{\mathrm{ph}}`$ in Figures 2, 3, & 4 indicate that $`\gamma `$-ray deposition after about day 50 in SN 1998bw is more efficient than predicted by CO138E50 and E30, although these models give an appropriate description of the early light curve and spectra. Higher deposition can be achieved if there exists a significant amount of low-velocity and high-density material. Indications of the presence of a low-velocity, high-density region suggest that the ejecta distribution is not spherically symmetric as will be discussed in §6. We note that such indications are also found for SN 1997ef (Iwamoto et al. 2000; Mazzali, Iwamoto, & Nomoto 2000a), a lower-energy analogue of SN 1998bw.

### 3.3. Late Phase

After day $``$ 200 the decline of the model light curve becomes slower, and it approaches the half-life of 56Co decay around day 400. At $`t\text{ }>400`$ days most $`\gamma `$-rays escape from the ejecta. The $`\gamma `$-ray deposition fraction at 400 days is 1%, 1.5%, 6%, and 13% for CO138E50, E30, E7, and E1, respectively. On the other hand, kinetic energies of positrons emitted from the decay of 56Co are supposed to be fully thermalized because of the postulated weak magnetic field (e.g., Colgate & Petchek 1979). Therefore, positron kinetic energy deposition determines the light curve at $`t\text{ }>400`$ days (dotted line in Figure 2). Here the luminosity by positron deposition is given as

$`L(\mathrm{Co},e^+)=1.4\times 10^{43}\mathrm{erg}/\mathrm{s}(M(\mathrm{Ni})/M_{})`$
$`\times \mathrm{exp}(t/111.26(\mathrm{day}))\times 0.035,`$ (1)
where the positron fraction in energy is 3.5 percent (e.g., Axelrod 1980). (See Cappellaro et al. 1997 and Milne et al. 1999 for the light curves including positron escape.) If the observed tail should follow the positron-powered light curve, the 56Co mass could be determined directly. Since positron deposition should occur almost on the spot, this determination does not depend much on any asphericity of the ejecta.

On June 11, 2000, which corresponds to a SN epoch of $`t=`$ 778 days, HST observations detected a point-like source at the position of SN1998bw (Fynbo et al. 2000). The observed magnitude (V-mag 25.41 $`\pm `$ 0.25) is consistent with the prediction of CO138E7 (Figure 2) but brighter than CO138E50 (Figure 3). The observed luminosity would even be higher than CO138E7, if the bolometric correction was significant. In any case, these comparisons suggest that the positron contribution is not dominant yet around day 800.

## 4. Synthetic spectra

In Figure 6 we show the synthetic spectra obtained for the same 3 epochs fitted in IMN98 (continuous lines) compared with spectra observed at ESO (bold lines). We used model CO138E50 and computed synthetic spectra with a Monte Carlo model (Mazzali & Lucy 1993), improved with the inclusion of photon branching and a new extended and improved line list (Lucy 1999; Mazzali 2000). The synthetic spectra were computed using the luminosity derived from the light curve, a distance modulus of $`\mu =32.89`$ mag and $`A_V=0.05`$. The assumption of low reddening is supported by the upper limit of 0.1 Å in the equivalent width of the Na i D line obtained from high-resolution spectra (Patat et al. 2000). The observed spectra used here (solid lines) are the ‘definitive’, fully reduced version of the same ESO spectra shown in IMN98, and are calibrated with respect to the V photometry. The residual correction factors for the other bands are usually very close to 1, but they are $`1.1`$ for B in the May 11 and 23 spectra. Therefore, the new models (i.e. 3 epochs of CO138E50) have somewhat different parameters than those of IMN98. Both the luminosities and the photospheric velocities are larger than in IMN98. The photospheric velocity is now in better agreement with the measured velocity of the Si ii line (IMN98, Figure 3).

The synthetic spectra clearly improve over those of IMN98, especially at the earliest epochs. Absorptions not caused by broad blends of many lines of moderate strength, such as the Si ii feature near 6000Å and, in particular, the O i+Ca ii feature between 7000 and 8000Å, are now much broader, in significantly better agreement with the data. Nevertheless, the blue sides of those absorptions are still too narrow, indicating that even the new model CO138E50 may not contain enough mass at the highest velocities.

Therefore we introduced an arbitrary change to the original CO138E50 density structure. Several possibilities were tested, and improved results were obtained when the density slope was reduced from $`\rho r^8`$ to $`\rho r^6`$ at $`v>30,000`$ km s-1. This does not introduce a significant change in $`M_{\mathrm{ej}}`$, and increases $`E_\mathrm{K}`$ by only about 10%, but it does increase the density at high velocities, leading to significant absorption at $`v60,000`$ km s-1 in the strongest lines, especially the Ca ii IR triplet, extending the absorption troughs to the blue. The corresponding synthetic spectra are shown as the dotted lines in Figure 6. The effect of the change is of course largest at the earliest epochs. Although the overall agreement with the observed spectra is better, several problems remain, the most severe of which is clearly the excessive strength of the O i line at 7200Å on May 11 and 23. The composition of the highest velocity ejecta is dominated by O in our 1D models, and it is difficult to make that line become weaker. On 23 May, the synthetic Ca ii IR triplet matches the weak feature at 8000Å, which is first seen on 11 May and which continues to grow until it finally causes the wavelength of the absorption minimum of the entire broad feature to shift to $`8200`$Å (Patat et al. 2000). This is rather a peculiar behavior, because on 3 May the O i and Ca ii lines had to blend much more to give rise to the observed broad feature, which then had a minimum at 7000Å. The synthetic O i line is too strong and too fast. The core of the broad feature, if it is interpreted as O i 7774Å, indicates a velocity of 10500 km s-1 on May 11 and 6000 km s-1 on May 23. This is significantly lower than the corresponding velocity of the model photosphere. Actually, the definition of a photosphere at very red wavelengths is not very accurate, because the density of spectral lines is low and so line opacity does not define a pseudo-continuum in that region. Therefore, the observations probably indicate that a large fraction of the O is located at low velocities.

A very flat ($`\rho r^2`$) density distribution was also used by Branch (2000) to fit the spectrum of SN 1998bw. This dependence is however too flat when we use our MC model, because the ionization of, e.g., Ca ii does not fall as steeply as he assumed. On the other hand, Branch’s value of $`E_\mathrm{K}`$ ($`5\times 10^{52}`$ erg) is similar to ours, but he quotes a mass of $`6M_{}`$ above 7000 km s-1, while in our case the mass above that velocity is as large as $`10M_{}`$. Such a flat density distribution at high velocities is also required to fit the spectrum of another hypernova, SN 1997ef (Mazzali et al. 2000a). This might indicate that the progenitors of these hypernovae underwent very extensive mass loss; the outer density structure of the ejecta is flatter than that of ordinary giants, and comparable to that of a mass-losing star.

Clearly, a definitive solution has not been found yet. It is quite possible that only by taking into account departures from spherical symmetry will it be possible to obtain a really accurate fit to the spectra. Nevertheless, considering the complexity of the problem, our fits at least demonstrate that a large $`E_\mathrm{K}`$ is necessary and that the O-dominated composition of the SN Ic model yields quite a reasonable reproduction of the observations.

## 5. Explosive Nucleosynthesis

We calculated explosive nucleosynthesis using a detailed nuclear reaction network (Thielemann, Nomoto & Hashimoto 1996; Nakamura et al. 1999b). Our calculations are performed in two steps. The first step is a hydrodynamical simulation of the explosion with a small nuclear reaction network containing only 13 alpha nuclei (4He, 12C, 16O, 20Ne, 24Mg, 28Si, 32S, 36Ar, 40Ca, 44Ti, 48Cr, 52Fe, and 56Ni), as described in §2. In the second step, post-processing calculations are performed at each mesh point of the hydrodynamical model with the extended reaction network (Hix & Thielemann 1996), which contains 211 isotopes up to 71Ge.

The top panel of Figure 7 shows the isotopic composition of the ejecta of the hypernova model CO138E50 as a function of the enclosed mass $`M_r`$ and of the expansion velocity. The nucleosynthesis in other hypernova and supernova models (CO138E30, CO138E7, CO138E1) is also shown in Figure 7 for comparison. The yields of the hypernova and supernova models are summarized in Table 2. Table 3 and Figure 8 give respectively more detailed yields and the abundances of stable isotopes relative to the solar values for model CO138E50. From these figures and tables, we note the following characteristics of nucleosynthesis of hyper-energetic explosions compared with normal energy explosions.

1. The complete Si-burning region where 56Ni is produced is extended further out (in mass coordinates) as the explosion energy increases. How much processed matter is ejected from this region depends on the mass cut. Compared with normal core-collapse supernovae, the much larger amount of 56Ni ($`0.4M_{}`$) observed in SN1998bw implies that the mass cut is deeper, so that the elements synthesized in this region, such as 59Cu, 63Zn, and 64Ge (which decay into 59Co, 63Cu, and 64Zn, respectively), are ejected more abundantly. Among hypernova models, more energetic models produce more 56Ni in the incomplete Si burning region (See (3) below), and thus $`M_{\mathrm{cut}}`$ is larger to eject $`0.4M_{}`$ 56Ni as constrained from the light curve of SN 1998bw (Figure 7).
2. In the complete Si-burning region of the hypernova models (CO138E50/E30), elements produced by $`\alpha `$-rich freezeout are enhanced because nucleosynthesis proceeds at lower densities than in CO138E1. Figure 7 clearly shows a trend that a larger amount of 4He is left in more energetic explosions. Hence, the mass fractions of the species synthesized through $`\alpha `$-particle capture, such as 44Ti and 48Cr (which decay into 44Ca and 48Ti, respectively) is larger in CO138E50/E30/E7 than CO138E1. The integrated mass of these species depends on $`M_{\mathrm{cut}}`$. For CO138E50, the ejected mass of 44Ti is smaller than in other models because of its larger $`M_{\mathrm{cut}}`$. Note that the 4He produced even in the most energetic models has a velocity significantly smaller than that of the He shell identified in SN 1998bw at a velocity of 18300km s-1 by Patat et al. (2000) on the basis of near-IR spectra. Spectral evidence for low velocity He in SN 1998bw is unclear, however the distribution of He depends on mixing in velocity space.
3. The incomplete Si-burning region is more massive in more energetic explosions. The main products in this region are 28Si, 32S, and 56Ni. CO138E50 produces 0.2$`M_{}`$ of 56Ni in this region. Other important species such as 52Fe, 55Co, and 51Mn (decaying into 52Cr, 55Mn, and 51V, respectively) are synthesized more abundantly in the more energetic explosions.
4. For the larger explosion energy, oxygen burning takes place in more extended, lower density regions. O, C, and Al are burned more efficiently in these cases, and the abundances of the elements in the ejecta are smaller, while a larger amount of ash products such as Si, S and Ar is synthesized by oxygen burning.

## 6. Conclusions and Discussion

In this paper, we have presented a model for SN 1998bw which is in better agreement with the early observations than the previous model in IMN98. Models with different $`E_\mathrm{K}`$ yield different synthetic spectra, and by comparing with the observed early-time spectra of SN 1998bw and trying to fit the very broad absorption features, we selected model CO138E50 with $`M_{\mathrm{ej}}=10M_{}`$ and $`E_\mathrm{K}=5\times 10^{52}`$ erg as the best match to the early data. The large value of $`E_\mathrm{K}`$ qualifies SN 1998bw as ‘the’ Type Ic Hypernova. (SN IIn 1997cy may be called a “Type IIn Hypernova”; Germany et al. 2000; Turatto et al. 2000). The mass of the progenitor C+O star is $`13.8M_{}`$, corresponding to a main sequence mass of $`40M_{}`$. All models require $`M`$(56Ni)$`0.4M_{}`$ to power the bright light curve peak 111 After submitting this paper, Sollerman et al.’s paper (2000) was published. They found that 56Ni of $`0.3`$ \- $`0.9M_{}`$ was necessary to power the late light curve based on the model CO138 (IMN98), which is consistent with our results.. This is about an order of magnitude larger than in typical core-collapse SNe. The compact remnant is probably a black hole, because its mass exceeds $`3M_{}`$ as constrained from the mass of 56Ni.

Although the early light curve of SN1998bw ($`t\text{ }<`$ 50 days) is reproduced well by the most energetic model CO138E50, the observed tail declines more slowly than this model does (§3). The lower energy model CO138E30 is in better agreement with observations at $`t\text{ }<`$ 100 days but it still declines is too fast. Model CO138E7 has a slower tail and can reproduce the observed light curve tail if a smaller 56Ni mass of 0.28 $`M_{}`$ is adopted. This suggests that there might be a high density core with low velocities in SN1998bw where the $`\gamma `$-rays deposit efficiently.

Such a dense core could be formed in spherically symmetric models if the exploding star had a massive He envelope so that a strong reverse shock was formed at the C+O/He interface and largely decelerated the inner core as found in SN Ib models (Hachisu et al. 1991, 1994). However, the apparent absence of the He I lines in the optical spectra in SN 1998bw (Patat et al. 2000) is not consistent with the presence of such a massive He envelope. Also late time energy input from 56Co bubbles might create a non-homologous structure. However, the energy from radioactive decays is much smaller than the large explosion energy of SN 1998bw ($`10^{52}`$ ergs), so that its effect is negligible.

We suggest that the peculiar density distribution is the result of a non-spherically symmetric explosion. If the outburst in SN1998bw took the form of a prolate spheroid, for example, the explosive shock was probably strong along the long axis, ejecting material with large velocities and producing abundant 56Ni, which might have caused the early bright light curve. In directions away from the long axis, on the other hand, oxygen would not be consumed, and the density could be high enough for $`\gamma `$-rays to be trapped even at advanced phases, thus giving rise to the slowly declining tail. These features can be seen in the hydrodynamical models of jet-like explosions (Maeda et al. 2000; MacFadyen & Woosley 1999; Khokhlov et al. 1999; Nagataki et al. 1997). The fact that the late light curve is fitted by model CO138E7 with 0.28 $`M_{}`$ of 56Ni, while the early light curve requires as much as 0.4 $`M_{}`$, could also support this suggestion, and point at the geometrical ratio between the extension of the high- and low-density regions.

There are several observations which support the above non-spherical explosion scenario. The observed polarization ($`0.5`$%) of the early optical light (Patat et al. 2000; IMN98) suggests that the ejecta of SN 1998bw is aspherical. The observed abnormal distribution of elements in velocity space as seen in the emission line profiles in the nebular phase, with significant amounts of O being located at lower velocity than Fe (Patat et al. 2000; Danziger et al. 1999; Nomoto et al. 2000) can better be explained if the explosion is non-spherically symmetric (Maeda et al. 2000).

The possible connection between SN 1998bw and GRB980425 also supports the conjecture that SN 1998bw was aspherical. The energy of the photons produced by synchrotron emission at the relativistic shock is approximately given by $`h\nu 160`$ $`\mathrm{keV}`$ $`(\mathrm{\Gamma }/100)^4n_1^{1/2}`$ (Piran 1999), where $`\mathrm{\Gamma }`$ is the Lorentz factor of the shock and $`n_1`$ (cm-3) is the density of the interstellar matter. In order to produce an observable GRB, $`\mathrm{\Gamma }`$ should be as large as $`\mathrm{\Gamma }100`$. However, even the most energetic model (CO138E50) has only a very small mass of relativistic ejecta ($`10^{10}M_{}`$ with $`\mathrm{\Gamma }100`$), although relativistic hydrodynamical calculations are necessary to obtain the accurate mass. Such a small amount of relativistic material, which is consistent with previous estimates (IMN98; Woosley et al. 1999), is not enough to produce GRB980425 in a spherically symmetric model. However, if the explosion is axi-symmetric, for instance, the energy can be carried by only a small fraction of the material, which might then attain a large Lorentz factor.

Regarding asphericity, we note that Danziger et al. (1999) and Mazzali et al. (2000b; see also Nomoto et al. 2000a) estimated a 56Ni mass of 0.35 - 0.65 $`M_{}`$ from the nebular lines of Fe of SN 1998bw. This estimate does not depend much on the asphericity, and is in good agreement with the 56Ni mass of the spherical models CO138. On the other hand, Höflich et al. (1999) suggested that the 56Ni mass can be as small as 0.2 $`M_{}`$ if aspherical effects are large. Because the difference between these results is not so large, aspherical effects might be modest in SN 1998bw. If the 56Ni mass could be determined more accurately from the late observations, it would provide a good measure of the degree of the asphericity.

We would like to thank Ferdinando Patat and Enrico Cappellaro for useful information and discussion on the observations of SN1998bw. This work has been supported by the grant-in-Aid for Scientific Research (12640233, 12740122) and COE research (07CE2002) of the Japanese Ministry of Education, Science, Culture, and Sports.