CONSEQUENCES OF MAGNETIC FIELD STRUCTURE FOR HEAT TRANSPORT IN MAGNETOHYDRODYNAMICS

Interfaces between hot and cold plasmas can occur in astrophysics where understanding the rate of thermal conduction may be an important part of the astrophysical phenomenology. One example occurs in wind blown bubbles (WBBs) of evolved stars where magnetized hot supersonic outflow shock heats the cooler ambient magnetized interstellar medium (ISM). For such WBBs, there are examples where the presumed shock-heated bubble is cooler than expected if only radiative cooling is considered (Zhekov & Park 2011). A possible explanation is that heat loss through the interface of hot bubble into the cold shell via thermal conduction reduces the temperature of the hot bubble (Zhekov & Myasnikov 1998, 2000). However, the source of heat into the cold side of the interface will continuously evaporate material there and potentially induce interface instabilities and mass mixing (Stone & Zweibel 2009) that could tangle the magnetic field. Understanding the thermal conduction and its dependence on magnetic structure is important for determining the thermal properties of the plasma on either side of the interface.

A second example is the unexpected slow mass deposition rate of the cooling flows in some galaxy cores, which might be inhibited by a restricted thermal conduction (Rosner & Tucker 1989; Balbus & Reynolds 2008; Mikellides et al. 2011). In the intracluster medium (ICM), the tangled magnetic field can potentially produce a strongly anistropic thermal conductivity that may significantly influence temperature and density profiles (Chandran & Maron 2004; Maron et al. 2004; Narayan & Medvedev 2001; Mikellides et al. 2011).

For the ISM and ICM, it is usually valid to assume that the electrons are totally inhibited from moving across field lines (McCourt et al. 2011), as the electron mean free path is much greater than the electron gyroradius. The underlying assumption is that the field configuration is ideal and there are no stochastic fluctuations existing. The magnetic field structure then plays a key role in controlling the rate of thermal conduction since electrons can move freely only along the field lines. It is worth pointing out that in reality, however, the stochastic field can change the cross field diffusion even if its amplitude is small. Using the conditions given by Rechester & Rosenbluth (1978), it can be shown that the ratio of ion gyroradius and the field length scale presented in our simulation is small. This means that even a small added stochastic field is likely to make a difference in the anisotropicity of the thermal diffusion. The result is a strong thermal conductivity parallel to the field lines and a weak conductivity across the field lines.

The quantitative subtleties of how a complicated magnetic field structure affects thermal conduction raises the open question of whether there is a simple measure of field tangling that allows a practical but reasonably accurate correction to the isotropic conduction coefficient for arbitrarily tangled fields. In this context, two classes of problems can be distinguished. The first is the conduction in a medium for which forced velocity flows drive turbulence, which in turn tangles the field into a statistically steady state turbulent spectrum (Tribble 1989; Tao 1995; Maron et al. 2004). The second is the case in which the flow is laminar and the level of conduction inhibition is compared when the field starts from initial states of different levels of tangling subject to an imposed temperature difference across an interface. This second problem is the focus of the preset paper. It should be stated that the conductivity in our simulations remains that associated with the micro-physical scale throughout the evolution of our simulation. That is, our flow remains laminar so we do not have a broad turbulent spectrum of magnetic fluctuations or a corresponding increase in the effective conductivity as in Narayan & Medvedev (2001).

Using the ASTROBEAR magnetohydrodynamics code with anisotropic thermal conduction, we investigate the influence of initial magnetic structure on thermal conduction in an otherwise laminar flow. The key questions we address are: (1) does the interface become unstable? and (2) how fast is the thermal conduction across the interface compared to the unmagnetized case?

We study these questions using different initial magnetic configurations imposed on a planar thermal interface to determine how the conduction depends on the amount of field tangling across the interface.

In Section 2, we review the basic equations of MHD with anisotropic thermal conduction. In Sections 3 and 4, we provide detailed description of the simulation setup. In Sections 5 and 6, we present the simulation results and analyses. In Section 7, we discuss the simulation results in the context of the WBB cooling problem and the cooling flow problem in cores of galaxy clusters. The appendix provides more detailed information on the testing of the ASTROBEAR code.

The MHD equations with anisotropic heat conduction that we will solve are given by

$$\begin{equation} \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \bf {v}) = 0, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial (\rho \textbf {v})}{\partial t}+\nabla \cdot \left[\rho \textbf {vv}+\left(p+\frac{B^2}{8\pi }\right)\textbf {I}-\frac{\textbf {BB}}{4\pi }\right] = 0, \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial \textbf {B}}{\partial t}+\nabla \times (\textbf {v}\times \textbf {B}) = 0, \end{equation} \tag{ 3 }$$

