format: text

text:
# Macroscopic entanglement jumps in model spin systems

## Abstract

In this paper, we consider some frustrated spin models for which the ground states are known exactly. The concurrence, a measure of the amount of entanglement can be calculated exactly for entangled spin pairs. Quantum phase transitions involving macroscopic magnetization changes at critical values of the magnetic field are accompanied by macroscopic jumps in the $`T=0`$ entanglement. A specific example is given in which magnetization plateaus give rise to a plateau structure in the amount of entanglement associated with nearest-neighbour bonds. We further show that macroscopic entanglement changes can occur in quantum phase transitions brought about by the tuning of exchange interaction strengths.
PACS Nos.: 03.65.Ud, 75.10.Jm

Physics Department, Bose Institute, 93/1, A.P.C. Road, Calcutta-700009, India

Entanglement is a characteristic feature of quantum mechanical systems which has no classical analogue . The state of a pair (or more than a pair) of quantum systems is entangled if the corresponding wave function does not factorize, i.e., is not a product of the wavefunctions of the individual systems. A well-known example of an entangled state is the singlet state of two spin$`\frac{1}{2}`$ particles, $`\frac{1}{\sqrt{2}}\left(\right||)`$, which cannot be written as a product of the spin states of individual spins. Measurement on one component of an entangled pair fixes the state of the other implying non-local correlations. Interest in quantum entanglement is extensive because of its fundamental role in quantum communication and information processing such as quantum teleportation , superdense coding , quantum cryptographic key distribution  etc. Experimental implementation of some of the protocols has so far been achieved in simple physical settings. Solid state devices, specially, spin systems have been proposed as possible candidates for large scale realizations . In particular, the Heisenberg spin-spin exchange interaction gives rise to entangled states in spin systems and has been shown to provide the basis for universal quantum computation . Examples of other interacting many body systems in which entanglement properties have been studied include the harmonic chain , the 1D Kondo Necklace model  and the BCS condensate .

Entanglement in a state like its energy is quantifiable and has been computed both at $`T=0`$ and at finite $`T`$ (thermal entanglement) for a variety of spin models in both zero and finite external magnetic fields. The models include the Heisenberg $`XX`$, $`XY`$, $`XXX`$, $`XXZ`$ and transverse Ising models in one dimension (1D) \[12-18\]. The computational studies show that the amount of entanglement between two spins in a multi-spin state can be modified by changing the temperature and/or the external magnetic field. Some recent studies have explored the relations between entanglement and quantum phase transitions in the ferromagnetic (FM) $`XY`$ model in a transverse magnetic field and in a special case of the model, namely, the transverse Ising model . A quantum phase transition (QPT) can take place at $`T=0`$ by changing some parameter of the system or an external variable like the magnetic field . The ground state wavefunction undergoes qualitative changes in a QPT and model studies indicate that entanglement develops special features in the vicinity of a quantum critical point.

The studies carried out so far have been confined to 1D systems with only nearest-neighbour (NN) exchange interactions. There are several quasi-1D and 2D antiferromagnetic (AFM) spin models which describe frustrated spin systems . Frustration arises if conflicting interactions are present in the system, i.e., when all the interactions between spins cannot be simultaneously satisfied. A good example of a frustrated spin system is the AFM Ising model on the triangular lattice. An elementary plaquette of the lattice is a triangle. The Ising spin variables have two possible values, $`\pm 1`$, corresponding to the up and down spin orientations. An antiparallel spin pair has the lowest interaction energy $`J`$. A parallel spin pair has the energy $`+J`$. In an elementary triangular plaquette, there are three interacting spin pairs. Due to the topology of the plaquette (odd number of NN bonds), all the three spin pairs cannot be simultaneously antiparallel in the ground state. There is bound to be one parallel spin pair giving rise to frustration in the system. Consider another example in which the spins in a linear chain interact via NN as well as next-nearest-neighbour (NNN) AFM interactions. In a triplet of successive spins, the two NN spin pairs and one NNN spin pair cannot be simultaneously antiparallel and the linear spin chain is frustrated. Frustration in a system may occur due to the topology of the lattice (triangular, kagomé etc.) or due to the inclusion of further-neighbour interactions leading in most cases to spin-disordered ground states. The models exhibit QPTs as the exchange interaction parameters are tuned to particular ratios or at critical values of the external magnetic field. Some of the models exhibit the phenomenon of magnetization plateaus in which plateaus appear in the magnetization/site $`m`$ versus the external magnetic field $`h`$ curve at quantized values of $`m`$ (*m* and $`h`$ are chosen to be dimensionless) . The condition for the appearance of plateaus is

$$S_Um_U=\text{integer}$$
(1)
where $`S_U`$ and $`m_U`$ are the total spin and magnetization in unit period of the ground state. The plateaus indicate that the spin excitation spectrum is gapped so that the magnetization remains constant. In this paper, we consider some spin models for which the ground states are known exactly in certain parameter regimes. An exact measure of the entanglement between two spins can be obtained in these states in a straightforward manner because of the simple structure of the ground states. We report two major results. In the presence of an external magnetic field, macroscopic magnetization jumps at critical values of the magnetic fields (first order QPTs) are accompanied by macroscopic changes in the amount of pairwise entanglement. Furthermore, some examples are given in which the entanglement structure is modified due to QPTs (again first order) brought about by the tuning of exchange interaction strengths.

A measure of entanglement between the spins $`A`$ and $`B`$ is given by a quantity called concurrence . To calculate this, a knowledge of the reduced density matrix $`\rho _{AB}`$ is required. This is obtained from the ground state wavefunction by tracing out all the spin degrees of freedom except those of the spins $`A`$ and $`B`$. Let $`\rho _{AB}`$ be defined as a matrix in the standard basis $`\left\{\right|,|,|,|\}`$. One can define the spin-reversed density matrix as $`\stackrel{~}{\rho }=\left(\sigma _y\sigma _y\right)\rho ^{}\left(\sigma _y\sigma _y\right)`$, where $`\sigma _y`$ is the Pauli matrix. The concurrence $`C`$ is given by $`C=\text{max}\{\lambda _1\lambda _2\lambda _3\lambda _4,0\}`$ where $`\lambda _i^{}s`$ are the square roots of the eigenvalues of the matrix $`\rho \stackrel{~}{\rho }`$ in descending order. The spins $`A`$ and $`B`$ are entangled if $`C`$ is nonzero, $`C=0`$ implies an unentangled state and $`C=1`$ corresponds to maximum entanglement. The models we consider in this paper belong to the Majumdar-Ghosh (MG)  and the Affleck-Kennedy-Lieb-Tasaki (AKLT)  families of models. The MG model is the simplest frustrated model in 1D. The spins have magnitude $`\frac{1}{2}`$ and the Hamiltonian is given by

$$H_{MG}=J\underset{i=1}{\overset{N}{}}\stackrel{}{S}_i.\stackrel{}{S}_{i+1}+\frac{J}{2}\underset{i=1}{\overset{N}{}}\stackrel{}{S}_i.\stackrel{}{S}_{i+2.}$$
(2)
$`N`$ is the total number of spins and the strength of the next-nearest-neighbour (NNN) exchange interaction is half that of the NN interaction. The boundary condition is periodic. The exact ground state is doubly degenerate with the wave functions

$`\varphi _1`$ $`=`$ $`[12][34]\mathrm{}\mathrm{}.[N1N]`$
$`\varphi _2`$ $`=`$ $`[23][45]\mathrm{}\mathrm{}..[N1]`$ (3)
where $`[lm]`$ denotes a singlet (valence bond (VB)) of two spins located at the lattice sites $`l`$ and $`m`$. Thus in the ground state NN spin pairs in singlet spin configurations are maximally entangled with concurrence $`C=1`$. The value of $`C`$ for all the other NN or further-neighbour pairs is zero. This is consistent with the special property of entanglement that in a set of three spins $`A`$, $`B`$ and $`C`$, if $`A`$ and $`B`$ are maximally entangled, the entanglement between $`A`$ and $`C`$ is zero . Translational invariance in the ground states can be restored by taking linear combinations of the ground states, $`\frac{\varphi _1\pm \varphi _2}{\sqrt{2}}`$. In this case, all NN pairs are entangled with $`C=0.5`$. The value of $`C`$ in the case of the $`S=\frac{1}{2}`$ Heisenberg AFM spin chain is $`C=0.386`$ . Frustration appears to increase the NN entanglement in a spin system.

We now give examples of some spin models in the presence of an external magnetic field for which macroscopic magnetization jumps are accompanied by macroscopic jumps in entanglement. The first model is a frustrated two-leg $`S=\frac{1}{2}`$ ladder model  (Fig. 1) the Hamiltonian of which is given by

$$H_{ladder}=\underset{ij}{}J_{ij}\stackrel{}{S}_i.\stackrel{}{S}_jh\underset{i=1}{\overset{N}{}}S_i^z.$$
(4)
The exchange interaction strength $`J_{ij}=J_R`$ along the rungs. The intra-chain NN and the diagonal exchange interactions are of equal strength $`J`$. When $`h=0`$ and $`\frac{J_R}{J}>\lambda `$ $`(\lambda 1.401)`$, the exact ground state consists of singlets (VBs) along the rungs. The spin pairs along the rungs are thus maximally entangled and all other spin pairs are unentangled. We define a quantity

$$f=\frac{2}{N}\underset{i=1}{\overset{\frac{N}{2}}{}}C(i)$$
(5)
which is the average concurrence per rung of the ladder. The summation is over the rung index with $`\frac{N}{2}`$ being the total number of rungs. Since all the other spin pairs are unentangled, the average concurrence per NN bond, $`f_{NN}=\frac{1}{3}`$ (the NN bonds include $`\frac{N}{2}`$ rung bonds and $`N`$ intra-chain NN bonds). The magnetization properties of the frustrated ladder model are simple .

