The Magnetohydrodynamics of Shock-Cloud Interaction in Three Dimensions

Although most previous studies of shock-cloud interaction used a density discontinuity at the cloud surface (including the 3D hydrodynamical calculations presented in XS), for small clouds it is more realistic to represent the density distribution with a smooth profile at the boundaries. Following previous authors (Kornreich & Scalo 2000; Nakamura et al. 2006), the density profile is given by

where $\rho_{0}=1$ is the density of the ambient medium, $\rho_{c}=10$ is the density at the center of the cloud, and the core radius $r_{c}=0.62$. Thus, the cloud is overdense compared to its surroundings by a factor $\chi=\rho_{c}/\rho_{0}=10$. Larger values of $n$ give steeper density profiles at the surface; we adopt $n=8$ for all calculations presented in this paper.

The density distribution defined in equation 7 extends to infinity. Thus, we must adopt an arbitrary definition for the boundary of the cloud. In this paper, we define the boundary at $r=r_{b}$, where $\rho(r_{b})=1.002\rho_{0}$. For the parameters given above, this gives $r_{b}=1.77$. Thus, the total column density of the cloud in these units is

The value of $\Sigma_{cl}$ is important to the dynamics because we expect significant evolution to occur once the cloud has swept up a comparable column of the postshock gas.

Initially the cloud is in pressure equilibrium with its surroundings. We set the gas pressure $P=1$ everywhere, giving a sound speed in the ambient gas of $C_{s}=\sqrt{5/3}$.

We study three different initial directions and strengths for the initial magnetic field. Depending on the orientation of the magnetic field in the preshock gas with respect to the unit normal perpendicular to the shock front, our simulations use either parallel, perpendicular, or oblique magnetic fields. The latter use fields inclined at $45^{\circ}$ to the shock normal. The strength of the magnetic field can be defined by a plasma $\beta$ calculated in the preshock medium. For each field geometry, we study fields with a strength given by $\beta=10$, 1 and 0.5

The final quantity required to specify the initial conditions is the sonic Mach number of the incident shock, $M_{s}$. We study the evolution of strong shocks with $M_{s}=10$ in this paper. Given $M_{s}$, and the magnetic field orientation and strength, it is straightforward to calculate the properties of the postshock gas from the shock jump conditions. The evolution of stronger shocks should be the same as presented here with the appropriate dimensionless scaling (see M94). The evolution of much weaker shocks may be different, but is not considered further here.

All of the numerical simulations presented in this paper use Athena, a newly developed Godunov scheme for astrophysical MHD. Details of the algorithms implemented in Athena are given in Gardiner & Stone (2005; 2007). These algorithms have been extended to integrate equation 2 along with the usual equations of ideal MHD.

The computational domain covers the region $-3.66\leq x\leq 6.34$, $-5\leq y\leq 5$, and $-2.5\leq z\leq 2.5$ for all of the simulations. The cloud is initially centered on the origin. The shock is located at $x=-2.66$, and propagates along the $x-$axis. The initial orientation of the magnetic field in the preshock gas is along the $x-$ axis for parallel shocks, along the $y-$axis for perpendicular shocks, and at $45^{\circ}$ in the $x-y$ plane for oblique shocks.

Once the shock impacts the cloud, the cloud is rapidly accelerated, and would leave the computational domain through the right-hand boundary in the $x-$direction in only a few dynamical times. To keep the shocked cloud centered in the grid, we compute the mass-averaged cloud velocity

every time-steps, and integrate the MHD equations in a frame of reference moving at $\langle v_{x}\rangle$. This keeps the center of the cloud at rest with respect to the grid. Essentially, our domain moves with the cloud, considerably reducing the computational volume needed to follow the evolution. We report all our results with respect to the initial frame of reference which is at rest with respect to the preshock medium.

All of the simulations presented in §4 use a grid of $N_{x}\times N_{y}\times N_{z}=680\times 680\times 340$, or about 160 million zones. This is thirty times larger than the 3D hydrodynamic shock-cloud interaction simulation first presented in Stone & Norman (1992). This resolution results in roughly 120 grid points covering the cloud radius $r_{b}$. To investigate the convergence of our results (§3), we have also computed the evolution on a grid with one half the resolution, or $N_{x}\times N_{y}\times N_{z}=340\times 340\times 170$ zones.

The time taken for the incident shock to propagate across the interior of the cloud defines the “cloud-crushing time” $t_{cc}$, a characteristic time scale for the shock-cloud interaction problem (KMC). For a cloud with the density profile equation 7, the cloud-crushing time is defined as

