\documentclass[12pt]{article}
\usepackage[latin1]{inputenc}
\usepackage{palatino,url}
\usepackage{graphicx}
\begin{document}

\title{Heat transport across hot - cold interface with magnetic field}
\author{Shule Li\\Department of Physics and Astronomy, University of Rochester}
\date{December 2010}
\maketitle

\begin{abstract}
The low temperature of the hot bubbles observed in the wind blown bubbles (WBB) remains a mystery for many years. It is assumed that the heat conduction plays an important role in the cooling of the hot bubbles. Meanwhile, in a stellar object such as WBB, the magnetic field can be strong and highly tangled which results in the anisotropic heat conduction which only allows the electrons to carry thermal energy along the magnetic field lines. A natural question is, in the situation such as a WBB, how would the heat transfer rate which relates to the hot bubble cooling rate be related with the magnetic field structure? In this study, we start from two simple cases: (1) field lines aligned with the hot bubble - cold shell interface; (2) field lines perpendicular to the hot bubble - cold shell interface. We explore the relation between the magnetic field topology and the heat transfer rate by adding tangled field to the above configuration and demonstrate that the degree of tanglement and the heat transfer rate has a simple mathematical relation which matches well with the analytical calculation. This relation can be useful in determining the heat transfer efficiency when there is a magnetic field in the environment of astrophysical context.\\\\
\textit{\textbf{keywords:} magneto-hydrodynamics, planetary nebula, magnetic reconnection, wind blown bubbles, anisotropic heat conduction}
\end{abstract}

\section{Introduction}
\paragraph{} The interface bounding hot and cold interstellar medium (ISM) exist in astrophysics problems of various scales. The thermal conduction through such interface is critical in determining the heating process of the materials behind the interface. However, the thermal conduction in ISM is usually controlled by the magnetic field topology, which allows the electrons to move freely only along the field lines. This results in a strong thermal conductivity along the field lines and weak conductivity across the field lines. When considering ISM, it is usually a valid assumption that the electrons are totally inhibited from moving across field lines (Braginskii 1965, McCourt et al. 2010), giving the electron mean free path is much greater comparing to the electron gyroradius. Under this assumption, the interface bounding the hot and cold materials will be significantly modified by the presence of tangled magnetic field, creating small structures at the field variation length scale. One typical example of such problems being the wind blown bubbles (WBB) which are observed to have hot bubble with observed temperature lower than expected (Zhekov et al. 2010). Many believe that heating through the hot bubble - cold shell interface via thermal conduction may have induced the cooling of the hot bubble (Zhekov et al. 1998, Zhekov et al. 2001). Meanwhile, the presence of the temperature gradient and global field will continuously cause materials of the cold side to evaporate, causing interface instabilities and mass mixing (Stone 2009), which would further modify the magnetic field topology. Another example is the unexpected slow mass deposition rate of the cooling flow in central galaxy cores, for which many suggest is a result of global thermal conduction (Rosner et al. 1989, Balbus et al. 2008, Mikellides 2010). In the situation of intracluster medium (ICM), the tangled magnetic field often renders a high anistropic thermal conductivity that may significantly change the temperature and density profile (Narayan et al. 2001, Mikellides et al. 2010). The key questions here are: (1) does the interface evolve into an instability driven by the anisotropic thermal conduction or a steady state stabilized by the anisotropic thermal conduction? (2) how fast is the thermal exchange between the hot and cold sides comparing to the isotropic thermal conduction case taking into account the dynamic feature of the interface? In this paper, we investigate the various situations with different magnetic field configurations imposed on a planar hot cold interface and study the effect of the tangled magnetic field on the evolution of the interface. Using ASTROBEAR code with anisotropic thermal conduction, we run simulations with different initial setup. In section 2, we review the basic equations of MHD with anisotropic thermal conduction. Sections 3 and 4 provide detailed description of the simualtion setup, section 5 and 6 present the simulation results and analysis. In section 7, we apply the simulation results to astrophysical problems at different scales: the WBB cooling problem and the cooling flow problem in cores of galaxy clusters. The appendix offers detailed information on the testing of ASTROBEAR code.