Fig. 2(a) shows the magnetization per rung $`M`$ as a function of the magnetic field $`h`$. For $`0<h<h_{c_1}=J_R`$, $`M=0`$. For $`h_{c_1}<h<h_{c_2}=J_R+2J`$, the rungs are alternately in singlet and $`S^z=1`$ triplet spin configurations in the ground state. The value of $`M`$ is now $`\frac{1}{2}`$. For $`h>h_{c_2}`$, saturation magnetization is obtained with $`M=1`$. One can verify that the quantization condition (1) is obeyed at each plateau. Fig. 2(b) shows the plot of the average concurrence per rung $`f`$ versus the magnetic field $`h`$. A similar plot is obtained for $`f_{NN}`$ with the difference that $`f_{NN}=\frac{1}{3}`$ for $`0<h<h_{c_1}`$ and $`f_{NN}=\frac{1}{6}`$ for $`h_{c_1}<h<h_{c_2}`$. For $`h_{c_1}<h<h_{c_2}`$, half the total number of rungs are in $`S^z=1`$ triplet spin configurations in which the two spins of a rung point up. The parallel spin pairs are unentangled with $`C=0`$. The other rungs are in singlet spin configurations so that $`f=\frac{1}{2}`$. In the fully polarised state $`f=0`$. Fig. 2(b) provides an example of a changing magnetic field giving rise to macroscopic changes in entanglement. The ground state of the frustrated spin ladder model for $`h=0`$ has a simple structure and can be expressed as a product over rung singlet states. The reduced density matrix $`\rho _{AB}`$, needed to calculate the entanglement between the spins $`A`$ and $`B`$, can be calculated exactly and analytically in a straightforward manner. Once the matrix elements of $`\rho _{AB}`$ are known, the concurrence $`C`$ can be calculated. Because of the uncomplicated structure of the ground state, $`C`$ can alternatively be calculated using simple arguments. Each spin in the ladder belongs to a rung singlet or VB so that it is maximally entangled with its partner spin. Thus, the entanglement between a spin *j* and any other spin, belonging to other rung singlets, is zero. The concurrence $`C`$ of a maximally entangled pair is 1 so that $`f,`$ the average concurrence per rung, is 1 for $`0<h<h_{c_1}.`$ Similar arguments hold true for $`h_{c_1}<h<h_{c_2}`$ when the ladder has a different ground state, again with a simple structure.

Fig. 3 shows another spin ladder model with modulated exchange interactions . The model consists of four-spin plaquettes coupled to two-spin rungs (solid lines) through NN and diagonal exchange interactions (dotted lines) of strength $`J_2`$. Within each plaquette, the NN and diagonal exchange interactions are of strength $`J_1`$ and $`J_3`$ respectively. The rung exchange interaction strength is $`J_R`$. In a wide parameter regime, the exact ground state has a simple product form: the plaquettes and the rungs are individually in their ground state spin configurations. The ground state of a plaquette is a resonating valence bond (RVB) state with wavefunction given by $`\psi _{RVB1}(\psi _{RVB2})`$ for $`J_3<J_1(J_3>J_1)`$. The ground state of a rung is a singlet. The RVB states are given by

(6)
The solid lines represent singlets and the arrow signs are drawn according to the phase convention that in a VB between the sites $`i`$ and $`j`$, if the arrow points away from the site $`i`$, then the spin configuration is $`\frac{1}{\sqrt{2}}\left(\right|(i)(j)|(i)(j))`$. It is easy to check that the NN spins in $`|\psi _{RVB1}`$ are entangled with concurrence $`C=0.5`$. The NNN spins along the diagonals are unentangled. On the other hand, the NN spins in $`|\psi _{RVB2}`$ are unentangled and the NNN spin pairs are entangled with $`C=1`$. At $`J_1=J_3`$, there is a QPT from the ground state in which the plaquette spin configurations are described by $`|\psi _{RVB1}`$ to the ground state in which the same are described by $`|\psi _{RVB2}`$. In both the phases, the rungs are in singlet spin configurations. The QPT is accompanied by macroscopic changes in the amount of NN and NNN entanglements. The average concurrence per NN bond, $`f_{NN}`$, in the full ladder is $`\frac{1}{3}`$ for $`J_3<J_1`$ and $`\frac{1}{9}`$ for $`J_3>J_1`$. The average concurrence per NNN bond, $`f_D`$, is $`0`$ for $`J_3<J_1`$ and $`\frac{1}{3}`$ for $`J_3>J_1`$. In a finite magnetic field $`h0`$, the exact ground state maintains its product form in an extended parameter regime. This gives rise to magnetization plateaus in the magnetization/site $`m`$ versus $`h`$ curves. Again, the jumps in the magnetization are accompanied by jumps in the amount of entanglement. To give one specific example, consider the case $`J_3<J_1`$. The average concurrence/NNN bond in the full ladder is zero for $`0<h<h_{c_1}=J_1`$ and $`\frac{1}{6}`$ for $`h_{c_1}<h<h_{c_2}=2J_2`$. At $`h_{c_1}`$, $`m`$ jumps from zero to the value $`\frac{1}{6}`$. Each plaquette is in the ground state spin configuration $`\frac{1}{2}\left(\right|||+|)`$ in the field range $`h_{c_1}<h<h_{c_2}`$.

Schulenberg et al.  have constructed exact eigenstates which consist of independent, localized one-magnon states in a class of frustrated spin lattices and have shown that these are the ground states in high magnetic fields. The eigenstates are obtained provided certain conditions are satisfied. In a kagomé lattice, the magnon states are localised in $`\frac{N}{9}`$ hexagons where $`N`$ is the total number of spins. The hexagons in which the magnon excitations occur are isolated from each other (Fig. 4). Above the saturation magnetic field, all the spins in the lattice are in a FM spin configuration $`|0`$ and each spin pair is unentangled. The state $`|0`$ is the vacuum for magnon excitations. At a critical field below the saturation magnetic field, the exact eigenstate consisting of $`\frac{N}{9}`$ localised one-magnon excitations becomes the ground state. The localised magnon excitation has the wave function

$$|1=\frac{1}{\sqrt{6}}\underset{l=1}{\overset{6}{}}(1)^lS_l^{}|0.$$
(7)
This gives rise to a macroscopic change in the magnetization curve. In each of the $`\frac{N}{9}`$ hexagons, a spin is equally entangled with all the other five spins and the magnitude of the concurrence is $`C=\frac{1}{3}.`$ This signifies a macroscopic change in the amount of entanglement. The magnetization jump occurring at a critical value of the magnetic field signifies a first order QPT. We have shown through specific examples that such QPTs give rise to macroscopic jumps in the amount of entanglement associated with NN and/or further-neighbour spin pairs.

QPTs can also be brought about by tunning exchange interaction strengths. We have already given an example of this in the case of the ladder model shown in Fig. 3. For $`h=0`$, a QPT occurs at $`\frac{J_3}{J_1}=1.`$ In the frustrated two-leg ladder model shown in Fig. 1, a QPT takes place at the critical value $`(\frac{J_R}{J})_c1.401.`$ Below the critical point, the ground state is that of an effective $`S=1`$ chain with the $`S=1`$ spins forming out of the pairs of $`S=\frac{1}{2}`$ rung spins. The spin ladder is now in the Haldane phase of a $`S=1`$ chain. Above the critical point, the ladder is in the rung-singlet (RS) phase in which the rung spins are in singlet spin configurations. Kolezhuk and Mikeska  have proposed a generalised frustrated ladder model in which the first order QPT between the RS and the Haldane phases can be studied in an exact manner. This model includes biquadratic interactions besides NN and NNN (diagonal) ones. The Hamiltonian is given by,

$$H=\underset{i}{}h_{i,i+1},$$

$$h_{i,i+1}=\frac{y_1}{2}(\stackrel{}{S}_{1,i}.\stackrel{}{S}_{2,i}+\stackrel{}{S}_{1,i+1}.\stackrel{}{S}_{2,i+1})+(\stackrel{}{S}_{1,i}.\stackrel{}{S}_{1,i+1}+\stackrel{}{S}_{2,i}.\stackrel{}{S}_{2,i+1})$$

$$+y_2(\stackrel{}{S}_{1,i}.\stackrel{}{S}_{2,i+1}+\stackrel{}{S}_{2,i}.\stackrel{}{S}_{1,i+1})+x_1(\stackrel{}{S}_{1,i}.\stackrel{}{S}_{1,i+1})(\stackrel{}{S}_{2,i}.\stackrel{}{S}_{2,i+1})$$

$$+x_2(\stackrel{}{S}_{1,i}.\stackrel{}{S}_{2,i+1})(\stackrel{}{S}_{2,i}.\stackrel{}{S}_{1,i+1})$$
(8)
$$\text{where},\text{ }x_1=\frac{4}{5}(32y_2),\text{ }x_2=\frac{4}{5}(3y_22).$$
(9)
The indices 1 and 2 correspond to the lower and upper legs of the ladder respectively and $`i`$ is the rung index. The exact phase boundary between the RS and the Haldane phases is given by $`y_1=\frac{4}{5}(1+y_2).`$ The ground state in each phase is known exactly. In the Haldane phase the ground state is the AKLT state . The ground state energy/rung in the RS and the Haldane phases are $`E_{RS}=\frac{3}{4}y_1+\frac{3}{20}(1+y_2)`$ and $`E_{AKLT}=\frac{1}{4}y_1\frac{13}{20}(1+y_2)`$ respectively.

As mentioned before, in the RS phase the spin pairs along the rungs are perfectly entangled with concurrence $`C=1.`$ All the other spin pairs are unentangled. The amount of entanglement in the AKLT phase can be computed in the following manner. For $`[H,S^z]=0,`$ where $`S^z`$ is the $`z`$ component of the total spin, the reduced density matrix of a spin pair located at the sites $`i`$ and $`j`$ has the form

$$\rho _{ij}=\left(\begin{array}{cccc}u_+& 0& 0& 0\\ 0& w_1& z& 0\\ 0& z& w_2& 0\\ 0& 0& 0& u_{}\end{array}\right).$$
(10)
The concurrence quantifying the entanglement is given by

$$C=2\text{max}[0,\left|z\right|\sqrt{u_+u_{}}].$$
(11)
Wang and Zanardi  have shown that the matrix element of $`\rho _{ij}`$ can be expressed in terms of the various correlation functions $`G_{\alpha \beta }=\sigma _{i\alpha }\sigma _{j\beta }=Tr(\sigma _{i\alpha }\sigma _{j\beta }\rho ),`$ $`(\alpha =x,y,z),`$ where $`\rho `$ is the density operator. The magnetization/site $`m=\frac{1}{N}Tr(_{i=1}^N\sigma _i^z\rho ).`$ In particular, the following relations hold true:

$$u_\pm =\frac{1}{4}(1\pm 2m+G_{zz})$$