$$\begin{equation} \frac{\partial E}{\partial t}+\nabla \cdot \left[\textbf {v}\left(E+p+\frac{B^2}{8\pi }\right)-\frac{\textbf {B}(\textbf {B}\cdot \textbf {v})}{8\pi }\right]+\nabla \cdot Q = 0, \end{equation} \tag{ 4 }$$

where ρ, $\textbf {v}$, $\textbf {B,}$ and p are the density, velocity, magnetic field, and pressure, and E denotes the total energy given by

$$\begin{equation} E = \epsilon +p\frac{\textbf {v}\cdot \textbf {v}}{2}+\frac{\textbf {B}\cdot \textbf {B}}{8\pi }, \end{equation} \tag{ 5 }$$

where the internal energy is given by

$$\begin{equation} \epsilon = \frac{p}{\gamma -1} \end{equation} \tag{ 6 }$$

and γ = 5/3. In our simulations, we will assume that the heat flux is confined to be parallel to the magnetic field lines. This assumption applies only when the ratio of electron gyroradius to field gradient scales is small. Under this assumption, the heat flux parallel to field lines can be written as

$$\begin{equation} Q = -\kappa _{\parallel }(\nabla T)_{\parallel }, \end{equation} \tag{ 7 }$$

where the subscript || indicates parallel to the magnetic field and κ_∥ is the classical Spitzer heat conductivity: κ_∥ = κ_c T^2.5, with κ_c = 2 × 10⁻¹⁸ cm s g⁻¹ K^−2.5. We take κ_∥ to be a constant throughout our simulations and so hereafter write it simply as κ.

The ASTROBEAR code uses an operator splitting method to solve these MHD equations with heat conduction. The viscosity and resistivity are ignored in our calculation so the dissipation is numerical only. The ideal MHD equations are then solved with the MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) primitive method with total variation diminishing preserving Runge–Kutta temporal interpolation. The result is then sent to the implicit linear solver utilizing the high performance preconditioners (Hypre) to solve the anisotropic heat conduction equation. The linear solver requires temporal sub-cycling technique to maintain its accuracy. It should be pointed out that the sub-cycling is implemented to improve the efficiency of simulations with extremely strong conductivity. However, in this particular study, the conductive diffusion rate is never fast enough to trigger the sub-cycling. The code runs in parallel with fixed grid domain.

Our initial set up involves hot and cold regions separated by a thin planar interface. We study how the magnetic field configuration alters the heat transfer rate between the hot and cold regions in presence of anisotropic heat conduction. We study the problem in two dimensions.

To guide subsequent interpretation of the results, we first compare two simple but illustrative limits of magnetic field orientation: (1) a uniform magnetic field aligned with the direction normal to the interface and (2) a uniform magnetic field perpendicular to the normal direction of the interface. In case (1), because the angle between the magnetic field and temperature gradient is zero everywhere, heat conduction across the interface is expected to take on the Spitzer value associated with isotropic heat conduction. In case (2) however, the angle between the magnetic field and the temperature gradient is always 90°, so with our approximations, heat cannot flow across the interface.

We define a heat transfer efficiency ζ equal to the magnetic-field-regulated heat transfer rate divided by the isotropic Spitzer rate, namely,

$$\begin{equation} \zeta = \frac{q}{q_i}, \end{equation} \tag{ 8 }$$

where q is defined as the amount of thermal energy transported through the interface per unit time. The average angle θ between the temperature gradient and the uniform magnetic field then plays an important role in determining ζ. At θ = 0, ζ = 1. At θ = π/2, ζ = 0.

We now address the influence of both a mean field and a tangled field on ζ. Consider there to be a strongly tangled local field that has no mean value in the direction normal to the interface, i.e., B_{0, x} whose total magnitude is B₀, and a global magnetic field B_d aligned with the normal of the interface of magnitude B_d. If B_d ≫ B₀, the magnetic field around the interface only slightly deviates from the normal direction and ζ should be close to 1. If B_d ≪ B₀, ζ should be close to zero. The asymptotic limit of the heat transfer efficiency is given by Chandran & Cowley (1998) that ζ ≃ 1/(ln(d/ρ_e)), where d is the scale length of the magnetic field fluctuation and ρ_e is the electron gyroradius. In our problem, we estimate this limit at about 4.7%. In the subsequent context, we will use ζ₀ to denote this asymptotic limit.