\section{MHD Equations with Anisotropic Heat Conduction}
The MHD equations with anisotropic heat conduction are:
\begin{equation}
\frac{\partial \rho}{\partial t}+\nabla \cdot (\rho \textbf{v}) = 0,
\end{equation}
\begin{equation}
\frac{\partial (\rho \textbf{v})}{\partial t}+\nabla \cdot [\rho \textbf{vv}+(p+\frac{B^2}{8\pi})\textbf{I}-\frac{\textbf{BB}}{4\pi}] = 0,
\end{equation}
\begin{equation}
\frac{\partial \textbf{B}}{\partial t}+\nabla \times (\textbf{v}\times \textbf{B}) = 0,
\end{equation}
\begin{equation}
\frac{\partial E}{\partial t}+\nabla \cdot [\textbf{v}(E+p+\frac{B^2}{8\pi})-\frac{\textbf{B}(\textbf{B}\cdot \textbf{v})}{8\pi}]+\nabla \cdot Q = 0,
\end{equation}
where $E$ denotes the total energy:
\begin{equation}
E = \epsilon+p\frac{\textbf{v}\cdot \textbf{v}}{2}+\frac{\textbf{B}\cdot \textbf{B}}{8\pi}
\end{equation}
and the internal energy:
\begin{equation}
\epsilon = \frac{p}{\gamma -1}
\end{equation}
with $\gamma = 5/3$. In our simulation, it is assumed that the heat flux is confined to be parallel to the magnetic field lines. This assumption is valid when the electron gyroradius is much smaller compared to the scale of the problem, i.e. the magnetic field strength has a threshold to let the anisotropicity to be effective. Under this assumption, the anisotropic heat flux can be written as:
\begin{equation}
Q = -\kappa_{\parallel}(\nabla T)_{\parallel}
\end{equation}
where $\kappa_{\parallel}$ is the classical Spitzer heat conductivity. We assume the heat conductivity to be a constant throughout the simulation and will be writing it simply as $\kappa$.
\paragraph{}The ASTROBEAR code uses operator splitting method to treat the MHD equations with heat conduction. The ideal MHD equations are solved with the MUSCL primitive method with TVD preserving Runge-Kutta temporal interpolation, the result is then sent to the implicit linear solver utilizing High Performance Preconditioners (HYPRE) to solve the anisotropic heat conduction equation. The linear solver requires temporal sub-cycling technique to maintain its accuracy. The time step taken by the linear solver can be written as:
\begin{equation}
t_{diff} = \frac{\rho l^2}{\kappa}
\end{equation}
where $l$ is the minimum scale length of the temperature variation: $l = min(T/|\nabla T|)$. The code runs in parallel with fixed grid domain.

\section{Problem Description and Analytical Model}
\paragraph{}The problem involves a hot region and a cold region which are separated by a thin planar interface. We will study how the magnetic field configuration alters the heat transfer rate between the hot and cold regions in presence of anisotropic heat conduction. \setlength{\parskip}{0cm}

\paragraph{}We begin by considering two simple but illustrative cases: (1) a uniform magnetic field aligned with the normal direction of the interface; (2) a uniform magnetic field perpendicular to the normal direction of the interface. In case (1), because of the angle between the magnetic field and the temperature gradient is zero everywhere, the anisotropic heat conduction acts exactly the same as isotropic heat conduction. In case (2) however, the angle between the magnetic field and the temperature gradient is always $90^{\circ}$, thus the interface would stay unchanged. If we define a heat transfer rate $\eta$ that is normalized by the transfer rate in the isotropic conduction case, namely:
\begin{equation}
\eta = \frac{q}{q_i}
\end{equation}
then the average angle $\theta$ between the temperature gradient and the uniform magnetic field plays an important role in determining $\eta$. At $\theta = 0$, $\eta = 1$. At $\theta = \pi/2$, $\eta = 0$. \setlength{\parskip}{0cm}
\begin{figure}[h]
\begin{center}
\includegraphics[width=.6\columnwidth]{ini_ste}
\caption{\small{The initial and steady state field configuration. (a): the initial field forms complete loops that only allows heat transfer within the interaction region. (b): the steady state field reconnects itself so that it allows heat transfer between regions deeply into the hot and cold areas. We call the interaction region in (b) the ``field membrane'' because it allows the heat transfer between the hot and cold areas but the transfer is slower than the isotropic case since the heat flux is confined on the field lines. }}
\end{center}
\end{figure}
\paragraph{}We now consider the following configuration: there is a strongly tangled local field $B_0$ whose average on the normal of the interface is zero. There is also a global magnetic field $B_d$ aligned with the normal of the interface. If $B_d \gg B_0$, the magnetic field around the interface merely deviates from the normal direction and $\eta$ should be close to 1. If $B_d \ll B_0$, one would expect that there will be local energy redistribution, but the global energy transfer should be slow. Thus the global temperature structure should remain unchanged and $\eta$ should be close to zero although an expansion of the interface will still occur. \setlength{\parskip}{0cm}

\paragraph{} If $B_d$ and $B_0$ are comparable, $\eta$ should be between zero and one. $\eta$ in this case will change throughout the evolution since the strucuture of the magnetic field will be modified by the heat transfer. One may ask whether the feedback from magnetic field topology evolution will amplify the energy channeling by creating more channels or shut it down eventually. If the feedback is positive, one would expect a continuous increase on the overal heat flux and expansion of interface in the weak global field case, resulting in the locally connected field loops to reconnect to the global field structure, making the heat flux be able to penetrate. Otherwise one would expect the isolated local field loops to stay for longer period of time. \setlength{\parskip}{0cm}

\paragraph{} From now on we will call the tangled field region the interaction region. Fig.1.(a) shows the initial and steady state field configuration in such a tangled field interface situation. From the figure we can see that the initial field configuration forms a ''wall'' which does not allow the energy transfer between the two ends of the interaction region. But the steady state sees an expansion of the interaction region and magnetic reconnection which allows the field to penetrate through the entire region so that the field wall is destructed. In this paper, we will call the destructed field wall as the ``field membrane'', which describes a stable region with tangled field distribution that allows a certain amount of heat to penetrate through in terms of anisotropic heat conduction, as in Fig.1.(b). In the following study, we will check whether our assumption of the formation of field membrane is valid or not in the numerical simulation and quantitatively discuss its effect on the energy transfer. \setlength{\parskip}{0cm}

