The anisotropic conduction equation looks like

\frac{\partial T}{\partial t} = \nabla \cdot \left [ n \hat{b} \left ( \chi_\parallel - \chi_\perp \right ) \left ( \hat{b} \cdot \nabla T \right ) + n \chi_\perp \nabla T \right ]

however, the conduction coefficients have a Temperature dependence.

\chi_\parallel = \kappa_\parallel T^{5/2}

\chi_\perp = \kappa_\perp \frac{n}{B^2 T^{1/2}}

Let's first just consider the \chi_\parallel term.

\frac{\partial T}{\partial t} = \nabla \cdot \left [ n \hat{b} \kappa_\parallel T^\lambda \left ( \hat{b} \cdot \nabla T \right )\right ]

We need a way to write this implicitly but we need it to also be linear in T_{*}

\frac{\partial T}{\partial t} = \nabla \cdot \left [ n \hat{b} \kappa_\parallel T_{*}^\lambda \left ( \hat{b} \cdot \nabla T_{*} \right )\right ]

We could Taylor expand T^\lambda_{*}=T^\lambda + \lambda T^{\lambda-1} \left ( T_{*}-T \right ) = \left ( 1-\lambda \right ) T^\lambda + \lambda T^{\lambda-1} T_{*}

but then we would still have a non-linear term like \lambda T^{\lambda-1} T_{*} \nabla T_{*}

We could write this as \lambda T^{\lambda-1} 1/2 \nabla T_{*}^2 and Taylor expand again to get \lambda T^{\lambda-1} 1/2 \nabla \left ( - T + 2 TT_{*} \right )

but we've now done a Taylor expansion on a Taylor expansion…

Alternatively, we can rewrite the diffusion equation

\frac{\partial T}{\partial t} = \nabla \cdot \left [ n \hat{b} \frac{\kappa_\parallel}{\lambda+1} \left ( \hat{b} \cdot \nabla T_*^{\lambda+1} \right )\right ]

and then perform a single Taylor expansion

\frac{\partial T}{\partial t} = \nabla \cdot \left [ n \hat{b} \frac{\kappa_\parallel}{\lambda+1} \left ( \hat{b} \cdot \nabla \left (-\lambda T^{\lambda+1} + \left ( \lambda + 1 \right ) T^{\lambda}T_{*} \right ) \right ) \right ]

Switching to Einstein notation, we have

\partial_t T = \frac{\kappa_\parallel}{\lambda+1} \partial_i n b_i b_j \partial j \left ( -\lambda T^{\lambda+1} + \left ( \lambda + 1 \right ) T^{\lambda}T_{*} \right )

\partial_t T = \frac{\kappa_\parallel}{\lambda+1} \partial_i n b_i b_j \partial j \left ( -\lambda T^{\lambda+1} + \left ( \lambda + 1 \right ) T^{\lambda}T_{*} \right )

Let's also take a moment to write T_* = \phi T + \psi T' where T' is the new temperature, and \phi + \psi = 1. Backward Euler would have \phi=0 and \psi=1 where Crank-Nicholson would have \phi=\psi=1/2

\partial_t T = \frac{\kappa_\parallel}{\lambda+1} \partial_i n b_i b_j \partial j \left ( -\lambda T^{\lambda+1} + \left ( \lambda + 1 \right ) T^{\lambda}\left ( \phi T + \psi T' \right ) \right )

\partial_t T = \frac{\kappa_\parallel}{\lambda+1} \partial_i n b_i b_j \partial j \left ( \left ( \phi - \psi \lambda \right ) T^{\lambda+1} + \psi \left ( \lambda + 1 \right ) T^{\lambda} T' \right )

\partial_t T = \frac{\kappa_\parallel}{\lambda+1} \left [ \left ( \partial_i n b_i b_j \right ) \partial j \left ( \left ( \phi - \psi \lambda \right ) T^{\lambda+1} + \psi \left ( \lambda + 1 \right ) T^{\lambda} T' \right ) + n b_i b_j \left ( \left ( \phi - \psi \lambda \right ) \partial_i \partial_j T^{\lambda+1} + \psi \left ( \lambda + 1 \right ) \partial_i \partial_j T^{\lambda} T' \right ) \right ]

Now if we write the equation as

\partial_t T = B_j \partial_j T^{\lambda+1} + C_j \partial_j T^\lambda T' + D_{ij} \partial_i\partial_j T^{\lambda+1} + E_{ij} \partial_i \partial_j T^\lambda T'

we get expressions for

A = \frac{1}{\Delta t}

B_j = \frac{\kappa_\parallel \left ( \phi-\psi\lambda \right )}{\lambda+1} \partial_i n b_i b_j

C_j = \kappa_\parallel \psi \partial_i n b_i b_j

D_{ij}=\frac{\kappa_\parallel \left ( \phi-\psi\lambda \right )}{\lambda+1} n b_i b_j

E_{ij} = \kappa_\parallel \psi n b_i b_j

We then need expressions for

\partial_t T = A \left ( T'_0 - T_0 \right )

\partial_j T^{\lambda+1} = \frac{T^{\lambda+1}_{\hat{j}} - T^{\lambda+1}_{-\hat{j}}}{2 \Delta x}

\partial_j T^{\lambda}T' = \frac{T^{\lambda}_{\hat{j}}T'_{\hat{j}} - T^{\lambda}_{-\hat{j}}T'_{-\hat{j}}}{2 \Delta x}

\partial_i\partial_j T^{\lambda+1} = \frac{T^{\lambda+1}_{\hat{i}+\hat{j}} - T^{\lambda+1}_{\hat{i}-\hat{j}}-T^{\lambda+1}_{-\hat{i}+\hat{j}} + T^{\lambda+1}_{-\hat{i}-\hat{j}}}{4 \Delta x^2}\left(1-\delta_{ij}\right ) + \frac{T^{\lambda+1}_{\hat{i}} - 2T^{\lambda+1}_{0} + T^{\lambda+1}_{-\hat{i}}}{\Delta x^2}\delta_{ij}

\partial_i\partial_j T^{\lambda}T' = \frac{T^{\lambda}_{\hat{i}+\hat{j}}T'_{\hat{i}+\hat{j}} - T^{\lambda}_{\hat{i}-\hat{j}}T'_{\hat{i}-\hat{j}}-T^{\lambda}_{-\hat{i}+\hat{j}}T'_{-\hat{i}+\hat{j}} + T^{\lambda}_{-\hat{i}-\hat{j}}T'_{-\hat{i}-\hat{j}}}{4 \Delta x^2}\left(1-\delta_{ij}\right ) + \frac{T^{\lambda}_{\hat{i}}T'_{\hat{i}} - 2T^{\lambda}_{0}T'_{0} + T^{\lambda}_{-\hat{i}}T'_{-\hat{i}}}{\Delta x^2}\delta_{ij}

We can also write the equation as

\beta = \alpha_0 T'_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T'_{\pm \hat{j}} + \sum_{\pm i, \pm j,|i| \ne |j|} \alpha_{\pm i, \pm j} T'_{\pm \hat{i} \pm \hat{j}}

and then work out the coefficients for the matrix equation

\alpha_0 = -\frac{1}{\Delta t} -\frac{2E_{jj}}{\Delta x^2}T_0^\lambda

\alpha_{\pm j} = \pm C_jT^\lambda_{\pm \hat{j}} + \frac{E_{jj}}{\Delta x^2}T_{\pm \hat{j}}^\lambda

\alpha_{\pm i \pm j} = \pm \pm \frac{E_{ij}}{4 \Delta x^2}T^\lambda_{\pm \hat{i} \pm \hat{j}}