Resistivity

The resistive equation can be written as:
frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B} The magnetic diffusivity eta is a function of temperature as we all know.
Now for the first equation, what if we want to include the temperature dependence into the code? From the next section, you can see that when a
field configuration in equilibrium is subject to strong diffusion, usually heating would occur and surppress the local resistivity and thus the
diffusivity. By expanding the first equation, we have the form:
frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 \textbf{B} + \nabla \eta \times (\nabla \times \textbf{B} Now the second term depends on all of the three components of B. So we end up with equations in which the time variance of Bx, By and Bz depend
on each other.
So the coefficient array is thus no longer a tri-diagonal matrix. For cases where resistive speed is slow (which is usually the case), we can use
explicit solver to treat the problem instead.
One may wonder if it is possible to throw away the second term on the right hand side of the diffusion equation to just let the diffusivity to vary
with position but ignore its own spacial variation. This would give us a form of resistive MHD similar to that of the thermal diffusion case but
with temperature dependence built in.
frac{\partial \textbf{B}}{\partial t}=\eta (T) \nabla^2 \textbf{B Unfortunately, this does not work because the diffusion equation itself has to be divergence free. Our roughest approximation, treating the resistivity
as a constant satisfies the requirement as long as the divergence and the Laplacian are commutable. If we try to do the same thing to the above equation,
we end up getting:
frac{\partial \textbf{B}}{\partial t}=\eta \nabla^2 (\nabla \cdot \textbf{B}) + (\nabla \eta \cdot \nabla)\textbf{B The first term on the right hand side is zero but the second term is not, especially at sharp temperature fronts where grad(T) is large.
So this approximation only works under slowly varying temperature situation.
Upon a closer look at the simulations conducted below, we can observe that the energy inside the domain does not conserve. There is a small increase of
total energy during the evolution, especially for the force-free field case. This phenomenon is explained below.
In AstoBEAR, we explicitly calculate the resistivity induced current on the cell edges, following equation:
textbf{J} = \eta \nabla \times \textbf{B
The stencil for this explicit solver is a 3 by 3 cube surrounding the cell we want to update. The magnetic field is represented by the aux field, which is
centered on the cell faces. Its curl therefore reside on the cell edges. Here is an example on calculating the diffusive current on the x direction, notice
that the red arrow is where we are calculating the diffusive current, the green arrows are where the magnetic field originally resides:

The actual code looks like the following in 2D (jy are initialized to 1):

    DO i=mjy(1,1),mjy(1,2); DO j=mjy(2,1),mjy(2,2); DO k=mjy(3,1),mjy(3,2)
       IF(nDim==2)THEN
          jy(i,j,k) = -jy(i,j,k)*(Info%q(i,j,k,iBz)-Info%q(i-1,j,k,iBz))/dx
       ELSE
          jy(i,j,k) = -jy(i,j,k)*((Info%aux(i,j,k,3)-Info%aux(i-1,j,k,3))+(Info%aux(i,j,k-1,1)-Info%aux(i,j,k,1)))/dx
       END IF
    END DO; END DO; END DO 

After calculating and storing jx, jy and jz of the given grid, the change of magnetic field due to diffusion can be calculated as：
frac{\partial \textbf{B}}{\partial t}=\nabla \times \textbf{J
In the case of AMR, we need to store the emf using function storefixupfluxes. In case of resistivity, it is simply the diffusive current. In 2D case, the emfs
are only in z. So we need to use f instead. The code looks like:

    IF(nDim==2)THEN
       ALLOCATE(fjx(mjy(1,1):mjy(1,2),mjy(2,1):mjy(2,2),mjy(3,1):mjy(3,2),1))
       ALLOCATE(fjy(mjx(1,1):mjx(1,2),mjx(2,1):mjx(2,2),mjx(3,1):mjx(3,2),1))
       ALLOCATE(bflux(1))
       bflux(:)=iBz
       fjx(:,:,:,1)=jy(:,:,:)*dt/dx
       fjy(:,:,:,1)=-jx(:,:,:)*dt/dx
       CALL storefixupfluxes(Info,mjy,1,fjx,bflux)
       CALL storefixupfluxes(Info,mjx,2,fjy,bflux)
       DEALLOCATE(bflux)
       DEALLOCATE(fjx); DEALLOCATE(fjy)
    ELSE
       CALL StoreEmfs(Info,mjx,1,jx*dt/dx)
       CALL StoreEmfs(Info,mjy,2,jy*dt/dx)
    END IF

In the case of resistive MHD, the energy can be dissipated in the form of Joule heat, comparing to the infinite conductivity case, where the voltage inside
the fluid is everywhere zero, and no heat is generated by the current.
If we dot the resistive induction equation with the magnetic field B, we obtain the time evolution equation for magnetic energy:

partial (B^2)/\partial t + \nabla \cdot \textbf{S} = - \textbf{J}\cdot\textbf{E} = -j^2/\et

where textbf{S}=\textbf{J} \times \textbf{B is the magnetic energy flux caused by resistive diffusion and =|\textbf{J} is the magnitude
of the diffusive current.
In this equation, the S term accounts for the redistribution of magnetic energy (and thus the redistribution of total energy), and the last term accounts for
the loss of magnetic energy due to reconnection.
The total energy change for the resistive step is therefore:

partial \epsilon/\partial t + \nabla \cdot \textbf{S} =

Here the ^2/\et dissipation term is absent because the dissipation of magnetic energy does not change the total energy: the loss of magnetic energy
is converted into thermal energy.
In the code, the energy flux as a result of magnetic diffusion needs to be calculated explicitly using:
textbf{S}=\textbf{J} \times \textbf{B
The energy fluxes reside on the cell faces while the diffusive current reside on the cell edges. We therefore need to compute a face average of the diffusive
current as well as the magnetic field components which are not normal to the face using the surrounding edges. In the last diagram, the blue arrows connected by
dashed lines are what used to compute the energy flux.
The magnetic field can be updated from the diffusive currents by:
frac{\partial \textbf{B}}{\partial t}=\nabla \times \textbf{J
The updating of magnetic field and the calculation of the energy fluxes can be done at the same time, as in the code:

    DO i=msx(1,1),msx(1,2); DO j=msx(2,1),msx(2,2); DO k=msx(3,1),msx(3,2)
       IF(nDim==2)THEN
          Info%aux(i,j,k,1) = Info%aux(i,j,k,1)+(jz(i,j,k)-jz(i,j+1,k))*dt/dx
          xflux(i,j,k,1)=-0.25*(jz(i,j+1,k)+jz(i,j,k))*(Info%q(i-1,j,k,iBy)+Info%q(i,j,k,iBy))
       ELSE
          Info%aux(i,j,k,1) = Info%aux(i,j,k,1)+((jy(i,j,k+1)-jy(i,j,k))+(jz(i,j,k)-jz(i,j+1,k)))*dt/dx
          xflux(i,j,k,1)=0.25*(jy(i,j,k+1)+jy(i,j,k))*(Info%q(i-1,j,k,iBz)+Info%q(i,j,k,iBz))-0.25*(jz(i,j+1,k)+jz(i,j,k))*(Info%q(i-1,j,k,iBy)+Info%q(i,j,k,iBy))
       END IF
    END DO; END DO; END DO

We then call storefixupfluxes for the three energy flux components (as xflux, yflux and zflux in the code).
Finally, the energy is updated using the divergence of the energy fluxes. This finishes the diffusive process.
The resistive time scale is explicitly calculated by:
t_{resistive} = \frac{c dx}{\eta
where c is the cfl number (set to be 0.5 in the code). To change the cfl number, the following line of code can be modified:

    explicit_maxspeed(level) = max(explicit_maxspeed(level),2d0*resistivity/dx)

where 2d0 can be replaced by any value above 1 (for example, putting this number to 5d0 will reduce the cfl to 0.2).

Resistivity Test Problems

Sweet-Parker Problem

To investigate the Sweet-Parker problem, we introduce the following magnetic field configuration called "sheer pinch":
_y(x) = b_0 tanh(x/a

where in computational units we choose _0 = and = 0..
The density profile is chosen so that the pressure equilibrium can be maintained with constant temperature.

rho(x) = \rho_0/cosh^2(x/a)+\rho_

where rho_0 = and = 0.. The temperature is set to be constant as 0.5.
The domain is set to be -6.4 < x < 6.4 and -12.8< y < 12.8, with fixed resolution 480 * 960. The boundaries are all open. The initial profile is plotted below:


The initial state is in pressure equilibrium though unstable. There are two ways to generate instabilities. The first way is to artificially increase the resistivity
at the center of the domain. This increase will result in a higher reconnectivity, which will eventually bend magnetic field. This creates an X point where field
lines continue to come in and annihilate because of the lower field pressure at the center. The reconnection heat will drive outflows out of the X point, parallel
to the direction of the sheer pinch. The box surrounding the X point where the outflows (Petschek shock) come out of is called the "Sweet-Parker Box". The following
diagrams show how increased resistivity at the center of a sheer pinch drives Petschek shock.


The Mach number in 2D case is plotted in pseudo-color:

To test Sweet-Parker problem in AstroBEAR, we construct the sheer pinch using the above setup, and modify .

The following figure shows the Petschek shock from a reconnection spot at the center. Colored variable is the kinetic energy in log scale, magnetic field is illustrated
by white lines.

A movie with 2-level AMR:
​http://www.pas.rochester.edu/~shuleli/multi_test/spmovie.gif open movie

The second way to perturb the sheer pinch is to put in a sinusoidal perturabtion on magnetic field:

_x = B_p sin(kx)cos(ky
_y = -B_p cos(kx)sin(ky

where Bp is the perturbation amplitude and the wave number:
= 2\pi/L_.
This perturbation creates periodical X points and O points, which leads to bright and dark spots of magnetic pressure. Material flow into the dark spots as a result
of pressure imbalance, which creates periodical dense "islands".The growth rate and the size of the "islands" depend on resistivity and the strength of the perturbation.

To watch a full movie, ​click here

Notice that in our island formation problem, the perturbation is anti-symmetric about the center. So the density pattern is also anti-symmetric. It is trivial to use
a symmetric perturbation in order to obtain a symmetric island pattern.

Microphysical Resistivity

In cgs units, the real magnetic diffusivity can be written as:

eta = c^2/4\pi \sigm

where sigma is the material conductivity which depends on the electron-ion collision rate mu which in term depends on the electron temperature to the three halves.
The cross field diffusivity can thus be expressed as function of electron temperature if we assume the plasma is fully ionized:

eta = 8.2243 \times 10^5 Z_{eff} F(Z_{eff}) \ln \Lambda T_{e}^{-3/2

where T is the electron temperature in eV, Zeff is the effective ion charge.The Coulomb logarithm for interested number density is plotted below:

Here, T is the same as before, n is the electron number density in the units of per cubic centimeter. F is a function of effective ion charge:

(Z) = \frac{1+1.198Z+0.222Z^2}{1+2.996Z+0.753Z^2

Since the electron's maxwellian will be distorted when there is a strong outer magnetic field applied, the realistic resistivity can depend on the local field orientation.
The anisotropicity is around 2 for the parallel and cross field resistivity. This requires the electron mean free path to be much longer than the gyroradius under such a field:
_B << \lambda_{mfp

or, the gyro frequency to be much greater than the mean electron ion collision frequency.
This effect is not considered currently.

How to Turn on Different Types of Resistivity

Currently in the code, the microphysical resistivity gets calculated at each cell edges by averaging the temperatures and density information surrounding that edge.
The resistive source term can be turned on in astrobear by putting the resistive flag to true:

lResistive     = .false.    ! Turns on resistivity [.false.]

The second parameter ResType determines which type of resistivity one requires: 1 for constant computational resistivity; 2 for artificial resistivity, 3 for microphysical resistivity.

ResType        = 2          ! Resistivity type. 1-constant, 2-user defined distribution, 3-Spitzer. [1]

The third parameter resistivity is the user defined constant resistivity. When ResType = 2, the resistivity is determined by the routine "realResistivity". The
distribution can be manually determined. In this routine, all variables are in computational units. If ResType = 3, the resistivity is microphysical and is automatically
calculated using the scaling parameters. the "resistivity" parameter will be ignored.

resistivity    = 1.2d-1     ! resistivity value [0d0]

The F function gets calculated automatically via and is a global const. If ResType = 3, The microphysical resistivity automatically scales with the input scaling parameter.