text:
# 1 Introduction

## 1 Introduction

Now it is well known that competition between magnetostatic and exchange energy in a very small ($``$20 nm) particle leads to a single - domain state. If the radius of a particle is sufficiently large, nonuniform distribution of magnetization has minimal energy. It is a vortex in isotropic magnetic disk. Such particle was referred to as circular nanomagnet . In the case of the individual nanoparticle the vortex has perpendicularly magnetized core . Nowadays the vortex distribution of magnetization in nanomagnets is under detailed investigation . Such systems are considered as perspective for use in RAM (Random Access Memory) device . It is obvious that distribution of magnetization in nanomagnets depends on the interaction between particles. The fundamental type of interaction, which can lead to the long-range ordering  and to collective behavior in the system of particles, is the magnetostatic interaction. On the other hand, the energy and a character of the magnetostatic interaction is determined by the magnetization state of the particles.

In this work we experimentally investigate magnetization curves of regular rectangular lattices of permalloy nanoparticles for different external field orientations. We demonstrate that the magnetization distribution within single particle depends on the magnetization process and external field orientation to the lattice axis. It is a result of the interplay of the interparticle interaction and the single particle state. The particles can be both at single-domain and vortex state at zero field. The appearance of the magnetization vortices leads to the appearance of the remanent perpendicular magnetization while magnetizing perpendicular to the easy axis of the system. The competition between the single-domain and vortex states in the system of two magnetostatically interacted magnetic nanodisks has been numerically studied.

## 2 Experiment

Arrays of magnetic particles were fabricated by the electron-beam lithography. The main feature of our method is the usage of fullerine as a resist for electron lithography. Small sizes of $`C_{60}`$ molecules and the ability of fullerine to modify its physical and chemical properties under exposure of electrons allow to use this material for a high-resolution nanofabrication. The capabilities of $`C_{60}`$ as a negative e-resist have been recently demonstrated by fabrication of 20-30 nm Si pillars . The main steps of the procedure of manufacturing permalloy nanoparticles are thin films deposition, exposing by e-beam, development and two-stage etching. We used a double-layer mask containing the $`C_{60}`$ film as a sensitive layer and Ti film as a transmitting layer. Permalloy and Ti films were prepared by pulse laser evaporation on the substrate at room temperature. Fullerine films were deposited by sublimating of a $`C_{60}`$ powder at temperature of 350 C in a vertical reactor with hot walls and cooling holder for the substrate. Transmission electron microscopy, selected area diffraction and X-ray diffraction of the metal and fullerine films were carried out to check the crystalline structure and the thickness of layers. Maximum sizes of the metal crystallites did not exceed 5 nm. The $`C_{60}`$ films have an amorphous structure. The thickness of the magnetic layers was varied from 25 to 45 nm. The thickness of the masking films was 20 nm for the $`C_{60}`$ layer and 30 nm for the Ti film.

The fullerine was patterned in the JEM-2000EX electron microscope with a scanning electron microscopy (SEM) mode by 200 kV e-beam. we had possibility to change the e-beam diameter from 10 nm and over. Utilizing of the high-energy electron beam decreases the amount of backscattering electrons and the shape of patterns becomes defined better. The doses of electron beam irradiation were 0.05-0.1 C/$`cm^2`$, because it assures the reproducibility and uniformity of patterns sizes. Electron beam irradiation of $`C_{60}`$ films reduces the solvability of fullerine in organic solvents. The most likely reason of the changes of the solvability is electron induced polymerization of $`C_{60}`$ molecules accompanied by partially graphitization . The exposed samples were developed in the toluene during 1 min and then patterns were transferred to the Ti layer by the plasma etching with $`CF_2Cl_2`$ atmosphere. The last step of the fabrication of the magnetic particles was the $`Ar^+`$ ion milling of permalloy films using this double-layered mask. Basically, the resistance of the $`C_{60}`$ films to the ion milling is sufficient to use this single mask with little bit greater thickness, but the double etch steps are necessary to ensure a uniformity and reproducibility of the sizes of the particles. By carefully monitoring the elemental composition of the samples by means of EDS (energy dispersion spectroscopy), qualitative microanalysis and checking up the morphology of particles by SEM we can better detect the final points of the plasma etching and ion milling processes. However, usually we did some overmilling at the last step, to prevent presence of any magnetic substance between the prepared particles. One of the SEM images of the arrays of the ferromagnetic particles is presented at Fig.1. The shape of the ferromagnetic particles is a disc, which height equals the thickness of the initial permalloy film.

The parameters of the investigated samples are summarized in the Table. The symbols a and b denote the rectangular lattice parameters (the distances between the centers of the particles), h is the height of the particles and d is their diameter. The total number of the particles is equal to $`10^5`$.

Table

| N | a(nm) | b(nm) | h(nm) | d(nm) |
| --- | --- | --- | --- | --- |
| 1 | 90 | 180 | 45 | 50 |
| 2 | 120 | 240 | 25 | 80 |

To carry out measurements of the magnetic properties we choose a differential Hall microsensor technique. Recently it was shown that Hall magnetometer is a very powerful tool for investigation of magnetic properties of 2D nanoparticles lattices . In our work we applied a commercial magnetometer based on the Hall response in a semiconductor (InSb) to investigate a collective behavior of the fabricated 2D permalloy nanoparticle arrays. The widths of the current and the voltage probes are 100$`\mu `$m and 50$`\mu `$m respectively and the thickness of the semiconductor layer is 10$`\mu `$m. The differential magnetometer consisted of two Hall crosses with series connection of the voltage probes. The lattice of the particles was produced in the active area of one of the Hall crosses. The difference in the Hall voltage between this sample cross and the closely spaced empty cross is measured using the bridge circuit . If the bridge is properly balanced, the resulting output voltage is proportional to the sample contribution to the perpendicular component of the magnetic field induction. The large Hall response in the combination with the good coupling of the small samples to the device results in the excellent spin sensitivity (ratio signal/noise is approximately 100). The sensor works over a large range of the magnetic field and temperature. We have investigated the magnetic properties of the samples by measuring the perpendicular magnetization as a function of the direction and magnitude of the applied field. We have carried out the investigation for three orientations of the external magnetic field: 1) the field is perpendicular to the sample plane (The polar angle $`\theta =0^{}`$); 2) the field is directed at $`45^{}`$ to the sample plane along the short side of the rectangle cell ($`\theta =45^{}`$, the azimuthal angle $`\varphi =0^{}`$); 3) the field is directed at $`45^{}`$ to the sample plane along the long side of the rectangle cell ($`\theta =45^{}`$, $`\varphi =90^{}`$). The direction $`\varphi =0^{}`$ corresponds to the direction along the chains formed by the particles. The chains are elongated along the short side of the elementary rectangle (Fig. 1A). The experimental results for the first sample for $`T=4.2K`$ are represented on the Figs. 2, 3, 4.

The difference of the magnetization curves indicates the collective behavior of the system. It is a result of the magnetostatic interaction between the particles. The hysteresis if the field is directed at $`\theta =45^{}`$, $`\varphi =0^{}`$ (Fig. 3) is the attribute of the easy axis of the magnetization which is directed along the short side of the rectangle cell. The perpendicular remanent magnetization is absent in this case. The existence of such anisotropy in the dipole system was discussed . The magnetization curves if the field is directed at $`\theta =0^{}`$ or $`\theta =45^{}`$, $`\varphi =90^{}`$ (Figs. 2, 4) are qualitatively similar. They have hysteresis in the weak magnetic field with the remanent magnetization. Remanent magnetization is approximately 5 % of the saturation magnetization. The magnetization curves for the second sample are qualitatively similar to those of the first sample, although the samples have the different shape of the particles (See Table). The particles have the polycrystal structure (this was determined by the X-ray diffraction) and do not have the anisotropy of the form in the plane of the system. In this case the difference of the magnetization curves for the different orientation of the external magnetic field with respect to the sample points to the collective behavior of the particles. The existence of the anisotropy axis in the plane of the rectangular lattice was discussed earlier  and was expected. As for the hysteresis and the remanent magnetization for the sample with the rectangular lattice if the external magnetic field direction is $`\theta =0^{}`$ or $`\theta =45^{}`$, $`\varphi =90^{}`$, their existence seems to be unexpected. The effect can not be explained by the properties of single particle. Really, for noninteracting particles the remanent magnetization does not depend on a direction of the external magnetic field.

The first sample was also investigated at T=77K. The hysteresis if the field was directed at $`\theta =45^{}`$, $`\varphi =0^{}`$ did not be observed. The hysteresis if the field was directed at $`\theta =0^{}`$ or $`\theta =45^{}`$, $`\varphi =90^{}`$ was qualitatively changed (Fig. 5).

We suppose that the observed behavior of the system is connected with the fact that the particles of the examined sizes can be in two states. First one is a single domain state, the second one is a vortex distribution of the magnetization. In the last case the core of the vortex is magnetized perpendicularly . The interplay of the lattice anisotropy and the anisotropy of dipole interaction between particles leads to the following behavior of the system. The particles turn to be in the single domain state if the external field has the component directed along the particle chains. In this case the interaction within chain has the ferromagnetic character and therefore stabilizes the single domain state. On the other hand, if the system is demagnetized within the field perpendicular to the chains the dipole interaction has antiferromagnetic character within the chain and the particles turn to be in the vortex state. As core magnetization has its own coercivity, all cores are ferromagnetically ordered and the system has the remanent perpendicular magnetization. We performed numerical simulation to approve these suggestions.

## 3 Numerical simulation

Let us consider two magnetostatically interacting nanoparticles. This model allows to investigate the magnetization distribution within a particle in the anisotropic system with interaction. We numerically solved the system of stochastic Landau-Lifshitz equations. The effective magnetic field takes in consideration the following components: the applied external magnetic field, the demagnetization field of the system, the field of anisotropy, is the exchange field within a particle and the random field defined by the thermal fluctuations. The explicit Euler method was used to solve such stochastic differential equation. Note that such an approach to solve stochastic LLH equation is well known and widely used in magnetic simulation . Besides we have used the approximation of uniform distribution of the magnetization in the perpendicular direction. Our numerical scheme is represented in detail in our article .