If B_d and B₀ are comparable, we expect ζ₀ < ζ < 1. We also expect ζ can change throughout the evolution if the structure of the magnetic field is modified by the dynamics of heat transfer. It is instructive to ask whether the feedback from the magnetic field structure evolution will amplify the heat transfer by creating more channels or shut it down. The answer depends on the influence of magnetic reconnection as we will see from the simulations. Only if magnetic reconnection acts to smooth out local small-scale structures and link the initially isolated structures to the global mean field across the interface would we expect the heat conductivity to increase.

In what follows, we refer to the initial tangled field region as "the interaction region." Figure 1(a) shows a schematic of initial and hypothetical evolved steady state field configurations for such a tangled field set up. From the figure, we can see that the initial field configuration forms a "wall" that restricts energy transfer across the two interaction regions. However, if the subsequent evolution evolves to the steady state shown in (b), then expansion of the interaction region and magnetic reconnection has allowed the field to penetrate through the entire region. Thus, the initial "wall" of tangled field is destroyed and thermal conduction will be less inhibited than initially. We will check how accurately this proposed picture of destruction of field wall is valid from analyzing our numerical simulations and quantitatively discuss the effects on the energy transfer.

Zoom In Zoom Out Reset image size

For our initial conditions, we set up an interface between hot and cold regions in mutual pressure equilibrium. The temperature distribution on the horizontal (x) axis is given by

$$\begin{equation} T(x) = T_0[1-4(x-x_0)^2/w^2]^{0.4} \end{equation} \tag{ 9 }$$

in the region [x₀, x₀ + w] with T₀ = 100 in computational units. This region is the interaction region as described in the previous section, with x₀ as the left end and w as the width. In the regions (0, x₀) and (x₀ + w, L_x), where L_x is the domain length, we simply assume that T(0 < x < x₀) = T(x = x₀) and T(x₀ + w < x < L_x) = T(x = x₀ + w): in other words, the temperature profile has a sharp gradient inside [x₀, x₀ + w], while outside this region, it remains flat. The obtained temperature distribution is plotted in Figure 2(a). We set x₀ = 0.4 and w = 0.1. The region 0.4 < x < 0.5 is therefore the interaction region. The temperature is constant and uniform across the regions of each respective side of the box connecting to that side of the interaction region. We are primarily interested in the region of the box where the heat transfer occurs and noticeably evolves during the simulation run time. This means that we will mainly focus on the interaction region. The horizontal length of the interaction region in the simulation domain is L_x = 0.8 in computational units.

Zoom In Zoom Out Reset image size

The thermal pressure is set to be in equilibrium over the entire box, that is

$$\begin{equation} P(x) = P_0 \end{equation} \tag{ 10 }$$

with P₀ = 100. The density distribution is set up by the ideal gas law, namely,

$$\begin{equation} \rho (x) = \frac{P(x)}{T(x)} \end{equation} \tag{ 11 }$$

in computational units.

For the Spitzer diffusion coefficient, we assume that the diffusion is linear as in Equation (7), and use the approximation: κ_∥ = κ_c T^2.5_mid, where κ_c is the classical conductivity and T_mid is taken to be the middle value of temperature across the interface, about 0.5 T₀.

We choose the initial field configuration:

$$\begin{equation} B_x = B_d+B_0\,\sin (n\,\pi \,y/\lambda), \end{equation} \tag{ 12 }$$

$$\begin{equation} B_y = B_0\sin (n\,\pi \,x/\lambda), \end{equation} \tag{ 13 }$$

where n and λ are the mode number and wavelength of the tangled field, respectively, B₀ = 10⁻³ in computational units, and B_d can assume various initial values that reflect the evolving global field as the result of reconnection. The magnetic field configuration is laminar, and there is no broad spectrum of magnetic fluctuations. The magnetic spectrum is concentrated at length scale λ. This initial field configuration is therefore one of a locally tangled field surrounding the interface with one measure of the tangle given by

$$\begin{equation} R = B_d/B_0. \end{equation} \tag{ 14 }$$

When R = 0, there are only locally confined field lines, whereas R = ∞ indicates a straight horizontal field without any tangling. As R increases, the relative fraction of field energy corresponding to lines that penetrate through the interaction region increases. In our simulations, we consider cases with R = 0.0, 0.2, 0.4, 0.6, 1, 2, 4, ∞. Figure 2(a), Figure 3(a), and Figure 4(a) show the magnetic field configuration for initial R values of 0.0, 0.4, and 1.0.