\section{Simulation Setup}
\paragraph{}To study the hot - cold interface energy transfer in circumstances such as the WBB, we construct the initial condition by setting up an interface at pressure equilibrium. The temperature follows the distribution:
\begin{equation}
T(x) = T_0(1-x^2)^{1/2.5}
\end{equation}
in the region of $0 < x < 1-10^{-5}$ with $T_0 = 100$ in computational units. Later we will see that our simulation can be used to describe certain real astrophysical situations once correct scaling parameters are chosen. The temperature distribution is plotted in Fig.2.(a). The region $0.4 < x < 0.5$ is the interaction region we defined in the previous section. At the two sides, the temperature is set to be constant and continuous from the ends of the interaction region. During the simulation, we are only interested in the region where the heat transfer can have a noticeable impact during the simulation run time. In other words, we will focus ourselves on the region surrounding the interaction region, not regions deeply into the hot material. The horizontal length of the simulation domain is $0.1$ in computational units. The thermal pressure is set to be in equilibrium:
\begin{equation}
P(x) = P_0
\end{equation}
with $P_0 = 100$. The density distribution is set up by the perfect gas law, namely:
\begin{equation}
\rho(x) = \frac{P(x)}{T(x)}
\end{equation}
For the Spitzer diffusion constant, we use the approximation of $\kappa = \kappa_c\,T_{mid}^{2.5}$, where $\kappa_c$ is the classical conductivity, and $T_{mid}$ is taken to be the middle value of temperature across the interface, which is about $0.6\,T_0$.\
\paragraph{} In our study, we choose the field configuration:
\begin{equation}
B_x = B_d+B_0\,\sin(n\,\pi\,y/\lambda),
\end{equation}
\begin{equation}    
B_y = B_0\sin(n\,\pi\,x/\lambda)
\end{equation}
where $n$ and $\lambda$ are the mode number and wavelength of the tangled field respectively, $B_0 = 10^{-3}$ in computational units, and $B_d$ can take various values to reflect the gradually changing global field comparing to $B_0$. This field configuration creates a locally tangled field surrouding the interface. We use ratio
\begin{equation}
R = B_d/B_0
\end{equation}
to describe the ''tangleness'' of the magnetic field. $R = 0$ indicates locally confined field lines, $R = \infty$ indicates straight field with all the field lines go from the hight temperature region to the low temperature region without any twist. One can imagine that with $0 < R < 1$, the field lines will be twisted and form sets of circles, but some of the loopss are incomplete so that some of the field lines still go from the hot region into the cold region which become channels for energy transfer. In our simulation, we will consider cases with $R = 0.0,\,0.2,\,0.4,\,0.6,\,1,\,2,\,4,\,\infty$. Fig.2.(a), Fig.4.(a) and Fig.5.(a) show the magnetic field configuration for initial $R$ values of $0.0$, $0.4$, $1.0$.\

\paragraph{} We run the simulation with standar resolution 2048 cells on the x axis in fixed grid. Runs with doubled resolution are also carried out, the results of which demonstrate no significant difference comparing to the standard resolution runs. We use fixed boundary condition on the x direction: the pressure, density and temperature at the two ends are fixed to their initial values so as the magnetic field. By this boundary condition, the two ends on the x axis can be seen as the regions deeply into the hot and cold regions, whose thermal dynamic variables stay constant. On the y axis, we use the periodic boundary condition.\setlength{\parskip}{0cm}

\paragraph{}There are five parameters involved in the simulation:\setlength{\parskip}{0cm}

\paragraph{}\textbf{1. Magnetic $\beta$.} $\beta$ has little effect on diffuson since the field lines only determine the direction of the diffusion. However, extremely large $\beta$ may affect the anisotropicity of the system. By contrast, extremely small $\beta$ may result in pressure imbalance in the initial condition. In our simulations, $\beta$ falls in the range that it does not cause a strong pressure imbalance (the magnetic pressure is $10^{-8}$ in magnitude comparing to the thermal pressure), and also does not affect the anisotropicity (effective gyroradius is about $10^{-5}$ of the scale of the problem considered).\setlength{\parskip}{0cm}

\paragraph{}\textbf{2. Tangle vs straight measure: $R = B_d/B_0$.} If $R$ is large, the local tangled field can mostly be ignored, we would get high energy transfer efficiency, and vice versa.\setlength{\parskip}{0cm}

\paragraph{}\textbf{3. Ratio of the diffusion time scale and the hydrodynamic time scale: $r = t_{diff}/t_{hy}$.} By simple physics consideration, we have:
\begin{equation}
r = \frac{\rho\,C_s\,l}{\kappa}
\end{equation}
where $l$ is the characteristic scale length of temperature: $l = min(\frac{T}{|\nabla\,T|})$, $C_s$ is the sound speed, $\kappa$ is the linear diffusion constant. The initial kick is driven by the diffusion, the later evolution of the system depends on the value of $r$. If $r$ is much smaller than 1, diffusion would dominate the initial phase and the pressure equilibrium may be broken by the fast energy transfer, if $r$ is large, then the pressure equilibrium is well maintained throughout the entire process and the energy transfer may be viewed as a slow relaxation process. In the cases when the two time scales are comparable, both effects are present. In our simulation, $r \approx 0.3$ initially, so that the diffusion would trigger the heat transfer and pressure imbalance fastly at first. Then the hydrodynamics process would catch up and become faster than the diffusion process as a result of the decreasing of the temperature gradient and the fast expansion of the interface. \setlength{\parskip}{0cm}