Let us now discuss some results of numerical experiments. First of all, for single magnetic disk we found that the vortex magnetization state becomes ground state only when the height and radius of disk exceeded some critical values. This fact was pointed out earlier in  for the square-shape particles. Fig. 6 demonstrates an example of vortex state in cylindrical particle of diameter 50nm and height 12.5nm at zero external magnetic field. Note that such state has nonzero perpendicular component of magnetic moment. This can lead to appearance of hysteresis loop in magnetization curve in systems of magnetic particles.

If there are two interacting particles, the magnetization curve is changed. Two different cases were considered: first, when the projection of the external magnetic field on disk plane lies along the pair of the particles (it corresponds to $`\theta =45^{}`$, $`\varphi =0^{}`$ in experiment), and second, when this projection lies perpendicularly to pair (it corresponds to $`\theta =45^{}`$, $`\varphi =90^{}`$ in experiment). Diameter of the particle was 50nm, its height was 18nm, the distance between the particle centers was 50nm. The corresponding curves for the perpendicular magnetization are shown in Fig. 7. The hysteresis loop (stars in Fig. 7) is the result of vortex penetrating into the particles at some value of the applied field. In contrast with this situation, the magnetization state remains single-domain at all values of external field when its projection lies along the pair of the particles, and as a result, the curve for perpendicular magnetization has no a hysteresis loop (Fig. 8).

## 4 Discussion

The following mechanism for the appearance of the remanent magnetization implies from the results of experimental and numerical simulation suggest with the external field orientation perpendicular to the particle chain. If the external magnetic field has component parallel to the chains, the magnetostatic interaction between particles decreases the total energy of the system and the remanent (H = 0) state is uniform and the perpendicular component of magnetic moment is absent. In the case when the external magnetic field is perpendicular to the chains, magnetostatic interaction increases the total energy of the system. In order to decrease the energy of the magnetostatic interaction, distribution of the magnetization of the particle takes the vortex form. As discussed earlier, a particle in the vortex state has a nonzero magnitude of the remanent perpendicular magnetization.

Besides there were fundamental changes in the magnetization curve of the sample observed when the sample temperature was raised to 77K (Fig. 5). It points to the fact that thermal fluctuations play a significant role in the system behavior, in spite of $`T_c`$ of bulk permalloy is 885K. It is possible that the absence of the remanent magnetization at 77 K is caused by thermally - induced switching between the vortex states with different polarization . If the coercivity of the vortex core become less than the antiferromagnetic magnetostatic interaction between cores the remanent magnetization of the system will be equal to zero. This hypothesis requires further theoretical investigation.

So it is shown that the interaction in the regular 2D rectangular lattice of nanoparticles plays a significant role. The interplay of anisotropy of magnetostatic interaction and lattice anisotropy determines the magnetization state of the particles. It can be both vortex and single-domain state. Due to possibility of the vortex state existence, the system can have the remanent magnetization directed perpendicularly to the system plane.

## Acknowledgment.

We are grateful to Prof. A.S.Arrott for helpful discussion. The work was supported by the RFBR (00-02-16485) and PSSNS Program grants.


format: text

text:
# Nonlinear ac response of anisotropic composites

## I Introduction

When a suspension of highly polarizable dielectric particles is subjected to an intense electric field, the induced dipole moments cause the particles to form chains along the field direction, resulting in complex anisotropic structures. These chain structures alter the apparent viscosity of the suspension by several orders of magnitude, and are the basis of the electrorheological (ER) effect. Under the action of the ER effect, the arrangement of particles, or the shape of the lattice changes and the local electric field deviates from the Lorentz cavity field . As a result, these structures have anisotropic physical properties, such as in the effective conductivity and permittivity , and in the optical nonlinerity enhancement . The rapid field-induced aggregation and the large anistropy in their properties render this field-induced-structured material potentially important for technological applications .

Moreover, the applied electric field used in most ER experiments is quite high, and important data on nonlinear ER effects induced by a strong electric field have been reported recently . A convenient method of probing the nonlinear characteristics of the field-induced structures is to measure the harmonics of the nonlinear polarization under the application of a sinusoidal (ac) electric field. Thus by measuring the ac responses both parallel and perpendicular to the uniaxial anisotropic axis of the structures, it is possible to perform a real-time monitoring of the field-induced aggregation process.

Recently, the effect of nonlinear characteristics on the interparticle force has been analyzed in an ER suspension of nonlinear particles  and further extended to a nonlinear host medium . The nonlinear characteristics are due to the field-dependence of the dielectric constants of the materials used in ER fluids, which have a constitutive relation $`𝐃=ϵ𝐄+\chi 𝐄^3`$. When a nonlinear composite with nonlinear dielectric particles embedded in a host medium, or with a nonlinear host medium is subjected to a sinusoidal field, the electrical response in the composite will in general be a superposition of many sinusoidal functions . It is natural to investigate the effects of a nonlinear characteristics on the ac response in an ER fluid which can be regarded as a nonlinear composite medium . The strength of the nonlinear polarization should be reflected in the magnitude of the harmonics.

In this work, we will develop a perturbation expansion  method and a self-consistent theory  to calculate the ac response of a nonlinear compsite. The paper is organized as follows. In the next section, we calculate the dipole moment of the polarized spherical particles in an anisotropic composite and extract the harmonic response. In section III, we perform a series expansion of the local field inside the particles from the self-consistent solution, in an attempt to obtain analytic expressions of the higher harmonics of the dipole moment. Numerical results are performed in section IV to validate the analytic results. Discussions on the results will be given.

## II Nonlinear polarization and its higher harmonics

We first examine the effect of a nonlinear characteristics on the induced dipole moment. We concentrate on the case where the suspended particles have a nonlinear dielectric constant, while the host medium has a linear dielectric constant $`ϵ_2`$. The nonlinear characteristics gives rise to a field-dependent dielectric coefficient . In which case, the electric displacement-electric field relation inside the spheres is given by:

$$𝐃_1=ϵ_1𝐄_1+\chi E_1^2𝐄_1=\stackrel{~}{ϵ}_1𝐄_1,$$
(1)
where $`ϵ_1`$ and $`\chi `$ are the linear coefficient and the nonlinear coefficient of the suspended particles respectively.

This constitutes an approximation: the local field inside the particles is assumed to be uniform and the assumption is called the decoupling approximation . It has been shown that such an approximation yields a lower bound for the accurate result for the local field . We further assumed that both $`ϵ`$ and $`\chi `$ are independent of frequency, which is a valid assumption for low-frequency processes in ER fluids. As a result, the induced dipole moment under an applied field $`𝐄=E(t)\widehat{𝐳}`$ is given by:

$$\stackrel{~}{p}=\stackrel{~}{ϵ}_ea^3\stackrel{~}{b}E(t),$$
(2)
where $`\stackrel{~}{b}`$ is the field-dependent dipolar factor and is given by:

$$\stackrel{~}{b}=\frac{\stackrel{~}{ϵ}_1ϵ_2}{\stackrel{~}{ϵ}_1+2ϵ_2}=\frac{ϵ_1+\chi E_1^2ϵ_2}{ϵ_1+\chi E_1^2+2ϵ_2}.$$
(3)
In order to obtain the effective dielectric constant $`\stackrel{~}{ϵ}_e`$, we invoke the Maxwell-Garnett approximation (MGA) for anisotropic composites . For the longitudinal field case when the ac field is applied along the uniaxial anisotropic axis, the MGA takes on the form

$$\frac{\stackrel{~}{ϵ}_eϵ_2}{\beta _L\stackrel{~}{ϵ}_e+(3\beta _L)ϵ_2}=f\frac{\stackrel{~}{ϵ}_1ϵ_2}{\stackrel{~}{ϵ}_1+2ϵ_2},$$
(4)
whereas for an transverse field case when the ac field is applied perpendicular to the uniaxial anisotropic axis, the MGA reads

$$\frac{\stackrel{~}{ϵ}_eϵ_2}{\beta _T\stackrel{~}{ϵ}_e+(3\beta _T)ϵ_2}=f\frac{\stackrel{~}{ϵ}_1ϵ_2}{\stackrel{~}{ϵ}_1+2ϵ_2},$$
(5)
where $`\beta _L`$ and $`\beta _T`$ denote the local field factors parallel and perpendicular to the uniaxial anisotropic axis respectively, and $`f`$ is the volume fraction of particles. These factors are defined as the ratio of the local field in the particles to the Lorentz cavity field . For isotropic composites, $`\beta _L=\beta _T=1`$, while both $`\beta _L`$ and $`\beta _T`$ will deviate from unity for an anisotropic distribution of particles in composites. The $`\beta `$ factors have been evaluated in a tetragonal lattice of dipole moments  and similar quantities were calculated in various field-induced-structured composites ; they satisfy the sum rule:

$$\beta _L+2\beta _T=3.$$

In what follows, we denote $`\beta _L`$ as $`\beta `$ for notation convenience.

On the other hand, the effective nonlinear dielectric constant can be given by

$$\stackrel{~}{ϵ}_e=\frac{1}{E^2(t)V}_Vϵ(𝐫)|𝐄(𝐫,t)|^2dV=\frac{f\stackrel{~}{ϵ}_1}{E^2(t)}E_1^2+\frac{(1f)ϵ_2}{E^2(t)}E_2^2,$$
(6)
where $`V`$ the volume of the composite and $`f`$ is the volume fraction. Accordingly, the local electric field inside the spheres can be expressed in terms of the derivative of $`\stackrel{~}{ϵ}_e`$ with respect to $`\stackrel{~}{ϵ}_1`$:

$$E_1^2=\frac{1}{f}E^2(t)\frac{\stackrel{~}{ϵ}_e}{\stackrel{~}{ϵ}_1}.$$
(7)
If we apply a sinusoidal electric field, i.e. $`E(t)=E_0\mathrm{sin}(\omega t)`$, the induced dipole moment $`\stackrel{~}{p}`$ will depend on time sinusoidally, too. By virtue of the inversion symmetry, $`\stackrel{~}{p}`$ is a superposition of odd-order harmonics such that