where $v_{s}$ is the shock velocity in the cloud interior. We will report all of our results with time given in units of $t_{cc}$. Generally, our simulations are run until $t\sim 20t_{cc}$.

For the cloud parameters given in §2.1, $t_{cc}\sim 0.15$. To convert to typical values of small clouds in the ISM, if we assume $r_{c}=4000~{\rm AU}$, $C_{s}=10^{5}~{\rm cm~s^{-1}}$ for $\chi=10$ and $M_{s}=10$, then $t_{cc}\sim 6,000$ years. The physical size of the computational domain in this case is $\sim 0.3~\rm{pc}\times 0.3~\rm{pc}\times 0.15~\rm{pc}$.

There are a variety of diagnostic quantities we will use to analyze our results. We define the cloud-mass-weighted average of any quantity $f$ as

For example, the average cloud velocity defined in equation 9 is just the cloud-mass-weighted average of the $x-$component of velocity.

To follow changes in the shape of the cloud, the mass-weighted moment along the $x-$axis is used (MacLow et al. 1994)

The expressions for the moments along the $y-$ and $z-$axis, $b$ and $c$ respectively, are similar. For the parallel shock simulations, we expect symmetry about the $x-$axis, so that $b=c$. This can be used as a code check.

The acceleration of the cloud can be studied using the mass-weighted velocity along each axis, e.g. $\langle v_{x}\rangle$ (already given by equation 9), $\langle v_{y}\rangle$ and $\langle v_{z}\rangle$.

In order to trace mixing between the cloud and ambient material, we use the mixing fraction $f_{mix}\equiv M_{mix}/M_{cl}$, where following XS we define $M_{mix}$ to be the total mass of cells in which the color variable $0.1\leq C\leq 0.9$. Such cells are those in which the cloud and ambient gas are well mixed.

Finally, the RMS of velocity along each orthogonal axis is useful to understand mixing and the generation of turbulence (Nakamura et al. 2006). The root-mean-square (RMS) of the $x-$component of the velocity is

with similar expressions for the RMS along the $y-$ and $z-$axes, $\delta v_{y}$ and $\delta v_{z}$ respectively. We also expect symmetry in $\delta v_{y}$ and $\delta v_{z}$ for the parallel shock case.

Figure 2 presents volumetric renderings of the cloud density $\rho C$ at four times during the evolution of the weak and strong parallel fields, runs PA10 and PA05 respectively. The image is computed viewing along the direction given by the unit normal $\hat{\bf n}=(0.5,0.866,0)$, that is at $30^{\circ}$ to the $y-$axis, and is centered on the origin. We expect the shock to have propagated across the diameter of the cloud by $t=2t_{cc}$. Thus, at $t=3.76t_{cc}$ shown in the first panel in figure 2, both the incident and transmitted shocks are well past the cloud, the initial compression phase of the cloud has ended, and the shocked cloud material is now re-expanding. Note there is little difference in the structure of the cloud between the weak and strong field cases at this time, since the mechanical energy of the shock rather than the magnetic energy of the cloud dominates. The only significant difference at this time is that the column of compressed material behind the head of the cloud is broader for $\beta=0.5$, since the strong field resists compression towards the axis by the incident shock as it is refracted around the cloud. However, at later times the structure of the cloud in the weak and strong field cases begins to show significant differences. For the weak field, the tip of the cloud is shredded into a complex network of filaments due to shear instabilities associated with the vorticity deposited at the cloud surface, similar to the hydrodynamic problem in 3D (SN, XS). On the other hand, for a strong field the cloud resists fragmentation. Shear instabilities cause wrinkles, fingers, and clumping at the cloud surface, however to a much smaller degree in comparison to the weak field. Direct comparison of figure 2 with, e.g. figure 2 in XS, shows that even for the weak field case, the destruction of the cloud is decreased in comparison to hydrodynamics.