We note that our MHD approximation a priori implies that the electron gydroradius is much smaller than the length scale of one grid cell. Thus, the dissipation scale and all field gradient scales are larger than the electron gyroradius by construction in our simulations.

We run simulations with typical resolution of 2048 cells on the horizontal axis in fixed grid mode. Runs with doubled resolution showed no significant differences compared to the standard resolution runs. We use fixed boundary conditions at the x-axis boundaries: the pressure, density, and temperature at the two ends are fixed to their initial values, as is the magnetic field. We use periodic boundary conditions for the y-axis boundaries.

There are five parameters whose influence determines the simulation behavior and guide interpretation of results.

• 1.  
  Plasma β. β ≡ 8πP/B² has little effect on diffusion because even with very high values of the plasma β used in the simulation, we are still in the MHD regime and the gyroradii of electrons are assumed small. Thus, the direction of thermal conduction is not locally affected by β. It is possible that instabilities could arise in the low β limit that affect pressure balance during the evolution of the simulations but that turns out not to be the case for the β range of 10⁵ ∼ 10⁸ that we use. The value of β in this range does not exhibit any influence on the simulation result as indicated by our numerical experiment.
• 2.  
  Initial tangle measure R = B_d/B₀. If R ≫ 1, the local small-scale field can mostly be ignored and Spitzer thermal conductivity is expected, whereas if R ≪ 1 a value much less than Spitzer is expected.
• 3.  
  Ratio of the diffusion timescale to the sound crossing timescale for one grid cell.

  $$\begin{equation} r = t_{{\rm diff}}/t_{hy}=\frac{\rho \,C_s\,l}{\kappa _{\parallel }}, \end{equation} \tag{ 15 }$$

  where ρ is the density, l is the characteristic gradient length scale of temperature: l = min (T/|∇ T|), and C_s is the sound speed. If r ≪ 1, thermal diffusion would initially dominate and the pressure equilibrium would be broken by this fast energy transfer. If r ≫ 1, then the pressure equilibrium would be well maintained throughout the entire evolution and the energy transfer may be viewed as a slow relaxation process. In our simulation, r ≈ 0.3 initially, so that diffusion induces a pressure imbalance. Eventually, as the heat transport slows, the pressure equilibrium catches up and is maintained.
• 4.  
  Ratio between the temperature gradient scale length and the wavelength of the tangled field: h = 2 π l/λ = k l. If h = 0 there is no tangled field, and no inhibition to heat transfer. As h increases, the field becomes more tangled, and the energy is harder to transfer. However, a large h value may also result in increased magnetic reconnection, because the Lundquist number of field confined in a smaller region is larger for the same field strength. This would then lower h.
• 5.  
  Mean global energy transfer rate: q = δE/t_bal, where t_bal is defined as the time needed for the hot region and cold region to reach a certain degree of temperature equilibrium by a transfer of heat energy δE across the interface.

A mathematical expression for the heat transfer rate can be derived by considering a slab with a planar interface aligned with the y-direction at the middle of the interaction region with a tangled magnetic field, and an average temperature gradient aligned in the x-direction. As in Section 4, we denote the region [x₀, x₀ + w] as the interaction region which contains the interface and the tangled field. Define the average temperature gradient inside the interaction region as |∇ T|_g = (T_hot − T_cold)/w, where the subscripts "hot" and "cold" denote the characteristic temperatures of the hot and cold sides at the two ends of the interaction region. We assume that the resulting effective heat flux depends on how much straight mean field can penetrate through the interaction region. We also assume that the heat flux depends on the average temperature gradient of the interaction region in the form: q∝|∇ T|_g. By integrating the strength of the straight mean field normalized by the total field strength over the volume of the interaction region (and since there is no z-dependence, the essential content is an area integral) we obtain the effective heat flux through this region:

$$\begin{equation} \overline{q} = D\,|\nabla \,T|_{g}\int {\frac{B_d}{|B|}\,dx\,dy}, \end{equation} \tag{ 16 }$$

where D is a constant that depends on neither the magnetic field nor the average temperature gradient and |B| is the local field strength. The two-dimensional integration is carried out over the interaction region: with x₀ < x < x₀ + w, 0 < y < L_y. Note that this expression is valid only when the magnetic field is varying at a length scale smaller than the interaction region length.

Using Equations (12), (13), and (14) in Equation (16) and the approximation that the areal average in the interaction region 〈B₀ · B_d〉 ∼ 0 so that 〈(B₀ + B_d)²〉 ∼ 〈B²₀ + B²_d〉, we obtain

