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}}

So we can rewrite the equations as

\frac{\partial T}{\partial t} = \nabla \cdot \left [ n \hat{b} \left ( \kappa_\parallel T^{\lambda_\parallel} - \frac{n \kappa_\perp}{B^2} T^{\lambda_\perp} \right ) \left ( \hat{b} \cdot \nabla T \right ) + \frac{n \kappa_\perp}{B^2} T^{\lambda_\perp} \nabla T \right ]

or

\frac{\partial T}{\partial t} = \nabla \cdot \left [ n \hat{b} \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( \hat{b} \cdot \nabla T^{\lambda_\parallel+1} \right ) - \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( \hat{b} \cdot \nabla T^{\lambda_\perp + 1} \right ) \right ) + \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp +1 \right )} \nabla T^{\lambda_\perp+1} \right ]

or in Einstein notation

\partial_t T = \partial_i \left [ n b_i \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( b_j \partial_j T^{\lambda_\parallel+1} \right ) - \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( b_j \partial_j T^{\lambda_\perp + 1} \right ) \right ) + \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp +1 \right )} \partial_i T^{\lambda_\perp+1} \right ]

or

\partial_t T = \partial_i \left [ n b_i \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( b_j \partial_j T^{\lambda_\parallel+1} \right ) - \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( b_j \partial_j T^{\lambda_\perp + 1} \right ) \right ) + \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp +1 \right )} \delta_{ij} \partial_j T^{\lambda_\perp+1} \right ]

or

\partial_t T = \partial_i n b_i \frac{\kappa_\parallel}{\lambda_\parallel+1} b_j \partial_j T^{\lambda_\parallel+1} + \partial_i n \left ( \delta_{ij} - b_i b_j \right ) \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp + 1\right )} \partial_j T^{\lambda_\perp+1}

or

\partial_t T = \frac{\kappa_\parallel}{\lambda_\parallel+1} \partial_i n b_i b_j \partial_j T^{\lambda_\parallel+1} + \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp + 1\right )} \partial_i n \left ( \delta_{ij} - b_i b_j \right ) \partial_j T^{\lambda_\perp +1}

Now both of these terms are of the form

A \partial_i B_{ij} \partial_j T^{\lambda+1}

Now to solve this implicitly, we need to replace T with T_* where

