=======================================

The 2.5-D Magnetized Cloud with Resistivity

=======================================

The 2.5-D magnetized cloud with resistivity is explored with the following simple configuration:
The Bz field is centralized inside a cylinder-shaped cavity. The density inside the cavity is 10 times smaller than that of the ambient, and the temperature is uniformly distributed.
This configuration can maintain the balanced state for a long time although the field is not force-free.
When there is resistivity, the Bz field will drift away from the center of the cylinder and carry flow away. Here, there are two competing effects on the fluid parcels:
(1)the fluid parcels are carried away by the field depending on the field strength and the resistivity. In the perfect conducting case, the fluid will move along when the field is drifting because of the "frozen in" phenomenon.
(2)when the field is leaving the cavity, it will create pressure imbalance between the cavity and the ambient medium. The pressure gradient will thus drive the flow inward towards the cavity.
(3)the drifting away of the field provides heating which drives up the cavity temperature. This increase in temperature along with the inflow of material will drive the cavity pressure to be higher than the ambient thus drive the flow backwards (outflow).

We use vr, va and cs to identify the speed of the resistive speed, Alfven speed and sound speed respectively. In the discussion we will see that the relative greatness of the three speeds can sometimes define very different evolution behaviors. For sound speed, we have:

and for the resistive speed which characterize how fast the field is drifting away from its center, we have:

where eta is the constant resistivity as defined in the previous section, l is the length scale of field variation. In the discussion, we will take l as the clump or cavity radius.

To illustrate the above effects, we have the following set of simulation results corresponding to the density, temperature, field strength and velocity vector on each row.
This simulation is done with beta = 0.5 inside the cavity. The relationship between vr, va and cs is vr > va > cs.
The field strength shows no surprise that the field drifts away from where it is concentrating, and the material flow in due to the push of pressure imbalance from the ambient (see the velocity plot, row 4, column 2). It is interesting that the fall-in velocity seems to always track the front of the field strength, when comparing row 1 and row 3.
Next let us look at the density. The shape of the density distribution roughly stays the same throughout the simulation, although we can see a material deposition inside the cavity which makes sense. The deposition rate is about 0.1241 at t = 2.0. The temperature rise up shortly after the simulation starts, and stays relatively stable throughout the simulation. In the long run, it is the inflow of material that makes up for the pressure drop.

As we stated before, the above run is done with beta = 0.5, in other words, the field is supporting the cavity initially. Next let us look at a run with high beta: beta = 100. The relation between vr, va and cs is vr > cs >> va. We compare the high beta run with the first run in terms of the density deposition and field diffusion.
As seen in the following picture, the field diffusion occurs at roughly the same rate since we are keeping vr in the two runs. But the mass deposition is much slower, almost unrecognizable in the high beta case, which makes sense: the energy carried away by the field only takes a tiny potion of the total energy. The heating is not plotted here but the same thing happens: the heating is negligible in the high beta case. This result supports our assumption that the heating of the cavity is triggered by the energy carried away by the field diffusion.

One of our intuition is that the mass deposition observed inside the cavity is caused by the pushing from the ambient material. The simulation under weak conduction case, can give us some insight on this topic. In the following simulation, we compare the case with resistivity ten times greater comparing to figure 1. This situation can be called vr > > cs > va. The following graph shows the comparison between field strength, density, and density at different time, for each row. The first row is comparing the field diffusion speed, which does not need more description than that the greater the resistivity, the faster the spreading of the field. The second row compares mass deposition into the cavity at the same time t = 0.1. As we see from the scale, the right plot has a greater deposition rate and more blurred edge. So the weak conduction does help material to deposit into the cavity. So how fast it is? By comparing the third row, we find that the mass deposited into the blue region of the cavity for weak conduction at t = 0.2, is almost identical to the one for strong conduction at t = 2.0. The deposited density does show how resistive the material is. And it seems to scale to the resistivity.

Another assumption one can make is that by looking at the velocity plot, we may think that rotation may happen due to the Lorentz force at the edge of the field front where the field is pointing at z direction and the speed is radial. The initial setup is not force-free so at the edge of the field front, the Lorentz force which is pointing on the azimuthal direction may not be negligible. So the real inflow might have a shape of spiral inward to the center of the cavity instead of straight fall-in. However the tangential velocity is small so the vector plot on figure 1 row 4 looks like mostly radial.

In the next section we examine a more complicated scenario: the force-free configuration.

===================================

The Force Free Field Subject to Resistivity

===================================

We still use the 2.5-D model for the force-free field, where the field has a B_z, B_phi component that only depends on Bessel functions of r. The following vector graph shows how the field looks like initially.