\paragraph{}\textbf{4. Ratio between the temperature scale length and the wavelength of the tangled field: $h = 2\,\pi\,l/\lambda = k\,l$} where l follows the same definition as that in the previous paragraph. $h = 0$ marks the situation of no tangled field, which means there is nothing inhibitting the energy transfer. As $h$ gets larger, the field becomes more confined, the energy is harder to transfer in the tangled field region, but a large $h$ value may also result in a large chance of magnetic reconnection. The number of modes may also affect the result. We currently only look at the single mode case in which the spectrum is a delta function.\setlength{\parskip}{0cm}

\paragraph{}\textbf{5. Mean global energy transfer rate: $\eta = 1/t_{bal}$}, where $t_{bal}$ can be defined as the time needed for the hot region and cold region to reach a certain degree of temperature equilibrium.\setlength{\parskip}{0cm}

\paragraph{}A mathematical expression for the heat transfer rate can be derived by considering a slab with a sharp temperature interface aligned with the y direction and an average temperature gradient aligned on the x direction. There is an interaction region surrounding the temperature interface with tangled magnetic field in it. Define the global temperature gradient as $|\nabla\,T|_{g} = (T_{hot}-T_{cold})/(T_0\,L)$, where the subindices ``hot'' and ``cold'' denote the characteristic temperatures of the hot and cold regions, $T_0$ is a normalization factor which makes $|\nabla\,T|_{g}$ having dimension of 1/length, $L$ is the width of the interaction region. We can then integrate over the interaction region to 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}
where D is a constant that depends on neither the magnetic field nor the temperature distribution, $|B|$ is the local field strength. Notice that this expression is valid only when the magnetic field is varying at a length scale smaller than the interaction region length.\\
Using Eqs.(13), (14) and (15), we find the approximated relation 
\begin{equation}
\overline{q} \approx D\,\frac{|\nabla\,T|_{g}\,R}{\sqrt{1+R^2}}
\end{equation}
Meanwhile, in the isotropic heat conduction case, the heat transfer rate purely depends on the temperature gradient:
\begin{equation}
\overline{q}_i = D\,|\nabla\,T|_{g}
\end{equation}
Deviding Eq.(18) by Eq.(19), we get the appoximated heat transfer efficiency over the interaction region to be:
\begin{equation}
\eta = \frac{R}{\sqrt{1+R^2}}
\end{equation}
Thus if the initial temperature profiles are identical for different field configurations, i.e. if $|\nabla\,T|_{g}$ are the same for the cases we are looking at, we can calculate the energy tranfer rate from situations with various field configuration and normalize them by the heat transfer rate of the isotropic heat conduction case to obtain the heat transfer efficiency. Later we will plot the heat transfer efficiency ($\eta$) against the field configuration ($R$) in our simulation results to verify the approximate relation Eq.(20).\setlength{\parskip}{0cm}

\paragraph{}If magnetic reconnection happens during the transfer process then channels can be opened up and the energy exchange process can be accelerated, we would expect the actual curve of $\eta$ vs $R$ to be higher than the value Eq.(20) predicts at regions where $R$ is low. Meanwhile, at regions of high $R$, the analytical prediction and the real physical outcome should both approach the horizontal line $\eta = 1$, which denotes perfect efficiency.

\section{Simulation Result}