$$\begin{equation} \overline{q} \approx D\,\frac{|\nabla \,T|_{g}\,R}{\sqrt{1+R^2}}. \end{equation} \tag{ 17 }$$

For the unmagnetized isotropic case, or for transfer with a field entirely aligned with the temperature gradient, we have instead

$$\begin{equation} \overline{q}_i = D\,|\nabla \,T|_{g}. \end{equation} \tag{ 18 }$$

Dividing Equation (18) by Equation (19), we obtain an approximation for the heat transfer efficiency over the interaction region:

$$\begin{equation} \zeta = \frac{R}{\sqrt{1+R^2}}. \end{equation} \tag{ 19 }$$

It should be pointed out that our approximation does not take into account the ζ₀ "leakage" from the magnetic field fluctuation as stated in Section 3 and Chandran & Cowley (1998). If the initial temperature profiles are identical for different field configurations, this formula can then be used to estimate the expected energy transfer rate from situations with various field configuration. By normalizing the heat transfer rate to that of the isotropic heat conduction case, we obtain the heat transfer efficiency, ζ. The accuracy of Equation (19) can be tested by plotting the heat transfer efficiency obtained from the simulations against measured values of R.

If magnetic reconnection occurs during the time evolution of the heat transfer process, then conduction channels can open up and the energy exchange can be accelerated. We would then expect the actual curve of ζ versus R to evolve to be higher than the value Equation (19) predicts in situations with low R values. Meanwhile, for high R, the analytical prediction and the real physical outcome should both approach the horizontal line ζ = 1, which denotes conductive efficiency consistent with the unmagnetized case. We emphasize that R as used in this paper is always calculated with the initial values of the magnetic field, not time evolved values, and that Equation (19) is valid when estimating a cold–hot interface with initial tangle measured as the ratio of the initial global straight field to the initially local tangled field. To follow a measure of the tangle that evolves with time, a generalized tangle measure should be calculated in a more sophisticated manner and the integral form (Equation (16)) should be applied.

We choose initial conditions with values R = 0.0, 0.2, 0.4, 0.6, 1.0, 2.0, 4.0 to run the simulations. The simulation run time is taken to be 1.2, which corresponds to 12, 000 years in real units for WBBs. The initial cuts of temperature and magnetic field lines for R = 0.0, 0.4, 1.0 are shown in Figure 2(a), Figure 3(a), and Figure 4(a), respectively. Figure 5(a) shows the initial cut of the density distribution in the R = 0.0 run. We also run simulations with purely horizontal magnetic field lines, equivalent to the R = ∞ case, and runs with purely vertical field lines. Frames (b)–(d) in Figures 2–5 are from the late stages of the evolution, and the final frames always display the steady state of the runs. A steady state is facilitated by the fact that the boundaries are kept at a fixed temperature throughout the simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 6, we plot the mean cuts of the temperature T_c, obtained by averaging the temperature along y-axis, against the x-position for selected evolution times. Since the anisotropic heat conduction is initially faster than the pressure equilibration rate, the energy distribution around the temperature interface changes rapidly until about t = 0.4. This energy transfer is mostly confined to the interaction region for the low R₀ runs, since in these cases only a few field lines can penetrate into the entire interaction region.

Zoom In Zoom Out Reset image size

During the initial heat exchange phase, the thermal energy and density quickly redistribute in the interaction region. As seen in Figure 2(b), islands at x = 0.48 are formed by material bounded by the magnetic field lines, since the field orientation blocks heat exchange with the surroundings. Around x = 0.4, there are also cavities formed where the thermal energy is inhibited from flowing. The magnetic field lines, which form complete sets of loops in the R = 0.0 case, begin to distort. It can be observed that the field lines are more strongly distorted in the low density part of the interaction region than in the high density part. This occurs because velocity gradients are driven by the early rapid redistribution of heat (pressure) by conduction.

At time t = 0.4 (see Figure 2(b)), the field lines surrounding the cavities at x = 0.4 reconnect, making thermal exchange possible. During the evolution, field lines begin to link the interaction region to the hot material on the left. This phenomenon is most apparent in Figure 2(d), which marks the final state of the thermal energy exchange. We also see that there is little difference between Figures 2(c) and (d), because at such a late stage of the process, the thermal diffusion gradually slows so that the magnetic field configuration approaches a steady state.

