1D & 2D Models

In a 1D model, the Fokker-Planck form of the transport is equivalent to the SDEs

\[\begin{split}dX & = \left(\frac{\partial\kappa}{\partial x} + V_x\right)ds + \sqrt{2\kappa} dW_\sigma(s),\\ dp & = -\frac{p}{3}\frac{\partial V_x}{\partial x}ds,\end{split}\]

Note

In a 1D or 2D model, the particle drift \(\boldsymbol{V}_d\) is out of the simulation plane and is not considered here.

where \(dW_\sigma\) is the normalized distributed random number with mean zero and variance \(\sqrt{\Delta t}\), and \(\Delta t\) is the time step for stochastic integration. This corresponds to a Wiener process. Numerical approximation is often used for the Wiener process to replace the normal distribution. We use a uniform distribution in \([-\sqrt{3}, \sqrt{3}]\) in the code.

Note

Future tests should be included to test this method thoroughly.

In a 2D model, the corresponding SDEs are

\[\begin{split}dX & = (\nabla\cdot\boldsymbol{\kappa} + \boldsymbol{V})ds + \sum_\sigma\boldsymbol{\alpha}_\sigma dW_\sigma(s),\\ dp & = -\frac{p}{3}(\nabla\cdot\boldsymbol{V})ds,\end{split}\]

where \(\sum_\sigma\alpha_\sigma^\mu\alpha_\sigma^\nu = 2\kappa^{\mu\nu}\).

\[\begin{split}\boldsymbol{\alpha}_1 = \begin{pmatrix} \sqrt{2\kappa_\perp} \\ 0 \end{pmatrix}, \quad \boldsymbol{\alpha}_2 = \begin{pmatrix} 0 \\ \sqrt{2\kappa_\perp} \end{pmatrix}, \quad \boldsymbol{\alpha}_3 = \sqrt{2(\kappa_\parallel - \kappa_\perp)} \begin{pmatrix} B_x/B \\ B_y/B \end{pmatrix}.\end{split}\]

The parameters used at particle locations are calculated from \(v_x\), \(v_y\), \(B_x\), \(B_y\), \(\nabla\cdot\boldsymbol{v}\), \(\partial B_x/\partial x\), \(\partial B_x/\partial y\), \(\partial B_y/\partial x\), and \(\partial B_y/\partial y\), which are all obtained from the MHD simulations. We interpolate these parameters to the particle positions and then calculate the other required parameters:

\[\begin{split}\begin{aligned} \frac{\partial\kappa_{xx}}{\partial x} & = \frac{\partial\kappa_\perp}{\partial x} - \frac{\partial(\kappa_\perp-\kappa_\parallel)}{\partial x}\frac{B_x^2}{B^2} - 2(\kappa_\perp-\kappa_\parallel)\frac{\frac{\partial B_x}{\partial x}B_xB- \frac{\partial B}{\partial x}B_x^2}{B^3}, \\ \frac{\partial\kappa_{yy}}{\partial y} & = \frac{\partial\kappa_\perp}{\partial y} - \frac{\partial(\kappa_\perp-\kappa_\parallel)}{\partial y}\frac{B_y^2}{B^2} - 2(\kappa_\perp-\kappa_\parallel)\frac{\frac{\partial B_y}{\partial y}B_yB- \frac{\partial B}{\partial y}B_y^2}{B^3}, \\ \frac{\partial\kappa_{xy}}{\partial x} & = -\frac{\partial(\kappa_\perp-\kappa_\parallel)}{\partial x} \frac{B_xB_y}{B^2} - (\kappa_\perp-\kappa_\parallel) \frac{\left(\frac{\partial B_x}{\partial x}B_y+ B_x\frac{\partial B_y}{\partial x}\right)B - 2B_xB_y\frac{\partial B}{\partial x}}{B^3}, \\ \frac{\partial\kappa_{xy}}{\partial y} & = -\frac{\partial(\kappa_\perp-\kappa_\parallel)}{\partial y} \frac{B_xB_y}{B^2} - (\kappa_\perp-\kappa_\parallel) \frac{\left(\frac{\partial B_x}{\partial y}B_y+ B_x\frac{\partial B_y}{\partial y}\right)B - 2B_xB_y\frac{\partial B}{\partial y}}{B^3}, \\ \frac{\partial B}{\partial x} & = \frac{1}{B}\left(B_x \frac{\partial B_x}{\partial x} + B_y\frac{\partial B_y}{\partial x}\right), \\ \frac{\partial B}{\partial y} & = \frac{1}{B}\left(B_x\frac{\partial B_x}{\partial y} + B_y\frac{\partial B_y}{\partial y}\right). \end{aligned}\end{split}\]