\paragraph{} We choose field set up 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 0.6 million years in real units for the problem of WBB. The initial cuts of temperature and magnetic field lines for $R = 0.0,\,0.4,\,1.0$ are shown in Fig.2.(a), Fig.4.(a) and Fig.5.(a) respectively. Fig.3.(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 which is equivalent to the $R=\infty$ case and purely vertical field lines. The rest frames in Fig.2 to Fig.5 are from the later stages of the evolution, and the final frames always display the steady state of the runs.\

\begin{figure}[h]
\begin{center}
\includegraphics[width=\columnwidth]{int1}
\caption{\small{Evolution of temperature distribution with $R = 0.0$. The cuts are at (a): $t=0.0$, the initial state, (b): $t=0.4$, (c): $t=0.8$, (d): $t=1.2$, the steady state.}}
\end{center}
\end{figure}

\paragraph{} In Fig.6, we plot the mean cut temperature $T_c$ against the $x$ position for selected evolution time. $T_c$ is obtained by averaging the temperature on the y direction, because the question we are interested in is how the heat transfer would perform along the x axis. Since initially the anisotropic heat conduction is faster than the hydrodynamic process, the energy distribution around the temperature interface will see a fast change until about $t = 0.4$. This energy transfer is mostly confined to the interaction region for the low $R$ runs, because in these cases there are little to no field lines that can penetrate into the entire interaction region. 

\begin{figure}
\begin{center}
\includegraphics[width=\columnwidth]{rho1}
\caption{\small{Evolution of density distribution with $R = 0.0$. The cuts are at (a): $t=0.0$, the initial state, (b): $t=0.4$, (c): $t=0.8$, (d): $t=1.2$, the steady state.}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=\columnwidth]{int2v2}
\caption{\small{Evolution of temperature distribution with $R = 0.4$. The cuts are at (a): $t=0.0$, the initial state, (b): $t=0.4$, (c): $t=0.8$, (d): $t=1.2$, the steady state.}}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[width=\columnwidth]{int3v2}
\caption{\small{Evolution of temperature distribution with $R = 1.0$. The cuts are at (a): $t=0.0$, the initial state, (b): $t=0.4$, (c): $t=0.8$, (d): $t=1.2$, the steady state.}}
\end{center}
\end{figure}

During this initial heat exchange phase, the thermal energy quickly redistribute itself in the interaction region as well as the density. As we can see from Fig.2.(b), islands at $x=0.48$ are formed by materials bounded by the magnetic field lines, since they cannnot exchange heat with the environment. \

\paragraph{} The mass flow from the cold material pushes the interface, being heated up and forming an evaporation which is subject to interface instability because of the materials being accelerated to super-Alfven speed at the interface. This instability does not grow further but is rather suppressed because of the tangled field lines. There are also cavities formed in between the spikes at around $x=0.4$ because the thermal energy is inhibitted from flowing out of them. The magnetic field lines, which forms complete sets of loops in the $R=0.0$ case, begin to distort. We can see that the field lines are distorted heavier in the low density part of the interaction region rather than in the high density part, because of the different speed they gain during the density redistribution. At time $t=0.8$, we observe that the field lines surrounding the holes at $x=0.4$ reconnect, making the thermal energy of the outter part of the holes begins to exchange. Field lines that run from the high temperature region on the left side to the low temperature region deeply into the interaction region begin to emerge, making possible the energy exchange between the hot side and the cold side. This phenomenon is more apparent in Fig.2.(d), which marks the final state of the thermal energy exchange. We also see that there are little difference between Fig.2.(c) and Fig.2.(d), because at the later stage of the process, the magnetic field configuration is approaching the steady state by which we mentioned as ``field membrane'' earlier in this paper. It is impossible for the temperature gradient to drive the material to move further at the steady state since the field membrane allows a certain degree of temperature jump to happen across it. This effect is also apparent by comparing Fig.6.(c) with Fig.6.(d), in which we observe that the mean cut temperature distribution changes little for all the different $R$ values investigated. The mean cut temperature $T_c$ sees a sharp jump in the region of $x=0.35\,\sim\,0.5$, but is relatively smooth on the two sides. Because except for the interaction region which later becomes the field membrane and induces a temperature jump, the two sides of the domain are mainly governed by the heat conduction under horizontal magnetic field, which is identical to the situation of isotropic heat conduction. \

\begin{figure}
\begin{center}
\includegraphics[width=\columnwidth]{tcut}
\caption{\small{Evolution of mean cut temperature averaged on y direction with different $R$ values labeled by different colors. The cuts are at (a): $t=0.0$, the initial state, (b): $t=0.4$, (c): $t=0.8$, (d): $t=1.2$, the steady state.}}
\end{center}
\end{figure}

\paragraph{} 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 find that there are still magnetic reconnection happenning, which causes the magnetic blobs to merge. The observed behavior resembles the process displayed by Fig.(1). When $R=1.0$, we can hardly find temperature islands that are bounded by magnetic blobs. The topological evolution of the field lines is better decribed as being straightened.\

\section{Discussion}
\begin{figure}[h]
\begin{center}
\includegraphics[width=.8\columnwidth]{meanflux}
\caption{\small{(a) top left: time evolution of mean heat flux at the interface, (b) top right: time evolution of average temperature difference between the hot and cold regions, (c) bottom left: time evolution of interface width, (d) bottom right: time evolution of the mean value of $|curl\textbf{B}|$}}
\end{center}
\end{figure}
\paragraph{}We first analyze the evolution of heat flux. The average heat flux per computation cell for different $R$ values is plotted as a function of time in Fig.7.(a). The two processes in effective are the temperature equilibration across the interface and the evolution of the magnetic field topology. From Fig.7.(a), we can see that the average heat flux for $R = \infty$ remains zero, which indicates that the field structure inhibits heat transfer completely across the interface. For cases with $R > 1$, heat flux decreases throughout the evolution because of the temperature equilibration. For the other $R$ values, especially those of $R < 0.5$, a phase of heat flux amplification is observed. This shows that the magnetic reconnection occurred in the early phase of the evolution which tends to open up channels for heat to transfer from hot region to cold region. At the late stage of the evolution when there are enough channels opened up, the temperature equilibration becomes the dominating effect which leads to a similar decreasing phase as observed in the $R > 1$ cases. The similarity between the $R > 2$ cases and the $R = \infty$ case is predicted by Eq.(20): as the global field getting strong, the inhibition imposed by the local tangled field tends to be ignorable.\setlength{\parskip}{0cm}
\begin{figure}[h]
\begin{center}
\includegraphics[width=.9\columnwidth]{fluxcompare}
\caption{\small{(a) Comparison of averaged heat flux for situation with field loops filling up the entire domain and situation with field loops only fill the interaction region. (b) Comparison of averaged heat flux for situations with different tangeld field length scale.}}
\end{center}
\end{figure}

\paragraph{} The influence of magnetic reconnection on the heat transfer can be demonstrated comparing the averaged heat flux of the following two cases shown in Fig.8.(a): (1) temperature interface with tangled field filling up the entire domain and (2) temperature interface with tangled field filling the region surrounding the interface, the rest of the domain is filled with straight field lines. Case (1) is much slower in terms of heat transfer comparing to case (2) since the reconnected field lines in case (2) are connected to the global straight field lines, comparing to those in case (1) reconnects only to form larger field loops. Another comparison is also made to demonstrate the impact of the scale length of field loops on the reconnection rate, which is shown in Fig.8.(b). Here we plot the averaged heat flux at time $t = 0.3$ for three different field loop length scale: $\lambda$ as used in the previous simulation result, $2\lambda$ and $\lambda/2$. We can clearly see from the plot that the run with the smallest field loop length scale has the greatest averaged heat flux which indicates the strongest reconnection rate.\setlength{\parskip}{0cm}

\paragraph{}We next analyze the two contributing physics processes in detail. The temperature equilibration across the interface is plotted in Fig.7.(b). It shows the difference between the averaged temperature at the hot side and the cold side. One significant feature in Fig.7.(b) is that the temperature difference decrease to a steady value $T_{end}$ in all cases. This resembles the percolation of a membrane which allows a density jump to happen when filtering two fluids. Fig.7.(c) shows the distance required for the temperature to drop $80$ percent 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 for all different $R$ values at a different rate. The expansion for all the cases of nonzero $R$ approach a steady value which is also a characteristic feature of the temperature equilibration evolution. \setlength{\parskip}{0cm}

\paragraph{}We now analyze the second aspect in the evolution: the modification on magnetic field configuration. 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 results in the distortion of these circular loops, which eventually induces the magnetic reconnection which oppens up channels connecting the hot and cold regions. In our analysis, we use $J_B = |curl \textbf{B}|$ as a characteristic value measuring the degree of magnetic field tanglement. Fig.7.(d) shows the evolution of the mean value of the strength of $curl \textbf{B}$ in the interaction region. We observe that in the vertical and straight field case, $|curl \textbf{B}|$ remains constant, but decreases to a fixed value for $R \geq 2$ cases. This means the field in high $R$ cases is straightened by the stretching of the interaction region as seen in Fig.7.(c). For the $R \leq 1$ cases, $|curl \textbf{B}|$ increase. This is due to the energy deposited into magnetic field in the interface instabilty with low $R$. \setlength{\parskip}{0cm}
\begin{figure}[h]
\begin{center}
\includegraphics[width=.9\columnwidth]{fig8}
\caption{\small{(a) Comparison on evolution of local field energy in terms of $B_x$ and $B_y$. Circles corresponds to the $B_x^2/2$ curve, stars corresponds to the $B_y^2/2$ curve. The different colors denote various $R$ values. (b) Eccentricity of the ellipses constructed by assigning the mean values of local $|B_x|$ and $|B_y|$ to the major and minor axes, respectively. The set of curves show different evolution patterns for different $R$ values.}}
\end{center}
\end{figure}
\paragraph{}The local field distortion can be clearly demonstrated by studying the energy evolution of magnetic energy stored in forms of $B_x$ and $B_y$. In Fig.9.(a), we plot the evolution of mean magnetic energy stored in the vertical field, namely $\bar{B}_y^2/2$, comparing to that stored in the mean horizontal field $\bar{B}_x^2/2$. Notice that to exclude the influence of the global field, here we only count in the tangled field part. From Fig.9.(a), we observe that the $B_y$ energy decrease while the $B_x$ energy either increase or stay stable for all cases. The magnetic energy evolution can thus be viewed as a conversion of field energy from vertical pointing to horizontal pointing, made possible by the magnetic reconnections. However, this conversion may not conserve the total magnetic energy of the local tangled field, because of the existence of the induced interface instability and magnetic reconnection. By comparison, in the $R > 1$ cases, the thermal energy and local magnetic energy would both decrease and add to the kinetic energy of the material surrounding the interface, because of the fast diffusion caused 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 Fig.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 at the steady state. This is due to the fast expansion of the interface induced by the strong global field. To summarize the characteristics of the magnetic field topological evolution, large $R$ value induces more distorted local field loops and less tangled total field with the fast interface expansion, while small $R$ values results in less eccentric local field loops but with more tangled total field and strong magnetic reconnection.\

\begin{figure}[h]
\begin{center}
\includegraphics[width=.7\columnwidth]{heatrate}
\caption{\small{heat transfer rate observed in the simulation compared with the analytic model}}
\end{center}
\end{figure}

\paragraph{}To compute the estimated heat transfer rate in the simulation, we calculate the averaged slope of the curve plotted in Fig.7.(b), and compare it to the analytic model in Section 4. Although the equilibration rate represented by the slope of the curves in Fig.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 and comptute the averaged heat transfer rate. By comparing the heat transfer rate under various situations with different $R$ values and the ideal heat transfer rate with magnetic field parallel to the gradient. We can then find out the heat transfer efficiency under different field topologies. From Fig.10, we can see that the two 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 interception on the y axis. This interception indicates that even if there are few to no channels for energy transfer initially, the magnetic reconnection can open up channels and allow the heat transfer to happen. Eq.(20) is valid for predicting the cooling rate of the hot material throughout the early phase of the heat equilibration process. It also offers us insights on the strength of the local field around the interface once we know the cooling rate and the global magnetic field strength. \setlength{\parskip}{0cm}

\paragraph{} As a quick summary of the discussion, we find that the average heat flux is in general a decreasing function for all $R$ values, the curve also has an inital increasing phase which denotes the stage of strong magnetic reconnection. The average temperature difference decreases to a constant value $T_{end}$ which is a function of initial $R$. Because of this feature, the magnetized interface can be viewed as an ``energy membrane'' that allows a certain energy jump to happen. The width of the interface expands to a fixed value during the simulation. The term that describes the structural change of the magnetic field is $J_B$. It decreases to a constant value for large $R$ cases while increases to a constant value for small $R$ values. We proved that Eq.(20) can be used to estimate the energy transfer rate in the case of complex field topology by considering the relative strength of the local field and the global field. For those cases with $R$, Eq.(20) becomes invalid since the energy transfer in these cases are mainly induced by a feedback from the magnetic field reconnection. By comparing cases with different field loop length scales, we demonstrated that the smaller the field loop length scale, the faster the reconnection rate. \\

\section{Astrophysical Applications}

\begin{tabular}{c|c|c}
Variables & Computional Units & WBB \\
Number Density & $1$ & $1\,cm^{-3}$ \\
Temperature & $100$ & $1\,kev$ \\
Domain Length & $0.1$ & $0.1\,pc$ \\
Local Field Strength  & $10^{-3}$ & $2^{-8}\,Gauss$ \\
Global Field Strength  & $10^{-4}$ & $2^{-9}\,Gauss$ \\
Evolution Time  & $1.2$ & $60,000\,yrs$ \\
Heat Conductivity  & $10^{-2}$ & $2\times10^{-18}\,cm\,s\,g^{-1}\,K^{-2.5}$ \\
\end{tabular}\\

The problem of magnetically modified temperature interface exist in many astrophysical contexts. One long-standing problem being the observed low temperature of the hot bubble in a Wind Blown Bubble (WBB) formed by a stellar wind. The WBB is predicted to have hot bubbles with a characteristic temperature greater than $2\,kev$, but the temperature observed from X ray emission detections is located at the $0.5\,kev$ to $1\,kev$ range (for details of observed NGC 6888, see Zhekov S.A. et al. 2010). This seemingly unnatural low temperature can be related to radiative cooling, or heat conduction. By scaling our computational data with the real physics data and putting into simulations, we find out two facts that are important in determining the actual temperature history of the WBB evolution: $(1)$ the heat conduction in the WBB is strong enough to alter the evolution and impose an important influence on the temperature of the expanding hot bubble and the cold shell bounding it; $(2)$ the magnetic field in the WBB, usually in the mili-Gauss range, is much stronger than the field strength in the simulations. Thus the magnetic field in a realistic WBB is usually strong enough to regulate the behavior of heat conduction. Since the heat transfer does not directly depend on the magnetic $beta$, we can thus apply our analysis to the WBB interface if we approximate the interface to be planar and stationary, which is accurate as the radius of WBB being much greater than the scale of the problem considered. We also assume that the global magnetic field is mostly radial. The computational parameters used in our simulations and the real physics parameters typical in a WBB is listed in the first two columns of Tab.1. We choose the domain length to be $0.1\,pc$, which is about $4$ percent of the radius of the actual WBB. Tab.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 analysing the simulation results and the analytical expressions, especially Eq.(20), can be used in analysing real WBB evolution. \\

\section{Conclusion}
\paragraph{}The problem of heat transfer efficiency of a region filled with tangled magnetic field lines is crucial to many astrophysical problems. In this paper, we investigate this problem using the simulation results of ASTROBEAR code with anisotropic heat conduction. There are three conclusions that are : (1) hot and cold regions separated by a tangled field region with only locally confined field loops and no field lines penetrating through may still have heat transfer because 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 equilibration of temperature through such a penetrated tangeld field region tends to reach a steady state dependent on the energy difference of the hot and cold reservoirs on the two ends, the steady state field configuration allows a certain jump of temperature which prevents an isothermal final state. (3) we can use Eq.(20) to determine the when both confined local field and straight global field are present to estimate the heat transfer efficiency across the interface.\setlength{\parskip}{0cm}

\section{Future Works}
\paragraph{} The future works of this study may include a multi-mode study, which investigates the impact of the spacial spectrum of the magnetic field distribution on the heat transfer efficiency and the topological evolution of the thermal dynamic variables. The interception on the y axis in Fig.13 indicates the R = 0 case heat transfer efficiency. When there are multiple modes or spectrum is continuous, it would be useful to predict how the efficiency would depend on the spectrum. One of the limitation of the study presented in this paper is that it rules out the regions where magnetic field becomes weak in the tangled field region, which creates a small area where isotropic heat conduction may become possible. This will create a leaking channel through the tangled field region even if the global field is zero. Thus it may be helpful to study the anisotropicity effect as a follow-up study. We also know that the radiative cooling would be effective once the temperature is above $10^4\,k$. Since in our simulation the temperature has an order of magnitude of $10^7\,k$, the cooling would have a noticeable influence on the evolution. Another possible direction is to run the simulation of the hot cold interface with anisotropic heat conduction, radiative cooling, turbulence and gravity present.\\\\\\

\appendix{\textbf{Appendix: Code Test}}
\paragraph{} We use the magneto-thermal instability (MTI) problem to test the accuracy of the ASTROBEAR code with anisotropic heat conduction. The problem involves setting up a 2D temperature profile with uniform gravity pointing on the y direction. The domain is square with length of $0.1$ in compuational units. The temperature and density profiles are:

\begin{equation}
T = T_0\,(1-y/y_0)
\end{equation}
\begin{equation}
\rho = \rho_0\,(1-y/y_0)^2
\end{equation}

with $y_0=3$. The pressure profile is set up so that the fluid can be balanced under uniform gravity with gravitational acceleration $g = 1$ in computational units. We also set $T_0=1$ and $\rho_0=1$ in computational units. There is a uniform magnetic field on the x direction with field strength $B_0=1.0\times10^{-3}$ in computational units. The anisotropic heat conductivity is set to be $\kappa=1\times10^{-4}$ in computational units. We use the pressure equilibrium condition for the top and bottom boundaries, that is, the pressure in the ghost cells are set up so that its gradient can balance the gravitational force. On  the x direction, we use the periodic boundary condition. \setlength{\parskip}{0cm}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=.9\columnwidth]{mti1}
\caption{\small{Field line evolution of magneto-thermal instability. (a): initial state. (b): $t=75\tau_s$. (c)$t=150\tau_s$. (d):$t=250\tau_s$.}}
\end{center}
\end{figure}