$$z=\frac{1}{4}(G_{xx}+G_{yy}).$$
(12)
The correlation functions in the AKLT phase can be calculated in an exact manner using the transfer matrix method in the matrix product (MP) formalism . The AKLT ground state can be written in the MP form as

$$\psi _{AKLT}=Tr\underset{i=1}{\overset{\frac{N}{2}}{}}g_i$$
(13)
$$\text{with},\text{ }g_i=\left[\begin{array}{cc}|t_0_i& \sqrt{2}|t_+_i\\ \sqrt{2}|t_{}_i& |t_0_i\end{array}\right].$$
(14)
The product in Eq. (13) is over the $`\frac{N}{2}`$ rungs of the ladder. The two spin-$`\frac{1}{2}`$’s of each rung are in a triplet spin configuration in the AKLT phase giving rise to an effective spin 1. In Eq. (14), the states $`|t_\mu `$ with $`\mu =+1,`$ 0 and $`1`$ represent a spin state with $`S^z=+1,`$ 0 and $`1`$ respectively. Calculation of the correlation functions $`G_{\alpha \alpha }\text{ }(\alpha =x,y,z)`$ of the spin $`\frac{1}{2}`$ pairs in the AKLT state can be carried out following standard procedure . We quote the final results. The various correlation functions are:

$$\sigma _{1,n\alpha }\sigma _{1,n+1\alpha }=\sigma _{2,n\alpha }\sigma _{2,n+1\alpha }=\sigma _{1,n\alpha }\sigma _{2,n+1\alpha }=\sigma _{1,n+1\alpha }\sigma _{2,n\alpha }$$

$$=\frac{4}{9}\text{ }(\alpha =x,y,z)$$
(15)
$$\sigma _{1,n\alpha }\sigma _{2,n\alpha }=\frac{1}{3}\text{ }(\alpha =x,y,z).$$
(16)
The correlation functions, in which the distance $`l`$ separating the rungs on which the spins are located is $`1`$ (in Eq. (15) $`l=1`$), involve the factor $`4(1)^l3^{l1}`$. With the knowledge of the correlation functions, the concurrence can be determined using the relations (11) and (12) ($`m`$ in the AKLT phase is zero). One finds that in the AKLT phase, the spin pairs along the rungs are unentangled. In the RS phase, the same spins are maximally entangled. The intra-leg NN and the NNN (diagonal) spin pairs are entangled in the AKLT phase with concurrence $`C=\frac{1}{6}`$ in each case. In the RS phase, these spin pairs are unentangled. Again, the first order QPT from the RS phase to the AKLT phase is accompanied by macroscopic changes in the entanglement structure.

The generalised ladder model studied by Kolezhuk and Mikeska has a rich phase diagram  describing both first order and second order QPTs. At the second order phase boundary, the gap in the excitation spectrum goes to zero. There are five possible phases: RS, AKLT, FM, D1 and D2. The phases D1 and D2 have spontaneously broken symmetry but the exact ground states are not known in these phases. The phase boundaries separating the FM-D1, AKLT-D1, AKLT-D2, RD-D1 and RD-D2 are known exactly and have been determined in the MP formalism. The corresponding QPTs are second order transitions. The two first order phase boundaries separating the RD-AKLT and AKLT-FM phase are also known exactly. Preliminary calculations  near the second order phase boundaries suggest that the pairwise entanglement does not extend beyond the NNN distance. The same is true in the case of first order QPTs. Osterloh et al.  have considered the range of entanglement in the vicinity of the quantum critical point of the transverse Ising model in 1D. They have found that even at the critical point, where the spin-spin correlations are long ranged, the concurrence is zero for spin pairs separated by more than NNN distance. There is, however, one crucial difference between this model and the AKLT-type models described by MP states. The MP states are finitely correlated, i.e., the spin-spin correlation function $`S_i^zS_{i+l}^z=A_Se^{\frac{l}{\xi }}`$ decays exponentially with the correlation length $`\xi `$ equal to a few lattice spacings. As the transition point $`\tau =\tau _c`$ ($`\tau `$ is some model parameter) is approached, the correlation length $`\xi `$ either does not exhibit any singularity or diverges in a power-law fashion as $`(\tau \tau _c)^2`$. In the latter case, however, the prefactor $`A_s(\tau \tau _c)`$ becomes zero at the transition point . Thus long range spin-spin correlations cannot develop in the system.

In summary, we have shown through specific examples that first order QPTs can bring about macroscopic changes in the amount of pairwise entanglement in spin systems. We have given the specific example of a spin model in which magnetization plateaus give rise to a plateau structure in the average concurrence per rung as well as NN bond. There are several low-dimensional AFM compounds in which magnetization plateaus have been observed experimentally . The appearance of magnetization plateaus has been explained in terms of metal-insulator transitions of magnetic excitations driven by a magnetic field . In the insulating (plateau) phase, the magnetic excitations give rise to crystalline order and in the metallic (non-plateau) phase they are itinerant. It will be of interest to explore the possibility of a plateau structure in the amount of $`T=0`$ as well as finite $`T`$ entanglements in such systems. Our study shows that an external magnetic field can be employed to give rise to large changes in the amount of entanglement and provides the basics for the construction of an entanglement ‘amplifier’ or ‘switch’. There is a large number of spin models which exhibit QPTs of significant interest. Some of these models describe 2D systems. Examples include the Shastry-Sutherland model , the $`S=\frac{1}{2}`$ AFM model on the $`\frac{1}{5}`$-depleted square lattice , a lattice of weakly-coupled two-leg ladders  etc. Some of these models are reviewed in Refs. . All these models exhibit second order QPTs from a spin-disordered gapped phase to a gapless phase with long range spin-spin correlations. Osterloh et al.  and Osborne and Nielsen  have found evidence of entanglement showing scaling behaviour in the vicinity of the quantum critical point. Similar studies should be carried out for the spin models mentioned to augment our knowledge of the relationship between QPTs and entanglement.

Acknowledgement

E. Chattopadhyay is supported by the Council of Scientific and Industrial Research, India under sanction No. 9/15(186)/97-EMR-I.


text:
# Even-odd effects in magnetoresistance of ferromagnetic domain walls

## Abstract

Difference in density of states for the spin’s majority and minority bands in a ferromagnet changes the electrostatic potential along the domains, introducing the discontinuities of the potential at domain boundaries. The value of discontinuity oscillates with number of domains. Discontinuity depends on the positions of domain walls, their motion or collapse of domain walls in applied magnetic field. Large values of magnetoresistance are explained in terms of spin-accumulation. We suggest a new type of domain walls in nanowires of itinerant ferromagnets, in which the magnetization vector changes without rotation. Absence of transverse magnetization components allows considerable spin accumulation assuming the spin relaxation length, $`L_S`$, is large enough.

As it was demonstrated first in Johnson ; Sohn , spin accumulation effects in the presence of a current flowing between ferromagnetic and normal metals may result in a considerable contribution to a contact’s resistivity. This phenomenon is the key element for the GMR devices with the CPP (current perpendicular to the plane) geometry. The CPP-GMR effect has been studied theoretically for the spin-valve systems for a simple triple layer and multi-layered structures in the pioneering work Valet . The subsequent experiments (see e.g. Pratt ; Piraux ; Dubois ) were in excellent agreement with the predictions made in Valet  concerning the dependence of the resistivity on the width of magnetic and non-magnetic components and the role played by spin-relaxation mechanisms. Most recently it was discovered that nanocontacts Garcia  and domain walls in magnetic nanowires Ebels  possess significant magnetoresistance.

Realization of different experimental configurations allows determination of parameters present in the expressions in Ref. \[Valet, \] such as the values of resistivity for each of the GMR components and what is most important the spin relaxation length, $`L_S`$, characterizing the width of a non-equilibrium distribution of spins near the contacts Piraux ; Dubois .

The formulas in Ref. \[Valet, \], however, have been derived in assumption that while the conductivities of the majority and minority spins are different, the corresponding densities of states (DOS), $`g_\alpha `$, remain equal. This assumption is not realistic. In the presentation below we address this issue to demonstrate that taking difference in the DOS into account, the changes in the expressions for the distribution of the electrostatic potential lead to some new observable effects. In Ref. \[Garcia, ; Ebels, \] it was speculated that the pronounced GMR effects are caused by the significant role of spin accumulation. We suggest, as we believe, for the first time, that large magnetoresistance observed in Ref. \[Ebels, \] is due to the non-rotational character of the domain walls Zhirnov  which are possible in itinerant ferromagnets Dzero .

Following Valet , we re-write the expression for the current $`j_\alpha =eD_\alpha n{}_{}{}^{}{}_{\alpha }{}^{}\sigma _\alpha U^{}(x)`$, ($`U(x)`$ \- the electrostatic potential, the index, $`\alpha =\pm `$, stands for the majority and minority spins correspondingly), into the form:

$$𝒋_\alpha =\frac{\sigma _\alpha }{e}\mathbf{}\mu _\alpha ,$$
(1)
where

$$\mu _\alpha =g_\alpha ^1n_\alpha eU(x)$$
(2)
is a non-equilibrium (in the presence of a total current, $`𝑱`$) electrochemical potential for each spin component, and the relation $`D=\sigma /(e^2g)`$ is used in (1, 2).

The equations for the current components:

$$\text{div}𝒋_\alpha =\frac{en_\alpha }{\tau _S}$$
(3)
together with (2) and the electro-neutrality condition:

$$\underset{\alpha }{}n_\alpha =n_++n_{}=0$$
(4)
present the complete system of equations for each side of an interface (in (3) $`\tau _S`$ is a spin relaxation time).

To simplify the analysis, we at first assume the ballistic regime for the interface, i.e. the width of the corresponding domain wall is small (contributions due to the spin scattering inside the barrier will be discussed at the end of the paper). Correspondingly $`\mu _\alpha `$ and $`j_\alpha `$ are taken continuous at the boundary.

In convenient notations of Ref. \[Rashba, \] with the help of (4) we obtain for $`\mu _S=\mu _+\mu _{}`$:

$$\mu _S=\left(g_+^1+g_{}^1\right)n_+.$$
(5)
Subtracting one of Eqs. (3) from another and making use of (5), one obtains the equation for the distribution of $`\mu _S`$ in the “bulk” on each side of the contact:

$$\mu {}_{}{}^{\prime \prime }{}_{S}{}^{}=\frac{\mu _S(x)}{L_S^2},$$
(6)
where

$$L_S^2=\frac{e^2}{\tau _S}\left(\frac{\sigma _+^1+\sigma _{}^1}{g_+^1+g_{}^1}\right).$$
(7)
Equation (6) coincides in its form with Eq. (14) from Ref. \[Valet, \]. For our purposes it is convenient to re-write (6) as:

$$\mu {}_{}{}^{\prime \prime }{}_{S}{}^{}\mu _S/L_S^2=\delta (x)\mathrm{\Delta }(\mu {}_{}{}^{}{}_{S}{}^{}).$$
(6’)
Equation (6) applies now everywhere in the sample. In (6’), $`\delta (x)`$ stands for the Dirac delta-function at the boundary ($`x=0`$) (equation ((6’) is valid when the boundary between the domains is abrupt). While $`\mu _S`$ is continuous, $`\mu _{S}^{}{}_{}{}^{}`$ has a jump at the interface:

$$\mathrm{\Delta }(\mu {}_{}{}^{}{}_{S}{}^{})=\mu {}_{}{}^{}{}_{S}{}^{}(+0)\mu {}_{}{}^{}{}_{S}{}^{}(0).$$
(8)
The jump can be expressed in terms of the total current, $`J`$. To do so, we write down the expressions for $`J`$, spin current $`j_S=j_+j_{}`$, to obtain with the help of (4,5):

$$\begin{array}{c}\hfill J=\frac{1}{e}\left(\frac{\sigma _+g_+^1\sigma _{}g_{}^1}{g_+^1+g_{}^1}\right)\mu _S^{}(x)(\sigma _++\sigma _{})U^{}(x),\\ \hfill j_S=\frac{1}{e}\left(\frac{\sigma _+g_+^1+\sigma _{}g_{}^1}{g_+^1+g_{}^1}\right)\mu _S^{}(x)(\sigma _+\sigma _{})U^{}(x).\end{array}$$
(9)
We specify that *on the left* from the domain wall spins up belong to the majority band. Magnetization changes sign across the wall and so does the band’s occupation so that one has to interchange between $`\sigma _\pm ,g_\pm ^1`$ in (9). Since the currents are continuous, for the jump of electrochemical potential gradient, $`\mathrm{\Delta }(\mu _S^{})`$, we get:

$$\mathrm{\Delta }(\mu _S^{})=e\left(\frac{\sigma _+\sigma _{}}{\sigma _+\sigma _{}}\right)JC.$$
(10)
For a multilayer system, Eq. (6’) reads:

$$\mu {}_{}{}^{\prime \prime }{}_{S}{}^{}\frac{\mu _S}{L_S^2}=C\underset{i=1}{\overset{N}{}}(1)^{i+1}\delta (xx_i),(i=1,\mathrm{},N).$$
(11)
It is also assumed below that the width of the very left and very right banks are larger then $`L_S`$. The solution of (11) is a superposition of solutions for a single domain wall:

$$\mu _S(x)=\frac{CL_S}{2}\underset{j=1}{\overset{N}{}}(1)^je^{xx_j/L_S}.$$
(12)
We consider first the behavior of the electrostatic potential, $`U(x)`$, close to a single domain wall. With the help of Eq. (5) one can write:

$$\mu _++\mu _{}=\left(\frac{g_+g_{}}{g_++g_{}}\right)\mu _S(x)2eU(x).$$
(13)
Making use of continuity of the electrochemical potential at the boundary, one immediately sees from (13) that potential $`U(x)`$ is discontinuous with the jump $`\mathrm{\Delta }U_1(x_i)=U(x=x_i,R)U(x=x_i,L)`$ at the $`i`$th domain wall:

$$\mathrm{\Delta }U_1(x_i)=\frac{(1)^i}{e}\left(\frac{g_+g_{}}{g_++g_{}}\right)\mu _S(x_i).$$
(14)
The spatial distribution of $`U(x)`$ can be found by integrating the first equation (9) along each side of the domains:

$$\begin{array}{cc}& U^{}(x)=\frac{J}{\sigma _++\sigma _{}}+\frac{1}{2e}[\left(\frac{\sigma _+\sigma _{}}{\sigma _++\sigma _{}}\right)\hfill \\ & \left(\frac{g_+g_{}}{g_++g_{}}\right)]\mu ^{}_S.\hfill \end{array}$$
(15)
The total potential drop across an isolated single domain wall in the chosen geometry is the sum of two terms:

$$\begin{array}{c}\hfill \mathrm{\Delta }U=\mathrm{\Delta }U_1+\frac{1}{e}\left(\frac{\sigma _+g_+^1\sigma _{}g_{}^1}{g_+^1+g_{}^1}\right)\times \\ \hfill 2\underset{0}{\overset{\mathrm{}}{}}\mu _S^{^{}}(x)𝑑x=\frac{L_S}{2\sigma _+\sigma _{}}\frac{(\sigma _+\sigma _{})^2}{(\sigma _++\sigma _{})}J\end{array}$$
(16)
(the second term in (16) comes about from the currents distribution (12) on the distances of the order of $`L_S`$ on the two sides of the wall). In the following notations

$$\sigma _\pm =\frac{\sigma }{2}(1\pm \beta ),g_\pm =\frac{g}{2}(1\pm \delta ),$$

one obtains:

$$\mathrm{\Delta }R=\frac{\mathrm{\Delta }U}{J}=2\frac{\beta ^2L_S}{\sigma (1\beta ^2)},$$
(17)
i.e. $`g_\pm `$ drop out from the *total* magnetoresistance. Equations (16-17) coincide with the results obtained in Valet . The differences in the DOS for the minority and majority spins lead to an appearance of the discontinuities (14) and changes in the dependence of the electrostatic potential $`U(x)`$ (see (15)) along the domains.

To make a numerical estimate for the total potential drop we will use the data from Ref. \[Ebels, \] obtained for Co nanowires: $`L_S60`$nm, $`\rho 1.310^5\mathrm{\Omega }cm`$. The typical values of $`\delta `$ are of the order of $`0.4÷0.5`$ Piraux ; Dubois , while $`\beta `$ is of the order of $`0.5÷0.7`$ Janak . After substituting these values in (17), the resistance drop $`\mathrm{\Delta }R`$ per unit area is $`510^{11}\mathrm{\Omega }`$, or $`\mathrm{\Delta }R112\mathrm{\Omega }`$ for the geometry used in Ebels . The value of the potential drop at domain boundary is approximately the same as the value of the total potential drop in (16):

$$\mathrm{\Delta }U_1/\mathrm{\Delta }U=\delta /\beta .$$
(17’)
We would like to emphasize, that the ratio in (17’) may also be negative.

As an example, let us consider in some more details, the drop of the potential, $`\mathrm{\Delta }U(x_i)`$, across the very left domain for a system of $`N`$ walls (we also take $`x_1=0`$). From (12) and (14) we have

$$\mathrm{\Delta }U_N=\frac{1}{2}\frac{(g_+g_{})(\sigma _+\sigma _{})}{\sigma _+\sigma _{}(g_++g_{})}JL_SS_N,$$
(18)
where

$$S_N=\underset{j=1}{\overset{N}{}}(1)^je^{x_j/L_S}.$$
(18’)
For the $`N`$ walls the value of $`\mathrm{\Delta }U_N`$ shows an interesting “even-odd” effect:

$$\begin{array}{c}\hfill \left(\delta U\right)_N\left(\mathrm{\Delta }U_{N+1}\mathrm{\Delta }U_N\right)=\\ \hfill =\frac{(1)^N}{2}\left(\frac{g_+g_{}}{g_++g_{}}\right)\left(\frac{\sigma _+\sigma _{}}{\sigma _+\sigma _{}}\right)Je^{Nd/L_S}.\end{array}$$
(19)
Below in Fig. 1 we plot $`\left(\delta R\right)_N/\mathrm{\Delta }R=\left(\delta U\right)_N/J\mathrm{\Delta }R`$ as a function of number of the domain walls. The “even-odd” effect is well-pronounced at $`NdL_S`$, where $`d`$ is a size of domain.

Another interesting feature is that changes in Eq. (18) could follow a motion of a single domain wall, say at $`x=x_i`$, through its contribution $`(1)^ie^{x_i/L_S}`$ into (18’), or even a collapse of a domain caused by the applied magnetic field (such a collapse has been experimentally seen in Ref. \[Ebels, \]). In order to explicitly demonstrate this effect, we first introduce the following notations:

$$\frac{\mathrm{\Delta }R_{coll}}{\mathrm{\Delta }R}=\underset{j=1}{\overset{i}{}}(1)^je^{x_j/L_S}+\underset{j=i+1}{\overset{N}{}}(1)^je^{x_j/L_S}.$$

The result of our calculations of $`\mathrm{\Delta }R_{coll}/\mathrm{\Delta }R`$ as a function of the distance between two neighboring domains $`L_c=x_{i+1}x_i`$ is plotted in Fig. 2. Experimentally, the motion of a domain wall can be detected by the STM technique.

Now we would like to briefly discuss the change in our results in case when one takes the finite width of the domain wall into account. In particular, we consider the importance of depolarizing effects in the Bloch or Neél type of domain wall for the results above. Since the width of the wall is usually much smaller then the spin diffusion length, $`L_S`$, the arguments that lead us to results given by (14-16) still hold. The only modification would come from the change in the boundary conditions for the spin current.

Electrons going through the Bloch or Neél domain wall loose part of its polarization because the transverse component of magnetization inside the wall creates a torque which causes spin’s re-orientation. Obviously this process will reduce the spin accumulation. If $`d_w`$ is the width of a domain wall, spins of electrons traveling through the wall are rotated in the transverse component of the exchange field, $`H_{exch}`$, in the wall by an angle $`\vartheta _\alpha (\alpha =)`$:

$$\vartheta _\alpha =(\pm )\frac{\mu _BH_{exch}d_w}{\mathrm{}v_\alpha }.$$
(20)
From (20) one obtains:

$$j_S(x=+0)j_S(x=0)=2en_+v$$
(21)
with $`v`$ given by:

$$v=\frac{2\mu _BH_{exch}d_w}{\mathrm{}}$$
(21’)
and after some algebra with help of (5):

$$j_S(x=+0)j_S(x=0)=2ev\frac{g_+g_{}}{g_++g_{}}\mu _S(0).$$
(22)
Taking into account (22), the expression for the jump of $`\mu _S^{}`$ at the interface is:

$$\begin{array}{cc}& \mathrm{\Delta }(\mu _S^{})=e\left(\frac{\sigma _+\sigma _{}}{\sigma _+\sigma _{}}\right)J\hfill \\ & +\frac{e}{2}\left[\frac{2evg_+g_{}(\sigma _++\sigma _{})}{(\sigma _+\sigma _{})(g_++g_{})}\right]\mu _S(0).\hfill \end{array}$$
(23)
As a result, the expression for the potential drop $`\mathrm{\Delta }U`$ across an isolated domain wall acquires the form:

$$\begin{array}{c}\hfill \mathrm{\Delta }U=\frac{\beta ^2}{(1\beta ^2)}\frac{2L_S\sigma ^1}{1+\alpha }J,\\ \hfill \alpha =\frac{e^2vgL_S\sigma ^1}{2}\frac{1\delta ^2}{1\beta ^2}.\end{array}$$
(24)
To estimate the value of $`\alpha `$, we first re-write the second expression in (24) as:

$$\alpha 3\left(\frac{v}{v_F}\right)\frac{L_S}{l},$$
(25)
where $`l`$ is a mean free path and $`v_F`$ is the Fermi velocity. Using the data, provided in Ebels , we have: $`\mu _BH_{exch}0.1`$eV, $`d_w10`$nm, $`L_S/l10`$, $`v4.510^7`$cm/s, $`v_F10^8`$cm/s. Thus our estimates provide $`\alpha 14`$, which according to (23) significantly reduces the spin accumulation effect. This result can also be explained slightly differently: from (20) we estimate the depolarization angle $`\vartheta 1`$, i.e. at crossing the Bloch or Neél type of domain wall the electronic spins would adiabatically follow magnetization. To reconcile the above estimates with the significant magnetoresistance experimentally observed in Ref. \[Piraux, ; Ebels, \], we suggest that the domain walls in these experiments were neither of the Bloch nor Neél type. Instead, in itinerant ferromagnets another type of domain walls, “linear walls”, is realized Dzero . In a linear wall the direction of magnetization does not change, while its absolute value goes through zero inside the wall. Theoretically, “linear” domain walls were considered in Zhirnov  for local spin ferromagnetic systems at temperatures $`T`$ close to critical temperature $`T_{Curie}`$.

To summarize, we have shown that taking into account difference in the density of states between the minority and majority spin bands drastically changes the distribution of the electrostatic potential along the domains. The discontinuities of the potential across each domain wall are of particular interest. The jumps in the values of potential can be measured directly by the STM technique. Jumps possess such characteristic features as the “even-odd” effects counting the total number of domain walls in magnetic nanowire in the presence of an external magnetic field. Our results can be directly extended to the GMR structures consisting of F/N layers which were studied theoretically in Valet  in the approximation of equal density of states of majority and minority bands. To ascribe large values of magnetoresistance observed in Piraux ; Ebels  to spin accumulation effects, it was also necessary to suggest that in nanowires made of itinerant ferromagnets the domain walls are of a linear type in which magnetization changes without creating a perpendicular component that would revert spins of polarized electrons.

The work of M.D. was supported by DARPA through the Naval Research Laboratory Grant No. N00173-00-1-6005, L.P.G. by the NHMFL through the NSF cooperative agreement DMR-9527035 and the State of Florida, A.K.Z. through the NISVP Supplement to the core NHMFL Grant No. 5024-599-22 and the State of Florida, and K. A. Z. by the RFBR grant No. 02-02-17389.


text:
# Radio Emission from SN 1988Z and Very Massive Star Evolution

## 1 Introduction

Supernova (SN) 1988Z was independently discovered in MCG +03-28-022 (Zw 095-049) near $`\mathrm{m}_\mathrm{B}16.4`$ by both G. Candeo on 1988 December 12 (Cappellaro, Turatto, & Candeo 1988) and C. Pollas on 1988 December 14 (Pollas 1988). The SN was classified as a Type II, based on the detection of hydrogen emission lines in the optical spectra by Cowley & Hartwich (Heathcote et al. 1988), which also indicated that the SN was quite distant, with redshift $`z0.022`$. Filippenko (1989) confirmed the Type II identification, with possible resemblance to SN 1987F. The spectra revealed at early times that SN 1988Z was peculiar, with a remarkably blue color, a lack of absorption lines and P-Cygni profiles, very narrow ($``$100 km s-1 FWHM) \[O III\] circumstellar emission lines, and a steadily growing, relatively narrow ($``$2000 km s-1 FWHM) component to the H I and He I lines (Stathakis & Sadler 1991; Filippenko 1991a,b). Stathakis & Sadler (1991), and Turatto et al. (1993) for later times, analyzed in detail the spectroscopic and photometric observations of SN 1988Z. Schlegel (1990) proposed that SN 1988Z, along with other similar SNe, constituted a new class of SNe, the Type II-“narrow,” or SNe IIn.

The properties of the optical spectra and light curves indicated strong SN shock-circumstellar shell interaction (Filippenko 1991a,b; Turatto et al. 1993) which, together with its resemblance optically to SN 1986J (Filippenko 1991a,b; see also Rupen et al. 1987, Leibundgut et al. 1991), strongly suggested that SN 1988Z should be a very luminous radio emitter, as SN 1986J is a strong radio source (see Weiler, Panagia, & Sramek 1990). Sramek et al. (1990) reported the radio detection at 6 cm wavelength of SN 1988Z with the Very Large Array (VLA)111The VLA telescope of the National Radio Astronomy Observatory is operated by Associated Universities, Inc. under a cooperative agreement with the National Science Foundation. on 1989 December 21. The supernova is located at RA(J2000) = $`10^h51^m50\stackrel{\mathrm{s}}{\mathrm{.}}138`$, Dec(J2000) = $`+16\mathrm{°}00\mathrm{}00\stackrel{}{\mathrm{.}}16`$, with an uncertainty of $`\pm `$ 0$`\stackrel{}{\mathrm{.}}`$2 in each coordinate, which is coincident, to within the uncertainties, with its optical position (Kirshner, Leibundgut, & Smith 1989). This initial announcement of the detection of radio emission was followed by a more thorough study and analysis of the first three years of multifrequency radio observations by Van Dyk et al. (1993a). They compared SN 1988Z in detail with SN 1986J and concluded that the two SNe are very similar in their radio properties. They suggested that the progenitor to SN 1988Z was a massive \[$`20M(M_{})30`$\] star which underwent a high mass-loss phase ($`\dot{M}10^4M_{}`$ yr-1) before explosion (see also Stathakis & Sadler 1991).

Luminous radio emission has also been an indicator of X-ray emission from SNe, with synchrotron radio emission being produced as the forward SN shock interacts with the CSM and X-rays being emitted from the corresponding reverse shock region interacting with the SN ejecta (Chevalier & Fransson 1994). For instance, SN 1986J was detected as a luminous X-ray source by Bregman & Pildis (1992) and Houck et al. (1998). Fabian & Terlevich (1996) reported the detection of X-rays from SN 1988Z with ROSAT. SN 1988Z is a very luminous X-ray emitter with a bolometric X-ray luminosity of $`10^{41}\mathrm{erg}\mathrm{s}^1`$.

Chugai & Danziger (1994) offered two models to explain the unusual characteristics of SN 1988Z: 1) shock interaction with a two component wind consisting of a tenuous, homogeneous medium with embedded higher density clumps, or 2) shock interaction with a similar tenuous, homogeneous medium and a higher-density, equatorial, wind-established, disk-like component, favoring the former over the latter. However, they unexpectedly conclude that the SN ejecta is of low mass ($`M<1M_{}`$), and that SN 1988Z may have originated from a relatively low-mass star of $`M_{\mathrm{ZAMS}}8`$–10 $`M_{}`$, in sharp contrast to the high-mass progenitor suggested by Van Dyk et al. (1993a) and Stathakis & Sadler (1991).

