Back to BearclawOutline

i_cosray

Location:

contrib/astro/physics/i_cosray.f90

Included In:

contrib/astro/problem.f90

Description:

This routine will handle the cosmic rays (CR) for AstroBEAR. The CR's are advected by AstroBEAR and then the diffusion, injection, and feedback physics are dealt with in this module. For details on the physics and solver read Jones & Kang, Astroparticle Physics 24 (2005) 75-91 (hereafter JK05) located ​here, a brief description of the physics and solver will be provided below. The code for this was taken from mhdode which was an application of this CR scheme to a 1D TVD mhd code. The code units for momentum are in units of m_p*c for both the protons and electrons.

Physics:

Diffusive Shock Acceleration:

Diffusive Shock Acceleration (DSA) is the current standard model for CR acceleration. The basic model is based on 1st Order Fermi Acceleration. The basic idea is equivalent to bouncing a ball off a moving wall. If one throws a ball at a stationary wall the ball comes back to them at the same speed as it left them just opposite direction. Now if the person is moving towards the wall, or equivalently the wall is moving towards to person, and throws the ball then in the person's frame the ball comes back faster than it left them. The ball actually gains energy from the wall.

If one now considers the ball to be a charged particle and the wall to be a magnetic field that is in motion towards the particle which is also moving towards the magnetic field then we can have the particle gain energy from reflecting/scattering off of the magnetic field. Now imagine that we are in a shock. If one has a particle in a shock as the particle moves upstream of the flow it will see material and magnetic field coming towards it and will back scatter. Then as it heads back downstream across the shock it will again see material and magnetic field coming towards it and will scatter again. Each time it back scatters and crosses the shock the particle gains energy from the shock. This turns out to follow a diffusion process which gives rise to the term DSA.

Naturally this is rather simplified but the beauty of DSA is that if you give a particle a reasonable diffusion coefficient and escape chance it will actually inherently produce a power law spectrum of particles if you send through a large number of particles. This is exactly what is seen of particle spectra in nature. This has led DSA to be the leading theory for particle acceleration.

The favored diffusion coefficient is what is called Bohm Diffusion. Remember that your standard diffusion coefficient has the form of kappa=1/3*v*lambda, where kappa is the diffusion coefficient, v is the speed of the particle and lambda is the characteristic length over which a particle will backscatter. We are dealing with relativistic particles so v=c, where c is the speed of light. We can make the reasonable assumption that a particle's characteristic scattering length will be equal to the gyroradius (r_g) of the particle. This then gives a diffusion coefficient of kappa=1/3*r_g*c=(m*c3)/(3*e*B)*p2/(p2+1)(½) where m is the particle mass, e is the particle charge, B is the magnetic field perpendicular to the particle's motion and p is the particle's momentum in units of m*c. Notice that this depends on the particle momentum which means that electrons and protons of the same momentum behave the same way in DSA.

Notice that it is the magnetic field that is perpendicular to the particle's motion that is important not the one that is parallel. Typically one assumes that for propagation parallel to the magnetic field that the CR's streaming along the magnetic field lines excite Alfven waves with characteristic amplitude of the parallel magnetic field strength and a characteristic wavelength of the particles gyroradius. This allows one to use the parallel magnetic field value for the diffusion parallel to the magnetic field. Diffusion perpendicular is difficult for the particles and does not excite Alfven waves which leads to a different diffusion coefficient for this direction as Bohm diffusion may not be appropriate anymore. This is where anisotropic diffusion comes in which can allow for the difference in diffusing along a magnetic field line as opposed to perpendicular to it.

Injection:

One may wonder how does a particle get started in this magical DSA process? This is where injection comes in. While injection of electrons isn't well known, the injection of protons is fairly well understood. For protons the particles from the suprathermal tail of the Maxwellian distribution downstream of the shock have a chance backscatter and travel upstream of the shock where they can enter into the DSA regime. This is because the suprathermal tail contains particles that are travelling faster than the post shock flow and therefore can propagate upstream and cross the shock. Injection is thought to be easier for a parallel shock, where the shock normal is parallel to the magnetic field as opposed to when it is perpendicular to the shock normal. This is because it is easier for the particle to propagate upstream along a field line than perpendicular to it.

Cosmic Ray Feedback:

Being relativistic particles CR's have a ton of energy. Even a small number of CR's can contain a non-trivial amount of energy. This energy manifests itself as an addition pressure that modifies the flow. This CR pressure effect is most prominent near shocks where the majority of CR's are. Shocks that have significant CR components can be modified by the additional pressure. What develops is what is called a shock precursor. The actual shock (called a sub-shock) is proceeded by a pressure ramp that causes the gas to compress. This compressed gas from the precursor then passes through the subshock and is compressed even further. Heavily modified shocks can have total compression ratios well in excess of the factor of 4 one gets from a normal strong gasdynamic shock.

Energy Losses:

While protons being very massive do not lose significant energy to radiation, electrons do via synchrotron and inverse Compton processes. These energy losses apply everywhere and therefore effect the acceleration dynamics as well as the spatial distribution of the electrons, particularly at high energy. This can heavily modify the electron spectrum creating a high energy cut off based on the severity of the losses and the location relative to the acceleration region.

Alfven Wave Transport:

As mentioned before the CR protons streaming ahead of the shock produce Alfven waves. These Alfven waves contribute two effects. The first is to shift the velocity of the scattering center for the CR's by the Alfven velocity. Second the Alfven waves heat the gas as they dissipate.

Solver:

Coarse-Grained Momentum finite Volume

The Coarse-Grained Momentum finite Volume (CGMV) method is at the core of the CR solver. This method assumes that the particle spectrum will take a powerlaw form, which is true for the vast majority of CR particle spectra in nature. Using this assumption one constructs two moments of the CR distribution, n(p)=p3*f(p) where f(p) is the normal distribution function and g(p)=p4*f(p). By tracking the evolution of these two moments and utilizing the powerlaw assumption one can solve the CR momentum convection problem (the acceleration of the particles to higher energy) using large momentum bins, in logarithmic space this would be a dy~1 where y=ln(p).

CR Variable Storage:

The relevant CR information is assumed to be stored in Info%q in indices above where the tracers are stored. The first CR indice is Pc which is only relevant when CR feedback is being used. Next are the proton and electron moments with each indice being a momentum bin. If CR feedback is off then the first set of moments are the CR tracer particles (usually electrons), if the CR feedback is on then the first set of moments are the protons and the second (if used) are the electrons. For example, let's say that we have 2 tracer variables and we want 2 species of CRs with 3 momentum bins. So we would set nplow to 11, np=3, and nj=2. The Info%q would look like this in structure:

Info%q(1-8)=MHD state variables

Info%q(9-10)=Tracer variables

Info%q(11)=CR Pressure

Info%q(12-14)=n moment of CR proton distribution

Info%q(15-17)=g moment of CR proton distribution

Info%q(18-20)=n moment of CR electron distribution

Info%q(21-23)=g moment of CR electron distribution

Info%q variables 9-23 are all considered as tracers to AstroBEAR, which means one would need to set NrTracerVars=15. One needs to make sure the iLTracer=2 for the CR's to work properly, this also means that any other tracers you have will also be treated with iLTracer=2.

Crank-Nicolson Scheme

Whereas the actual particle acceleration (momentum convection) is handled by the CGMV scheme the full DSA problem is really a solution to the diffusion-convection equation. As such we need to solve the particle diffusion problem as well. On a fixed grid with no subdomains or subgrids the Crank-Nicholson scheme using a tridiagonal solver, which will be referred to as just the CN scheme, is very fast and is unconditionally stable. However to maintain a converged solution dx < 0.1*x_d, where x_d is the relevant particle's diffusion length.

