A Simulation Study of the Interaction between a Fast Shock and the Heliopuase

In a one-dimensional hybrid simulation of a fast shock impinging on the heliopause (HP), it is found that two regions of increased plasma density are formed between the transmitted and reflected shocks. The region downstream of the shock, on the interstellar side of the heliopause, has typically twice the original plasma density. The density increase is dependent on the impinging shock speed. This region has a greatly enhanced ion temperature anisotropy T(subscript ┴)/T(subscript ∥)-100 and a plasma beta β~0.1-1.0, which may lead to the generation of ion cyclotron waves and mirror waves. The region downstream of the reflected shock, on the solar wind side of the heliopause, has typically a 25% increase in plasma density. The enhanced temperature anisotropy is T(subscript ┴)/T(subscript ∥)~1.2, but the plasma beta is very high, β-20, which may also lead to the generation of mirror waves and ion cyclotron waves.


INTRODUCTION
Voyager spacecraft observations of very low frequency radio emissions(2-3.6 kHz) in the outer heliosphere [Kurth et al. 1984, Macek et al. 1991, Gurnett et al. 1993 have recently stirred interest in the termination shock, the outer solar wind, and the heliopause. The termina tion shock is a reverse fast shock where the solar wind jumps f r om supersonic to subsonic speed. The heliopause separates the high beta plasma of the solar wind f r om the low beta plasma of local interstellar space. Its existence was first proposed by Davis [1955]. Over the next decade, Voyager 1 and 2 and Pioneer 11 spacecraft may provide answers about the loca tion, shape, and extent of the termination shock.
Gurnett et al. [ 1993] proposed that the fast leading shock of a global merged interaction region (GMIR) triggered the observed radio emissions at or near the heliopause, and then estimated the distance to the heliopause to be between 116 and 177 AU. A GMIR is an ex-panding spherical shell of high-speed solar wind plasma having a thickness of several AU and is believed to be formed from the merging of a series of solar flares [Burlaga et al., 1993]. Whang and Burlaga [1994] and Whang et al. [1995] studied the interaction of an interplan etary MHD shock with the heliopause using an MHD formulation. They identified dense plasma regions which could account for the observed radio frequencies and also calculated possible locations for the termination shock and heliopause. Figure 1 shows the interaction of the 1991 GMIR shock with the termination shock and heliopause using the most likely loca tions for the termination shock and heliopause. The fast GMIR shock slows down after pass ing the termination shock and then crosses the subsonic outer heliosphere. Notice that the heliocentric distance to the termination shock increases after the passage of the fast GMIR shock. After contact with the heliopause, the fast GMIR shock splits into reflected and trans mitted shocks. Whang et al. [1995] used stagnation conditions downstream of the termination shock to estimate the average solar parameters just inside the heliopause. The requirement of pressure balance across the heliopause puts constrains on interstellar value of plasma tempera ture, density, and magnetic field.
In this paper we present numerical simulation results of the interaction of a fast shock with the heliopause. The simulation uses a hybrid code that combines a particle treatment of the ions with a massless fluid treatment of the electrons. The heliopause is modeled as a tan gential discontinuity (TD). First we discuss the simulation and the initial set-up. Next we present and explain the simulation results for a typical case. A large temperature anisotropy 7;_ I JZ = 100 is found downstream of the transmitted shock on the interstellar side of the heliopause. The effect of different shock speeds on the results is then examined by a parameter search. Finally we predict that the mirror and ion-cyclotron waves would be generated by anisotropic ion heating associated with the transmitted and reflected shocks.
Previous simulations of quasi-perpendicular shocks have shown that a substantial perpen dicular heating does occur [Lee et al., 1988;Schopke et al., 1990]. The generation of ion cyclotron waves in a plasma with a significant pressure anisotropy has been studied by many authors [e.g., Cornwall, 1965]. The generation of mirror waves by an anisotropic temperature has been studied by Hasagawa [1969], Price et al. [1986], Lee et al. [1988], and McKean et al. [1994]. Burgess and Schwarz [1988] and Wu et al. [1993] have simulated the interaction of a fast perpendicular shock with a tangential discontinuity in the solar wind. With a hybrid code, Mandt and Lee [1991] have simulated the interaction of a weak fast shock with a tangential discontinuity, representing Earth's magnetopause. In this paper, we simulate the interaction of a fast (GMIR) and the heliopause.

SIMULATION MODEL
The hybrid code used in our simulation is described by Swift and Lee [1983]. The code treats all variables as functions of one spatial component (x) and time. The ions are treated as particles and their dynamics are followed exactly (one spatial and all three velocity compo nents), while the electrons are treated as a massless fluid of finite pressure. The initial electron to ion pressure ratio is set to 1. The displacement current is dropped from the equations since we are interested in low frequency waves. Charge neutrality is assumed in solving the field The leading edge of the fast GMIR shock penetrated through the temination shock at -66 AU before it interacts with the heliopause at -150 AU. The impinge ment of the fast GMIR shock at the heliopause produces a reflected shock (R) and a transmitted shock (T). The two regions bounded by these shocks have increased densities and the corresponding local plasma frequencies that are in tile 2-3.6 kHz range. 347 equations. We do not use anomalous resistivity in the code. The numerical procedure inte grates two sets of equations, one for the vector field and one for the particles. The second order differential equation for the field is implicitly integrated using a finite differencing method that is accurate to second order. The particle velocity equation is integrated every half time step using a leapfrog finite difference scheme. Buffer zones and inactive reservoirs are set up at each end of the simulation domain to ensure that the interruption of particle orbits does not introduce spurious boundary currents. Each buffer zone has a thickness of an ion gyrodiameter. This prevents the disruption of the ion orbits that pass into the buffer zone. This procedure is accurate as long as any perturbations coming from the ends of the simulation domain do not reach the central third, where we will examine the numerical results. This hybrid code has been used successfully by others. Swift Figure 2 illustrates the various plasma regions (using density) in the simulation domain.
We model the heliopause (HP) as a TD which is centered in the simulation domain. The solar wind (region 1) and the downstream region of the fast GMIR shock (region 3) are on the right side of the HP and the interstellar plasma (region 2) is on the left. The resulting transmitted shock (TS) and reflected shock (RS) are the leading edges of region 4 and 5, respectively. The positive x-axis points toward the Sun along the shock normal and perpendicular to the mag netic field. The variation of plasma parameters across a discontinuity or shock is described by the Rankine-Hugoniot jump conditions. There are two conditions relating the plasma param eters across a TD: the normal net flow is zero and the total pressure (plasmas and magnetic) on either side remains equal. In our simulation, the plasmas on the two sides of TD are related by two ratios for density and magnetic field strength. The plasma and field parameter in region 3 are related to those in region 1 by the Rankine-Hugoniot conditions across the fast GMIR shock, which is determined by the upstream Alfven Mach number M A of the shock and the plasma /3.
The simulation is initialized with the fast GMIR shock just making contact with the HP, as in Fig. 2b. Thus there are initially only two regions in the simulation. First the interstellar plasma (region 2) is initialized in the code. Then the downstream plasma (region 3) is initial ized using the plasma density ratio, magnetic field ratio across the HP and the upstream Alfven Mach number of the fast GMIR shock. The plasma parameters used in our simulation for the case in which the fast GMIR shock Alfven number MA = 8.5 are listed in Table 1    shows the perpendicular ion pressure � ( e) the ion perpendicular plasma beta /3, 1-, and (f) the ion pressure anisotropy ratio � I �1• the tangential discontinuity's transition region is not yet discernible at t = 22Q3-1. The flow velocity in plot (c) is continuous across the HP and there is no net mass flux across the HP. In plot (f) it is seen that the temperature anisotropy is very large in region 4 since Ef._ I £!t = 'J;_ I� = 100. In region 5 the temperature anisotropy is found to be 1.2 (average). These two values of anisotropy can be expected from estimates using the Rankine-Hugoniot equations [Mandt and Lee, 1991). A large anisotropy can be generated in a low beta plasma even with a low shock Mach number, but not in a high beta plasma. Figure 4 is a time series plot of the density profile for the case shown in Figure 3. The heliopause ( H P) is centered by using a frame moving to the left at -275 km/s, slightly faster than the speed of the HP, which is convecting away from the sun. The transmitted shock (TS) and reflected shock (RS) fronts are clearly seen. From Table 1, the speed of a fast GMIR shock with M A = 8.5 is 465 km/s in the solar frame. The speed of the transmitted and reflected shocks are -350 km/s and -170 km/s, respectively. The heliopause moves from its original location at -180 km/s. The effects of GMIR shock speed are presented in Figure 5. Plot (a) shows the temperature anisotropy plotted vs. the resultant perpendicular beta /3 �for both downstream regions of the transmitted and reflected shocks. The results from five different shock speeds are shown.
For the transmitted shocks, the lowest /3 � is obtained for the case with GMIR shock speeds V s = 7.0V A which corresponds to the fast magnetosonic Mach number M Ms = 2.0. Here V A is the Alfven speed in region 1, the solar wind plasma just inside of the heliopause. Plot (b) shows the plasma density ratio across the transmitted and reflected shocks. The reflected shock in creases the density only slightly more than what the fast GMIR shock has already done. The downstream densities in regions 4 and 5 are obtained by averaging the density over the entire width of each of the two regions. The strength of a shock is best seen in the change of density. The transmitted shock is weaker than the original fast GMIR shock in all cases. For V s = 8.5V A , the density change across the GMIR shock is -2.18 and the density change across the transmitted shock is -1.9.  The greatly enhanced temperature anisotropy in region 4 is from nonadiabatic ion heating of the low f3 interstellar plasma across the transmitted shock. The region 4 plasma can be unstable to mirror and ion cyclotron wave mode due to the large temperature anisotropy. The region 5 plasma has a much lower temperature anisotropy of -1.2 but a higher plasma f3 -20, and this region is also unstable to mirror waves and ion cyclotron waves.

GENERATION OF MIRROR WAVES AND ION CYCLOTRON WAVES
The criterion for the mirror instability is J] I;; > 1 + 1/ /3 1-, which is plotted in Fig. 5a.
The stable and unstable regions are below and above the solid line, respectively. Fig. 5a shows that the plasmas in regions 4 and 5 for all cases simulated are unstable to mirror waves. The wave number k of the most unstable mode is typically k"" 0.5p7 1 where p,. is the local ion gyradius. The growth rate y of the most unstable mirror waves can be estimated as y ""0. 3Q ; , in region 4, and y :::: O.OlQ ; , in region 5, where Q ; is the local ion cyclotron frequency [e. g. , Price et al. , 1986;Lee et al. , 1988]. The temperature anisotropy in regions 4 and 5 can also excite ion cyclotron waves. The growth of unstable ion cyclotron wave growth driven by a temperature anisotropy has been previously explored in hybrid simulations [e. g. , Price et al. , 1986;Mckean et al. , 1994]. The wavenumber of the most unstable mode is k"" 0.5p71• The growth rate can be estimated as y ""O. lQ i ' However, the presence of minor ions may reduce the growth rate of the ion cyclo tron waves [Price at al., 1986].
In summary, we have used a hybrid code to simulate the interaction of the fast GMIR shock and the heliopause (HP). As a result of the interaction, a transmitted and reflected shock are present on the two sides of the heliopause. For the simulated parameters, the downstream region of the transmitted shock (interstellar side of H P) has an ion temperature anisotropy of l] I;; -100 and a plasma beta /3 -0.1 -1.0. In the downstream region of the weaker reflected shock, l] I;; -1.2 and /3 -20. In both region, it is expected that both mirror waves and ion cyclotron waves are excited as a result of temperature anisotropy.