Aretxaga et al. (1999) collect the observations from X-ray to radio and attempt to estimate the integrated electromagnetic energy radiated by SN 1988Z in the first 8.5 years after discovery. They obtain a value of $`2\times 10^{51}`$ erg, perhaps as high as $`10^{52}`$ erg, which they consider is sufficiently high to suggest that SN 1988Z could be classified as a “hypernova,” approaching the 2–$`5\times 10^{52}`$ ergs estimated to have been released in SN 1998bw, the possible counterpart of GRB 980425 (Iwamoto et al. 1998), and perhaps indicative of the collapse of the stellar progenitor core into a black hole. (It is interesting to note here that two SNe IIn, 1997cy and 1999E, may also have been associated with $`\gamma `$-ray bursts; see Pastorello et al. 2002 and references therein). Aretxaga et al. (1999) suggest an ejecta mass of $`15M_{}`$, and, therefore, a very high-mass progenitor.

Obviously, SN 1988Z is an extremely interesting and, in many ways, unusual object. Fortunately, due to its intrinsic brightness it has been relatively well-studied in many wavelength bands. Here we report radio observations of SN 1988Z at multiple wavelengths, including new observations which add another six years of monitoring, more than doubling the coverage reported by Van Dyk et al. (1993a). We further interpret its radio emission using a clumpy wind model. We conclude that SN 1988Z, like SNe 1986J and, possibly, other SNe IIn, as well as the unusual SN 1998bw, arise from the explosions of very massive stars surrounded by highly filamentary CSM.

## 2 Radio Observations

New radio observations of SN 1988Z have been made with the VLA at 20 cm (1.425 GHz), 6 cm (4.860 GHz), 3.6 cm (8.440 GHz) and 2.0 cm (14.940 GHz) from 1993 May 4 through 2001 January 25 and are presented in Table 1 and Figure 1. The previously published results from Van Dyk et al. (1993a), along with a few additional and reanalyzed flux density values, are also included in Table 1 and Figure 1 for ease of reference.

Note that the estimated explosion date of 1988 January 23 obtained by Van Dyk et al. (1993a) from their model fit was very poorly constrained by their fitting procedure. Here, we adopt the explosion date of 1988 December 1 estimated from optical data (Stathakis & Sadler 1991). We have correspondingly modified the SN age at each of the epochs in Table 1 and Figure 1 relative to those reported in Van Dyk et al. (1993a). This change in assumed explosion date, although large, is not critical for the model fitting and does not significantly affect the quality of the fit or the conclusions made by Van Dyk et al. (1993a).