$$\stackrel{~}{p}=p_\omega \mathrm{sin}\omega t+p_{3\omega }\mathrm{sin}3\omega t+\mathrm{}.$$
(8)
Also, the local electric field contains similar harmonics

$$\sqrt{E_1^2}=E_\omega \mathrm{sin}\omega t+E_{3\omega }\mathrm{sin}3\omega t+\mathrm{}.$$
(9)
These harmonic coefficients can be extracted from the time dependence of the solution of $`\stackrel{~}{p}`$ and $`E_1(t)`$.

## III Analytic Solutions

In what follows, we will apply two methods to extract the harmonics of the induced dipole moment and the eletric field: the perturbation expansion method (PEM)  and the self-consistent (SC) theory . The self-consistent theory can deal with the case of strong nonlinearity, while the perturbation expansion method is applicable to weak nonlinearity only, i.e. $`\chi E_1^21`$, limited by the convergence of the series expansion. In the case of weakly nonlinear response, the SC theory should agree with PEM.

### A Perturbation expansion method

We expand $`\stackrel{~}{p}`$ and $`\sqrt{\chi E_1^2}`$ into a Taylor expansion

$`\stackrel{~}{p}`$ $`=`$ $`a^3E(t){\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}a_n(\chi E_1^2)^n,`$ (10)
$`\sqrt{\chi E_1^2}`$ $`=`$ $`{\displaystyle \frac{\sqrt{\chi E^2(t)}}{f^{1/2}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}d_m(\chi E_1^2)^m,`$ (11)
where

$$a_n=\frac{1}{n!}\frac{^n}{\stackrel{~}{ϵ}_1^n}(\stackrel{~}{ϵ}_e\stackrel{~}{b})|_{\stackrel{~}{ϵ}_1=ϵ_1},d_m=\frac{1}{m!}\frac{^m}{\stackrel{~}{ϵ}_1^m}\sqrt{\frac{ϵ_e}{\stackrel{~}{ϵ}_1}}|_{\stackrel{~}{ϵ}_1=ϵ_1}.$$

For weak nonlinearity, we can rewrite Eqs.(10) and (11), keeping the lowest orders of $`\chi E^2(t)`$ and $`\chi E_1^2`$:

$`\sqrt{\chi E_1^2}`$ $`=`$ $`{\displaystyle \frac{\sqrt{\chi E^2(t)}}{f^{1/2}}}(d_0+d_1\chi E_1^2),`$ (12)
$`\stackrel{~}{p}`$ $`=`$ $`a^3a_0E(t)+{\displaystyle \frac{1}{f}}a^3a_1d_0^2\chi E^3(t)h_1E(t)+h_3\chi E^3(t).`$ (13)
In the case of a sinusoidal field, we can expand $`E^3(t)`$ in terms of the first and the third harmonics. The comparison with Eq.(8) yields the harmonics of the induced dipole moment:

$$p_\omega =h_1E_0+\frac{3}{4}h_3\chi E_0^3,p_{3\omega }=\frac{1}{4}h_3\chi E_0^3.$$

Similarly, we find the harmonics of the local electric field

$$\chi ^{1/2}E_\omega =j_1\chi ^{1/2}E_0+\frac{3}{4}j_3(\chi ^{1/2}E_0)^3,\chi ^{1/2}E_{3\omega }=\frac{1}{4}j_3(\chi ^{1/2}E_0)^3,$$

where

$$j_1=d_0/f^{1/2},j_3=d_0^2d_1/f^{3/2}.$$

In the above analysis, we have used the identity $`\mathrm{sin}^3\omega t=\frac{3}{4}\mathrm{sin}\omega t\frac{1}{4}\mathrm{sin}3\omega t`$ to obtain the first and the third harmonics. Similar analysis can be used to extract the higher-order harmonics, by retaining more terms in the series expansion.

### B Self-consistent theory

Eq.(7) is actually a self-consistent equation since its right-hand side depends on the local field itself. It must be solved self-consistently, and the local electric field reads:

$$\sqrt{\chi E_1^2}=\frac{2(3)^{1/3}R+2^{1/3}(9Q\sqrt{T}+\sqrt{3}\sqrt{4R^3+27Q^2T})^{2/3}}{6^{2/3}\sqrt{T}(9Q\sqrt{T}+\sqrt{3}\sqrt{4R^3+27Q^2T})^{1/3}}$$
(14)
where

$$Q=3ϵ_2\chi ^{1/2}E(t),R=ϵ_1T+ϵ_2(2+f\beta ),T=1f\beta .$$

For a sinusoidal applied electric field, we determine numerically the harmonics of the induced dipole moment and the local electric field, which is a subject for the next section.

## IV Numerical Results

In this section, we perform numerical calculations to investigate the effects of a nonlinear characteristics on the harmonics of the induced dipole moment and the local electric field. As shown in Section III, the harmonics of the induced dipole moment and the local electric field are both related to the field-induced anisotropy parameter $`\beta `$. In order to investigate the effect of field-induced anisotropy on the harmonics, we perform the numerical calculations. Without loss of generality, we let $`f=0.09`$, $`a=1`$, $`ϵ_1=10`$, and $`ϵ_2=1`$ for model calculations. We let $`p_0=ϵ_ea^3bE_0`$ to normalize the induced dipole moment.

In Fig.1, we plot the normalized harmonics versus $`\beta `$ by using the SC theory for the longitudinal electric field. It is evident that, for increasing $`\beta `$, the harmonics of the induced dipole moment and the local electric field are both increased. Similarly, in Fig.2 the same quantities are plotted but for the transverse electric field. We find that the harmonics are decreased as $`\beta `$ increases.

Both Fig.1 and Fig.2 also show that the harmonics are strongly dependent on the nonlinear response of the suspended particles. Moreover, increasing the nonlinear response leads to an increase in the harmonics. This is consistent with the results of our recent work , in which a pair of nonlinear particles suspended in a linear host was considered.

In Figs.3 and 4, we compare the SC theory and PEM. For weak nonlinearity, $`\chi E_0^2=0.09`$, the PEM results coincide with the SC results, while for moderate nonlinearity $`\chi E_0^2=0.9`$, there are some deviations. In both cases, we find that the PEM results are in good agreement with the SC results for weak nonlinearity. For strong nonlinearity, there are large deviations (not shown here) and the SC theory must be used.

The $`\beta `$ dependence of the local field can qualitatively be understood with the spectral representation . From Eq.(7), denoting $`s=(1ϵ_1/ϵ_2)^1`$, we find that $`E_1=sE_0/(ss_1)`$, where $`s_1=(1\beta f)/3`$. Using the numerical values of our model calculations, we find $`s=1/9`$, and $`E_1=(1/9)E_0/(1/9+s_1)`$. When $`\beta `$ increases, $`s_1`$ decreases, leading to an increase in the local electric field $`E_1`$.

## Discussion and conclusion

Here a few comments on our results are in order. In the present study, we have examined the case of nonlinear particles suspending in a linear host. We may extend our considerations to a nonlinear host medium . In this case, preliminary results show that qualitatively similar yet more complex behaviors in the ac response have been observed.

So far, we have not considered the frequency dependence of the particle dielectric constant. In a realistic situation, the dielectric constant of the particles can decrease with the increase of the frequency. For simplicity, we may adopt the Debye relaxation expression for $`ϵ_1`$. Preliminary results show that the ratio of the third to first harmonic decreases with frequency, results that are in accord with recent experimental data of Ref..

To our knowledge, an accurate evaluation of the local field factor $`\beta `$ during field-induced aggregation is lacking. We suggest that the Ewald method  can be extended to compute the $`\beta `$ factor in computer simulation of ER fluids. This is a formidable task for the future. Our present theory may be applied to the electrorheology of milk chocolate . By measuring the nonlinear ac response, it may be possible to monitor the food production processes.

In conclusion, a self-consistent formalism has been employed to compute the induced dipole moment for suspensions in which the suspended particles have a nonlinear characteristics, in an attempt to investigate the anisotropy in the ac response.

## Acknowledgments

This work was supported by the Research Grants Council of the Hong Kong SAR Government under project number CUHK 4284/00P. K. W. Y. acknowledges useful discussion with Prof. G. Q. Gu.

Huang, Wan, Lo, Yu

Huang, Wan, Lo, Yu

Huang, Wan, Lo, Yu

Huang, Wan, Lo, Yu


text:
# Onset of Wave Drag due to Generation of Capillary-Gravity Waves by a Moving Object as a Critical Phenomenon

## Abstract

The onset of the wave resistance, via generation of capillary gravity waves, of a small object moving with velocity $`V`$, is investigated experimentally. Due to the existence of a minimum phase velocity $`V_c`$ for surface waves, the problem is similar to the generation of rotons in superfluid helium near their minimum. In both cases waves or rotons are produced at $`V>V_c`$ due to Cherenkov radiation. We find that the transition to the wave drag state is continuous: in the vicinity of the bifurcation the wave resistance force is proportional to $`\sqrt{VV_c}`$ for various fluids.

An object, moving uniformly on a free fluid surface, is subjected to a drag force, that can be of different physical nature. The most common one is a viscous drag which at low Reynolds numbers, $`Re1`$, is the Stokes drag, $`F_v=3\pi \eta dV`$, that is proportional to an object velocity, $`V`$, and at higher $`Re`$ is the drag which originates from either laminar or turbulent wakes. (Here $`d`$ is the sphere diameter, and $`\eta `$ is the viscosity.) However, there exists another drag force, which originates from generation of capillary-gravity waves by an uniformly moving object. These waves remove momentum from the object into infinity, that produces the wave resistance force, acting on the object.

The dispersion relation for the surface waves is

$$\omega ^2=gk+\sigma k^3/\rho ,$$
$`(1)`$
where $`\omega `$ is the circular frequency, $`k`$ is the wave number, $`\rho `$ is the fluid density, $`g`$ is the gravity acceleration, and $`\sigma `$ is the surface tension. According to Eq. (1), a phase velocity of the waves $`c=\omega /k`$ has a mimimum $`V_c=(4g\sigma /\rho )^{1/4}`$ at the capillary wave number $`k_c=\sqrt{\rho g/\sigma }`$. Below $`V_c`$ the surface waves cannot be emitted, so no wave resistance force acts on the object. The stationarity condition of the wave pattern in the frame moving with the object, leads to