where \(\kappa_\parallel\) and \(\kappa_\perp\) can be functions of \(B_x\), \(B_y\) and \(B\), so \(\partial \kappa_\parallel/\partial x\), \(\partial \kappa_\parallel/\partial y\), \(\partial \kappa_\perp/\partial x\), and \(\partial \kappa_\perp/\partial y\) still depend on the derivatives \(\partial B_x/\partial x\), \(\partial B_x/\partial y\), \(\partial B_y/\partial x\), and \(\partial B_y/\partial y\). The detailed expressions depend on the diffusion model to choose. Considering the expression of \(\kappa_\parallel\), we get

\[\begin{aligned} \frac{\partial\kappa}{\partial x}\sim\kappa\left( \frac{2}{3L_c}\frac{\partial L_c}{\partial x} - \frac{1}{3B}\frac{\partial B}{\partial x} - \frac{1}{\sigma^2}\frac{\partial(\sigma^2)}{\partial x} \right) \end{aligned}\]

Note

We typically assume \(L_c\) and \(\sigma^2\) are spatially independent. However, they could vary spatially in a large system.

Time step criteria

For a 1D problem, the particle moves a distance satisfying \(l_x^2=\text{max}\left(\left<\Delta x\right>^2, \left<\Delta x^2\right>\right)\) [Strauss17], where

\[\begin{aligned} \left<\Delta x\right> = \left(V_x + \frac{d\kappa(x)}{dx}\right)\Delta t, \quad \left<\Delta x^2\right> = 2\kappa(x)\Delta t, \end{aligned}\]

and \(l_x\) should be much smaller than the spatial variation scale of the fields. In this code, we require \(\left<\Delta x\right>^2 < \left<\Delta x^2\right>\) and choose \(\Delta t\) so that \(l_x < \delta_x\), where \(\delta_x\) is the grid size. For the 2D problems, we choose the following criteria to determine the time step:

\[\begin{split}\begin{aligned} \Delta t_x & = \text{min}\left[\frac{(0.5\delta x)^2}{2\kappa_\parallel}, \frac{2\kappa_\perp} {(V_x + \partial_x\kappa_{xx} + \partial_y\kappa_{xy})^2}\right], \\ \Delta t_y & = \text{min}\left[\frac{(0.5\delta y)^2}{2\kappa_\parallel}, \frac{2\kappa_\perp}{(V_y + \partial_y\kappa_{yy} + \partial_x\kappa_{xy})^2}\right],\\ \Delta t & = \text{min}(\Delta t_x, \Delta t_y). \end{aligned}\end{split}\]

Higher-order method

To get a higher-order solution, we can use a derivative-free Milstein method [Burrage04] to solve the SDEs. It is different from the usual method due to one more term, which makes it a higher-order method.

Note

This method was used by [Li18] but has since not been supported in default.

\[\begin{split}\begin{aligned} dX_t & = f(X_t,t)dt + g(X_t,t)dW_t, \\ X_{n+1} & = X_n + f_n h + g_n\Delta W_n + \frac{1}{2\sqrt{h}}[g(\bar{X}_n)-g_n][(\Delta W_n)^2-h], \\ \bar{X}_n & = X_n + f_n h + g_n\sqrt{h}, \\ \Delta W_n & = [W_{t+h}-W_t] \sim \sqrt{h}N(0,1), \end{aligned}\end{split}\]

where \(X\) corresponds to spatial positions \(x\), \(y\) and particle momentum \(p\) in our simulation. Here \(f(X_t,t)\) is the deterministic term, \(g(X_t,t)\) is the probabilistic term, \(h\) is the time step, and \(N(0,1)\) indicates a normal distribution, which is substituted with a uniform distribution \([-\sqrt{3}, \sqrt{3}]\) in our simulations to speed up the computation.

[Burrage04]

Burrage, K., Burrage, P.M. and Tian, T., 2004. Numerical methods for strong solutions of stochastic differential equations: an overview. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 460(2041), pp.373-402.

[Li18]

Large-scale Compression Acceleration during Magnetic Reconnection in a Low-β Plasma, Xiaocan Li, Fan Guo, Hui Li, and Shengtai Li, The Astrophysical Journal 866, no. 1 (2018): 4.

[Strauss17]

Strauss, R. and Effenberger, F., 2017. A hitch-hiker’s guide to stochastic differential equations. Space Science Reviews, 212(1), pp.151-192.