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= \left ( 1 - \phi \left ( \lambda + 1 \right ) \right ) and D = \phi \left ( \lambda + 1 \right )

Now we can expand the derivatives and get

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

and the whole equation is then

\partial_t T = \displaystyle \sum_{\parallel,\perp}{A \left ( \partial_i B_{ij} \right ) \left ( C \partial_j T^{\lambda + 1} + D \partial_j T^{\lambda} T' \right ) + A B_{ij} \left ( C \partial_i \partial_j T^{\lambda + 1} + D \partial_i \partial_j T^{\lambda} T' \right )}

where A, B, C, D, and \lambda are functions of \perp,\parallel

We can also write this as

\partial_t T = \displaystyle \sum_{\parallel,\perp}{CE \partial_j T^{\lambda + 1} + DE \partial_j T^{\lambda} T' + CF \partial_i \partial_j T^{\lambda + 1} + DF \partial_i \partial_j T^{\lambda} T'}

where

E_j = A \left ( \partial_i B_{ij} \right )

F_{ij} = A B_{ij}

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 - \displaystyle \sum_{\parallel, \perp} C \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 + \displaystyle \sum_{\parallel, \perp} D \left [ \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}} \right ]

where

\alpha_0 = -\frac{2F_{ij}\delta_{ij} \Delta t}{\Delta x^2}

\alpha_{\pm j} = \pm \frac{E_j \Delta t}{2 \Delta x} + \frac{F_{jj}\Delta t}{\Delta x^2}

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