Model 1D Vlasov Equation

Question

Vlasov Equation

The non-relativistic form of the Vlasov equation is given by: $$ \partial_{t} f\left( \mathbf{x}, \mathbf{v}, t \right) + \mathbf{v} \cdot \nabla f\left( \mathbf{x}, \mathbf{v}, t \right) + \frac{ q }{ m } \left[ \mathbf{E} + \mathbf{v} \times \mathbf{B} \right] \cdot \partial_{\mathbf{v}} f\left( \mathbf{x}, \mathbf{v}, t \right) = 0 \tag{0} $$ where $f\left( \mathbf{x}, \mathbf{v}, t \right)$ is the one-particle velocity distribution function (VDF) (e.g., see discussion at https://physics.stackexchange.com/q/216819/59023) at $\mathcal{R}^{3}$ spatial ordinates $\mathbf{x}$, $\mathcal{R}^{3}$ velocity ordinates $\mathbf{v}$, and time $t$, $q$ is the particle charge, $m$ is the particle mass, $\mathbf{E}$ and $\mathbf{B}$ are the electric and magnetic fields, and $\partial_{k} = \tfrac{ \partial }{ \partial k }$ (just a shorthand, not covariant versus contravariant). Note that $\mathbf{x}$, $\mathbf{v}$, and $t$ are all abscissa of $f$, so at any given $\mathbf{x} = \mathbf{x}_{o}$ and $t = t_{o}$, there are a defined or arbitrary number of $\mathbf{v}$ ordinates describing either a discrete or continuous, respectively, VDF.

Velocity Distribution Models

Generally $f$ is modeled as some continuous function for mathematical convenience, among other things. As an example, the Maxwellian VDF at some given point and time is given by: $$ f\left( v_{x}, v_{y}, v_{z} \right) = \frac{ n_{o} }{ \pi^{3/2} \ V_{T x} \ V_{T y} \ V_{T z} } \ exp\left[ - \sum_{k} \left( \frac{ v_{k} - v_{o, k} }{ V_{T k} } \right)^{2} \right] \tag{1} $$ where $x$, $y$, $z$ refer to the coordinate basis components, $n_{o}$ is the particle number density [cm^-3], $V_{T_{j}}$ is the $j^{th}$ thermal speed (actually the most probable speed) [km s^-1], and $v_{o, j}$ is the $j^{th}$ component of the bulk drift velocity of the distribution (i.e., from the 1st velocity moment) [km s^-1].

General Statement of Problem

I am trying to numerically model a 1D Vlasov equation using a series of differential equations and initial conditions/boundary conditions using NDSolve (if something else is more efficient/appropriate, I am not opposed to changing the function).

The purpose is to perform something called Liouville mapping (e.g., see https://physics.stackexchange.com/a/177972/59023 for details on Liouville's theorem) across a magnetic field gradient, whereby one provides an initial VDF and then transports it using the Vlasov equation to some other region of space with known spatially dependent fields and some known conservation quatities.

Assumptions and ICs/BCs

We can assume the following hold for the system in question:

• $\partial_{t} f\left( \mathbf{x}, \mathbf{v}, t \right) = 0$, i.e., time stationary;
• assume all electric and magnetic fields are quasi-static such that $\mathbf{E} = -\nabla \ \phi$, where $\phi$ is some scalar potential, and $\partial_{t} \mathbf{E} = \partial_{t} \mathbf{B} = 0$;
• $\mu = \tfrac{ m \ v_{\perp}^{2} }{ 2 \ B_{o} } = C_{0} \equiv$ the single-particle magnetic moment, where $\perp$($\parallel$) is the direction with respect to the quasi-static magnetic field vector, $\mathbf{B}_{o}$, and $C_{0}$ is some constant;
• $KE = \tfrac{ m }{ 2 } \left( \sum_{k} v_{k}^{2} \right) + q \ \phi\left( \mathbf{x} \right) = C_{1} \equiv$ the single-particle kinetic energy, where $v_{k}$ is the $k^{th}$ velocity component, $q$ is the particle charge, $m$ is the particle mass, $\phi\left( \mathbf{x} \right)$ is an electrostatic potential with a spatial dependence, and $C_{1}$ is some constant;
• for simplicity sake, let $x \ = \ \parallel$ and $y \ = \ z \ = \ \perp$ such that Equation 1 reduces to a bi-Maxwellian VDF and one can choose any appropriate values for $n_{o}$, $V_{T_{\parallel}}$, and $V_{T_{\perp}}$ as is necessary;
  • I am not sure but I think it may be useful to keep $y$ and $z$ as two separate perpendicular components during the differentiation/numerical analysis then relax them to being equal afterwards.
• let us also assume that all spatial gradients exist only along the $x \ = \ \parallel$ direction;
• to further increase simplicity, assume the particles are all electrons so that $q \rightarrow -e$ and $m \rightarrow m_{e}$, where $e$ is the fundamental charge [C] and $m_{e}$ is the electron rest mass [kg];and
• for a final simplifying assumption, let us assume $\phi\left( \mathbf{x} \right) \approx \psi \ B\left( x \right)$, where $B\left( x \right)$ is the magnetic field magnitude as a function of $x$ (or the parallel direction) and $\psi$ is some constant (assume ~4 V nT^-1 as an initial value).

From observations, let us assume that the magnetic field varies as a hyperbolic tangent as: $$ B\left( x \right) = B_{i} + \left( \frac{ B_{f} - B_{i} }{ 2 } \right) \left[ 1 + tanh\left( \frac{ x }{ x_{o} } \right) \right] \tag{2} $$ where $B_{i}$($B_{f}$) is the intial(final) magnetic field magnitude and $x_{o}$ is the gradient scale length half-width of the hyperbolic tangent (assume magnetic ramp centered at zero). If necessary, one can assume that $B_{i}$ = 5 nT and $B_{f}$ = 15 nT. We can also assume that the magnetic field is only along the $x$ component, i.e., why $x$ = $\parallel$ here.

General Questions

Note that I am quite certain at least one of my issues is with syntax in using NDSolve, so please let me know if I have done something simplistically wrong in that regard. I am also inclined to think that I am confusing the abscissa/ordinate usage in the code, which is causing me headaches. I am not sure if the initial VDF should be a 2D numerical list of values at some known list of $v_{x}$ and $v_{y}$ and $v_{z}$ or if I keep it analytical and allow Mathematica to find a numerical integrating solution within the bounds of the velocity range I chose.

In short, I am trying find $f\left( v_{x}, v_{y}, v_{z} \right)$ at any arbitrary position $x$ (equivalent to the parallel direction for simplicity) given an intial form/function for $f$ at some initial position $x_{i}$.

My 1st Attempt

First, I define some physical constants and numerical values for initial and final conditions for later use:

(* Define Constants [2014 CODATA/NIST values] *)
mme = 9.1093835600 10^-31;  (* electron mass [kg] *)
cc = 2.9979245800 10^8;  (* speed of light in vacuum [m/s] *)
qq = 1.6021766208 10^-19;  (* fundamental unit of charge [C] *)
(* Define numerical velocity grid [m s^-1] *)
nva = 300;
vamax = 2 10^4 10^3;
dva = (2 vamax)/(nva - 1);
va1d = Table[i, {i, -vamax, vamax, dva}];
(* Define parallel drift speed or initial bulk flow speed of particles [m s^-1] *)
voa0 = 400 10^3;
(* Define asymptotic upstream number density [m^-3] *)
nb0 = 10 10^6;
(* Define asymptotic upstream thermal speeds [m s^-1] *)
vtab0 = 1500 10^3;
(* Define asymptotic B-field magnitudes [T] *)
bup0 =  5 10^-9;
bdn0 = 15 10^-9;
(* Define shock ramp half-width [m] *)
xo = 10^2 10^3;
(* Define potential factor [V T^-1] *)
psi0 = 4*10^9;
(* Define maximum delta x [m] to use in numerical solutions *)
dxmax = xo*0.005;
(* Define maximum spatial distance from ramp [m] to consider *)
xmax = xo*10;

Now I define the relevant equations to numerically solve:

(* Define magnetic moment = constant *)
mumom = (m*(vy^2 + vz^2))/(2 (bup + ((bdn - bup)/2) (1 + Tanh[x/x0]))) == C[0];
(* Define kinetic energy = constant *)
kinen = m/2/e (vx^2 + (vy^2 + vz^2)) - psi*(bu + ((bd - bu)/2) (1 + Tanh[x/x0])) == C[1];
(* Define stationary Vlasov equation *)
statvlaseq = vx D[f[x, vx, vy, vz], x] - e/m (({-(((bd - bu) psi Sech[x/x0]^2)/(2 x0)), 0, 0} + {0, vz (bup + ((bdn - bup)/2) (1 + Tanh[x/x0])), -vy (bup + ((bdn - bup)/2) (1 + Tanh[x/x0]))}).{D[f[x, vx, vy, vz], vx], D[f[x, vx, vy, vz], vy], D[f[x, vx, vy, vz], vz]}) == 0;
(* Define initial VDF condition *)
finit = f[-xmax, vx, vy, vz] == nb0/(\[Pi]^(3/2) vtab0^3) Exp[-3*((vx - voa0)/vtab0)^2];

When I try to numerically solve this, I know I am missing several things and screwing up the syntax but am failing to see how, but here's my latest attempt:

(* Try to solve numerically *)
NDSolve[({statvlaseq, mumom, kinen, finit} /. {bu -> bup0, bd -> bdn0, x0 -> xo, psi -> psi0, e -> qq, m -> mme}), f[x, vx, vy, vz], {x, -xmax, xmax}, {vx, -vamax, vamax}, {vy, -vamax, vamax}, {vz, -vamax, vamax}, MaxStepSize -> dxmax, MaxSteps -> 10^8, AccuracyGoal -> 8, PrecisionGoal -> 8, Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 8} ]

I had originally got an error about vy appears in the head of the expression vy[x] because I had made vi = vi[x] (i.e., functions of space), but this was wrong. Since removing the spatial dependence from the velocity coordinates and now treating them as independent variables, the new error message is The equations {lots of stuff} are not differential equations or initial conditions in the dependent variables {f}. My guess is this complaint derives from the magnetic moment and total kinetic energy being constant constraint equations, mumom and kinen, respectively. Though I have not figured out how to use such constraints on independent variables while solving such a PDE. I have managed to properly wrangle NDSolve in the past when all abscissa and functions were simply functions of only time and space (space being a function of time too), e.g., at https://mathematica.stackexchange.com/a/63526/19732.

Note: The $v_{k}$ are not the total time derivatives of the spatial positions $x_{k}$, so this is not a matter of defining something like vx[t] == x'[t] or anything like that. The $v_{k}$ are part of the 6D phase space for which $f\left( \mathbf{x}, \mathbf{v}, t \right)$ is defined at any given time. In this specific example, we ignore implicit time variations.

My 2nd Attempt

I started to realize that perhaps my issues were in some of the definitions or code and perhaps I should walk before I run. As @bbgodfrey suggested, I tried to do this for a 1D electrostatic case (i.e., magnetic field dependence disappears). Again I ran into odd issues with NDSolve so I tried to solve the equations in parametric space (kind of grabbing at straws at the moment hoping for a more useful error message) but now I get a very puzzling error from ParametricNDSolveValue.

So now f[x,vx,vy,vz] reduces to just f[x,vx] and assume the same constants as above, then I tried:

(* Try to solve in parametric space *)
ClearAll[fvlas]
fvlas = ParametricNDSolveValue[
  ({
    vx (D[f[x, vx], x]) - 
      qq/mme (({-(((bdn0 - bup0) psi0 Sech[x/xo]^2)/(2 xo)), 0, 
            0} + {0, 0, 0}).{(D[f[x, vx], vx]), 0, 0}) == 0,
     f[-xmax, vx] == 
     nb0/(\[Pi]^(3/2) vtab0^3) Exp[-3*((vx - voa0)/vtab0)^2],
    f[x, -\[Infinity]] == 0, f[x, \[Infinity]] == 0,
    (mme vx^2)/qq - 
      psi0 (bup0 + 1/2 (bdn0 - bup0) (1 + Tanh[x/xo])) == c1
    }),
  f,
  {x, -xmax, xmax}, {vx, -vamax, vamax}, {c1}
  ]

I see a legitimate looking ParametricFunction output box with no errors that contains the following:

• Expression: f
• Parameters: {c1}
• Generator: ParametricNDSolveValue
• Dependent variables: {f}
• Independent variables: {x, vx}

However, when I try to follow some of the examples from the ParametricNDSolveValue documentation on plotting or showing results, e.g., I enter fvlas[-500], I get the error message A value for {c1} has already been defined as part of the equations. I am at a complete loss here as c1 is a constant and I defined it properly as a parameter within the ParametricNDSolveValue call, I think (the output seems to suggest as much). Any suggestions would be greatly appreciated.