Previous studies (Hawley & Zabusky 1989) have shown that the vortensity $\omega^{2}\equiv(\nabla\times{\bf v})^{2}$ is a useful diagnostic for the hydrodynamic shock-cloud interaction problem. Figure 3 plots volumetric renderings of $\omega^{2}$ for the weak and strong parallel fields at the same four times and using the same viewing angle, as shown in figure 2. Again, the evolution with a weak field is similar, but not identical, to the hydrodynamic case (SN, XS). At the earliest time, a thin layer of vorticity has been deposited at the cloud surface, while the reflection of the incident shock on-axis behind the cloud has generated a very strong vortex ring downstream of the cloud. This ring is completely absent in the strong field case. The thin sheet of vorticity is associated with this ring marks the location of the incident shock, which has been curved by passage over the cloud. (By comparing figures 2 and 3, note the disk of cloud material behind the head of the cloud is downstream of the location of the incident shock, marked by this thin sheet of vorticity.) At later times, KH instability in the sheet of vorticity at the head of the cloud results in the formation of a complex network of tangled vortex filaments that produce the complex density structures at the corresponding times in figure 2. In contrast, for the strong field case the vortex sheet remains largely intact, indicating the KH instability is suppressed.

The stabilizing effect of a strong magnetic field for parallel shocks is evident in the time evolution of each component of the velocity RMS. Figure 4 shows the evolution of the parallel and transverse components, $\delta v_{x}$ and $\delta v_{\perp}\equiv(\delta v_{y}^{2}+\delta v_{z}^{2})^{1/2}$ respectively, normalized by the preshock sound speed for runs PA10, PA1, and PA05. The largest values of $\delta v_{x}$ occurs very early, at about one $t_{cc}$. At this time, the transmitted shock is only half way through the cloud, giving a large dispersion between the shocked and unshocked cloud material. The largest values of the transverse component of the dispersion occurs slightly later, at slightly less than $2t_{cc}$. This corresponds to the period of maximum transverse compression by the transmitted shock. Thereafter, the longitudinal component decreases uniformally, with only a small difference between different values of $\beta$ (with the largest $\beta$ having the largest dispersion). The transverse component also decreases to a constant value, however this value depends strongly on the field strength. For $\beta=10$, $\delta v_{\perp}/C_{s}\approx 0.5$, whereas for $\beta=0.5$, $\delta v_{\perp}/C_{s}\approx 0.2$. The longitudinal component of the velocity RMS is indicate of uniform acceleration of the cloud, whereas the transverse component is more indicative of postshock turbulence. This turbulence is subsonic with respect to the preshock gas (and therefore highly subsonic with respect to the postshock gas). The large decrease in the transverse velocity RMS with decreasing $\beta$ is indicative that strong fields suppress turbulence.

We expect the evolution of the velocity RMS and magnetic field to be tightly coupled. Figure 5 plots the time evolution of the energy associated with amplification of both the longitudinal field $(B_{x}^{2}-B_{x,0}^{2})/B_{x,0}^{2}$, and the transverse field $(B_{y}^{2}+B_{z}^{2})/B_{x,0}^{2}$ in a simulation box (where $B_{x,0}^{2}$ is the energy associated with the initial field), for each value of $\beta$. The weak field $\beta=10$ shows the largest amplifications in both plots. Moreover the longitudinal field is amplified significantly more than the transverse. This reflects the importance of compression of the longitudinal field on-axis to form a flux rope behind the cloud (M94) in comparison to turbulent amplification of the transverse field. The field geometries produced by the shock-cloud interaction will be discussed further in §5.

Finally, the relative changes in the time evolution of the mixed fraction $f_{mix}$ for different values of $\beta$ is shown in figure 6. The largest mixing occurs for $\beta=10$, while both of the stronger field cases ($\beta=1$ and 0.5) are comparable. Since $f_{mix}$ is not converged (as discussed in §3), we assign no relevance to the particular values of this parameter in any run. However, the decrease in the relative values with increasing field strength is meaningful, and reflects the decrease in fragmentation and mixing of the cloud evident in figure 2.

Figure 7 shows 3D renderings of the cloud in runs PE10, PE1, and PE05 at three times during the evolution. We have found that it is difficult to interpret the structure of the cloud in perpendicular shock runs with strong fields from volumetric renderings as shown in figure 2 (similarly, we find plotting vortensity is not useful in this case). Instead, figure 7 shows a combination of the $\rho C=6$ isosurface, a slice of $\rho C$ plotted on the horizontal plane $y=0$, and magnetic field vectors plotted on the vertical plane $z=0$. The viewing angle and focus point in each panel is identical to that in figure 2.

At the earliest time, $t=3.76t_{cc}$, the structure of the cloud is similar in all three runs (and remarkably similar to the parallel shock runs shown in figure 2). There is a slight asymmetry in the cloud shape along the directions perpendicular and parallel to the field, but this is clearly evident only in the $\beta=0.5$ case. At later times, however, the cloud structure with different field strengths is significantly different.