T_* = T + \phi ( T' - T ) = (1-\phi) T + \phi T'

Note for Backward Euler, \phi = 1 and for Crank Nicholson, \phi = 1/2

A \partial_i B_{ij} \partial_j T_*^{\lambda+1}

so we have

A \partial_i B_{ij} \partial_j \left ((1-\phi) T + \phi T' \right )^{\lambda+1}

Now to solve this using a linear system, we need to linearize terms involving T'

So we need to Taylor expand about T

\left ((1-\phi) T + \phi T' \right )^{\lambda+1} \approx T^{\lambda + 1} + \left ( \lambda + 1 \right ) T^{\lambda} \phi \left ( T' - T \right )

or

\left ((1-\phi) T + \phi T' \right )^{\lambda+1} \approx \left ( 1 - \phi \left ( \lambda + 1 \right ) \right )T^{\lambda + 1} + \phi \left ( \lambda + 1 \right ) T^{\lambda} T'

Now we could also have expanded T_*^{\lambda+1} and then plug in T_*=\left ( 1 - \phi \right ) T + \phi T'. In that case we also get

T_*^{\lambda + 1} \approx T^{\lambda + 1} + \left ( \lambda + 1 \right ) T^{\lambda} \left ( T_{*} - T \right ) = T^{\lambda + 1} + \left ( \lambda + 1 \right ) T^{\lambda} \phi \left ( T' - T \right ) = \left(1 - \phi \left ( \lambda + 1\right ) \right ) T^{\lambda + 1}+ \phi \left ( \lambda + 1 \right ) T^{\lambda} T'

So we have

A \partial_i B_{ij} \partial_j \left [ \left ( 1 - \phi \left ( \lambda + 1 \right ) \right )T^{\lambda + 1} + \phi \left ( \lambda + 1 \right ) T^{\lambda} T' \right ]

A \partial_i B_{ij} \left [ \left ( 1 - \phi \left ( \lambda + 1 \right ) \right ) \partial_j T^{\lambda + 1} + \phi \left ( \lambda + 1 \right ) \partial_j T^{\lambda} T' \right ]

A \partial_i B_{ij} \left [ C_\lambda \partial_j T^{\lambda + 1} + D_\lambda \partial_j T^{\lambda} T' \right ]

where C_\lambda= \left ( 1 - \phi \left ( \lambda + 1 \right ) \right ) and D_\lambda = \phi \left ( \lambda + 1 \right )

Now we can expand the derivatives and get

A \left ( \partial_i B_{ij} \right ) \left ( C_\lambda \partial_j T^{\lambda + 1} + D_\lambda \partial_j T^{\lambda} T' \right ) + A B_{ij} \left ( C_\lambda \partial_i \partial_j T^{\lambda + 1} + D_\lambda \partial_i \partial_j T^{\lambda} T' \right )

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

\phi' = \frac{\left ( \phi-\psi\lambda \right )}{\psi \left (\lambda+1 \right )}

B_j = \phi' C_j

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

D_{ij}=\phi'E_{ij}

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

We then need expressions for

\partial_t T = \frac{ T'_0 - T_0}{\Delta t}

\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

-T_0 - \phi' \left [ \alpha_0 T^{\lambda+1}_0 - \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda+1}_{\pm \hat{j}} - \sum_{\pm i, \pm j,|i| \ne |j|} \alpha_{\pm i, \pm j} T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}} \right ] = -T'_0 + \alpha_0 T^{\lambda}_0T'_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda}_{\pm \hat{j}} T'_{\pm \hat{j}} + \sum_{\pm i, \pm j,|i| \ne |j|} \alpha_{\pm i, \pm j} T^{\lambda}_{\pm \hat{i} \pm \hat{j}} T'_{\pm \hat{i} \pm \hat{j}}

where

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

\alpha_{\pm j} = \pm \frac{C_j \Delta t}{2 \Delta x} + \frac{E_{jj}\Delta t}{\Delta x^2}

\alpha_{\pm i, \pm j} = \pm \pm \frac{E_{ij}\Delta t}{4 \Delta x^2}

Now in the limit where we ignore spatial derivatives in density and magnetic field orientation, we can throw away the C terms - and if we also assume backward Euler, we can set

\psi=1

\phi' = \frac{-\lambda}{\lambda+1}

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

If we are using Crank-Nicholson, then

\psi=1/2

\phi' = \frac{1-\lambda}{\lambda+1}

E_{ij}=\frac{1}{2}\kappa_\parallel n b_i b_j

Now for the perpendicular term, we have

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

And performing a Taylor expansion gives

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

And in Einstein notation gives

\partial_t T = - \partial_i n b_i \frac{n \kappa_\perp}{B^2 \left ( \Lambda+1 \right ) } b_j \partial_j \left ( -\Lambda T^{\Lambda+1} + \left ( \Lambda + 1 \right ) T^{\Lambda}T_{*} \right ) + \partial_i n \frac{n \kappa_\perp}{B^2 \left ( \Lambda+1 \right ) } \partial_i \left ( - \Lambda T^{ \Lambda+1} + \left ( \Lambda + 1 \right ) T^{\Lambda} T_{*} \right )

The first term is very similar to the parallel case - and the second term can be made similar by replacing b_i b_j with \delta_{ij}

\phi' = \frac{\phi-\psi \Lambda}{\psi \left ( \Lambda + 1 \right )}

B_j = \phi' C_j

C_j = -\kappa_\perp \psi \partial_i n^2 B^{-2} \left ( b_i b_j - \delta_{ij} \right )

D_{ij}=\phi'E_{ij}

E_{ij} = -\kappa_\perp \psi n^2 B^{-2} \left ( b_i b_j - \delta_{ij} \right )