By comparing Figure 6(c) with Figure 6(d), we see that the mean cuts of temperature show little difference for all values of R. The mean cuts of temperature T_c exhibit a jump in the region of x = 0.35 ∼ 0.5, but are relatively smooth on either side of this region. This shows that even though the tangled field "wall" has been broken and allows channels of thermal conduction through it, the temperature profiles is not as smooth as in the purely straight field case.

For the cases of R = 0.4, there are field lines that penetrate the entire interaction region from the start. By observing the evolution of the magnetic field lines at about x = 0.38, we see that magnetic reconnection is still happening and causes the field loops to merge. The observed behavior resembles the process displayed by Figure 1. When R = 1, there are hardly any temperature islands that bounded by magnetic field loops. The evolution of the field lines shows less dramatic reconnection–the field lines evolve in what appears to be a more gentle straightening.

We begin our analysis with the evolution of the heat flux. The average heat flux per computation cell for different values of R is plotted as a function of time in Figure 7(a). Note that in the vertical field case (${\bf B} = B_y {\hat{y}}$), the heat flux remains zero as the field entirely inhibits electron motion across the interface. For cases with R > 1, the heat flux decreases throughout the evolution. Recall that R > 1 implies cases where the "tangled" portion of the field is relatively weak and heat is quickly transported from one side of the interface to the other. Thus, the trend we see for R > 1 occurs as the temperature distribution approaches its equilibrium value. For lower R values, especially those of R < 0.5, an initial phase of heat flux amplification is observed as magnetic reconnection in the early evolution opens up channels for heat to transfer from hot to cold regions. At the late stage of the evolution when reconnection has established pathways from deeper within the hot region to deeper within the cold region, temperature equilibration dominates, leading to a decreasing heat flux phase as observed in the R > 1 cases. Note that the similarity between the R > 2 cases and the R = ∞ case is predicted by Equation (19): as the global field comes to dominate, the heat flux inhibition imposed by anisotropic heat conduction in the local tangled field can be ignored.

Zoom In Zoom Out Reset image size

In order to understand the influence of magnetic reconnection on heat transfer rates we compare simulations with different filling fractions of the tangled field. Two cases are shown in Figure 8(a): (1) a temperature interface with a "volume filling" tangled field and (2) temperature interface with the tangled field filling only the region surrounding the interface. In case (2) the rest of the domain is filled with straight field lines connecting the hot and cold regions. From Figure 8(a), we see that case (1) shows much slower heat transfer rates compared with what is seen in case (2). This is because reconnected field lines in case (2) are linked to the globally imposed background field that in turn link the hot and cold reservoirs. In case (1) reconnection only leads to larger field loops but cannot provide pathways between the reservoirs. The effect of different scale lengths on the evolving field loops is shown in Figure 8(b), in which we plot the result from three simulations wavelengths for the tangled field component (tangled field "loops"). Note that λ is defined in Equations (12) and (13). We use a sequence of values for wavelengths: 2λ, λ, and λ/2. Figure 8(b) clearly shows that smaller field loop λ leads to the largest average heat flux, since smaller scale loops will reconnect before large loops for a given magnetic resistivity. This result demonstrates the link between the number of reconnection sites of the field and heat flux.

Zoom In Zoom Out Reset image size

We next analyze the temperature equilibration in detail. The averaged temperature difference across the interface is plotted in Figure 7(b). It shows the difference between the averaged temperature at the hot side and the cold side. One significant feature in Figure 7(b) is that the temperature difference decreases to a steady value T_end in all cases. This resembles percolation across a membrane, which allows a density jump to happen when filtering two fluids. Figure 7(c) shows the distance required for the temperature to drop 80% at the interface. This distance characterizes the length of the interaction region. Except for the vertical field case where no heat transfer is allowed, the interface is expanding at different rates for different R values. The expansion for all the cases of nonzero R approaches a steady value, which is also a characteristic feature of the temperature equilibration evolution.

We now analyze the modification of magnetic field configuration during the evolution. Throughout our simulations, the local magnetic field is initially a set of complete loops surrounding the interaction region. Once the energy transfer begins, the interaction region tends to expand as discussed previously. This expansion stretches the field lines on the x-direction and distorts these circular loops, eventually inducing magnetic reconnection that opens up channels connecting the hot and cold regions. From the current $J_B = |\nabla \times \textbf {B}|$, we can get information on how tangled the field is. Figure 7(d) shows the evolution of the mean value of the strength of $\nabla \times \textbf {B}$ in the interaction region. We observe that in the vertical and straight field case, $|\nabla \times \textbf {B}|$ remains constant, but decreases to a fixed value for R ⩾ 2 cases. This means the field in high R cases is straightened by the stretching of the interaction region as seen in Figure 7(c). For the R ⩽ 1 cases, we see that $|\nabla \times \textbf {B}|$ increases. This rise is due to magnetic energy brought in via the cold mass flow and the creation of fine field structures that amplify J_B faster than dissipation caused by interface expansion.