In the infinite conduction case where the field is frozen in, the force-free field configuration can sit literally forever because the Lorentz force is self-canceling everywhere. Now with the the resistivity, we would expect the field to drift and break the force-free status. We conduct simulation with strong resistivity and high beta to let the field diffusion play a prominent role in the evolution. The setup has a beta of beta = 0.2 and with vr > va > cs. One could say that this setup resembles figure 1, the major differences are that (1) the field is force-free; (2) there is no cavity initially. The density and temperature is uniform distributed. The following graph shows the field strength, velocity, density, temperature and the B_x component, respectively on each row.

Let us look at the field diffusion first. One of the drastic difference between the FF and non-FF case is that in the FF case although the field dissipates and reduce in strength, it does not spread in space very much. One of the reasons is that in the FF configuration, the field strength at the edge are largely composed by the B_x and B_y field because they are first order Bessel comparing to B_z which is zero order. B_z contributes mostly to the central area around r = 0. The B_phi field is reconnecting itself during the diffusion so that it reduces the strength in a way that it does not have to spread over. Let us look at row 4 for B_x evolution. The blue and red lobes in the B_x field correspond to opposite directions they point to. And the field diffusion happens so that the positive lobe would spread into the negative lobe's territory and vice versa. If the positive field and the negative field move into one zone, they inevitably cancel each other which at the macroscopic level, appears as field reconnection. The reconnection rate thus connects to how fast the colored lobes expand which in turn relates to how strong the resistivity is. The B_x evolution on row 4 shows that the positive and negative lobes are diminishing due to the reconnection and the zone with green (denotes zero field) increases. This dissipation happens without significant expansion since it happens at the inside and edge regions at the same rate, which is slightly different from traditional diffusion where edge always diffuse faster because it has a larger gradient.

From the density plot, we see a dense ring is formed at the edge of the total field strength plot, where the pressure gradient is the greatest. And the material flow is the fastest at the ring.
The material flow is interesting in itself that it reverts its direction during the evolution: during early phase of the evolution, the material is flowing into the strong field region, which is identical to figure 1: the inflow is there to compensate the pressure loss. Then, however, at time t = 0.1, it reverts its direction to point outwards. The material flow in this situation is subtle since it depends on which region is losing pressure faster. In the force-free configuration, the field cancellation caused by the reconnection is happening not just inside the strong field region but also at the edge. The reversed flow may indicate that the region outside the ring is losing pressure faster comparing to the region inside the ring at time t = 0.1. The temperature plot has no surprise that we also get heating where field pressure is dropping.

We next compare the former case to a weak field case: beta = 0.2 and beta = 20. In the beta = 20 situation, the velocity relation is vr > cs > va. The following plot shows comparison of field distribution, temperature and density respectively. Obviously the high beta case does not have a dense ring and the density does not redistribute much during the evolution. Another difference is that the central heating is not obvious in the high beta case, which is similar to the comparison we did before: the pressure loss in the high beta case is not significant enough to trigger the heating.

Finally we compare the first simulation to a case where the sound speed is comparable to the diffusion speed: vr = cs > va. This is achieved by raising the uniform temperature initially but keep everything else the same. The following graph shows the result. We can see that there is little difference between the hot case and the high beta case in terms of ring formation and central heating: both cases does not have ring formation and significant central heating. But in the hot case, the field is dissipating faster comparing to the original situation (See row 1).

===================================

3-D Tangled Field Cloud Evolution

===================================

In the 3-D tangled field case, the field will exert force on the cloud material because it is not force-free, and at the same time, it will diminish due to magnetic reconnection as a result of field diffusion.

===================================

Viscosity in AstroBEAR

===================================

The viscosity term in the Navier-Stokes Equation can be separated as:

dv/dt = \nabla \cdot T
For the applications we are considering, the fluid is isotropic. Under this assumption, the tensor can be expressed as the following combination of two scalar dynamic viscosities \mu and \mu':

T = 2 \mu \nabla v + \mu' \nabla \cdot v I

where I is the unit tensor.
It can be proven that the stress tensor should have zero trace, which leads to the following relation between \mu and \mu':

\mu' = -2/3\mu

Combining these relations, we can write down the final form of the separated viscous equation:

dv/dt=\nabla \cdot (2 \mu \nabla v) - 2/3 \nabla (\mu \nabla \cdot v)

Notice that this equation already counts the effect of all the elements in a viscosity tensor, assuming the fluid itself is isotropic.

===================================

Viscosity Implementation

===================================

The viscous source term is calculated rigorously by first finding the tensor elements at the surrounding face centers of a cell center considered: The \mu sources 2\mu \nabla v are naturally located at the cell centers, the \mu' sources however needs an averaging to obtain the cell centered velocity divergence.
Therefore the scheme requires a 20 point stencil in 3-D configuration, which makes the implicit scheme a bit complicated. The following graph shows such a stencil in 2-D setting:



Using this stencil setting, one can get face centered velocity divergence for the cell considered, and then calculate the gradient components located at the cell center.
In the viscous plasma context, the physical viscosity can be calculated rigorously. The plasma viscosity is proportional to the product of ion mean free path and the ion thermal energy:

\mu \sim \lambda_{ifp}v_{thermal}

In our code, we offer the option to calculate the realistic viscosity using the following equation:

\mu = \frac{m_{i}}{3 \sqrt{2} \pi d_{i}^2} (\frac{8kT}{\pi m_{i}})^{\frac{1}{2}}

The ion mass and local temperature should be readily available, but a good estimate of ion diameter d_{i} usually requires look up a ionic radius table.
As an estimate, one can use that of hydrogen, about 0.1~0.2 nm.
One can see that the viscosity itself is anisotropic and depends on the local temperature.
That is why we have both the \mu and \mu' terms in our viscosity equation, and have to calculate face centered source terms before applying the final divergence.

The viscosity interface is as follows:
The viscosity source term can be turned on in astrobear by putting the viscosity flag to true:
lViscous = .true.
The second parameter VisType determines which type of viscosity one requires: 1 for constant viscosity; 2 for real physical viscosity.
The third parameter viscosity is the user defined constant viscosity.
The viscosity module is tested to be working compatibly with implicit thermal conduction module and uniform gravity module.

=======================

Code Tests

=======================

The 2-D Harris current sheet test can be set up as follows:

B_y(x) = b_0 tanh(x/a)

where in computational units we choose b_0 = 1 and a = 0.5.
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_c

where \rho_0 = 1 and a = 0.5. 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 state is in pressure equilibrium though unstable. We put in a perturbation:

B_x = B_p sin(kx)cos(ky)
B_y = -B_p cos(kx)sin(ky)

where \B_p is the perturbation amplitude and k = 2\pi/L_y. The magnetic field and the density profiles are plotted below:



The magnetic energy is diminished at the X point which will convert to bulk ram pressure that drives the y direction outflow around the x = 0 line.
For the outflow pattern, the leading edge is bounded by intermediate shocks and the trailing edge is bounded by slow shocks(Petschek-type). There are totally four slow mode shocks stemming from the Sweet-Parker region, as shown in the following figure.



The time evolution model can be given by the Semenov solution, which is obtained by assuming a large resistivity box at the stagnation point. See the following figure.



The flow pattern for our setup is shown in the following plot:



The 2-D mach number plot is shown below:

The following plots show the AstroBEAR test of various settings in the current sheet outflow problem.



The Workman setting has central symmetry while the reconnection sites locate at the two edges. The positive radius of curvature of field lines at the reconnection site results in an inflow to the X point.
The first panel shows its kinetic energy distribution.
The Semenov setting is achieved by applying a resistivity function that varies with position. This resistivity gradient results in bending of magnetic field lines with negative radius of curvature. A vertical outflow from the X point is thus created. The resulted Sweet-Parker region is expanding throughout the simulation till the flow speed at its boundary reaches the Alfven speed. Panels 2 and 3 shows the kinetic energy and vertical flow momentum. Panel 4 shows the kinetic energy of Semenov outflow with an enhanced resistivity contrary.

Other possible tests:


(1) Flux Explusion Problem (Analytic 2D, resistivity only)


This problem studies the field redistribution when there is a rotating rigid cylinder with finite conductivity is placed in the flow. Can be found in text books.

(2) Hartman Flow (Analytic 2D, involves viscosity and resistivity)


This problem solves the state of a stationary flow of viscose and resistive fluid between two plates with superimposed transverse magnetic field. Can be found in many MHD text books and papers. See for instance: ​https://www.pas.rochester.edu/~shuleli/hartmann_apj08.pdf

(3) Hydromagnetic Rayleigh-Bernard Convection (2D or 3D. Involves thermal diffusion, viscosity, resistivity)


The stabilizing effect of vertical magnetic field on a classical Rayleigh-Bernard convection problem. This problem involves find the critical Chandrasekhar number with a given Rayleigh number. For a thermal convecting flow at fixed Rayleigh number, the applied vertical non-dimensional magnetic field exceeds the critical Chandrasekhar number. For details, see for instance, Chandrasekhar's 1961 book. Also experiments. see for instance: ​https://www.pas.rochester.edu/~shuleli/RB_pre00.pdf

The following LLNL journal contains useful tests on resistive MHD solver:

​https://e-reports-ext.llnl.gov/pdf/374198.pdf

Tests The stabilization effect of RT by constant viscosity and constant thermal conduction.

The Harris Sheet Test Simulation:

Click the following links for animations of kinetic energy and magnetic mach.
​GIF:Kinetic Energy​GIF:Magnetic Mach