$$\omega =kV\mathrm{cos}\theta ,$$
$`(2)`$
where $`\theta `$ is the angle between $`𝐕`$ and $`𝐤`$.

Eliminating $`\omega `$ from equations (1) and (2), one finds

$$\mathrm{cos}\theta (k)=c(k)/V.$$
$`(3)`$
This equation has evidently no solutions for $`V<V_c`$ and describes the opening of the Cherenkov radiation cone at $`V>V_c`$. Generation of surface waves is analogous to the Cherenkov emission of electromagnetic waves by a charge uniformly moving with a superluminal velocity in a medium. Due to the linear dispersion relation the Cherenkov radiation of the electromagnetic waves occurs at every wave number into a cone with an opening $`\theta `$ defined by Eq. (3).

For a given $`VV_c`$ the system of equations (1) - (2) has solutions in the range $`0\theta <\chi `$, where the angle $`\chi `$ is defined by $`\mathrm{cos}\chi =V_c/V`$. This range of $`\theta `$ corresponds to the interval $`k_1<k<k_2`$ of wave numbers, resulting from the intersection of Eqs. (1) with a straight line of slope $`V`$.

Therefore, unlike the electromagnetic waves, the capillary-gravity waves are radiated within a finite range $`\mathrm{\Delta }k=k_2k_1`$. It tends to zero as $`\mathrm{\Delta }k=4k_c\sqrt{ϵ}`$, when $`V`$ approaches $`V_c`$ from above, where $`ϵ=(VV_c)/V_c`$. The presence of the minimum on the dispersion curve of the capillary-gravity waves (eq.(1)) results in a gap in the spectrum $`\omega (k)`$, similar to the gap in the energy spectrum of a superfluid. The existence of the minimum phase velocity for the dispersion curve (1) is analogous to the existence of the Landau critical velocity for the phonon-roton energy spectrum in a superfluid helium . The problem of the drag onset at the Landau critical velocity $`V_0`$ due to the radiation of rotons was recently considered theoretically in Ref. . Pomeau and Rica observed the onset of the roton generation at a certain critical velocity, in numerical simulations of the generalized nonlinear Shr$`\ddot{o}`$dinger equation, and pointed out but without proof, that the roton drag force close to the onset depends on the control parameter $`ϵ=(VV_0)/V_0`$ as $`\sqrt{ϵ}`$ in 3D and as $`\mathrm{log}ϵ`$ in 2D.

The onset of the wave resistance and its behavior in the supercritical region for the capillary-gravity surface waves was theoretically considered in recent papers . The main result of the calculations is the prediction of a discontinuous transition to a nonzero wave resistance state at $`V=V_c`$, if the object size is much smaller than the capillary wavelength $`\lambda _c=2\pi /k_c`$.
In this Letter we investigate the dependence of the wave drag force, on the object velocity, $`V`$, in the vicinity of the transition for various fluids and object sizes. The results of the experiments clearly demonstrate that the transion to the wave drag state is continuous one in contradiction with the theoretical prediction.

The experiments were conducted in a circular rotating channel of 12.4 cm outer radius and 3.6 cm width, made of plexiglass, with a stainless steel ball, as an object, suspended on a wire, which is glued to an elastic fiber. The channel was driven by a stepping motor via a belt transmission with ratio 0.138 and with a velocity controlled better than 0.02%. Balls of various diameters, $`d`$, (1.57, 2.35, and 3.14 mm) and wires (0.3 and 0.7 mm dia) half-immersed into a fluid, were used. Each of them is much smaller than $`\lambda _c`$ for any fluid used in the experiments. A drag force as a function of the velocity was measured in two ways. The first method was to measure the drag force, applied to an object, when the object deviates from its equilibrium vertical position with increasing driving velocity in accordance with the drag force. The second method was to measure the same force but with an object, hold in an initial vertical position by using an appropriate feedback loop. The latter was actuated via eddy current commercial displacement gauge (EMD 1050 Electro Corporation), that measures a deviation of brass cylinder(2 mm dia), mounted on the same wire as the object, from its initial position. A signal from the displacement detector is used via an instrumentational amplifier and a power amplifier for the feedback loop to drive electromagnets which hold the object in its initial position (see Fig.1). A current supplied to the electromagnets is calibrated as a force. In another setup a ball deviation from its initial position was measured optically by CCD camera with spatial resolution of about 10 $`\mu `$m. Light refraction visualization was also used to define the onset of the surface waves generation, the angle of the Cherenkov cone opening, and the wavenumbers as a function of velocity. As working fluids water, silicone oils DC200 of the kinematic viscosities $`\nu =\eta /\rho =`$10 and 50 cS and of the same surface tension $`\sigma =`$21.2 mN/m, water-glycerol mixtures of $`\nu =`$18, 30, and 46 cS and almost the same $`\sigma =`$66 mN/m, were used. The depth of each fluid was sufficiently large to ensure the validity of an infinite depth approximation for the surface waves.

The raw data on the full drag force, $`F`$, as a function of the channel velocity for silicone oil(10 cS)-air interface and $`d=`$3.14 mm ball are presented in Fig.2. The data taken with increasing and decreasing values of velocities are reproducible and show no hysteresis. Analogous data for $`d=`$2.35 and $`d=`$1.57 mm balls were also obtained. Since the critical velocity of the transition to the wave generation state is rather low, the critical Reynolds numbers at the transition for different fluids and ball sizes are in the range between 2.2 (for 1.57 mm ball and silicone oil 50 cS) and about 700 (for 3.14 mm ball and water). It means that below the transition the drag consists of a viscous Stokes drag at $`Re1`$, which changes linearly with $`V`$, and the drag, proportional to $`V^2`$, when a boundary layer still remains laminar up to $`Re10^3`$. One fits the data below the transition by a second order polynomial and gets a reasonable value of the Stokes drag (see the lower inset in Fig.2, where the fit gives $`\nu =`$8.4 cS instead of 10 cS). Then by subtraction the viscous drag one obtains the reduced wave drag force, $`(FF_v)/F_c`$, as a function of the velocity. This plot clearly exhibits the continuous bifurcation to the wave drag state at $`V_c=16.5`$ cm/sec, that is rather close to the theoretical value of 17.1 cm/sec. Here $`F_v`$, $`F_c`$ are the viscous drag force and the critical value of the full drag force at the onset, respectively.

To emphasize a role of a contact line relocation due to the ball displacement, two set of data-one with a feedback as in Fig.2, and another without-are shown in Fig.3 for silicone oil of 50 cS. Since the capillary force exceeds the wave drag force near the onset, the relocation of the contact line due to a ball motion and variation of capillary forces in this way diminish partially or completely the result of the wave drag force in the measurements. So, the transition is almost smeared out in the measurements without the feedback. One should mention that in our set-up it was impossible to actuate the feedback with a ball less than 1.57 mm dia. Indeed, the viscous drag reduces with a radius between linear and square dependence but the stabilizing magnetic forces reduces in a third power. Fig.4 shows the data of the reduced drag force vs the reduced velocity for three different fluids with 3.14 mm ball: silicone oil(50 cS), glycerol-water(30 cS), and glycerol-water(46 cS). It is obvious from the plots that the increase in the wave drag force strongly depends on fluid properties ($`\nu `$ and $`\sigma `$). We would like also to point out that the wave drag for water(see the inset in Fig.4) increases even more dramatically. The inset presents the data for water with 3.14 mm ball in the narrower range of the control parameter, $`V/V_c`$, together with a fit by the stationary Ginzburg-Landau equation with a field: $`ϵ\xi a\xi ^3+h=o`$, where $`\xi (FF_v)/F_c`$ is the order parameter, $`ϵV/V_c1`$ is the control parameter, $`a`$ is the nonlinear coefficient, and $`h`$ is the field, which plays a role of the smearing factor. As we found from the fits both $`a`$ and $`h`$ depend on $`\nu `$ and $`\sigma `$ of fluids, and ball sizes. The scaling found experimentally, is $`a\eta \sigma /d^2`$. Indeed, using this scaling all the data for various fluids and ball sizes in the vicinity of the transition collapse onto one curve(see Fig.5 for 5 different sets of the data). Thus, the transition to the wave drag state is continuous, and the drag force behaves above the threshold as $`(FF_v)\sqrt{VV_c}`$.

Fig.6 presents an image of the water surface waves just above the transition that provides information about the opening angle of the Cherenkov cone and the wavelength. The inset in Fig.7 shows the dependence of $`\lambda `$ in water on the velocity together with the fit. The main plot presents the data for $`\mathrm{cos}\theta `$ as a function of the velocity. Due to dispersion the phase velocity of the waves depends on the wavenumber. We used the fit of the wavenumbers vs velocity from the inset in Fig.7 to obtain the corresponding values of $`c(k)`$ in eq.(3). The result of the calculations of $`\mathrm{cos}\theta `$ according to eq.(3), is shown by a solid line in Fig.7. If one takes into account that no fitting parameters were used in this procedure, the agreement with the data is rather surprising. Particularly, the critical velocity, $`V_c`$, obtained from extrapolation of $`\theta 0`$, is sufficiently close to the value, found from the force measurements.

Thus, the experiments presented shows unambiguously that the transition to the wave resistance state is continuous one, contrary to the theoretical prediction. What can be the reason for this discrepancy? The theory is based on the Kelvin model, which is a reasonable approach to long wavelength gravitational waves. Instead of a real ship-like object Kelvin considered a moving pressure point applied along its course on a water surface. In a case of an object much larger than $`\lambda _c`$, one can neglect a pressure redistribution due to capillary effects, and the Kelvin model works extremely well for long gravity waves. In the opposite case, a relocation of a contact line causes redistribution of the pressure. That results in a bump behind the moving object observed in experiments. This factor makes the applicability of the Kelvin model questionable. The role of capillarity and wetting is clearly demonstrated experimentally by switching off the feedback control, as we discussed above. Moreover, by using a thin wire as an object one can increase drastically the role of the capillary (wetting) forces in a force balance with the viscous and wave drag forces. Then strong fluctuations of the contact line position due to hysteretic behavior of wetting front, result in strong scatter of the data in the vicinity of the transition which we actually observed in the experiments with wires of 0.3 and 0.7 mm. This effect leads the authors of Ref. to an erronous conclusion about the type of the transition to the wave drag state.

