Plasma environment of magnetized asteroids : a 3-D hybrid simulation study

The interaction of a magnetized asteroid with the solar wind is studied by using a three-dimensional hybrid simulation code (fluid electrons, kinetic ions). When the obstacle’s intrinsic magnetic moment is sufficiently strong, the interaction region develops signs of magnetospheric structures. On the one hand, an area from which the solar wind is excluded forms downstream of the obstacle. On the other hand, the interaction region is surrounded by a boundary layer which indicates the presence of a bow shock. By analyzing the trajectories of individual ions, it is demonstrated that kinetic effects have global consequences for the structure of the interaction region.


Introduction
When the Galileo spacecraft flew by the asteroid Gaspra in 1991, the magnetometer detected significant perturbations in the direction of the interplanetary magnetic field which have been interpreted as a result of solar wind interaction with the obstacle's intrinsic magnetic field (Kivelson et al., 1993). In the following years, various studies on the subject of solar wind interaction with an asteroid have been carried out. Baumgärtel et al. (1994) examined the plasma environment of a magnetized asteroid by using a twodimensional Hall-MHD code. In these simulations, the asteroid's intrinsic magnetic field has been represented by a dipole. By chosing an adequate orientation of the dipole moment with respect to the direction of the interplanetary magnetic field, the authors were able to generate a plasma wake which qualitatively resembles the structures observed near Gaspra. Based on theoretical studies by Greenstadt (1971), Gurnett (1995a;1995b) analyzed the propagation character-Correspondence to: S. Simon (sven.simon@tu-bs.de) istics of the whistler waves originating from the interaction of a weakly magnetized asteroid with the solar wind. Further studies trying to identify the general characteristics of the interaction have been conducted by Wang and Kivelson (1996) who used one-dimensional full particle simulations as well as a Hall-MHD code to study the wave patterns in the immediate vicinity of an asteroid. A linear theory of solar wind interaction with asteroids has been developed by Baumgärtel et al. (1997). This study concentrates on the so-called sub-magnetospheric regime in whichnthe asteroid's intrinsic magnetic field is too weak to establish either a cavity, from which the solar wind is excluded, or a magnetic tail.
However, studying details of the interaction between the collisionless solar wind flow and an obstacle whose characteristic spatial scales are comparable to ion scale lengths (ion inertial length, ion gyroradius) is beyond the scope of a magnetohydrodynamic model. In order to include ion kinetic effects, one has to use a hybrid scheme treating the electrons as a fluid, whereas a completely kinetic description is used for the ions. The fist who analyzed the general characteristics of an asteroid's plasma environment by using an electromagnetic hybrid code were Omidi et al. (2002). This study clearly illustrates the importance of ion kinetic effects for the structure of the interaction region. Based on this work, Blanco-Cano et al. (2003) interpreted the magnetic field signatures observed near asteroids Gaspra and Ida. Further hybrid simulation studies performed by Blanco-Cano et al. (2004) and Omidi et al. (2004) show that if the obstacle's intrinsic magnetic moment is sufficiently strong, a complete mini-version of an Earth-like magnetosphere is formed around the obstacle.
The hybrid simulation studies mentioned above have been carried out in only two spatial dimensions. Specifically, Omidi et al. as well as Blanco-Cano et al. are using a twodimensional code retaining all three components of the electromagnetic fields. Such a model can describe neither particle deflection nor wave propagation in the third spatial dimension. Hence, one cannot expect such a model to cover the entire structure of the interaction region. For this reason, we have developed a three-dimensional hybrid code to study the interaction of a weakly magnetized obstacle with the solar wind. The purpose of this work is to identify the general charcteristics of the interaction region depending on the relative orientation of the solar wind magnetic field and the obstacle's intrinsic magnetic moment. Emphasis is placed on the analysis of the influence of ion kinetic effects on the structure of the interaction region.
After discussing the major features of the hybrid code that has been used for the simulations in Sect. 2, the results for two different geometries are presented in Sect. 3. The discussion includes an illustration of individual ion trajectories in the vicinity of the obstacle. It is shown the substructure of the boundary layer surrounding the interaction region as well as the plasma structures in the immediate vicinity of the obstacle cannot be fully understood without taking these kinetic effects into account.