The local field distortion can be clearly demonstrated by studying the energy evolution of magnetic energy stored in different field components. In Figure 9(a), we plot the evolution of mean magnetic energy stored in the vertical field $\bar{B}_y^2/2$, compared with $\bar{B}_x^2/2$. We note that the latter includes only the fluctuating contribution to the energy in the x field—that is, the contribution to the horizontal field that does not come from the global mean x-component.

Zoom In Zoom Out Reset image size

From Figure 9(a), we observe that the B²_y energy decreases while the B²_x energy either increases or remains the same for all cases. The magnetic energy evolution can thus be viewed as a conversion of vertical field to horizontal field. This conversion need not conserve the total magnetic energy of the local tangled field because of magnetic reconnection and because material advecting magnetic field can flow in and out of the interaction region. By comparison, in the R > 1 cases, the thermal energy and local magnetic energy can both decrease and add to the kinetic energy of the material surrounding the interface, because of the fast thermal diffusion enabled by the strong global field.

The distortion of the local field loops can also be demonstrated by plotting the mean eccentricity of the field loops. In Figure 9(b), we plot the mean eccentricity evolution. For all cases, the mean eccentricity is zero initially because of the circular shape of the field loops. Later in the evolution, large R cases tend to evolve into a state of large eccentricity in the steady state. This is caused by a rapid expansion of the interface induced by the strong global field. In short, large R induces more distorted local field loops and less tangled total field due to fast interface expansion, while small R values result in less eccentric local field loops but with more tangled total field and strong magnetic reconnection.

To compute the estimated heat transfer rate in the simulation, we calculate the averaged slope of the curve plotted in Figure 7(b) and compare it to the analytic model in Section 4. Although the equilibration rate represented by the slope of the curves in Figure 7(b) is changing throughout the evolution, an early phase of the evolution can be chosen when the field configuration has not been modified significantly for which we can then compute the averaged heat transfer rate. By normalizing the resulting heat transfer rate to the isotropic value, we can determine the heat transfer efficiency for different magnetic structures. From Figure 10, we can see that the analytic prediction and the simulation results agree quite well except for the situation when R is below 0.2. The simulation result does not converge to point (0,0) but ends at an intercept on the y-axis. This intercept, which is much larger compared with both the approximated model and the asymptotic limit ζ₀, indicates that even if there are initially negligibly few channels for energy transfer, the magnetic reconnection can open up channels and allow heat transfer. Equation (19) is valid for predicting the cooling rate of the hot material throughout the early phase of the heat equilibration process. It also provides insight on the strength of the local field in the vicinity of the interface once we know the cooling rate and global magnetic field strength.

Zoom In Zoom Out Reset image size

It should be pointed out that in our case the electron gyroradius is assumed small compared to the numerical resolution. If we had used an explicit resistivity, then the equivalent assumption would be that the gyroperiod is longer than the resistive time on a field gradient scale of order of the gyroradius. The numerical resistivity, which results in numerical reconnection is always present in our simulations and its effect does not seem to depend on resolution: simulations with double resolution show no significant difference in overall heat transfer efficiency. The existence of numerical resistivity allows the topology of the field to change when scales are approaching the grid scale. As long as this scale is very small compared to global scales, the overall heat transfer rates are not strongly sensitive to this value.

To summarize our results, we find first that the average heat flux at the end of our simulations is lower than at the beginning for all R values. Thus, we see an approach to thermal equilibrium. In some cases, we also see that the heat shows an initially increasing phase denoting a period of active magnetic reconnection.

In the simulations, we see the average temperature difference decreases to a constant value T_end, which is related to R. We also see the width of the initial interface expands to a fixed value during the simulation.

Analysis of the simulation behavior shows that J_B is an accurate measure of structural change in the magnetic field. Current decreases to a constant value for large R cases and increases to a constant value for small R values.

Finally, we have shown that Equation (19) can be used to estimate the energy transfer rate for an initially complicated field structure by considering the relative strength of the local field and the global field. For those cases for which R approaches 0, Equation (19) becomes invalid since the energy transfer in is mainly induced by a feedback from the magnetic field reconnection. By comparing cases with different field loop length scales, we demonstrate that the smaller the field loop length scale is, the faster the reconnection rate will be.