With subdomains and subgrids (AMR patches), the CN scheme develops unphysical answers owing to the size of the subgrid to the diffusion length and the lack of information about the rest of the grid. To solve this two fixes were applied to the CN scheme to minimize these pathologies. The first fix was to introduce adaptive subcycling into the CN scheme. An initial guess is made by the solver as to how many subcylces the problem needs, which is assumed to be twice the maximum diffusion length on the subgrid for the momentum being considered divided by the length of the subgrid. The solver then subcycles the grid that many times by dividing the timestep by the number of subcycles. During the subcycling the distribution is checked for unphysical answers, if they are found the solver aborts the current subcycling, multiples the current number of subcycles by a factor (the code currently has 10), then it restarts the subcylcing with the new number of subcycles. This continues till a good answer is reached.

The other fix is to reconstruct the subgrid boundaries such that they are continuous with the distribution on the subgrid but can still propagate boundary information to the subgrid. During each subcycle the solver takes the boundary zone furthest from the physical grid and uses that plus the first physical zone on the grid to construct a linear interpolation for the other zones on the grid. This ensures that the function is at least 0th order continuous across the boundary into the physical grid. For example let's say that we have 4 boundary zones and we are at the left boundary of an x-pass, which gives us -3, -2, -1, and 0 for the boundary zones. The boundary information at -3 is left alone and used in conjunction with the information at 1, the first zone on the grid, to construct a linear fit to the boundary. Then -2, -1, and 0 are interpolated based on this fit.

Both of these fixes allow the scheme to handle AMR patches that are smaller than the diffusion length. Caution must be used though as the subcycling can make the CR solver take quite a bit of time. Getting maximum performance requires a balance between the size of the AMR patches and the amount of subcylcing required for that patch.

Subroutines:

afterstepCR:

Main driver for CR solver. Call to afterstepCR located in afterfixup. This routine works in a dimensionally split way by taking 1-D stencils of the grid and solving the CR problem. It handles the Alfven transport, CR feedback and diffusion coefficient calculations. It also applies protection routines to make sure that the CR distribution is not below or equal to zero. This routine calls diffus, pconv and inject.

diffus:

Solves spatial diffusion using the Crank-Nicholson (CN) scheme with a tridiagonal solver. Subcycling is momentum dependant. This subroutine also includes energy losses for electrons and the application of the results of pconv and inject subroutines. diffus is called by afterstepCR.

pconv:

Solves momentum advection in compressive flow using Coarse-Grained Momentum finite Volume method. Whereas the actual particle advection is handled by AstroBEAR the compressive effects on the particles are handled here. Compressive effects include adiabatic losses, energy gains due to compressions, and density increases/decreases due to compressions. This routine is called by afterstepCR.

inject:

Include the effects of the injection of particles at the shock. The shock detector works along 1-D slices. This is the "old", circa 1990 injection scheme, meaning that this does not include a thermal injection model but rather injects a set fraction of thermal particles into the CR bins. This routine is called by afterstepCR

cosrayErrFlag:

Sets up the AMR flagging needed for the CR diffusion scheme to provide accurate solutions. This step is called by problemErrFlag. Flagging is done via finding shocks and gradients in the CR flow. To reduce number of loops needed to flag the grid a list is built up of rectangular regions that need to be refined. These rectangular regions are defined by a center location of the rectangle and the x and y zone length of the rectangle. The flgs array stores the following in flgs(i,j), i is the number of the rectangular region, j=1 is the x location of the rectangular region center, j=2 is y location, j=3 is half the length of the rectangular region in the x direction, j=4 is half the length in the y direction. The length of a region is characterized by the largest relevant diffusion length in that region. After this list is built up, a search routine is used to figure out if a particular zone lies in one of these regions. If it does then it is refined.

Questions?

Contact Paul Edmon (​pedmon@astro.umn.edu).