The techniques of observation, editing, calibration, and error estimation are described in previous publications on the radio emission from SNe (see, e.g., Weiler et al. 1986, 1990, 1991). The “primary” calibrator was 3C 286, which is assumed to be constant in time with flux densities of 14.45, 7.42, 5.20, and 3.45 Jy at 20, 6, 3.6, and 2 cm, respectively. The “secondary” calibrators222Secondary calibrators are chosen to be compact and unresolved by the longest VLA baselines. While compact and serving as good phase references, such objects are usually variable, so that their flux density must be recalibrated from the primary calibrators for each observing session. were J1051+213, used through 1991 January 17, and J1125+261, for all epochs after that, with defined positions of RA(J2000) = $`10^\mathrm{h}51^\mathrm{m}48\stackrel{\mathrm{s}}{\mathrm{.}}789`$, Dec(J2000) = +$`21^{}19\mathrm{}52\stackrel{}{\mathrm{.}}31`$ and RA(J2000) = $`11^\mathrm{h}25^\mathrm{m}53\stackrel{\mathrm{s}}{\mathrm{.}}712`$, Dec(J2000) = +$`26^{}10\mathrm{}19\stackrel{}{\mathrm{.}}98`$, respectively. After flux density calibration by 3C 286, they served as the actual gain and phase calibrators for SN 1988Z. As expected for secondary calibrators, the flux densities of J1051+213 and J1125+261 have been varying over the years, as can be seen in Table 2 and Figure 2.

The flux density measurement errors for SN 1988Z are a combination of the rms map error, which measures the contribution of small unresolved fluctuations in the background emission and random map fluctuations due to receiver noise, and a basic fractional error $`ϵ`$, included to account for the normal inaccuracy of VLA flux density calibration (see, e.g., Weiler et al. 1986) and possible deviations of the primary calibrator from an absolute flux density scale. The final errors ($`\sigma _f`$) given for the measurements of SN 1988Z are taken as

$$\sigma _f^2=(ϵS_0)^2+\sigma _0^2$$
(1)
where $`S_0`$ is the measured flux density, $`\sigma _0`$ is the map rms for each observation, and $`ϵ=0.05`$ for 20, 6, and 3.6 cm, and $`ϵ=0.075`$ for 2 cm.

## 3 Parameterized Radio Light Curves

Following Weiler et al. (1986, 1990) and Montes, Weiler & Panagia (1997), we adopt a parameterized model (N.B.: The notation is extended and rationalized here from previous publications. However, the “old” notation of $`\tau `$, $`\tau ^{}`$, and $`\tau ^{\prime \prime }`$, which has been used previously, is noted, where appropriate, for continuity.):

$`S(\text{mJy})=K_1\left({\displaystyle \frac{\nu }{\text{5 GHz}}}\right)^\alpha \left({\displaystyle \frac{tt_0}{\text{1 day}}}\right)^\beta e^{\tau _{\mathrm{external}}}`$
$`\times \left({\displaystyle \frac{1e^{\tau _{\mathrm{CSM}_{\mathrm{clumps}}}}}{\tau _{\mathrm{CSM}_{\mathrm{clumps}}}}}\right)\left({\displaystyle \frac{1e^{\tau _{\mathrm{internal}}}}{\tau _{\mathrm{internal}}}}\right),`$ (2)
with

$$\tau _{\mathrm{external}}=\tau _{\mathrm{CSM}_{\mathrm{homogeneous}}}+\tau _{\mathrm{distant}}=\tau +\tau ^{\prime \prime },$$
(3)
where

$$\tau _{\mathrm{CSM}_{\mathrm{homogeneous}}}=\tau =K_2\left(\frac{\nu }{\text{5 GHz}}\right)^{2.1}\left(\frac{tt_0}{\text{1 day}}\right)^\delta ,$$
(4)
$$\tau _{\mathrm{distant}}=\tau ^{\prime \prime }=K_4\left(\frac{\nu }{\text{5 GHz}}\right)^{2.1},$$
(5)
and

$$\tau _{\mathrm{CSM}_{\mathrm{clumps}}}=\tau ^{}=K_3\left(\frac{\nu }{\text{5 GHz}}\right)^{2.1}\left(\frac{tt_0}{\text{1 day}}\right)^\delta ^{},$$
(6)
with $`K_1`$, $`K_2`$, $`K_3`$, and $`K_4`$ corresponding, formally, to the flux density ($`K_1`$), homogeneous ($`K_2`$, $`K_4`$), and clumpy or filamentary ($`K_3`$) absorption at 5 GHz one day after the explosion date, $`t_0`$. The terms $`\tau _{\mathrm{CSM}_{\mathrm{homogeneous}}}`$ and $`\tau _{\mathrm{CSM}_{\mathrm{clumps}}}`$ describe the attenuation of local, homogeneous CSM and clumpy CSM that are near enough to the SN progenitor that they are altered by the rapidly expanding SN blastwave. $`\tau _{\mathrm{CSM}_{\mathrm{homogeneous}}}`$ is produced by an ionized medium that homogeneously covers the emitting source (“homogeneous external absorption”), and the $`(1e^{\tau _{\mathrm{CSM}_{\mathrm{clumps}}}})\tau _{\mathrm{CSM}_{\mathrm{clumps}}}^1`$ term describes the attenuation produced by an inhomogeneous medium (“clumpy absorption”; see Natta & Panagia 1984 for a more detailed discussion of attenuation in inhomogeneous media). Both terms have a radial dependence which, under a constant mass-loss rate, constant velocity wind assumption, is r-2 (see, e.g., Van Dyk et al. 1994 for an example where, for SN 1993J, the radial dependence of the CSM density is found to be flatter than r-2). The values of $`\delta `$ and $`\delta ^{}`$ determine the actual radial density profile if a constant shock velocity is assumed. The $`\tau _{\mathrm{distant}}`$ term describes the attenuation produced by a homogeneous medium which completely covers the source, but is so far from the SN progenitor that it is not affected by the expanding SN blastwave and is constant in time. All absorbing media are assumed to be purely thermal, singly ionized gas, which absorbs via free-free (f-f) transitions with frequency dependence $`\nu ^{2.1}`$ in the radio. The parameters $`\delta `$ and $`\delta ^{}`$ describe the time dependence of the optical depths for the local homogeneous and clumpy or filamentary media, respectively.

The f-f optical depth outside the emitting region is proportional to the integral of the square of the CSM density over the radius. Since in the simple model (Chevalier 1982a,b) the CSM density decreases as $`r^2`$, the external optical depth will be proportional to $`r^3`$, and since the radius increases as a power of time, $`rt^m`$ (with $`m1`$; i.e., $`m=1`$ for undecelerated blastwave expansion), it follows that the deceleration parameter, $`m`$, is

$$m=\delta /3.$$
(7)
The Chevalier model relates $`\beta `$ and $`\delta `$ to the energy spectrum of the relativistic particles $`\gamma `$ ($`\gamma =2\alpha 1`$) by $`\delta =\alpha \beta 3`$, so that, for cases where $`K_2=0`$ and $`\delta `$ is, therefore, indeterminate, we can use

$$m=(\alpha \beta 3)/3.$$
(8)
Since it is physically realistic and may be needed in some RSNe where radio observations have been obtained at early times and high frequencies, we have also included in equation (3) the possibility for an internal absorption term333Note that, for simplicity, we use an internal absorber attenuation of the form $`\left(\frac{1e^{\tau _{\mathrm{CSM}_{\mathrm{internal}}}}}{\tau _{\mathrm{CSM}_{\mathrm{internal}}}}\right)`$, which is appropriate for a plane-parallel geometry, instead of the more complicated expression (e.g., Osterbrock 1974) valid for the spherical case. The assumption does not affect the quality of our analysis because, to within 5% accuracy, the optical depth obtained with the spherical case formula is simply three-fourths of that obtained with the plane-parallel slab formula.. This internal absorption ($`\tau _{\mathrm{internal}}`$) term may consist of two parts – synchrotron self-absorption (SSA; $`\tau _{\mathrm{internal}_{\mathrm{SSA}}}`$), and mixed, thermal f-f absorption/non-thermal emission ($`\tau _{\mathrm{internal}_{\mathrm{ff}}}`$), so that

$$\tau _{\mathrm{internal}}=\tau _{\mathrm{internal}_{\mathrm{SSA}}}+\tau _{\mathrm{internal}_{\mathrm{ff}}},$$
(9)
where

$$\tau _{\mathrm{internal}_{\mathrm{SSA}}}=K_5\left(\frac{\nu }{\text{5 GHz}}\right)^{\alpha 2.5}\left(\frac{tt_0}{\text{1 day}}\right)^{\delta ^{\prime \prime }}$$
(10)
and

$$\tau _{\mathrm{internal}_{\mathrm{ff}}}=K_6\left(\frac{\nu }{\text{5 GHz}}\right)^{2.1}\left(\frac{tt_0}{\text{1 day}}\right)^{\delta ^{\prime \prime \prime }},$$
(11)
with $`K_5`$ corresponding, formally, to the internal, non-thermal ($`\nu ^{\alpha 2.5}`$) SSA and $`K_6`$ corresponding, formally, to the internal thermal ($`\nu ^{2.1}`$) f-f absorption mixed with nonthermal emission at 5 GHz one day after the explosion date, $`t_0`$. The parameters $`\delta ^{\prime \prime }`$ and $`\delta ^{\prime \prime \prime }`$ describe the time dependence of the optical depths for the SSA and f-f internal absorption components, respectively.

A cartoon of the expected structure of an SN and its surrounding media is presented by Weiler, Panagia, & Montes (2001; their Figure 1). The radio emission is expected to arise near the blastwave (Chevalier & Fransson 1994).

The success of the basic parameterization and modeling has been shown by the good correspondence between the model fits and the data for all subtypes of RSNe, e.g., the Type Ib SN 1983N (Sramek, Panagia, & Weiler 1984), the Type Ic SN 1990B (Van Dyk et al. 1993b), the Type II SNe 1979C (Weiler et al. 1991, 1992a; Montes et al. 2000) and 1980K (Weiler et al. 1992b, Montes et al. 1998), and SN 1988Z (Van Dyk et al. 1993a). (Note that, after day $`4000`$, the evolution of the radio emission from both SNe 1979C and 1980K deviates from the expected model evolution, and SN 1979C shows a sinusoidal modulation in its flux density prior to day $``$4000.) A more detailed discussion of SN radio observations and of modeling results is given in Weiler et al. (2002).