Model description
The plasma environment of a magnetized asteroid has been studied by using a three-dimensional hybrid simulation code. Since the basic principles of the code used here have been discussed extensively by Bagdonat andMotschmann (2001, 2002a, b), as well as by Bagdonat (2004), only a brief overview will be given. The hybrid model treats the electrons as a mass-less, charge-neutralizing (n e =n i ) fluid, whereas the ions are treated as individual particles according to the equations of motion In these expressions, x s and v s denote the position and the velocity of an individual particle, respectively. The electromagnetic field quantities E and B are obtained from and where u e and P e denote the mean electron velocity and the electron pressure, respectively. The electron fluid is assumed to be adiabatic, i.e. P e = P e,0 n e n e,0 κ , Due to the thermodynamic coupling in a collision-free plasma being realized by the electromagnetic fields and hence being effective in only two spatial dimensions,a value of κ=2 has been chosen for the adiabatic exponent (Bößwetter et al., 2004). The asteroid's intrinsic magnetic field is described by a magnetic dipole, where r is the radial distance to the center of the asteroid and M is the obstacle's intrinsic magnetic moment. In general, the hybrid code is also capable of handling an additional resistive term in the electromagnetic field equations, corresponding to the incorporation of a finite magnetic diffusion into the model (Bagdonat, 2004). However, for the simulations described in this paper, the resistive term has been neglected, i.e. the plasma conductivity is assumed to be infinite. Two different scenarios have been examined; the simulation frames are displayed in Figs. 1 and 3, respectively. In the first simulation, the solar wind magnetic field is orientated perpendicular to the obstacle's magnetic moment, whereas in the second simulation geometry, both vectors are antiparallel. The basic numerical parameters are the same for both simulations: The magnetic field strength and the mean plasma density of the ambient solar wind are B SW =2.7 nT and n 0 =1.4 · 10 6 m −3 , respectively. The plasma enters the simulation box with an Alfvénic Mach number of M A =8, i.e. the velocity of the undisturbed plasma is The initial flow velocity points in (+x) direction. Inflow boundary conditions have been chosen only at the left face of the cubic simulation domain, i.e. particles are continuously injected from this wall into the simulation box. At the other five walls, the plasma can leave the simulation domain (outflow boundaries). The obstacle's magnetic moment is given by M=9.5 · 10 18 Am 2 . Hence, the distance D p from the center of the obstacle to the position where the solar wind ram pressure is compensated by the obstacle's intrinsic magnetic pressure is given by D P ≈16.4c/ω p,i , where ω p,i is the ion plasma frequency in the undisturbed solar wind. This parameter has shown to be adequate to characterize the global structure of the interaction region as a function of the obstacle's magnetic moment (Omidi et al., 2002Blanco-Cano et al., 2003, 2004. In our simulations, the plasma betas of electrons and ions are supposed to be equal and are set to β e =β i =0.5. The spatial discretization is realized by using an equidistant Cartesian grid with 100×100×100 cells and a step size of g =1.8c/ω p,i in each direction. At the beginning of a simulation, each cell is filled with 10 ion macroparticles. The radius of the obstacle is r g =7c/ω p,i . Particles hitting the surface of the obstacle are removed from the simulation representing absorption by the asteroid. However, the mean particle density at the grid points inside the obstacle has to be kept on a constant, non-vanishing value. Due to the electric field depending on the pressure gradient, setting n e =0 at the grid nodes inside the obstacle would cause a strong electric field at the grid points immediately above the surface which is extremely critical for numerical stability. In order to obtain a completely self-consistent solution, the obstacle's finite conductivity would have to be incorporated into the simulation model. Hence, the entire set of Maxwell's equations would have to be solved at the grid points inside the asteroid. However, in order to restrict the simulation code's numericl complexity, the inner region is kept out of the simulation domain, i.e. the electromagnetic field quantities are not computed at the grid points inside the obstacle. This treatment is also mandatory since the strong magnetic field gradients near the center of the obstacle could compromise numerical stability. To sum up the basic idea of these boundary conditions, the electromagnetic field quantities as well as the moments of the particle distribution function are fixed on constant, timeindependent values in the inner region. Of course, the divergence of the magnetic field vanishes at the grid points inside as well as outside the obstacle. The set of hybrid equations is solved by means of a Particle-in-Cell technique as described by Bagdonat (2004). This scheme is based on a four-step computational cycle: In the first step, the moments of the distribution function (particle densities and currents) are calculated for each grid point by means of a weighting technique. In the second step, the electromagnetic field equations are solved by using these values for the particle densities and currents, yielding an updated electric and magnetic field vector for each grid point. In the third step, the electromagnetic field quantities are interpolated from the grid points to the individual particle positions. Finally, given the electromagnetic forces acting on the individual particles, the positions and velocities can be updated. Each computational cycle represents a certain interval t of real time. In our simulations, a time step of t=0.05 −1 has been chosen, where =eB SW /m p is the proton gyrofrequency. This value fulfills the Courant Friedrichs Lewy condition for numerical stability. The simulations have been carried out sufficiently long to reach a quasi-stationary state. The total simulation time of 120 −1 corresponds to a duration in which the undisturbed solar wind protons would cross the entire simulation domain five to six times. In order to guarantee the simulation's numerical stability, a 27 grid point smoothing procedure has been used.

Simulation results
We consider two different situations: Perpendicular and antiparallel orientation of M and B SW .

Perpendicular orientation of magnetic moment and interplanetary magnetic field
The simulation geometry is shown in Fig. 1, whereas Fig. 2 displays the simulation results in the (x, y) plane which we call the equatorial plane of the obstacle. In the vicinity of the dipole, the interplanetary magnetic field exhibits a strong draping pattern. Besides, a region where the mean plasma density lies significantly below the upstream value occurs downstream of the obstacle. Due to the solar wind plasma being excluded from this cavity, this region indicates the beginning formation of magnetospheric structures. As can be seen from Fig. 2b, the interaction region is surrounded by a boundary layer denoting a sudden increase of plasma density. A similar structure can be seen in the electric field as well as in the mean particle velocity. Because at the dayside of the obstacle, the plasma flow is penetrating the boundary layer almost perpendicularly, a strong deceleration of the plasma occurs. In contrast to this, the plasma flows almost tangential to the flanks of the boundary layer, and hence, only a slight reduction of plasma velocity can be seen in these regions. However, the crucial question is whether these simulation results indicate the presence of a complete mini version of a standard magnetosphere, especially the formation of a bow shock. On the one hand, the density signatures occuring in the downstream region indicate the formation of a dipole-dominated region which is not accessible for the solar wind. The presence of a shock-like structure is consistent with the formation of a boundary layer denoting a discontinuity of plasma density as well as a sudden deceleration of the plasma. On the other hand, the magnetic field signature exhibits no such boundary layer, i.e. the results show only a slight increase of magnetic field strength instead of a clear shock formation. Hence, one can only state that the interaction region exhibits first signs of magnetospheric structures, but the global structure of the interaction region does not coincide with a complete magnetosphere. Blanco-Cano et al. (2004) as well as Omidi et al. (2004) have been conducting numerical studies on the subject of solar wind interaction with weakly magnetized obstacles by using a two-dimensional hybrid code. The authors have proposed that the distance D P (in proton inertial lengths) to the center of the obstacle where the solar wind ram pressure is compensated by the obstacle's intrinsic magnetic pressure is an adequate parameter for characterizing the global structure of the interaction region as a function of the obstacle's intrinsic magnetic moment. In our simulations, this parameter is of the order of D P ≈16.4c/ω p,i . Although a direct comparison of 2-D and 3-D simulations is not possible, both studies correspond at least in so far that Blanco-Cano et al. (2004) have found out that only if D P >20c/ω p,i , the obstacle possesses an Earth-like magnetosphere. This is consistent with our results showing that the interaction region exhibits only certain sub-magnetospheric structures, but not a fully-developed magnetosphere.

Antiparallel orientation of magnetic moment and interplanetary magnetic field
In our second simulation geometry, the asteroid's intrinsic magnetic moment is orientated antiparallel to the solar wind magnetic field (cf. Fig. 3). The simulations clearly show that ion kinetic effects are of major importance for the structure of the interaction region and hence, using a fluid model instead of a hybrid code is not adequate when the scales of the interaction region are comparable to ion scale lengths. The results of a cut through the (x, y) plane (equatorial plane of the obstacle) are shown in Fig. 4. Certain features of the interaction region resemble the results of our first simulation: the draping of the interplanetary magnetic field around the obstacle, the formation of a cavity with very low plasma density from which the solar wind is excluded. Nevertheless, one of the most remarkable features of the interaction region is the substructure of the boundary layer that can be seen in Fig. 4b. The outer layer separating the interaction region from the undisturbed solar wind is followed by two other sharp curves denoting sudden increases of plasma density. Besides, the obstacle is surrounded by a belt denoting a significant increase of plasma density. The density plots for the (x, z) plane are shown in Figs. 5 and 6, respectively. The density plot at the r.h.s. of Fig. 5 shows the ring intersecting a plane that includes the vector of magnetic moment, whereas its structure in a plane perpendicular to the magnetic moment is shown at the l.h.s. of Fig. 5. These results indicate that this ring is orientated perpendicular to the obstacle's magnetic moment. The figure shows the projections of the particle trajectories on the planes z=− g /2 and y=− g /2 as well as the normalized plasma density n/n 0 in these planes for the quasi-stationary state of the simulation. It is important to know that both figures display the same particle trajectories, starting with z=0 at the left hand boundary of the simulation domain. It is illustrated that the scales of the interaction region are comparable to the characteristic length scales of ion motion. It is obvious that the substructure of the boundary layer -several sharp curves denoting an increase of plasma density -emerges from the finite gyroradius of individual solar wind particles.
The particle trajectories shown in Figs. 5 and 6 clearly illustrate the importance of ion kinetic effects. From several points of the intersection line of the (x, y) plane and the left face of the simulation domain, several particles have been followed through the simulation box. The figures show the projections of their trajectories on the planes z=− g /2 and y=− g /2, which almost coincide with the (x, y) and (x, z) plane of the coordinate system, respectively. As shown in Fig. 5, the formation of the boundary layer's substructure is correlated to the gyration of the individual ions. One one hand, as can be seen in Fig. 5a, the curves denoting an increase of plasma density are almost tangential to the gyration trajectories. On the other hand, Fig. 5b shows that the turning points of the gyration trajectories are also located in these regions with increased plasma density. A similar substructure of the bow shock (shocklet structure) has been observed by Shimazu (2001) who performed 3-D hybrid simulation studies on the subject of solar wind interaction with the iono-spheres of unmagnetized planets. Of course, solar wind interaction with an ionosphere cannot be compared directly to the interaction mechanism studied in our simulations. However, in agreement with our results, the author proposed that these shocklet structures originate from the finite gyroradius of the solar wind protons. Similar shocklet structures have been described near comets (Bagdonat and Motschmann, 2002b) as well as in the bow shock region of Mars (Bößwetter et al., 2004).
Finally, the particle trajectories shown in Fig. 6 allow to explain the formation of the belt surrounding the obstacle. The belt is generated by temporarily trapped particles. Ions entering the vicinity of the obstacle can be trapped in the magnetic dipole field, but the dipole is not strong enough to capture them permanently. The particle trajectories clearly indicate that the characteristic scales of ion motion are comparable to the size of the obstacle. The belt originates from protons being influenced by the the dipole field lines. The background color plots display the normalized plasma density in the cutting planes, whereas the black lines denote the projections of ten individual particle trajectories on the cutting planes. The belt surrounding the obstacle in the equatorial plane is generated by individual particles that are influenced by the magnetic dipole field. However, because the dipole is not strong enough to capture them permanently, some particles are also able to leave the vicinity of the obstacle. The particles are gyrating around the dipole field lines that are perpendicular to the equatorial plane of the obstacle.
It should be noted that the formation of such a belt did not occur in our first simulation with perpendicular orientation of M and B SW . In this situation, the proton gyration plane in the solar wind is almost perpendicular to the gyration plane in the dipole field. Thus, the particles have to alter their gyration plane when entering the vicinity of the obstacle. For this reason, the particle trajectories are considerable more complex than in the case of antiparallel orientation where the gyration planes in the dipole field and in the solar wind are almost parallel.

Summary
The plasma environment of a magnetized asteroid has been studied by using a 3-D hybrid simulation code. The obstacle's intrinsic magnetic field was assumed to be dipolar. Two different situations have been studied: On the one hand, the obstacle's intrinsic magnetic moment M was assumed to be perpendicular to the undisturbed interplanetary magnetic field B SW . On the other hand, antiparallel directions haven been chosen for M and B SW in a second simulation scenario. The studies refer to the sub-magnetospheric regime where the intrinsic magnetic field is not strong enough to form an Earth-like magnetosphere around the obstacle. However, the simulations have shown that the interaction region exhibits first signs of magnetosphere-like structures: A cavity where the plasma density lies significantly below the upstream value is formed in the downstream region. This region is dominated by the obstacle's intrinsic magnetic field and it is not accessible for the solar wind. The interaction region is surrounded by a sharp boundary layer denoting a significant decrease of plasma velocity at the dayside of the obstacle, but only a slight decrease at the flanks. These structures can be interpreted as the beginning formation of a bow shock. Nevertheless, the simulation results show no clear separation of bow shock and magnetopause.
The simulation results clearly demonstrate the importance of ion kinetic effects. In the analyzed upstream parameter case, the results show that the characteristic length scales of the interaction region are comparable to the scales of ion motion and hence, kinetic effects cannot be neglected. The substructure of the boundary layer surrounding the interaction region is a result of finite gyroradius effects and can only be understood by analyzing specific particle trajectories. Besides, gyration of individual ions results in the formation of a belt of temporarily trapped particles surrounding the obstacle in the equatorial plane. Since these effects play a central role for the structure of the interaction region, using a hybrid code instead of a (multi) fluid model allows to analyze significantly more details. Moreover, as most of the structures described above are not confined to a single plane, the results presented in this paper clearly emphasize the necessity of using a three-dimensional simulation code.