For the weak ($\beta=10$) field, the cloud is fragmented, with a clear asymmetry in the density structures between the $x-y$ and $x-z$ planes. Since initially the magnetic field lies in the $x-y$ plane, vortical motion in that plane twists the field, and thus the field resists such motion. However, in the $x-z$ plane, vortical motion simply interchanges field lines, and so the field does not affect such motion. Thus, strong vortices form in the $x-z$ plane, with cloud material entrained at the center of the vortices. The magnetic field vectors shown in the $x-y$ plane demonstrate that the field lines drape over the surface of the cloud as expected, and that the field is strongly amplified (as shown by the color and length of the arrows) as it is stretched downstream of the cloud.

For both of the stronger ($\beta=1$ and 0.5) fields, the cloud no longer becomes fragmented at late times, but instead resembles a sheet aligned with the field. In these two cases, stretching of the field over the cloud is evident at early times ($t=3.76t_{cc}$ and $t=7.97t_{cc}$). However, at late times the field is vertical everywhere. This occurs because acceleration of the cloud to the postshock velocity is rapid (e.g. for $\beta=0.5$, the postshock density and $v_{x}$ is about 3.54 and 9.58 respectively; this means the cloud sweeps up its own column density by $t\approx 0.5t_{cc}$, and by a few $t_{cc}$ the cloud has accelerated to the postshock flow speed). Once the cloud is at rest with respect to the postshock gas, the longitudinal perturbations in the field propagate away from the cloud as shear Alfvén waves, restoring the field to purely vertical near the cloud. These waves are clearly evident at $t=7.97t_{cc}$ for the $\beta=1$ case. Thus, unless the column density of the cloud is significantly higher (so that the cloud takes much longer to accelerate), there will not be strong longitudinal fields near the cloud surface at late times. Instead, strong longitudinal fields will be associated with the shear Alfvén waves far from the cloud surface.

Previous 3D studies of wind-cloud interaction using perpendicular fields (Gregori et al. 2000) discussed the importance of field compression at the leading surface of the cloud, forming a “magnetic bumper” which accelerates fragmentation. At $t=3.76t_{cc}$ there are small amplitude vertical striations in the leading surface of the cloud in both the $\beta=1$ and 0.5 runs, associated with MHD RT instabilities in the compressed field upstream of the cloud. The growth rate of the MHD RT instability is given by (Chandrasekhar 1961)

where $a$ is the effective acceleration of the cloud, $\rho_{1}$ and $\rho_{2}$ are the densities of the ambient gas and cloud respectively, $B$ is the magnetic field strength, and $k$ is the wave number of the perturbation. By measuring from the simulations the effective acceleration of the cloud, the densities across the surface of the cloud, the magnetic field strength, and the wavelength of the perturbations shown in figure 7, we find the growth time from the above equation is about 10 times shorter than the acceleration time. Thus, we expect the MHD RT instability to grow rapidly. By $t=7.97t_{cc}$ these striations have grown to large amplitudes, fragmenting the cloud into vertical columns. Even the $\beta=10$ begins to show the striations at this time. However, as discussed in §5, at late times the strongest amplified field is still associated with filaments downstream of the cloud, rather than the compressed field upstream.

The anisotropic shape of the cloud produced by a perpendicular field is best illustrated by the time evolution of the ratio of the moments of the cloud axes, $b/a$ and $c/a$ respectively. Figure 8 plots these moments for all three field strengths, runs PE10, PE1, and PE05 respectively. Initially both ratios show a peak between $1-2t_{cc}$ associated with the decrease in $a$ as the cloud is crushed along the direction of propagation of the shock. Thereafter the ratios decrease to constant values. Note for strong fields $b/a>c/a$, indicating the cloud is extended along the direction of the background field. For the weak field ($\beta=10$) the ratios are nearly equal and very similar to the values for the parallel shock shown in figure 1.