\paragraph{} Initially, the domain is in pressure equilibrium. We then seed a small velocity perturbation:

\begin{equation}
v_{per} = v_0\,\sin(n\,\pi\,x/\lambda)
\end{equation}

with $v_0=1\times10^{-6}$ and $\lambda=0.5$. This perturbation will cause the fluid elements to have a tiny oscillation on y axis as well as the field lines. Once the field lines are slightly bent, they open up channels for heat to transfer on the y direction thus allowing the heat on the lower half of the domain to flow to the upper half. It can be shown that this process has a positive feedback so that once the heat exchange happens, more channels will be openned up for heat conduction. Therefore this process forms an instability whose growth rate can be verified according to the linear theory on the growth rate of weak field stratified atmosphere. We use $\tau_s$ to denote the sound crossing time for the initial state. Fig.11 shows the time evolution of the field lines at various stages in our MTI simulation.

\begin{figure}[h]
\begin{center}
\includegraphics[width=.9\columnwidth]{lnvt}
\caption{\small{(a): $\ln\,v_y$ against evolution time in $\tau_s$. (b): calculated growth rate against evolution time in computational units. Initially the growth rate is stable around the theoretical value $0.4$ and then decreases sharply after $t=200$, which indicates the evolution has entered the nonlinear regime.}}
\end{center}
\end{figure}

