3D Model
The relationship
\(\sum_\sigma\alpha_\sigma^\mu\alpha_\sigma^\nu = 2\kappa^{\mu\nu}\)
is actually a matrix decomposition. We need to decompose
\(2\kappa=PP^T\), where
\(P=(\boldsymbol{\alpha}_1, \boldsymbol{\alpha}_2, \boldsymbol{\alpha}_3)\).
In a 2D problem, the third component of \(\boldsymbol{\alpha}_i\) is
essentially 0. In a 3D problem, we need to find all three components of
\(\boldsymbol{\alpha}_i\). We need some linear algebra for that.
Every real symmetric matrix can be written in the form
(https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix#Real_symmetric_matrices)
\[A=Q\Lambda Q^T\]
where \(Q\) is an orthogonal matrix whose columns are the
eigenvectors of \(A\), and \(\Lambda\) is a diagonal matrix
whose entries are the eigenvalues of \(A\). If the eigenvalues are
non-negative, then the real matrix \(P=Q\Lambda^{1/2}\), and
\[A=Q\Lambda^{1/2}\Lambda^{1/2}Q^T = \frac{PP^T}{2}\]
According to WolframAlpha, the eigenvalue of \(\kappa\) is
\(k_\parallel\), \(k_\perp\), and \(k_\perp\), and the
corresponding eigenvectors are
\[\begin{split}\begin{aligned}
v_1 = \left(\frac{b_x}{b_z}, \frac{b_y}{b_z}, 1\right), \\
v_2 = \left(-\frac{b_z}{b_x}, 0, 1\right), \\
v_3 = \left(-\frac{b_y}{b_x}, 1, 0\right).
\end{aligned}\end{split}\]
where \(b_x=B_x/B\), \(b_y=B_y/B\), and \(b_z=B_z/B\).
\(v_1\), \(v_2\), and \(v_3\) are not unit vectors, and
\(v_2\) and \(v_3\) are not orthogonal to \(v_1\), so we
need to re-organize \(v_2\) and \(v_3\) and normalize the
vectors.
\[\begin{split}\begin{aligned}
v_1 & = \left(b_x, b_y, b_z\right), \\
v_2 & = \left(-\frac{b_xb_z}{\sqrt{b_x^2+b_y^2}},
-\frac{b_yb_z}{\sqrt{b_x^2+b_y^2}}, \sqrt{b_x^2+b_y^2}\right),\\
v_3 & = \left(-\frac{b_y}{\sqrt{b_x^2+b_y^2}}, \frac{b_x}{\sqrt{b_x^2+b_y^2}}, 0\right)
\end{aligned}\end{split}\]
where \(v_2\) is calculated from the perpendicular component of the
original \(v_2\) w.r.t. \(v_3\). Then,
\[\begin{split}Q =
\begin{pmatrix}
b_x & -b_xb_z/\sqrt{b_x^2+b_y^2} & -b_y/\sqrt{b_x^2+b_y^2}\\
b_y & -b_yb_z/\sqrt{b_x^2+b_y^2} & b_x/\sqrt{b_x^2+b_y^2}\\
b_z & \sqrt{b_x^2+b_y^2} & 0
\end{pmatrix}\end{split}\]
\[\begin{split}\Lambda =
\begin{pmatrix}
\kappa_\parallel & 0 & 0\\
0 & \kappa_\perp & 0 \\
0 & 0 & \kappa_\perp
\end{pmatrix}\end{split}\]
\[\begin{split}P = \sqrt{2}Q\Lambda^{1/2} =
\begin{pmatrix}
b_x\sqrt{2\kappa_\parallel} & -b_xb_z\sqrt{2\kappa_\perp}/\sqrt{b_x^2+b_y^2} &
-b_y\sqrt{2\kappa_\perp}/\sqrt{b_x^2+b_y^2}\\
b_y\sqrt{2\kappa_\parallel} & -b_yb_z\sqrt{2\kappa_\perp}/\sqrt{b_x^2+b_y^2} &
b_x\sqrt{2\kappa_\perp}/\sqrt{b_x^2+b_y^2}\\
b_z\sqrt{2\kappa_\parallel} & \sqrt{b_x^2+b_y^2}\sqrt{2\kappa_\perp} & 0
\end{pmatrix}\end{split}\]
We can verify that \(PP^T=2\kappa\). For 3D simulation, we need to
calculate more terms of the gradient of the diffusion tensor. The
parameters used at particle locations are calculated from \(V_x\),
\(V_y\), \(V_z\), \(b_x\), \(b_y\), \(b_z\),
\(\nabla\cdot\boldsymbol{V}\), \(\partial_x b_x\),
\(\partial_y b_x\), \(\partial_z b_x\), \(\partial_x b_y\),
\(\partial_y b_y\), \(\partial_z b_y\), \(\partial_x b_z\),
\(\partial_y b_z\), \(\partial_z b_z\).
\[\begin{split}\begin{aligned}
\partial_x\kappa_{xx} & = \partial_x\kappa_\perp +
\partial_x(\kappa_\parallel-\kappa_\perp)b_x^2 +
2(\kappa_\parallel-\kappa_\perp)b_x\partial_xb_x, \\
\partial_y\kappa_{yy} & = \partial_y\kappa_\perp +
\partial_y(\kappa_\parallel-\kappa_\perp)b_y^2 +
2(\kappa_\parallel-\kappa_\perp)b_y\partial_yb_y, \\
\partial_z\kappa_{zz} & = \partial_z\kappa_\perp +
\partial_z(\kappa_\parallel-\kappa_\perp)b_z^2 +
2(\kappa_\parallel-\kappa_\perp)b_z\partial_zb_z, \\
\partial_x\kappa_{xy} & =
\partial_x(\kappa_\parallel-\kappa_\perp)b_xb_y +
(\kappa_\parallel-\kappa_\perp)(\partial_xb_xb_y + b_x\partial_xb_y), \\
\partial_y\kappa_{xy} & =
\partial_y(\kappa_\parallel-\kappa_\perp)b_xb_y +
(\kappa_\parallel-\kappa_\perp)(\partial_yb_xb_y + b_x\partial_yb_y), \\
\partial_x\kappa_{xz} & =
\partial_x(\kappa_\parallel-\kappa_\perp)b_xb_z +
(\kappa_\parallel-\kappa_\perp)(\partial_xb_xb_z + b_x\partial_xb_z), \\
\partial_z\kappa_{xz} & =
\partial_z(\kappa_\parallel-\kappa_\perp)b_xb_z +
(\kappa_\parallel-\kappa_\perp)(\partial_zb_xb_z + b_x\partial_zb_z), \\
\partial_y\kappa_{yz} & =
\partial_y(\kappa_\parallel-\kappa_\perp)b_yb_z +
(\kappa_\parallel-\kappa_\perp)(\partial_yb_yb_z + b_y\partial_yb_z), \\
\partial_z\kappa_{yz} & =
\partial_z(\kappa_\parallel-\kappa_\perp)b_yb_z +
(\kappa_\parallel-\kappa_\perp)(\partial_zb_yb_z + b_y\partial_zb_z)
\end{aligned}\end{split}\]
Or we may prefer to use current code structure that calculates
\(\partial_x B_x\), \(\partial_y B_x\), \(\partial_z B_x\),
\(\partial_x B_y\), \(\partial_y B_y\), \(\partial_z B_y\),
\(\partial_x B_z\), \(\partial_y B_z\), \(\partial_z B_z\).
Then, the derivatives are calculated as
\[\begin{split}\begin{aligned}
\partial_xB & = b_x\partial_xB_x + b_y\partial_xB_y + b_z\partial_xB_z, \\
\partial_yB & = b_x\partial_yB_x + b_y\partial_yB_y + b_z\partial_yB_z, \\
\partial_zB & = b_x\partial_zB_x + b_y\partial_zB_y + b_z\partial_zB_z, \\
\partial_x\kappa_{xx} & = \partial_x\kappa_\perp +
\partial_x(\kappa_\parallel-\kappa_\perp)b_x^2 +
2(\kappa_\parallel-\kappa_\perp)(B_xB\partial_xB_x - B_x^2\partial_x B)/B^3, \\
\partial_y\kappa_{yy} & = \partial_y\kappa_\perp +
\partial_y(\kappa_\parallel-\kappa_\perp)b_y^2 +
2(\kappa_\parallel-\kappa_\perp)(B_yB\partial_yB_y - B_y^2\partial_y B)/B^3, \\
\partial_z\kappa_{zz} & = \partial_z\kappa_\perp +
\partial_z(\kappa_\parallel-\kappa_\perp)b_z^2 +
2(\kappa_\parallel-\kappa_\perp)(B_zB\partial_zB_z - B_z^2\partial_z B)/B^3, \\
\partial_x\kappa_{xy} & = \partial_x(\kappa_\parallel-\kappa_\perp)b_xb_y +
(\kappa_\parallel-\kappa_\perp)[(B_y\partial_xB_x + B_x\partial_xB_y)B -
2B_xB_y\partial_xB] / B^3, \\
\partial_y\kappa_{xy} & = \partial_y(\kappa_\parallel-\kappa_\perp)b_xb_y +
(\kappa_\parallel-\kappa_\perp)[(B_y\partial_yB_x + B_x\partial_yB_y)B -
2B_xB_y\partial_yB] / B^3, \\
\partial_x\kappa_{xz} & = \partial_x(\kappa_\parallel-\kappa_\perp)b_xb_z +
(\kappa_\parallel-\kappa_\perp)[(B_z\partial_xB_x + B_x\partial_xB_z)B -
2B_xB_z\partial_xB] / B^3, \\
\partial_z\kappa_{xz} & = \partial_z(\kappa_\parallel-\kappa_\perp)b_xb_z +
(\kappa_\parallel-\kappa_\perp)[(B_z\partial_zB_x + B_x\partial_zB_z)B -
2B_xB_z\partial_zB] / B^3, \\
\partial_y\kappa_{yz} & = \partial_y(\kappa_\parallel-\kappa_\perp)b_yb_z +
(\kappa_\parallel-\kappa_\perp)[(B_z\partial_yB_y + B_y\partial_yB_z)B -
2B_yB_z\partial_yB] / B^3, \\
\partial_z\kappa_{yz} & = \partial_z(\kappa_\parallel-\kappa_\perp)b_yb_z +
(\kappa_\parallel-\kappa_\perp)[(B_z\partial_zB_y + B_y\partial_zB_z)B -
2B_yB_z\partial_zB] / B^3.
\end{aligned}\end{split}\]
Particle drift velocity
In the 3D model, we need the drift velocity, which is given by
\[\begin{split}\begin{aligned}
& \boldsymbol{V}_d = \frac{pcw}{3q}\nabla\times\left(\frac{\boldsymbol{B}}{B^2}\right)
= \frac{1}{3q}\frac{p^2c}{\sqrt{p^2+m^2c^2}}
\left(\frac{1}{B^2}\nabla\times\boldsymbol{B} -
\frac{2}{B^3}\nabla B\times\boldsymbol{B}\right) \\
& \nabla\times\boldsymbol{B} =
(\partial_y B_z - \partial_z B_y)\hat{i} +
(\partial_z B_x - \partial_x B_z)\hat{j} +
(\partial_x B_y - \partial_y B_x)\hat{k} \\
& \nabla B\times\boldsymbol{B} =
(B_z\partial_yB - B_y\partial_zB)\hat{i} +
(B_x\partial_zB - B_z\partial_xB)\hat{j} +
(B_y\partial_xB - B_x\partial_yB)\hat{k}
\end{aligned}\end{split}\]
where \(p=\gamma m v\) is particle momentum, \(c\) is the speed
of light, \(w=v/c\) is the normalized particle speed, and \(q\)
is particle charge. Using normalized quantities, we have
\[\begin{split}\begin{aligned}
\tilde{\boldsymbol{V}}_d & = \frac{1}{v_A}\frac{1}{3\tilde{q}e}\frac{\tilde{p}^2p_0^2c}{\sqrt{\tilde{p}^2p_0^2+m^2c^2}}\frac{1}{B_0L_0}
\left(\frac{1}{\tilde{B}^2}\tilde{\nabla}\times\tilde{\boldsymbol{B}} -
\frac{2}{\tilde{B}^3}\tilde{\nabla}\tilde{B}\times\tilde{\boldsymbol{B}}\right) \\
& = \frac{1}{\sqrt{d_1^2\tilde{p}^{-2}+d_2^2\tilde{p}^{-4}}}
\frac{1}{3\tilde{q}}\left(\frac{1}{\tilde{B}^2}\tilde{\nabla}\times\tilde{\boldsymbol{B}} -
\frac{2}{\tilde{B}^3}\tilde{\nabla}\tilde{B}\times\tilde{\boldsymbol{B}}\right)
\end{aligned}\end{split}\]
where \(\tilde{\boldsymbol{V}}_d=\boldsymbol{V}_d/v_A\),
\(\tilde{q}=q/e\), \(\tilde{\nabla}=L_0\nabla\),
\(\tilde{\boldsymbol{B}}=\boldsymbol{B}/B_0\),
\(\tilde{p}=p/p_0\), \(d_1=eB_0v_AL_0/(p_0c)\), and
\(d_2=emB_0v_AL_0/p_0^2\). Note that in the code, \(\tilde{p}\)
will be re-normalized. For example, \(\tilde{p}_0=1\) might
correspond to \(\tilde{p}_{n0}=0.1\) in simulations. The
re-normalized numerical momentum
\(\tilde{p}_n=\tilde{p}\tilde{p}_{n0}\). Thus,
\(\tilde{p} = \tilde{p}_n/\tilde{p}_{n0}\) in simulations, and we
need provide \(d_1\) and \(d_2\) based on the normalization.
Note
The velocity normalization \(v_A\) should be changed to \(v_0\) if \(v_0\neq v_A\)