Thus, the radio emission from SNe appears to be relatively well understood in terms of blastwave interaction with a structured CSM, as described by the Chevalier (1982a,b) model and its modifications by Weiler et al. (1986, 1990) and Montes et al. (1997). For example, the fact that the homogeneous external absorption exponent $`\delta 3`$, or somewhat less, for most RSNe is evidence that the absorbing medium is generally a wind with density $`\rho r^2`$, as expected from a massive stellar progenitor which explodes during the red supergiant (RSG) phase.

Additionally, in their study of the radio emission from SN 1986J, Weiler et al. (1990) found that the simple Chevalier (1982a,b) model could not describe the relatively slow turn-on. They therefore included terms described mathematically by $`\tau _{\mathrm{CSM}_{\mathrm{clumps}}}`$ in equations (3) and (6). This modification greatly improved the quality of the fit and was interpreted by Weiler et al. (1990) to represent the possible presence of filaments or clumps in the CSM. Such a clumpiness in the wind material was again required for modeling the radio data from SN 1993J (Van Dyk et al. 1994) and SN 1988Z (Van Dyk et al. 1993a). Since that time, evidence for filamentation in the envelopes of SNe has also been found from optical and UV observations (e.g., Filippenko, Matheson, & Barth 1994: Spyromilio 1994).

Through use of this modeling a number of physical properties of SNe can be determined from the radio observations. One of these is the mass-loss rate from the SN progenitor prior to explosion. From the Chevalier (1982a,b) model, the turn-on of the radio emission for RSNe provides a measure of the presupernova mass-loss rate to wind velocity ratio ($`\dot{M}/w_{\mathrm{wind}}`$). Weiler et al. (1986) derived this ratio for the case of pure, external absorption by a homogeneous medium. However, we now recognize that several possible origins for absorption exist and generalize equation (16) of Weiler et al. (1986) to

$`{\displaystyle \frac{\dot{M}(M_{}\mathrm{yr}^1)}{(w_{\mathrm{wind}}/10\mathrm{km}\mathrm{s}^1)}}=3.0\times 10^6<\tau _{\mathrm{eff}}^{0.5}>m^{1.5}\times `$
$`\left({\displaystyle \frac{v_\mathrm{i}}{10^4\mathrm{km}\mathrm{s}^1}}\right)^{1.5}\times \left({\displaystyle \frac{t_\mathrm{i}}{45\mathrm{days}}}\right)^{1.5}\left({\displaystyle \frac{t}{t_\mathrm{i}}}\right)^{1.5m}\left({\displaystyle \frac{T}{10^4\mathrm{K}}}\right)^{0.68}.`$ (12)
Since the appearance of optical lines for measuring SN ejecta velocities is often delayed relative to the time of explosion, we normally adopt $`t_\mathrm{i}`$ = 45 days. Because our observations generally have shown that $`0.8m1.0`$, and from equation (3) $`\dot{M}t_\mathrm{i}^{1.5(1m)}`$, the dependence of the calculated mass-loss rate on the date $`t_\mathrm{i}`$ of the initial ejecta velocity measurement is weak, $`\dot{M}t_\mathrm{i}^{<0.3}`$. Thus, we generally adopt the best optical or VLBI velocity measurements available, without worrying about the deviation of their exact measurement epoch from the assumed 45 days after explosion. For SN 1988Z we follow Turatto et al. (1993) and assume $`v_\mathrm{i}=20,000`$ km s-1 from the highest velocity measured for the broad components of the emission lines. One has to keep in mind that the actual shock velocity for SN 1988Z may be even higher, because it is very hard to derive accurate velocities from fitting very broad lines where the extreme wings may be lost in the continuum.

We also normally adopt values of $`T=20,000`$ K, $`w_{\mathrm{wind}}=10`$ km s-1 (which is appropriate for a RSG wind), $`t=(t_{6\mathrm{cm}\mathrm{peak}}t_0)`$ days from our best fits to the radio data, and $`m`$ from equation (7) or (8), as appropriate. The optical depth absorption term, $`<\tau _{\mathrm{eff}}^{0.5}>`$, however, is not as simple as that used by Weiler et al. (1986).

Weiler et al. (2001) were able to identify at least three possible absorption regimes: 1) absorption by a homogeneous external medium, 2) absorption by a clumpy or filamentary external medium with a statistically large number of clumps, and 3) absorption by a clumpy or filamentary medium with a statistically small number of clumps. Each of the three cases requires a different formulation for $`<\tau _{\mathrm{eff}}^{0.5}>`$. From our consideration of the radio light curves for SN 1988Z, we conclude that the Case 3 of Weiler et al. (2001, 2002) is most appropriate (see $`\mathrm{\S }4`$).

Case 3 occurs when the clump number density is small, the probability that the line-of-sight from a given clump intersects another clump is low, and both emission and absorption effectively occur within each clump. Such a situation will yield a range of optical depths, from zero for clumps on the far side of the blastwave-CSM interaction region, to a maximum corresponding to the optical depth through a clump for clumps on the near side of the blastwave-CSM interaction region. We expect an attenuation of the form, $`(1e^{\tau _{\mathrm{CSM}_{\mathrm{clumps}}}})\tau _{\mathrm{CSM}_{\mathrm{clumps}}}^1`$, as given in equation (3), but now $`\tau _{\mathrm{CSM}_{\mathrm{clumps}}}`$ represents the optical depth along a clump diameter. Since the clumps occupy only a small fraction of the volume, they have a volume filling factor $`\varphi 1`$. Additionally, the probability that the line-of-sight from a given clump intersects another clump is low, so that a relation between the size of a clump, the number density of clumps, and the radial coordinate can be written as

$$\eta \pi r^2RN<1,$$
(13)
where $`\eta `$ is the volume number density of clumps, $`r`$ is the radius of a clump, $`R`$ is the distance from the center of the SN to the blastwave-CSM interaction region, and $`N`$ is the average number of clumps along the line-of-sight, with $`N`$ appreciably lower than unity by definition. Finally, it is easy to verify that there is a relation between the volume filling factor $`\varphi `$, $`r`$, $`R`$ and $`N`$, of the form

$$\varphi =\frac{4}{3}\frac{r}{R}N.$$
(14)
We can then express the effective optical depth $`<\tau _{\mathrm{eff}}^{0.5}>`$ as

$$<\tau _{\mathrm{eff}}^{0.5}>=\frac{\sqrt{2}}{3}\tau _{\mathrm{CSM}_{\mathrm{clumps}}}^{0.5}\varphi ^{0.5}N^{0.5},$$
(15)
where, for initial estimates, we shall take $`N0.5`$ and a constant ratio $`rR^10.33`$, so that, from equation (14), $`\varphi 0.22`$. Before we can fully estimate the mass-loss rate, $`\dot{M}`$, for the progenitor star of SN 1988Z before explosion, we must first fit our parameterized radio light curve model to the radio dataset for the SN.

## 4 Radio Light Curve Description

Examination of Figure 1 and Table 1 shows that the multi-frequency radio data are best described by two time intervals: “early” data, which extends roughly from explosion through day 1479, and “late” data, which roughly extends from day 2129 through the end of the dataset. Figure 1 shows clear steepening of the light curves sometime between these two measurement epochs, but the actual “break” date is somewhat arbitrary, due to uncertainties in the flux density measurements for this relatively faint source and the likely smoothness of any transition region. We have therefore chosen to describe the flux density evolution separately for these two time intervals, with the period from day $`1500`$ to day $`2000`$ as a transition. With the explosion date assumed to be 1998 December 1 from optical estimates (Stathakis & Sadler 1991), and extensive trial fitting showing no evidence for distant, homogeneous thermal absorption ($`K_4=0`$), for SSA ($`K_5=0`$), for mixed free-free absorption/nonthermal emission ($`K_6=0`$), or for a local homogeneous thermal absorption component ($`K_2=0`$), the model fitting process needs to determine only five parameters ($`K_1`$, $`\alpha `$, $`\beta `$, $`K_3`$, and $`\delta ^{}`$). Applying our fitting procedures separately to the early and late periods, with the data points between day 1500 and day 2000 included in both fits to provide a smooth transition, we find that both the spectral index, $`\alpha `$, and the clumpy absorption parameters, K3 and $`\delta ^{}`$, are the same in the two time intervals, to within the fitting errors. However, the emission decay rate parameter, $`\beta `$, steepens significantly, from $`1.2`$ for the early period to $`2.7`$ for the late period. The results of this model fitting are listed in Table 3. Note that $`K_1`$ is a flux density scaling factor, so the fact that it has very different values for the early and the late periods is not physically significant.

For a purely clumpy CSM ($`K_2=0`$), the sharp steepening of $`\beta `$ around day 1750, while $`K_3`$ and $`\delta ^{}`$ remain unchanged, implies that: 1) the number density of clumps per unit volume, $`\eta `$, starts decreasing more rapidly with radius by approximately $`R^{1.5/m}`$ (i.e., $`\eta _{\mathrm{day}>1750}=\eta _{\mathrm{day}<1750}\left(\frac{R}{R_{\mathrm{day1750}}}\right)^{1.5/m}`$), with the average characteristics of the individual clumps remaining unchanged, and 2) most of the absorption occurs within the emitting clumps themselves. In other words, the spatial distribution is so sparse that the average number of clumps along the line-of-sight is less than one ($`N<1`$). This second condition is consistent with Case 3 from Weiler et al. (2001, 2002) and is the basis for our selection of Case 3 for SN 1988Z. It is perhaps significant that Weiler et al. (2001, 2002) found that Case 3 also applies to the radio light curves for the unusual SN 1998bw.