The evolution of the magnetic energy associated with the $x-$ and $z-$components of the magnetic field for runs PE10, PE1, and PE05 is shown in figure 9. The evolution of $B_{y}^{2}$ is dominated by linear growth due to compression as the shock propagates across the grid, and therefore is not shown. Once again the weakest field shows the largest growth. Moreover, the energy associated with the longitudinal field $B_{x}^{2}$ shows much larger growth than that associated with the transverse field $B_{z}^{2}$. This latter asymmetry simply reflects the fact the field is initially orientated in the $x-y$ plane, thus $B_{x}$ can be created by bending the field over the cloud surface (vortical motion in the $x-y$ plane), while $B_{z}$ cannot. The large decrease in both $B_{y}^{2}$ and $B_{z}^{2}$ at late times for $\beta=1$ and 0.5 is caused by the propagation of the shear Alfvén waves through the transverse boundaries of the grid. For example, it is evident in figure 9 that the largest values of $B_{x}$ are associated with these waves, and once they leave the grid the field becomes almost vertical everywhere, so $B_{x}$ drops to nearly zero. Note that no decrease is seen for the weak field $\beta=10$ run. In this case, field amplification is caused by winding up in vortices, and in this case the tension in the amplified field is never large enough to cause unwinding.

Finally, figure 10 plots $f_{mix}$ for runs PE10, PE1, and PE05. Once again, a strong decrease in the relative amounts of mixing is observed as the field strength increases, a factor of two between the $\beta=10$ and $\beta=0.5$ cases. While the absolute values of $f_{mix}$ are not converged, the relative decrease is real, and reflects the decrease in turbulence and fragmentation with stronger fields.

Figure 11 plots the $\rho C=6$ isosurface, a slice of $\rho C$ plotted on the horizontal plane $y=0$, and magnetic field vectors plotted on the vertical plane $z=0$ for the oblique field simulations, runs OB10, OB1, and OB05 respectively. The viewing angle and focus point for the $\beta=10$ run OB10 is identical to that used in previous figures, while for the $\beta=1$ and 0.5 plots the plots are centered on the point $(x,y,z)=(0,0,-1)$.

Once again, we find that at early times ($t=3.76t_{cc}$), there is little difference in the structure of the cloud with different field strengths (or by comparison to figures 2 and 7, for different field geometries). However, at later times, the cloud in the weak field run becomes highly fragmented, whereas with stronger oblique fields the cloud is more sheet-like. As with perpendicular shocks, small amplitude striations are evident in the leading surface of the cloud with strong fields; they indicate the presence of the MHD RT instability associated with compressed upstream field. At later times, the instability grows to large amplitude, and fragments the cloud into vertical columns. However, the most striking aspect of the evolution is that the cloud becomes significantly displaced from the $y=0$ plane with strong fields. For oblique shocks, $v_{y}$ is no longer zero in the postshock gas, instead it varies from -0.022 ($\beta=10$) to -0.434 ($\beta=0.5$). Thus, the cloud is simply dragged in the negative $y-$direction by the postshock flow, while retaining a sheet-like structure aligned with the field.

Fragmentation of the cloud in the weak field ($\beta=10$) oblique shock case is more complex than for perpendicular shocks. In addition to the suppression of vortical motion in the $x-y$ plane, oblique fields also affect vortical motion in the $x-z$ plane. Thus, while the $\beta=10$ perpendicular shock tended to fragment the cloud into ribbons and sheets in the $x-z$ plane, the oblique shock tends to fragment the cloud into isolated structures.

For the strong field runs, magnetic tension tends to erase longitudinal perturbations, and near the cloud the field is aligned with the direction in the postshock gas (about $15^{\circ}$ from the $y-$axis for both $\beta=1$ and $\beta=0.5$). Strong longitudinal perturbations in the field propagate away from the cloud as shear Alfvén waves. Figure 12 plots the time evolution of the magnetic energy associated with the $x-$ and $z-$components of the perturbed magnetic field for all three oblique field runs. Again, we do not show $B_{y}^{2}$ since it is dominated by shock compression. The evolution of the magnetic energies is similar to the perpendicular shock case as shown in figure 9. The strongest amplification is for the weakest fields, and the energy associated with the longitudinal field $B_{x}^{2}$ is much larger than that associated with the transverse field $B_{z}^{2}$.

Other diagnostics of the cloud evolution for oblique shocks are very similar to that shown for perpendicular shocks, and so are not shown here. For example, the evolution of the shape of the cloud as measured by the axis ratios $b/a$ and $c/a$ is similar to that shown in figure 8; the cloud is elongated in the $y-$direction and so $b/a$ tends to be largest. In addition, the mixing fraction $f_{mix}$ evolves in a similar manner to the perpendicular field case (shown in figure 10). In particular, the relative amount of mixing decreases strongly with increased field strength, indicating the suppression of turbulent fragmentation and mixing with strong oblique fields.