The issue of magnetized conduction fronts and their mediation of temperature distributions occurs in many astrophysical contexts. One long-standing problem that may involve anisotropic heat conduction are hot bubble temperatures in WBBs. WBBs occur in a number of setting including the planetary nebula (PN), luminious blue variables, and environments of Wolf–Rayet stars. When a central source drives a fast wind (V_wind ∼ 500 km s⁻¹), temperatures in the shocked wind material are expected to be of order 10⁷ K, which is greater than 2 keV. The temperatures observed in many WBB hot bubbles from X ray emission are, however, in the range of 0.5–1 keV range. NGC 6888 is a particularly well known and well studied example for a WR star (Zhekov & Park 2010). For PNe, Chandra X-ray observations have found a number of WBB hot bubbles with temperatures lower than expected based on fast wind speeds (Montez et al. 2005; Kastner et al. 2008). The role of wind properties and heat conduction in reducing hot bubble temperatures has been discussed by a number of authors (Steffen et al. 2008; Akashi et al. 2007; Stute & Sahai 2006). The role of magnetic fields and heat conduction was discussed in Soker (1994).

While our simulations herein were meant to be idealized experiments aimed at identifying basic principles of anisotropic heat conduction fronts, we can apply physical scales to the simulations in order to make contact with WBB evolution. Table 1 shows the results of such scaling. Upon doing so, we infer that (1) given field strengths expected for WBBs, heat conduction is likely to be strong enough to influence the temperature of the expanding hot bubble and the cold shell bounding it. We also note that magnetic fields in WBBs (for PN field strengths, see Vlemmings et al. 2006) are usually in the milli-gauss range and are relatively much stronger than the field strength that can be scaled to our simulations. Thus, the magnetic field in realistic WBBs is highly likely to result in anisotropicity and regulate the behavior of heat conduction. Since the heat transfer does not directly depend on the magnetic β, we can thus apply our analysis to the WBB interface if we approximate the interface to be planar and stationary, which is reasonable as the radius of curvature of WBBs is much greater than the interface scale of relevance. We must also assume that the global magnetic field is primarily radial.

The computational parameters used in our simulations and the real physics parameters typical in a WBB are listed in the first two columns of Table 1. We choose the domain length to be 0.025 pc, which is about 1% of the radius of the actual WBB. Table 1 shows that by choosing the proper scaling, our simulation fits well with the data observed in a typical WBB. Therefore, the conclusions we draw by analyzing the simulation results and the analytical expressions, especially Equation (19), can be helpful in analyzing WBB evolution.

We have investigated the problem of heat transfer in regions of initially arbitrarily tangled magnetic fields in laminar high β MHD flows using simulation results of ASTROBEAR code with anisotropic heat conduction. Three conclusions stand out as follows.

(1) Hot and cold regions initially separated by a tangled field region with locally confined field loops may still evolve to incur heat transfer. The local redistribution of fluid elements bend the field lines and lead to magnetic reconnection that can eventually connect the hot and cold regions on the two sides. (2) The temperature gradient through such a penetrated tangled field region tends to reach a steady state that depends on the energy difference between the hot and cold reservoirs on the two ends. (3) Equation (19), a measure of the initial field tangle, is a good predictor of the ultimate heat transfer efficiencies across the interface for a wide range of R.

A basic limitation of our simulations is that they are two dimensional. A three-dimensional version of this study would be of interest as the field would then have finite scales in the third dimension possibly allowing channels for heat transfer that are excluded in two dimensions. We have also not considered the effects of cooling in our simulations. The absence of cross field diffusion is also not realistic in our parameter regime. Future simulations should include both the diffusion parallel and perpendicular to the field.

Future directions of analysis could also include a multi-mode study, which investigates the effect of the spatial spectrum of the magnetic field distribution on the heat transfer efficiency. When there are multiple modes or a spectrum is continuous, it would be useful to predict how the efficiency would depend on the spectrum. In this context, a more detailed comparison of heat transfer in initially laminar versus initially turbulent systems would be of interest.

Financial support for this project was provided by the Space Telescope Science Institute grants HST-AR-11251.01-A and HST-AR-12128.01-A, by the National Science Foundation under award AST-0807363 and NSF PHY0903797, by the Department of Energy under award DE-SC0001063, and by Cornell University grant 41843-7012. We thank Jonathan Carroll, Kris Yirak, and Brandon Shroyer for useful discussions.