This experiment was initiated in numerous discussions with a generation of graduate students in the group. One of us (V.S.) is gratefult to M. Assenheimer, H. Davidowitz, A. Groisman, E. Kaplan, and D. Rinberg for valuable and helpful advices, particularly on initial stages of the research, and to M. Assenheimer and J. Groshaus for participation in early versions of the experiment. This work was supported by the Minerva Center for Nonlinear Physics of Complex Systems.


text:
# Rayleigh-Taylor Instabilities in Young Supernova Remnants Undergoing Efficient Particle Acceleration

## 1 INTRODUCTION

The presence of convective or Rayleigh-Taylor (R-T) instabilities in young, supernova remnants (SNRs) has been argued for a long time (e.g., Gull, 1973; Shirkey, 1978; Dickel et al., 1989) and examined in detail with numerical simulations of idealized models (e.g., Chevalier et al., 1992; Jun & Norman, 1996a; Wang & Chevalier, 2000). These instabilities are expected to arise in young, ejecta-dominated SNRs where a shell of dense, shocked ejecta is gradually decelerated by lower-density shocked circumstellar material. This situation is subject to the familiar R-T instability, which results in fingers of dense ejecta gas protruding into, and mixing with the shocked circumstellar material. These instabilities are relevant to remnants from both Type Ia and Type II explosions, and are expected to exist as long as there is a reverse shock present in the remnant.

These instabilities may be important for mixing ejecta out to, and possibly past, the forward shock, and for modifying the morphology of an otherwise spherical outer blast wave and/or reverse shock. One of the best studied young SNRs, Cas A exhibits ejecta material out to and beyond the nominal radius of the forward shock, as evidenced by optical knots of enriched gas (Fesen & Gunderson, 1996), as well as spectral analysis of X-ray knots and filaments (Hughes et al., 2000a). The outer blastwave of Cas A, as seen in X-ray emission, is roughly spherical, but possesses many bumps and protrusions that may be the result of the convective instability. Other examples of irregular blastwaves include the more extreme cases of SN1986J and SN1993J. VLBI observations of these SNRs (Bartel et al., 1991; Bietenholz, Bartel, & Rupen, 2001) show shells of radio emission with large radial protrusions, suggesting that very young SNRs can exhibit large deviations from spherical symmetry.

Another characteristic of young SNRs that may originate with the convective instability associated with the contact discontinuity is the slight dominance of a radial magnetic field as deduced from radio polarization measurements (Dickel et al., 1991). Jun & Norman (1996b) suggested that the polarization observations could be explained by convective instabilities dragging out an ambient magnetic field as the R-T fingers pushed from the contact discontinuity out toward the forward shock.

Despite the belief that convective instabilities alone should play an important role in mixing ejecta out to the forward shock and generating the observed radial magnetic fields, attempts to demonstrate this have failed. Chevalier et al. (1992) used a linear perturbation analysis and two-dimensional hydrodynamic simulations to study this instability, and found that the mixing region did not reach the forward shock. This result was confirmed by Jun & Norman (1996a) and found to hold true even in the presence of rapid cooling of the shocked ejecta (Chevalier & Blondin, 1995). Jun & Norman (1996b) showed that this limited mixing region was inconsistent with observations, for it produced a SNR with a two-shell structure when viewed in radio-synchrotron emission: an inner shell of strong emission and a dominant radial magnetic field, and an outer shell of weaker emission with a dominant tangential field. Jun & Jones (1999) extended this work by including a simplified scheme for injecting and accelerating test-particle electrons and calculated the radio synchrotron emission self-consistently in the turbulent fields. Their work emphasized the role of the reverse shock in accelerating electrons and producing radio emission.

One possible solution to these problems is the introduction of small-scale structure into the density of either the circumstellar medium (CSM) or the supernova ejecta. Jun, Jones, & Norman (1996) found that including small cloudlets in the surrounding CSM could enhance the growth of the Rayleigh-Taylor fingers enough so they reached all the way to the shock front. Alternatively, Blondin, Borkowski, & Reynolds (2001) examined the effects of including low-density bubbles in the supernova ejecta. They also found an enhanced turbulence in the region between the forward and reverse shocks, with some clumps of ejecta reaching out to the forward shock.

All of this previous work concerning R-T instabilities and driven waves has assumed the supernova ejecta and CSM can be treated as an ideal gas, typically modeled with an adiabatic index, $`\gamma =5/3`$.111 For discussions of other types of instabilities in shocks accelerating cosmic rays see, for example, Drury & Falle (1986) and McClements et al. (1996). However, these assumptions may not be appropriate for young SNRs, where the strong shock waves are expected to be efficient accelerators of energetic particles (i.e., cosmic rays). 222 In the heliosphere, where shocks are observed directly with spacecraft, there is clear evidence that the quasi-parallel earth bow shock (with a sonic Mach number, $`M_{\mathrm{S0}}<10`$) can place 10-30% of the solar wind kinetic energy flux into superthermal particles (e.g., Ellison et al., 1990). Quasi-perpendicular interplanetary shocks (IPSs) also accelerate ambient particles but generally with lower efficiencies due to their lower Mach numbers (generally $`M_{\mathrm{S0}}\begin{array}{c}<\hfill \\ \hfill \end{array}3`$ for IPSs) (e.g., Baring et al., 1997). In a few cases, however, IPSs with higher than average Mach numbers have been observed to accelerate particles with efficiencies comparable to the quasi-parallel Earth bow shock (Eichler, 1981; Terasawa et al., 1999). Hybrid plasma simulations showing direct injection and acceleration of ambient particles are reasonably consistent with these observations (e.g., Scholer et al., 1992), as are convection-diffusion models (e.g., Kang & Jones, 1997). Roughly, quasi-parallel and quasi-perpendicular shocks are those with an angle between the upstream magnetic field and the shock normal of less than $`45^{}`$ and greater than $`45^{}`$, respectively.  It has long been believed that SNRs are the primary sources of galactic cosmic rays below $`10^{15}`$ eV (e.g., Axford, 1981). However, for this to be the case, the acceleration mechanism must be efficient and place on the order of $`10\%`$ or more of the total ejecta kinetic energy into relativistic particles (e.g., Blandford & Eichler, 1987; Drury, Markiewicz, & Völk, 1989; Dorfi, 2000). Radio observations have long been proof that SNRs produce GeV electrons and recent evidence of X-ray synchrotron emission in young SNRs suggests they can produce TeV electrons as well (e.g., Koyama et al., 1995; Mastichiadis & de Jager, 1996; Allen et al., 1997; Tanimori et al., 1998; Reynolds & Keohane, 1999). While direct proof of ion acceleration in SNRs remains elusive, indirect evidence exists (e.g., Ellison, Slane & Gaensler, 2001) and there is a growing body of evidence linking the efficient production of cosmic-ray ions by the forward and reverse shocks with the thermal properties of the shock-heated X-ray emitting gas (e.g., Dorfi & Böhringer, 1993; Decourchelle, Ellison, & Ballet, 2000; Hughes, Rakowski, & Decourchelle, 2000b).

One of the important structural effects of efficient particle acceleration is a shock compression ratio greater than 4 (e.g., Ellison & Eichler, 1984; Berezhko & Ellison, 1999). For an ideal gas in the absence of particle acceleration, a high Mach number shock produces a compression ratio,

$$\sigma \frac{\gamma +1}{\gamma 1}.$$
(1)
For the nominal value of $`\gamma =5/3`$, $`\sigma =4`$. However, when efficient particle acceleration occurs, relativistic particles can be created and substantially contribute to the pressure, and superthermal particles can escape from the shock. These two effects combine to produce a shocked plasma which acts to a large extent like a gas with an effective adiabatic index, $`\gamma _{\mathrm{eff}}<5/3`$, allowing arbitrarily large compression ratios. In fact, it has been shown (Berezhko & Ellison, 1999) that, when injection of particles from the background occurs easily, nonlinear diffusive shock acceleration yields