We study the MTI growth rate by considering the acceleration of the fluid elements. The mean speed on the y direction for the fluid should follow the exponential growth:

\begin{equation} 
v_y = v_{per}e^{\gamma\,t}
\end{equation}

where $v_{per}$ is the strength of the initial velocity perturbation applied, $\gamma$ denotes the growth rate in the linear regime. We obtain the growth rate $\gamma$ by plotting $\ln\,v_y$ against the evolution time and then measuring the local slope through a certain time span. The $\ln\,v_y$ vs $t$ curve is plotted in Fig.12.(a), which shows a nice linear relation. We plot the growth rate against evolution time. It should be stable around the theoretical value $0.4$ initially and then decrease sharply due to the nonlinear effect. Fig.12.(b). shows that the simulation meets our expectation fairly well.\

\paragraph{} Next we look at the energy evolution in the linear regime. The mean kinetic energy should first stay stable and then enter into an exponential growing phase until it hits a cap at around $t=200$ which denotes the starting of the nonlinear phase. The evolution of magnetic energy should follow similar pattern as to the kinetic energy evolution, but lagged behind. In Fig.13, we plot the time evolution of the mean kinetic and magnetic energy evolutions. The results confirms the physical intuition quite well.\

\paragraph{} The dispersion relation of MTI can be expressed as a function of 