Considering other wavelengths, we note that the H$`\alpha `$ light curve from Aretxaga et al. (1999; see their Figure 2) deviates significantly from their models after day $``$1000. To investigate this further, we have plotted in Figure 3 the H$`\alpha `$ light curve after day 300, the interval for which we have radio data. Examination of Figure 3 shows that the H$`\alpha `$ light curve is consistent with a slope steepening from $`\beta {}_{1}{}^{}=1.84`$ to $`\beta {}_{2}{}^{}=2.94`$ after day $``$1250. The H$`\alpha `$ emission is proportional to $`\rho ^2`$, while for “typical” SNe the radio emission is proportional to $`\rho ^{1.4}`$ (i.e., if $`\beta `$ is the H$`\alpha `$ light curve decline rate, $`\beta {}_{1}{}^{}=1.84`$ should be the same as $`2/1.4\times \beta _{\mathrm{radio}}=1.42\times \beta _{\mathrm{early}}=1.73`$, where $`\beta _{\mathrm{early}}=1.22`$ from Table 3). Thus, very good agreement is found between the behavior of the decline for both the early optical H$`\alpha `$ emission-line and early radio light curves, to within the uncertainties. Also, the change in the slope ($`\mathrm{\Delta }\beta `$) at both H$`\alpha `$ and in the radio should be the same, because the break is due to a change in the number of clumps, not in the properties of the clumps and, therefore, should be wavelength independent: $`\mathrm{\Delta }\beta _{\mathrm{H}\alpha }=(2.94)(1.84)=1.10`$; $`\mathrm{\Delta }\beta _{\mathrm{radio}}=(2.73)(1.22)=1.51`$, which is also in fair agreement, considering the uncertainties. Thus, it appears that the steepening in the decline for the nonthermal radio light curves is consistent with the steepening of the thermal H$`\alpha `$ light curve.

### 4.1 The Mass-Loss Rate for the SN Progenitor

We can now use equations (3) and (15) to estimate from the radio absorption a mass-loss rate for the SN 1988Z progenitor star. With the assumptions for the blastwave and CSM properties discussed in §3 and the results for the best-fit parameters listed in Table 3, we obtain an estimated presupernova mass-loss rate, $`\dot{M}=1.2\times 10^4M_{}`$ yr-1.

This high mass-loss rate for SN 1988Z is only appropriate for the last $`10,000`$ years before explosion (see §5). At earlier epochs of the progenitor’s evolution the mass-loss rate was considerably lower, as indicated by the much more steeply declining $`\beta `$ with unchanged $`K_3`$ and $`\delta ^{}`$ after day $``$1750.

## 5 Discussion

Figure 1, Table 3, and $`\mathrm{\S }4`$ show that the radio light curves for SN 1988Z can be described by standard RSN models (Weiler et al. 1986, 1990, 2001; Montes et al. 1997), with only one fitting parameter varying with time — the index of the time evolution of the radio emission, $`\beta `$. At an age of $``$1750 days $`\beta `$ clearly steepened from $`1.2`$ to $`2.7`$. This change in radio emission characteristic, $`\beta `$, was not accompanied by a corresponding change in the absorption characteristics, $`K_3`$ and $`\delta ^{}`$. Thus, the change in the radio light curves for SN 1988Z cannot be described as only due to a change in the average CSM density, which so effectively described the changes in radio evolution for SN 1980K (Montes et al. 1998) and SN 1979C (Montes et al. 2000).

For SN 1988Z the observations are best interpreted as a decrease in the number density of clumps in the wind-established CSM after day $`1750`$, but with all clumps having essentially the same average internal characteristics as at earlier times. Since the shock velocity for SN 1988Z is assumed to be $`20,000`$ km s-1 and the pre-supernova wind velocity for an RSG is typically $`10`$ km s-1, an interval of 1,750 days after explosion implies that this change took place in the presupernova stellar wind $`9,600`$ years before explosion. With an estimated mass-loss rate in the interval prior to explosion of $`1.2\times 10^4M_{}`$ yr-1, it follows that during this last 10,000 years, the progenitor star shed $`1.1M_{}`$ in a massive “superwind.” The abrupt steepening of $`\beta `$ around day 1750 implies that the mass-loss rate was appreciably lower before that time, decreasing by about a factor of four by day 4,438. Even such a decreasing mass-loss rate, however, still accounts for an additional $`0.8M_{}`$ lost in the stellar wind by our last measurement on day 4,438 ($`25,000`$ years before explosion, using the same wind and blastwave velocity assumptions). Since a relatively high mass-loss rate must have been maintained over a much longer time, such as the $`10^5`$ years appropriate for a massive star’s evolution, an additional several $`M_{}`$ must have been lost earlier so that the presupernova mass loss for SN 1988Z over the entire RSG progenitor lifetime is expected to have been $`>2M_{}`$, and perhaps much greater.

As discussed in $`\mathrm{\S }1`$, Chugai & Danziger (1994) studied SN 1988Z in detail and developed a model from which they concluded that the ejecta mass is small ($`M_{\mathrm{ej}}<1M_{}`$), and that SN 1988Z, therefore, must have originated from a relatively low-mass star, $`M_{\mathrm{ZAMS}}8`$–10 $`M_{}`$, in sharp contrast with the high-mass (20–40 $`M_{}`$) progenitor proposed by Van Dyk et al. (1993a) and Stathakis & Sadler (1991). However, it should be noted that Chugai & Danziger derived a very high mass-loss rate for the SN 1988Z progenitor, $`7\times 10^4M_{}`$ yr-1, yielding what they call their “wind density parameter” (our $`\dot{M}/w_{\mathrm{wind}}`$) of $`5\times 10^{16}\mathrm{g}\mathrm{cm}^1`$. We believe that our mass-loss rate derived from the direct measurement of f-f absorption in the CSM, $`1.2\times 10^4M_{}`$ yr-1 (wind density parameter, $`\dot{M}/w_{\mathrm{wind}}=7.3\times 10^{15}\mathrm{g}\mathrm{cm}^1`$), is likely more appropriate. Since the critical factors in the modeling by Chugai & Danziger are the product of the ejecta mass and wind density parameter related to the energy of the explosion, for a given explosion energy our 6–7 times smaller wind density parameter implies, for their model, a factor of 6–7 times greater ejecta mass, $`M_{\mathrm{ej}}5`$–10 $`M_{}`$, rather than their estimate of $`M_{\mathrm{ej}}1M_{}`$. Even our value for the ejecta mass may still be an underestimate because, by the Chugai & Danziger model, it represents only the mass that was observed to interact with the pre-supernova wind and may not involve the entire ejecta mass. With an ejecta mass of as much as $`10M_{}`$ and a high radio luminosity, which places SN 1988Z among the brightest radio supernovae, we conclude that the stellar progenitor of SN 1988Z must have been a very massive star, perhaps as high as $`20`$ – 30 $`M_{}`$, as previously suggested by Van Dyk et al. (1993a).

It is worth noting that all of the well-studied cases of ultra-bright RSNe (say, having $``$10 times the radio luminosity of SN 1979C at 6 cm peak), such as SNe 1988Z, 1986J, and 1998bw, provide evidence for a highly-clumped CSM with presupernova wind mass-loss rates up to $`10^4M_{}`$ yr-1 and unusually high shock velocities $`20,000`$ km s-1. From optical and ultraviolet spectroscopy Fransson et al. (2002) also find for the SN IIn 1995N evidence for high shock velocities, a “superwind” phase for the progenitor characterized by a high mass-loss rate prior to explosion, and a clumpy CSM. Such high-rate mass loss, clumping of the CSM, and consequent high blastwave velocities may be the signature of particularly massive star progenitors (see also the discussion in Fransson et al. 2002). It is yet to be determined how the progenitors of SNe IIn, which must retain at least some of their hydrogen envelopes until the end of their lifetimes, and the progenitors of possibly extreme Type Ib/c SNe, such as SN 1998bw, which must have lost their entire hydrogen envelopes before explosion, are related, if at all.

## 6 Conclusions

The radio emission from SN 1988Z followed an evolution well described by standard models up to an age of $`1750`$ days ($`4.8`$ yrs), after which its behavior evolved into a much more rapid decline in radio emission without a corresponding change in radio absorption parameters. This is the first time we have witnessed this in any well-studied RSN, and it implies that highly radio-luminous RSNe, such as SN 1988Z and, possibly, SN 1986J, may follow a different evolutionary path in their presupernova mass-loss than that for more “normal” RSNe, such as SNe 1979C and 1980K. However, it is interesting to note that, if one assumes a blastwave velocity $`2,000`$ times faster than the RSG wind-established CSM ($`20,000`$ km s-1 vs. $`10`$ km s-1), the timescale for such a variation in SN 1988Z of $`10^4`$ yr, is similar to that seen for SNe 1979C and 1980K (see, e.g., for SN 1979C, Montes et al. 1998; for SN 1980K, Montes et al. 2000), and also for SN 1998bw (Weiler et al. 2001).

In addition, further analysis of the results from Chugai & Danziger (1994), which implied an extremely high mass-loss rate of $`7\times 10^4M_{}`$ yr-1 and a relatively low-mass progenitor with $`M_{\mathrm{ZAMS}}8`$–10 $`M_{}`$ for SN 1988Z, indicates that their estimates are probably unrealistic. Using their modeling, our lower mass-loss rate, $`1.2\times 10^4M_{}`$ yr-1, implies a significantly higher ejecta mass and, therefore, a higher ZAMS mass of $`20`$–30 $`M_{}`$ for the progenitor star.

We have noted that the H$`\alpha `$ data from Aretxaga et al. (1999), which are also indicative of the SN shock-CSM interaction (Chevalier & Fransson 1994), also begin to deviate significantly from their model by day $``$1250. Although a detailed comparison is beyond the scope of this paper, both the slope before the break and the magnitude of the break are roughly consistent with the radio results.

From our analysis we propose that SNe, such as SNe 1988Z, 1986J, and 1998bw (possibly the counterpart to GRB 980425), with possible very massive star progenitors ($`M_{\mathrm{ZAMS}}>20M_{}`$) and associated massive winds ($`\dot{M}10^4M_{}`$ yr-1), have very highly-clumped, wind-established CSM and unusually high blastwave velocities ($`>20,000`$ km s-1).

KWW thanks the Office of Naval Research (ONR) for the 6.1 funding supporting this research. This work was primarily done while CLW was at Montgomery Blair High School, Silver Spring, MD, serving as a summer student at NRL. CLW thanks the NRL SEAP and STEP programs for supporting his summer work. Additional information and data on radio supernovae can be found on http://rsd-www.nrl.navy.mil/7214/weiler/ and linked pages.