$$\sigma \{\begin{array}{cc}1.3M_{\mathrm{S0}}^{3/4}\hfill & \text{if }1M_{\mathrm{S0}}^2M_{\mathrm{A0}}\text{ ;}\hfill \\ \multicolumn{2}{c}{}\\ 1.5M_{A0}^{3/8}\hfill & \text{if }1M_{\mathrm{A0}}<M_{\mathrm{S0}}^2\text{ .}\hfill \end{array}$$
(2)
Here, $`M_{\mathrm{S0}}`$ ($`M_{\mathrm{A0}}`$) is the upstream sonic (Alfvén) Mach number (see also Kazanas & Ellison, 1986; Malkov, 1997). Furthermore, the increase in compression ratio is accompanied by a decrease in the post-shock temperature which links the efficient production of superthermal particles with the thermal properties of the shock-heated, X-ray emitting gas.

Of particular importance for the work reported here, the increased compression ratio also implies that the region between the forward and reverse shocks in young SNRs will be significantly narrower and denser when acceleration is efficient compared to the case where little energy is placed in relativistic particles. The increased density and narrower spatial extent of the interaction region leads to large density gradients suggesting that R-T instabilities may have faster growth rates and larger amplitudes at saturation than is commonly assumed (Decourchelle, Ellison, & Ballet, 2000). The dependence of the compression ratio on the sonic and Alfvén Mach numbers implies that these effects will be most pronounced in SNRs having high shock speeds and relatively low magnetic fields.

In this paper we explore the possibility that the familiar convective instability in young SNRs may be altered by high compression ratios. We use a well-tested hydrodynamic simulation of evolving SNRs, only replacing the standard $`\gamma =5/3`$ with an effective adiabatic index, $`\gamma _{\mathrm{eff}}<5/3`$. A softer equation of state (a low $`\gamma _{\mathrm{eff}}`$) causes the shocked plasma to be more compressible, yielding shock compression ratios $`\sigma >4`$, consistent with those expected when particle acceleration is efficient.

Somewhat unexpectedly, we find that the convective instability growth rate and mixing length are almost independent of the compression ratio. However, since the region between the forward and reverse shocks narrows as $`\sigma `$ increases, the likelihood that Rayleigh-Taylor fingers reach and perturb the forward shock is strongly correlated with $`\sigma `$ and, therefore, with the efficiency of cosmic-ray production. If $`\sigma `$ is large enough, the convective instabilities do reach and perturb the forward shock making it likely that clumps and filaments of dense ejecta material will be found there.

In section 2 we discuss the impact of a low value of $`\gamma _{\mathrm{eff}}`$ on a one-dimensional dynamical model of young SNRs. In section 3 we present hydrodynamic simulations in one, two, and three dimensions, and in section 4 we discuss the approximations and implications of our numerical results.

## 2 A DYNAMICAL MODEL

The dynamical evolution of young SNRs is determined by the interaction of the stellar ejecta, one to several solar masses of material thrown out by the SN explosion, with the surrounding circumstellar medium. In the simplest case, this interaction is characterized by a two shock structure: a forward shock driven ahead of the ejecta into the CSM, and a reverse shock that decelerates the supersonically expanding ejecta in the vicinity of the interaction region. To model this ejecta-dominated evolution, we begin with the self-similar driven wave (SSDW) model (Chevalier, 1982a), although we have found results similar to those reported here for the exponential ejecta model (Dwarkadas & Chevalier, 1998).

Chevalier (1982a) provided spherically symmetric, self-similar solutions to the structure of ejecta-dominated SNRs under the assumption that the mass density in both the ejecta and the CSM could be described by power laws. Since the steep power law density profiles for the ejecta lead to an infinite ejecta mass, we follow Chevalier and institute a cutoff in the ejecta density profile at small radii. Therefore, at radii less than some critical value, $`r_c=v_ct`$, the ejecta density is taken to be a spatially constant plateau, while beyond this radius it follows a power law, i.e.,

$$\rho _{\mathrm{ej}}(r,t)=\{\begin{array}{cc}\rho _cv_c^nr^nt^{n3}\text{for}\hfill & r>v_ct\hfill \\ \rho _ct^3\text{for}\hfill & r<v_ct\hfill \end{array}$$
(3)
where $`t`$ is time, the constant $`\rho _c`$ is given by

$$\rho _c=\frac{5n25}{2\pi n}E_{\mathrm{sn}}v_c^5,$$
(4)
the velocity at the intersection of the ejecta density plateau and the power law is given by

$$v_c=\left(\frac{10n50}{3n9}\frac{E_{\mathrm{sn}}}{M_{\mathrm{ej}}}\right)^{1/2},$$
(5)
and $`n`$ is a constant index. These quantities, $`\rho _c`$ and $`v_c`$, are related to the constant $`g`$ in Chevalier (1982a) by $`\rho _cv_c^n=g^n`$. Here, $`E_{\mathrm{sn}}`$ and $`M_{\mathrm{ej}}`$ are the kinetic energy and mass of the ejecta, respectively, and it’s clear that fixing the plateau allows us to relate the self-similar parameters to physical quantities normally associated with SNRs. Furthermore, we restrict our examples to times such that $`r_c`$ is less than the radius of the reverse shock to avoid complications from multiply reflected shocks etc., and to make direct comparisons between our numerical results and self-similar solutions straightforward.

The density of the ambient circumstellar medium is also described by a power law in radius in the Chevalier (1982a, b) solution, i.e.,

$$\rho _a=qr^s,$$
(6)
where $`s=0`$ corresponds to a uniform density (in which case $`q`$ is the ambient density) and $`s=2`$ represents the density profile of a steady-state stellar wind from the supernova progenitor of speed, $`v_w`$, and mass loss rate, $`dM/dt`$, such that $`q=(dM/dt)/(4\pi v_w)`$.

Thus, provided that the edge of the ejecta plateau region has not yet reached the reverse shock, the radius and expansion velocity of the forward shock are:

$`R_1`$ $`=`$ $`b\left({\displaystyle \frac{\rho _cv_c^n}{q}}\right)^{1/(ns)}t^{(n3)/(ns)},`$ (7)
$`V_1`$ $`=`$ $`\left({\displaystyle \frac{n3}{ns}}\right){\displaystyle \frac{R_1}{t}},`$ (8)
and the fluid variables in the shocked region immediately behind the forward shock are:

$`\rho _1`$ $`=`$ $`{\displaystyle \frac{\gamma _{\mathrm{eff}}+1}{\gamma _{\mathrm{eff}}1}}qR_1^s,`$ (9)
$`u_1`$ $`=`$ $`{\displaystyle \frac{2}{\gamma _{\mathrm{eff}}+1}}\left({\displaystyle \frac{n3}{ns}}\right){\displaystyle \frac{R_1}{t}},`$ (10)
$`p_1`$ $`=`$ $`{\displaystyle \frac{2q}{\gamma _{\mathrm{eff}}+1}}\left({\displaystyle \frac{n3}{ns}}\right)^2{\displaystyle \frac{R_1^{2s}}{t^2}},`$ (11)
where $`\rho `$ is the density, $`u`$ is the fluid speed, $`p`$ is the pressure, and the subscript $`1`$ indicates the forward shock. Note that all of the dependence on the adiabatic index in eq. (7) is contained in the constant $`b`$, which is related to parameters defined in Chevalier (1982a) through the expression $`b=(R_1/R_c)A^{1/(ns)}`$. The calculated values of the constant $`A`$ and the ratio of the forward shock radius to the radius of the contact discontinuity for $`\gamma =5/3`$ for various $`n`$ and $`s`$ are given in Table 1 of Chevalier (1982a).

We note that in an attempt to predict the $`\gamma `$-ray flux from pion-decay, Chevalier (1983) considered two-fluid, self-similar solutions where one fluid was a thermal gas with $`\gamma =5/3`$ and the other was a relativistic gas (i.e., cosmic rays) with $`\gamma =4/3`$. The effective adiabatic index was determined by

$$\gamma _{\mathrm{eff}}=\frac{5+3w}{3(1+w)},$$
(12)
where $`w`$ was the (constant) fraction of total pressure made up by relativistic particles at the shock. Chevalier did not address the possible effects $`\gamma _{\mathrm{eff}}`$ would have on R-T instabilities nor attempt to model the effects of particle escape during acceleration and assumed $`4/3\gamma _{\mathrm{eff}}5/3`$ with a maximum $`\sigma =7`$ in his calculations. Chevalier did note, however, that cosmic rays diffusing from the shocked gas might produce a dense shell with $`\sigma >7`$ similar to that which occurs in the later radiative phase.

The structure of a SSDW is illustrated in Figure 1, which shows the radial profiles of density, pressure, and velocity for values of $`n=7`$ and $`s=2`$. The interaction region is bounded by the forward and reverse shocks, which for this example have compression ratios $`\sigma =4`$ corresponding to a high Mach number adiabatic shock in an ideal fluid with $`\gamma =5/3`$. The pressure is monotonically increasing from the reverse shock to the forward shock, as all of the shocked gas is gradually decelerating.

An important characteristic of these solutions for $`s=2`$ is the convectively unstable region ahead of the contact discontinuity (Chevalier et al., 1992). This region is marked by opposing signs of the pressure and density gradients. The positive pressure gradient is gradually decelerating the shocked ejecta and shocked CSM, but the negative density gradient implies that the denser gas is being decelerated by lower density gas; a situation that is subject to the familiar R-T instability. A similar situation holds for $`s=0`$ solutions, but in this case the opposing gradients are behind the contact discontinuity rather than in front.

While the idealized self-similar hydrodynamical model has been very successful in describing general aspects of SNR evolution, it is clearly incomplete in some important ways. For instance, the hydrodynamic model ignores any effects from a magnetic field, the one-dimensional approximation eliminates the ability to describe irregularities in the ambient conditions or Rayleigh-Taylor generated instabilities, and particle acceleration at the collisionless forward and reverse shocks is neglected. Here we focus on how efficient cosmic-ray production influences R-T instabilities by setting $`\gamma _{\mathrm{eff}}<5/3`$. We find that lowering $`\gamma _{\mathrm{eff}}`$ causes the blast wave to expand slightly slower, but the radius of the forward shock still follows the self-similar result $`R_1t^{(n3)/(ns)}`$. More importantly, lowering $`\gamma _{\mathrm{eff}}`$ causes the width of the region between the forward and reverse shocks to shrink considerably and alters the internal structure of the self-similar driven wave.

We list some of the parameters describing SSDWs for a range of $`\gamma _{\mathrm{eff}}`$ in Table 1 and, for $`\gamma =5/3`$, these values are identical to those in Table 1 of Chevalier (1982a). These parameters, as well as the radial profiles shown in Figure 2, were generated following the procedure outlined in Chevalier (1982a). Note in particular that the value of $`b`$, which contains the $`\gamma _{\mathrm{eff}}`$ dependence, changes by less than 30% as the compression ratio is varied from 4 to 21. Thus the velocity of the blast wave is not strongly affected by the value of the compression ratio. In contrast, between $`\sigma =4`$ and $`21`$, the width of the interaction region shrinks by a factor of 4 in the case of $`n=7`$ and $`s=0`$.

The radial profiles of SSDWs for different values of $`\gamma _{\mathrm{eff}}`$ are shown in Figure 2. While the increase in the compression ratio and decrease in the width of the interaction region are clearly evident, there is also a change in the density gradient of the shocked CSM. As gas moves from the forward shock to the contact discontinuity it compresses/expands adiabatically to match the relatively constant pressure of the interaction region. A lower $`\gamma _{\mathrm{eff}}`$ leads to a slower drop in pressure due to radial expansion, so the gas does not need to be compressed as much as with a higher $`\gamma _{\mathrm{eff}}`$. The result is a more positive density gradient in the shocked CSM for a smaller $`\gamma _{\mathrm{eff}}`$. From Figure 2 we see that this can even lead to a change in the sign of the density gradient, which may in turn affect the stability properties of the waves. Nonetheless, we expect the low $`\gamma _{\mathrm{eff}}`$ driven wave to remain unstable, as there is still a shell of dense shocked ejecta being decelerated by lower-density shocked CSM.

## 3 HYDRODYNAMIC SIMULATIONS

To examine the effects of different compression ratios on the convective stability of driven waves, we use the Virginia Hydrodynamics (VH-1) time-dependent hydrodynamic numerical code in one, two, and three dimensions. All of the simulations were computed in spherical geometry on a numerical grid with 500 zones in each dimension. This resolution is comparable to the highest resolution simulations of Chevalier et al. (1992). Based on their results, we expect that a higher resolution simulation would show more small scale mixing, but the large-scale dynamics of the problem would not change. The angular extent of the simulation domain was chosen to produce roughly square numerical zones, i.e., $`R\mathrm{\Delta }\theta \mathrm{\Delta }R`$. Periodic boundary conditions were applied in both the $`\theta `$ and $`\varphi `$ directions for the multidimensional simulations. The radial boundary conditions were set to match the relevant power laws. The numerical grid was expanded to follow the forward blast wave so the evolution could be tracked for many expansion times.

We ran one-dimensional (1D) and two-dimensional (2D) simulations with the four combinations of $`s=0,2`$ and $`n=7,12`$. For each $`(n,s)`$ pair we ran four simulations with $`\gamma _{\mathrm{eff}}=5/3,4/3,1.2`$, and 1.1. In the following discussions we will focus on the set of simulations with $`n=7`$ and $`s=2`$ as a specific example, but the results are qualitatively the same for all of these parameter choices.

### 3.1 1-D SIMULATIONS

We ran 1D simulations for all the cases listed in Table 1 as a check on our ability to numerically recreate the analytic, self-similar solutions. The 1D hydrodynamic simulations were initialized with appropriate density power laws, and allowed to evolve until a self-similar state was reached. The gas pressure in the ejecta and CSM was kept sufficiently low to ensure that the Mach numbers of the forward and reverse shocks were always greater than 100.

In Figure 3 we show the results of our most extreme case with $`n=12`$ and $`s=0`$ and $`\gamma =1.1`$ overlayed on the semi-analytic solution obtained by numerically integrating an ordinary differential equation (Chevalier, 1982a). This particular SSDW is very thin and possesses steep gradients, making it a difficult test case for the hydrodynamic code. Nonetheless, our simulation does a reasonably good job of matching the analytic solution; note in particular the relative sharpness of the shock fronts and contact discontinuity. In the less extreme cases, the match between the simulation and the analytic solution is better than shown in Figure 3. However, the very sharp density peak at the contact discontinuity in the $`s=2`$ solutions, and the density trough in the $`s=0`$ solutions, are typically smoothed out over a few zones by this hydrodynamic method.

### 3.2 2-D SIMULATIONS

The multi-dimensional simulations were initiated with the analytic SSDW solutions normalized to $`R_1=1`$, $`V_1=1`$, and a preshock density of unity. The simulation thus begins at an initial time of $`t=(n3)/(ns)`$ in these normalized units. The initial solution is mapped onto a spherical grid with a radial span roughly twice the width of the interaction region. The convective instability was seeded by adding acoustic noise to the shocked interaction region. These simulations were evolved long enough for the convective instability to reach saturation, such that the mixing region is in a quasi-self-similar state with the growth of R-T fingers matched by their destruction due to shearing and advection. This typically required $`6`$–8 doubling times of the blast wave radius, and beyond this time, the results are qualitatively independent of the seed noise.

In all of the cases listed in Table 1, the layer of shocked ejecta was found to be unstable with an evolution similar to that seen in previous studies (i.e., Chevalier et al., 1992; Jun & Norman, 1996a). As an example, we show the density for driven waves with $`n=7`$, $`s=2`$, and different values of $`\gamma _{\mathrm{eff}}`$ in Figure 4. The results for other combinations of $`n`$ and $`s`$ are similar. These images are taken at the end of the simulations, once the instability has reached saturation. In addition to shrinking the interaction region, a lower value of $`\gamma _{\mathrm{eff}}`$ produces a narrower, denser shell of shocked ejecta and sharper, denser R-T fingers. Aside from the narrower interaction region, these results for small $`\gamma _{\mathrm{eff}}`$ are similar to those found by Chevalier & Blondin (1995) for driven waves with radiative cooling in the shocked ejecta. In this latter case the ejecta behaved as an ideal gas with $`\gamma 1`$, while the CSM was modeled as an ideal gas with $`\gamma =5/3`$.

We can obtain a rough estimate of the instability growth rate by following the turbulent energy density in the intershock region. Rather than attempt to subtract off the bulk radial velocity to find the contribution to the turbulent velocity, and to avoid confusion in comparing 2D and 3D simulations, we consider only a single angular component of the turbulent motion. Furthermore, to remove the effects of radial expansion, deceleration, and shock compression, we normalize this turbulent energy to the bulk kinetic energy density associated with the shock front. We thus define the ratio

$$\chi =\frac{\rho u_\theta ^2𝑑\tau }{\rho _1V_1^2𝑑\tau },$$
(13)
where the integration is only over the volume of shocked gas. The evolution of this turbulent energy parameter is shown in Figure 5, from which one can see that while the saturated levels of turbulent energy density differ by about a factor of two, the growth rates are not noticeably affected by changes in $`\gamma _{\mathrm{eff}}`$, despite the fact that the radial profile from the spherical solution has changed dramatically (e.g., the density gradient has changed sign).

The saturation of the instability, at least in terms of the width of the mixing region, is also relatively unaffected by changes in $`\gamma _{\mathrm{eff}}`$. In Figure 6 we plot an angle-averaged radial profile of the ejecta mass fraction. In a spherically-symmetric model this would be a step function from unity dropping to zero at the contact discontinuity. To provide an easy comparison, we have normalized the radii such that the reverse shock is located at the same radius for all models. From this we can see that the length of the R-T fingers is comparable for all four values of $`\gamma _{\mathrm{eff}}`$.

While the instability itself is insensitive to changes in the compression ratio, the overall width of the interaction region decreases with increasing compression ratio. This has the important effect that, for small $`\gamma _{\mathrm{eff}}`$, the R-T fingers can reach as far as the forward shock front. To measure the maximum extent of the R-T fingers, we identify an average mass fraction of 0.001 as the leading edge of the mixing region. Since the radial gradient of the mass fraction is quite steep at these small values, the exact choice of a cutoff value does not affect the results. The time dependence of the maximum radial extent (relative to the forward shock) of the R-T fingers is illustrated in Figure 7. In the run with the highest compression ratio ($`\gamma _{\mathrm{eff}}=1.1`$), the R-T fingers quickly reach and pass the average shock radius. We note that when the curves in Figure 5 become approximately flat, the system is in a near self-similar state and the ratio of the forward shock radius to the contact discontinuity radius remains approximately constant. If we allowed the system to evolve to the point where the reverse shock entered the plateau region in the ejecta density (which nominally happens when the mass of swept up CSM is comparable to the ejecta mass), the ratio of the forward shock radius to the contact discontinuity radius would start to increase and the R-T fingers would drop behind the forward shock.

Once the fingers of dense, shocked ejecta reach the forward shock, they can begin to distort the outer blast wave by pushing small regions out ahead of the average shock radius. The resulting bumps in the outer shock can be seen in Figure 4 for the simulation with $`\gamma _{\mathrm{eff}}=1.1`$. The shock front is nearly spherical for large values of $`\gamma _{\mathrm{eff}}`$, while values close to unity produce significant deviations from spherical symmetry.

We note that even in the most dramatic case of $`n=12`$, $`s=0`$, and $`\gamma _{\mathrm{eff}}=1.1`$, the deviation from a spherical outer shock remains relatively small; the penetrating fingers are not able to dramatically alter the forward shock. This can be seen in Figure 7, where the maximum radial extent of the fingers appears capped at $`R_1`$. Since the instability is relatively unaffected by the value of $`\gamma _{\mathrm{eff}}`$, we expect Figure 7 to show a series of parallel curves as the fingers grow until saturation occurs (at several 100 times the initial simulation time, although this depends on the magnitude of the initial perturbation). If the fingers reach the forward shock before saturation, as in the $`\gamma _{\mathrm{eff}}=1.1`$ case, they cannot push out much beyond the shock and their growth stalls. This limited ability to affect the shock front is due to the strong shearing flow created when the shock front is distorted. If a clump of dense ejecta has sufficient radial momentum to push out the forward shock in a local protrusion, the deformed shock front generates a substantial tangential post shock flow around the ejecta clump. In all of our simulations, this shear flow quickly disrupts the clump through the Kelvin-Helmholtz instability, with the remnants of the clump quickly advected back into the interaction region. Thus we never see a local protrusion of the shock front stick out more than a few percent of the blast wave radius.

### 3.3 3-D SIMULATIONS

We repeated two simulations in 3D to check that these results are not dramatically affected by the assumption of axisymmetry, using $`n=12`$ and $`s=0`$ for $`\gamma _{\mathrm{eff}}=5/3`$ and 1.1. In comparing the 2D and 3D simulations, we stress that the numerical code remained virtually identical. The only changes involved using a 3D grid and repeating the hydrodynamic updates in a third direction.

For the more familiar case of $`\gamma _{\mathrm{eff}}=5/3`$ our results are quite similar to previous 3D simulations of driven waves (e.g., Jun & Norman, 1996b), showing only minor differences between 2D and 3D simulations. This similarity is exhibited in Figure 8. The only significant difference visible in this Figure is the increased amount of small-scale structure in the 3D simulation, although this difference does not appear to affect any of the global properties of the driven wave. A quantitative comparison of these two simulations finds that the growth and saturation of the turbulent energy is virtually identical, the radial extent of the fingers is comparable (the 3D fingers reached slightly further than the 2D fingers), the forward shock remains spherical and the width of the interaction region grows slightly in both cases as the instability approaches saturation.

The situation is somewhat more complicated for the case of high shock compression, as shown in Figure 9. Here one sees much more small-scale structure in the 3D simulation, both within the interaction region and in the forward and reverse shocks. This more dramatic difference between 2D and 3D in the high-compression case is consistent with the fact that the R-T fingers do not reach the shock front when the shock compression is low, but they do reach - and perturb - the forward shock when the compression is large. For both values of $`\gamma _{\mathrm{eff}}`$ we see more small-scale structure in the R-T fingers in 3D, but in the high-compression case this small-scale structure can modify the forward shock to produce many small wavelength, but large amplitude, perturbations.

A consequence of the more oblique shocks created by the deformation of the forward shock in the 3D simulation is an overall drop in the average compression ratio and a corresponding increase in the width of the interaction region. Despite these differences, our statistical measures of the convective instability remain relatively unaffected by the dimensionality of the simulation. Furthermore, our primary conclusion from the 2D simulations still holds in 3D; the R-T fingers are able to deform the forward shock when the compression ratio is high.

## 4 DISCUSSION

If young supernova remnant shocks accelerate cosmic rays with high efficiency, the acceleration process can cause these high Mach number shocks to have compression ratios considerably greater than four. We have investigated the effects of this high compression on the hydrodynamic stability of young remnants with a simple and direct method; we perform hydrodynamic simulations with the effective adiabatic index, $`\gamma _{\mathrm{eff}}`$, set to values less than $`5/3`$, yielding compression ratios, $`\sigma (\gamma _{\mathrm{eff}}+1)/(\gamma _{\mathrm{eff}}1)>4`$. While this fluid approach is clearly an approximation to the real situation where the collisionless shocks accelerate particles to relativistic energies, we believe it accounts for the major effects of efficient particle acceleration on the overall dynamics and, most importantly, on the Rayleigh-Taylor instabilities that develop in the interaction region between the forward and reverse shocks.

Somewhat surprisingly, we find that a large compression ratio has little effect on the growth of R-T instabilities. However, since a large $`\sigma `$ dramatically shrinks the width of the interaction region between the forward and reverse shocks, several important effects emerge from of our high-compression simulations:

(1) In contrast to shocks with $`\sigma 4`$, R-T fingers reach closer to the forward shock front, and if $`\sigma `$ is large enough, ejecta material can be found at, and even slightly ahead of the average shock radius. This can occur fairly quickly after the explosion (depending on how quickly a reverse shock forms and the magnitude of any inhomogeneities seeding this instability) and the fingers should stay near the forward shock front throughout the time the reverse shock is in the power law portion of the ejecta density profile. In general this should apply to SNRs from an age of only months up to several 1000 years.

(2) The forward shock is perturbed on short wavelengths and with relatively small amplitudes. In addition to slightly altering the spherical shape of the shock front, this will lead to some small spread in the postshock temperatures. Note also that the overall temperature of the shocked gas will be substantially lower if particle acceleration is efficient and compression ratios are large than if little acceleration occurs (e.g., Ellison et al., 2000). This effect provides a coupling between the shock morphology and thermal X-ray emission (Decourchelle, Ellison, & Ballet, 2000; Hughes, Rakowski, & Decourchelle, 2000b).

(3) If the mixing region reaches the forward shock, the morphology of the magnetic field as seen through radio synchrotron polarization observations may be affected. Jun & Norman (1996b) followed the evolution of the ambient magnetic field in driven wave simulations, looking for an explanation for the origin of observed polarization in young SNRs in the elongation of the field by R-T instabilities. Our work offers a ready explanation for why this polarization can extend all the way to the forward shock.

While we do not include magnetic fields in the hydrodynamic simulations we perform here, we note that any ambient magnetic fields will also be compressed at the shock front. In addition, the increased turbulence associated with the high compression driven waves may be expected to further amplify the magnetic fields and show stronger radio emission and possibly more intense TeV emission from inverse-Compton and/or pion-decay than the more quiescent driven waves found for $`\gamma _{\mathrm{eff}}=5/3`$. It can be expected that the effects we see here from increased compression will add to those reported by Jun & Jones (1999).

The SNR 1E 0102.2–7219 in the Small Magellanic Cloud may represent a young remnant for which this model of a high compression driven wave is applicable. Recently, Hughes, Rakowski, & Decourchelle (2000b) have determined postshock electron temperatures in SNR 1E 0102.2–7219 using Chandra X-Ray Observatory observations. Using the measured forward shock speed of $`V_{\mathrm{FS}}6000`$ km s-1, they find that the electron temperature of 0.4–1 keV is at least 2.5 times lower than can be explained with standard (i.e., test-particle) shock heating even if only Coulomb electron heating occurred. They conclude that the forward shock is placing at least 50% of the shock kinetic energy flux in cosmic ray ions.

If one applies a SSDW model to the observed radius and expansion velocity of 1E 0102.2–7219 (e.g., Gaetz et al., 2000; Hughes, Rakowski, & Decourchelle, 2000b) assuming $`n=9`$ and $`s=0`$, one finds an age of 700 years, independent of the compression ratio. Choosing reasonable parameters of $`E_{\mathrm{sn}}=10^{51}`$ erg and $`M_{\mathrm{ej}}=1M_{}`$, one finds an ambient density of $`\rho _a0.1m_p`$ cm-3for $`\sigma =4`$. If instead we assume a large compression ratio of $`\sigma =11`$, the derived ambient density is lower by only a factor of 2.3. Despite the relative similarity in fitting the observed SNR parameters, these models, shown in Figure 10 are dramatically different in terms of the position of the reverse shock and the postshock temperature. If the shock compression is high, the reverse shock will be significantly closer to the forward shock, as shown in Figure 10. This is not consistent with X-ray observations which show the presence of shocked ejecta at a radius of only $`70\%`$ of the forward shock (Gaetz et al., 2000). However, if the ejecta mass is relatively small, this SNR may be evolving away from the SSDW phase, with the reverse shock now propagating in toward the center of the SNR.

With these values and one additional parameter, the ambient magnetic field strength $`B_0`$, we can estimate the particle acceleration using the model of Berezhko & Ellison (1999). We find that with $`B_03`$ $`\mu `$G, the forward shock obtains $`\sigma 11`$ (consistent with $`\gamma _{\mathrm{eff}}=1.2`$) with a shocked proton temperature of $`T_{p2}4\times 10^7`$ K, i.e., 10 times lower than obtained in a test-particle shock with the same parameters (for the acceleration calculation we assume a constant $`V_{\mathrm{FS}}=6000`$ km s-1 for $`700`$ yr). The low postshock electron temperature deduced from X-ray observations is now easily explained if the electron temperature is $`1/3`$ $`T_{p2}`$.333If we had used a lower $`\gamma _{\mathrm{eff}}`$, we could have obtained a consistent fit with a smaller $`T_{p2}`$. Furthermore, the forward shock accelerates electrons to $`70`$ TeV in $`700`$ yr, consistent with radio synchrotron emission. At the current age, $`M_{\mathrm{S0}}3700`$ and $`M_{\mathrm{A0}}150`$, and the acceleration is extremely efficient, i.e., more than 80% of the kinetic ram energy flux is placed in relativistic ions. While this model is not unique, it does show that efficient particle acceleration is consistent with reasonable supernova and ambient medium parameters, as well as with deduced values of radius, speed, age, and electron temperature.

While there are clearly differences between the actual situation in SNRs where some fraction of swept-up material is shock accelerated to relativistic energies, and simply lowering $`\gamma _{\mathrm{eff}}`$, we believe our results are qualitatively correct. The important differences include the fact that lowering $`\gamma _{\mathrm{eff}}`$ implies that the effects of relativistic particle pressure and particle loss occur everywhere rather than just in the shocked gas and in the precursor regions in front of the forward and reverse shocks. We note that if the finite size of the shock is limiting the maximum particle energy, there will be a significant fraction of the total pressure in high energy particles with upstream diffusion lengths on the order of $`1/10`$ of the shock radius. This difference in pressure distribution will have some effect on the shock evolution since, instead of a uniform $`\gamma _{\mathrm{eff}}`$, the actual remnant has a gas with a soft equation of state pushing one with a harder equation of state. While we have not tested this difference in detail, we expect the effects we found will actually be enhanced if $`\gamma _{\mathrm{eff}}`$ varies spatially.

Another difference comes about because we use eq. (1) with a constant $`\gamma _{\mathrm{eff}}`$ to determine $`\sigma `$ instead of eqs. (2). With (1), there is little variation in $`\sigma `$ as the SNR evolves as long as $`M_{\mathrm{S0}}1`$, while eqs. (2) can give a much larger variation depending on the parameters. 444 We note that it is not always the case that lowering the Mach numbers causes $`\sigma `$ to decrease. Equations (2) apply when injection into the shock acceleration mechanism is efficient. If injection is weak enough, high Mach number, test-particle solutions with $`\sigma 4`$ can result (see Berezhko & Ellison, 1999, for a full discussion).  In any case, since in all probability $`\sigma 4`$ during the time we consider before the reverse shock enters the ejecta plateau region, any differences resulting from $`\sigma `$ varying with time are likely to be small and go in the direction of increasing the effects we report.

Finally, since the shrinking of the interaction region between the forward and reverse shocks depends totally on the pressure in that region, effects other than particle acceleration which influence the pressure may be important. The most likely effect which we have neglected comes from the compression of the magnetic field, $`B`$. Since the magnetic pressure scales as $`B^2`$, increasing the compression ratio could produce magnetic pressures large enough to prevent the interaction region from becoming narrow enough to allow the R-T fingers to reach the forward shock. This, of course, will depend on the Alfvén Mach number and the angle the upstream field makes with the shock. The pressure effects of the magnetic field will be offset somewhat by the fact that, for a given compression ratio, the shocked thermal pressure is considerably less in a shock undergoing efficient acceleration compared to one with a low $`\gamma _{\mathrm{eff}}`$ and no acceleration.

We thank A. Decourchelle, J. P. Chièze, and S. Reynolds for helpful discussions. The numerical simulations reported here were performed at the North Carolina Supercomputing Center using 100 processors of an IBM SP2. We thank NCSC and IBM for their generous support of computing resources. Support for this work was provided by NASA under grant NAG-7153.