\begin{equation}
\phi=(\gamma-1)\,\frac{\kappa\,k^2}{N\rho_{mid}}
\end{equation}

where $k=2\pi/\lambda$, N is the $Brunt-V\ddot{a}s\ddot{a}l\ddot{a}$ frequency. Here we use the mid value of $\rho$. We can plot the MTI growth rate against $\phi$ to obtain the dispersion relation shown in Fig.14.(a).\

\begin{figure}
\begin{center}
\includegraphics[width=.9\columnwidth]{keme}
\caption{\small{(a): evolution of mean kinetic energy. (b): evolution of mean magnetic energy.}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=.9\columnwidth]{geta}
\caption{\small{(a): the dispersion relation of MTI linear growth. (b): impact of field tension on the growth rate}}
\end{center}
\end{figure}

\paragraph{} Finally we study the impact of magnetic tension on the MTI growth rate. A strong initial field would provide a tension force that tends to stabilize the exponential growth. Thus intuitively the growth rate
would be negatively correlated to the initial $Alf\acute{v}en$ speed. When the $Alf\acute{v}en$ speed approaches the value threshold value $\gamma_{max}/k$, we would expect a strong stabilizing tension to kill the instability growth. We use the stabilizing parameter as a function of

\begin{equation}
\eta=\frac{k^2v_A^2}{\gamma_{max}^2}
\end{equation}

Here $v_A$ is the initial maximum $Alf\acute{v}en$ speed, $\gamma_{max}$ is the maximum possible growth rate which can be derived from the mid values of the initial temperature and pressure profiles. We plot $\gamma$ against $\ln{\eta}$ in Fig.14.(b). We see that the field tension becomes increasingly effective as $\eta$ approaching 1.

\end{document}