# Non-Exemplar Online Class-incremental Continual Learning via
Dual-prototype Self-augment and Refinement
——Appendix——

Fushuo Huo1, Wenchao Xu1, Jingcai Guo1, 2, Haozhao Wang3111Corresponding
author, Yunfeng Fan1

## Overview

The appendix presents more experimental settings, results, and analyses as
follows:

Appendix A More implementation details

Appendix B More results on different dataset partitions.

Appendix C Results and analysis on different base session training strategies.

Appendix D Hyper parameter analysis

Appendix E Computation overhead analysis

Appendix A: More Implementation Details

Dataset Overview. We conduct experiments on three widely used datasets,
including CORE-50 (Lomonaco and Maltoni 2017), CIFAR 100 (Krizhevsky and
Hinton 2009), and Mini-ImageNet (Vinyals et al. 2016). Here we give brief
introductions. CORE-50 (Lomonaco and Maltoni 2017) is a benchmark designed for
class incremental learning with 50 classes. Each class has around 2,398
training images and 900 testing images, with the size of $3\times 128\times
128$. CIFAR 100 (Krizhevsky and Hinton 2009) contains 60000 images of
$32\times 32$ size from 100 classes, and each class includes 500 training
images and 100 test images. Mini-ImageNet (Vinyals et al. 2016) contains 100
classes and is divided into 10 sub-datasets for 10 disjoint tasks, and each
task contains 10 classes. Each task comprises 5,000 training images and 1,000
testing images, all with the size of $3\times 84\times 84$.

Training Details. For OCL methods, we employ the same dataset partitions and
training protocols of NO-CL, i.e., pre-training on the base classes and then
online class-incremental learning _with example buffers_. Other
hyperparameters are adopted as default. The example buffers are restored and
retrieved during the whole training procedure with the default updating
pipeline. MIR (Aljundi et al. 2019), GD (Prabhu, Torr, and Dokania 2020), ASER
(Shim et al. 2021), SCR (Mai et al. 2021), and DVC (Gu et al. 2022) are based
on the OCL codebase 222https://github.com/RaptorMai/online-continual-learning.
Other methods are implemented with the public released codes. For FS-CL
methods, FACT (Zhou et al. 2022a) and ALICE (Peng et al. 2022), we also adopt
the same protocols as NO-CL. The prototypes of novel classes are computed by
all data samples rather than few-shot samples. During the inference phase,
FACT and ALICE directly infer incremental data samples via computed prototypes
without finetuning the network. All methods are implemented with the public
released codes. For NE-CL methods (Zhu et al. 2021b, 2022), as (Zhu et al.
2022) does not provide training scripts, we adopt the three-party codes (Zhou
et al. 2022b)333https://github.com/G-U-N/PyCIL on CIFAR100 dataset. For (Zhu
et al. 2021b), we employ the same dataset partitions and training and testing
protocols of NO-CL. Note that all methods are employed the _same_ reduced
ResNet-18 as the feature extractor for fair comparisons. All experiments are
conducted with NVIDIA RTX3090 GPU on CUDA 11.4 using PyTorch framework.

Base Session Training Details. Here, we give the details of our base session
training strategy. We employ base training functions on the outputs of the
feature extractor and projection module to obtain vanilla and high-dimensional
prototypes for sequentially online sessions:
$L^{base}=L^{base}_{vp}+L^{base}_{hp}$.
$L^{base}_{vp}=Loss(Proj_{vp}(\theta_{1}(x)),y)$ and
$L^{base}_{vp}=Loss(Proj_{hp}(\theta{{}_{2}}(\theta_{1}(x))),y)$. $x$, $y$,
$\theta_{1}$, and $\theta_{2}$ denote input samples, labels, feature
extractor, and projection module. $Proj_{vp/hp}$ are linear layers to align
vanilla- and high-dimensional prototypes for loss calculations. For cross-
entropy (CE) loss functions, $Proj_{vp/hp}$ are one-layer MLP with the output
dimension of base class. For supervised contrastive (SC) loss (Khosla et al.
2020), we follow SCR (Mai et al. 2021) and adopt the _same_ hyper-parameters
of SCR. $Proj_{vp/hp}$ are two-layer MLP with the dimension of 160 and 128 and
the temperature is set to 0.1.

Appendix B: Results on Different Dataset Partitions

Due to space constraints in the main paper, in this subsection, we report the
additional results of different dataset partitions. Concretely, as shown in
Tables 1 and 2, we conduct experiments on the configuration of $40\%+6\%\times
10$ and $80\%+2\%\times 10$, where $40\%$ and $60\%$ classes are selected as
the base classes, and the rest classes are continually fed to the network in
10 sessions. Moreover, more incremental sessions (i.e., 20 sessions) with the
configuration of $60\%+2\%\times 20$ are also compared in Table 3. The five
representative state-of-the-art methods are compared with the same training
and inference protocols as Non-exemplar Online Class-incremental continual
Learning (NO-CL), including Online Class-incremental continual Learning (OCL)
methods (i.e., SCR (Mai et al. 2021), OCM (Guo, Liu, and Zhao 2022), and
DVC(Gu et al. 2022)), Non-Exemplar Class-incremental continual Learning (NE-
CL) method (i.e., PASS (Zhu et al. 2021b)), and Few-Shot Class-incremental
Learning (FS-CL) method (i.e., ALICE (Peng et al. 2022)). As we can see from
Tables 1 and 2, fewer base classes result in poor performance, both in base
and novel classes, because the network lacks enough pre-trained information to
generalize to novel classes. Meanwhile, as our method depends on the inner-
prototype computed by the pre-trained backbone, the well-trained backbone
benefits us much for the prototype refinement. Notably, even with $40\%$ base
classes, our method also achieves the best performance in Acc and HM metrics,
which validates the robustness of prototype refinement strategies. Moreover,
as for more incremental sessions in Table 3, the performance only drops
slightly. Overall, experiments on different dataset partitions validate the
effectiveness and robustness of our method.

Methods | CORE-50 | CIFAR100 | Mini-ImageNet
---|---|---|---
Metrics |  Acc(base/novel)$|$HM |  Acc(base/novel)$|$HM |  Acc(base/novel)$|$HM
ALICE | 35.0(49.3/25.5)$|$33.6 | 37.1(59.2/22.3)$|$32.4 | 36.5(58.6/21.8)$|$31.8
PASS | 25.2(62.8/0.2)$|$0.4 | 26.4(65.1/0.6)$|$1.2 | 25.8(63.4/0.8)$|$1.6
MS | 1000 | 2000 | 1000 | 2000 | 1000 | 2000
SCR | 37.3(36.9/37.5)$|$37.2 | 38.7(38.9/38.6)$|$38.7 | 32.7(34.2/31.7)$|$32.9 | 34.9(36.9/33.6)$|$35.2 | 29.8(29.1/30.3)$|$29.7 | 36.3(39.2/34.4)$|$36.6
SCRft | 34.2(45.2/26.8)$|$33.6 | 37.7(49.8/29.6)$|$37.1 | 30.9(45.2/21.3)$|$21.3 | 36.4(50.8/26.8)$|$35.1 | 32.9(43.2/26.1)$|$32.5 | 36.8(50.8/27.5)$|$35.7
OCM | 37.6(37.5/37.6)$|$37.5 | 39.5(40.5/38.9)$|$39.7 | 31.6(32.9/30.8)$|$31.8 | 37.1(36.1/37.8)$|$36.9 | 31.3(30.8/31.6)$|$31.2 | 31.7(32.9/30.9)$|$31.9
OCMft | 35.8(43.1/30.9)$|$36.0 | 36.7(46.1/30.6)$|$36.7 | 32.4(48.2/21.8)$|$30.1 | 35.9(50.6/26.1)$|$36.8 | 30.2(38.1/24.9)$|$31.4 | 37.7(51.2/28.7)$|$36.8
DVC | 37.5(36.9/37.9)$|$37.3 | 38.7(39.8/38.0)$|$38.9 | 29.9(30.0/29.8)$|$29.9 | 37.2(34.6/39.0)$|$36.6 | 29.7(29.9/29.6)$|$29.7 | 33.3(35.6/31.8)$|$33.6
DVCft | 36.7(45.8/30.6)$|$36.7 | 37.9(45.1/33.1)$|$38.2 | 32.3(43.6/24.7)$|$31.5 | 36.5(46.5/29.8)$|$36.3 | 31.4(39.4/26.1)$|$31.3 | 32.6(38.5/28.6)$|$32.8
Ours | 45.5(44.2/46.4)$|$45.3 | 38.6(43.8/35.2)$|$39.0 | 40.4(55.4/30.4)$|$39.3
Table 1: The quantitative analysis of dataset partition $\textbf{40\%}+\textbf{6\%}\times\textbf{10}$. Class-wise accuracy (Acc) by end of the training in terms of all classes, base classes, and novel classes and Harmonic accuracy (HM) are illustrated. MS and ft mean the example memory size and finetuning versions. The best results are marked in bold. Methods | CORE-50 | CIFAR100 | Mini-ImageNet
---|---|---|---
Metrics |  Acc(base/novel)$|$HM |  Acc(base/novel)$|$HM |  Acc(base/novel)$|$HM
ALICE | 44.7(48.2/30.6)$|$37.4 | 50.5(57.2/23.8)$|$33.6 | 50.3(56.4/26.1)$|$35.7
PASS | 48.4(60.3/0.6)$|$1.2 | 50.4(62.8/0.8)$|$1.6 | 49.8(62.0/0.9)$|$1.8
MS | 1000 | 2000 | 1000 | 2000 | 1000 | 2000
SCR | 39.9(39.2/42.9)$|$40.9 | 42.9(43.5/40.8)$|$42.1 | 40.7(42.2/34.7)$|$38.1 | 46.3(47.2/42.8)$|$44.9 | 40.4(40.6/39.7)$|$40.1 | 44.5(43.8/47.3)$|$45.5
SCRft | 45.4(48.9/31.2))$|$38.1 | 49.7(53.7/33.7)$|$54.6 | 47.4(50.6/34.8)$|$41.2 | 49.2(53.2/33.4)$|$41.0 | 40.2(45.2/32.8)$|$38.0 | 47.0(49.2/38.5)$|$43.5
OCM | 43.4(43.8/42.2)$|$42.9 | 43.2(43.6/41.8)$|$42.7 | 39.9(39.8/40.6)$|$40.2 | 43.6(43.9/42.3)$|$43.1 | 39.5(38.9/41.8)$|$40.3 | 43.5(42.9/46.0)$|$44.3
OCMft | 46.8(49.8/34.8)$|$40.9 | 47.8(49.8/39.7)$|$44.2 | 46.1(48.8/35.1)$|$40.8 | 48.7(51.8/36.3)$|$42.7 | 44.1(46.1/36.3)$|$40.6 | 46.3(48.6/37.2)$|$42.1
DVC | 43.0(42.8/43.8)$|$43.3 | 45.2(45.8/42.8)$|$44.2 | 39.4(38.6/42.5)$|$40.5 | 40.1(39.7/41.9)$|$40.8 | 41.6(42.6/37.4)$|$39.8 | 45.4(45.6/44.8)$|$45.2
DVCft | 48.0(50.3/39.0)$|$43.9 | 50.3(53.7/36.9)$|$43.7 | 41.1(41.6/38.9)$|$40.2 | 42.4(43.9/36.4)$|$39.8 | 45.3(48.6/31.9)$|$38.5 | 47.8(50.8/36.2)$|$42.2
Ours | 55.2(55.6/53.7)$|$54.6 | 54.3(56.2/46.8)$|$51.1 | 52.2(52.6/50.8)$|$51.7
Table 2: The quantitative analysis of dataset partition $\textbf{80\%}+\textbf{2\%}\times\textbf{10}$. Class-wise accuracy (Acc) by end of the training in terms of all classes, base classes, and novel classes and Harmonic accuracy (HM) are illustrated. MS and ft mean the example memory size and finetuning versions. The best results are marked in bold. Methods | CORE-50 | CIFAR100 | Mini-ImageNet
---|---|---|---
Metrics |  Acc(base/novel)$|$HM |  Acc(base/novel)$|$HM |  Acc(base/novel)$|$HM
ALICE | 39.5(46.2/29.5)$|$36.0 | 42.5(53.5/25.9)$|$34.9 | 41.1(51.4/25.7)$|$34.3
PASS | 35.2(58.3/0.8)$|$1.6 | 37.9(62.6/1.0)$|$2.0 | 37.2(61.2/1.1)$|$2.2
MS | 1000 | 2000 | 1000 | 2000 | 1000 | 2000
SCR | 38.6(37.2/40.6)$|$38.8 | 38.6(39.2/37.6)$|$38.4 | 36.9(38.8/34.1)$|$36.3 | 40.7(42.8/37.6)$|$40.0 | 34.4(34.6/34.2)$|$34.4 | 35.1(38.6/29.8)$|$33.6
SCRft | 36.6(42.4/27.9)$|$33.7 | 41.3(49.8/28.6)$|$36.3 | 38.6(49.2/22.8)$|$31.2 | 41.1(52.8/23.7)$|$32.7 | 37.9(43.2/30.0)$|$35.4 | 39.5(43.7/33.2)$|$37.7
OCM | 38.3(38.8/37.6)$|$38.2 | 41.1(42.1/39.6)$|$40.8 | 36.3(36.5/35.9)$|$36.2 | 41.4(42.7/39.6)$|$41.1 | 35.9(35.2/37.1)$|$36.1 | 39.5(43.7/33.2)$|$37.7
OCMft | 37.7(43.1/29.7)$|$35.2 | 43.3(47.9/36.4)$|$41.4 | 39.7(45.7/30.7)$|$36.7 | 40.9(45.9/33.4)$|$38.7 | 36.8(39.8/32.2)$|$35.4 | 41.5(43.5/38.4)$|$40.7
DVC | 38.1(37.8/38.6)$|$38.2 | 40.5(41.8/38.6)$|$40.1 | 37.5(37.2/38.1)$|$37.7 | 40.2(41.6/38.2)$|$39.8 | 34.4(32.9/36.7)$|$34.7 | 37.1(35.4/39.6)$|$37.4
DVCft | 39.8(45.9/30.7)$|$36.8 | 41.8(46.3/35.0)$|$39.9 | 37.8(42.1/31.4)$|$35.9 | 40.6(44.8/34.3)$|$38.9 | 35.1(38.9/29.4)$|$33.4 | 37.7(38.7/36.1)$|$37.4
Ours | 49.5(49.1/50.2)49.6 | 47.6(51.6/41.7)$|$46.1 | 49.1(53.8/42.1)$|$47.2
Table 3: The quantitative analysis of dataset partition
$\textbf{60\%}+\textbf{2\%}\times\textbf{20}$. Class-wise accuracy (Acc) by
end of the training in terms of all classes, base classes, and novel classes
and Harmonic accuracy (HM) are illustrated. MS and ft mean the example memory
size and finetuning versions. The best results are marked in bold.

Appendix C: Results and Analysis on Different the Base Session Training
Strategies

The stability and plasticity dilemma is a thorny problem in the area of
continual learning. To deal with this dilemma, previous NE-CL (Zhu et al.
2021b, a) and FS-CL (Peng et al. 2022; Kalla and Biswas 2022) methods employ
self-supervised learning (Jing and Tian 2021) and class and data augmentation
to learn task-agnostic and transferable representations. For the problem of
NO-CL, the base session training strategies also matter for the stability and
plasticity dilemma. For fair comparisons, similar to (Mai et al. 2021; Gu et
al. 2022; Guo, Liu, and Zhao 2022), we also employ supervised contrastive (SC)
learning. Here, we provide two training strategies. Concretely, we add the
extra self-supervised learning loss (Lee, Hwang, and Shin 2020) (+SSL) like
(Zhu et al. 2021b; Kalla and Biswas 2022) and use the data augmentation
strategy (+DA) proposed by (Zhu et al. 2021a; Peng et al. 2022). The results
in Table 4 show that the elaborately designed pre-training strategies improve
the accuracy both in the base and novel classes. Therefore, developing a more
robust pre-training strategies is a promising way for the proposed NO-CL
problem.

Ablations | CIFAR100 | Mini-ImageNet
---|---|---
Metrics |  Acc(base/novel) |  Acc(base/novel)
Ours(+CE) | 45.8(50.0/39.6) | 47.7(52.6/40.3)
Ours(+SC) | 48.6(52.4/42.9) | 50.7(56.1/42.6)
+SSL | 50.4(53.8/45.2) | 52.2(57.3/44.6)
+DA | 51.2(55.7/44.6) | 53.2(58.2/45.8)
Table 4: The results of Different training strategies on CIFAR100 and Mini-ImageNet. SSL and DA mean self-supervised learning and data augmentation, respectively. Figure 1: Quantitative results of varying $\lambda$. Figure 2: Quantitative results of varying $T$. Figure 3: Quantitative results of varying $K_{0}$, $K_{1}$, and $K_{2}$ | CIFAR100 | Mini-ImageNet
---|---|---
Ablations |  Acc(base/novel)/HM |  Acc(bse/novel)/HM
w/ 256 | 44.8(49.5/37.9)/42.9 | 47.4(53.6/38.1)/44.5
w/ 1024 | 47.4(51.8/40.8)/45.7 | 49.7(55.4/41.1)/47.2
w/ 2048 | 48.6(52.4/42.9)/47.2 | 50.7(56.1/42.6)/48.4
w/ 3074 | 48.4(52.3/42.7)/47.0 | 50.8(56.3/42.6)/48.4
Table 5: Quantitative results of varying dimension of hyperdimensional
embedding.

(a) CORE-50

(b) CIFAR100

(c) Mini-ImageNet

Figure 4: Training losses of bi-level optimization procedure | CIFAR100 | Mini-ImageNet
---|---|---
Metrics | ALICE | SSRE | SCR | DVC | OCM | Ours | ALICE | SSRE | SCR | DVC | OCM | Ours
Time(s) | 512 | 291 | 165 | 126 | 561 | 35 | 793 | 457 | 254 | 194 | 831 | 61
Memory(GB) | 1.9 | 3.2 | 1.8 | 1.6 | 12.8 | 1.4 | 4.4 | 6.8 | 4.1 | 2.8 | 21.4 | 1.9
Table 6: Quantitative comparisons of computation overhead in terms of online
training time and memory footprint.

Appendix D: Hyper Parameter Analysis

Here, we analyse hyperparameters including feature transform coefficient
$\lambda$, online iteration $T$, and the number of sampled prototypes $K$. We
provide quantitative results on the Mini-ImageNet dataset in Figures 1, 2, and
3. Also, the experiments of the varying dimension of hyperdimensional
embedding are conducted in Table 5.

In Figure 1, we can see that $\lambda>1$ leads to degraded performance as the
feature distribution is more concentrated close to 0. Meanwhile, decreasing
too much $\lambda$ makes the distribution scattered and less aligned to the
calibrated Gaussian distribution. Therefore, we set $\lambda$ as 0.5 in our
experiment.

In Figure 2, we vary online iteration $T$. Less $T$ iterations harm the
network to accommodate online novel classes by refining hyperdimensional
prototypes and aligning the projection module. Meanwhile, more online
iterations only lead to slight degradation, which validates the plasticity of
our method. Therefore, we set $T$ to 20 to achieve the stability-plasticity
trade-off.

In Figure 3, we vary the number of sampled prototypes $K$. Concretely, we vary
the number of sampled prototypes of base classes $K_{base}$, novel classes
$K_{novel}$, and all classes $K$. We can see that the imbalance sampling of
classes leads to performance degradation, which is similar to the class
imbalance problem (Hou et al. 2019; Wu et al. 2019). Also, increasing $K$ does
not bring many gains while inducing computation overheads. Therefore, we set
$K=20$.

From Table 5, we can learn that increasing the dimension of hyperdimensional
embedding benefits the proposed method while too large dimension brings little
gain. Therefore, we set the dimension of hyperdimensional embedding as 2048 in
our experiments.

Appendix E: Computation Overhead Analysis

For computation overheads during online learning, which is usually considered
in OCL scenarios (Fini et al. 2020), we provide more quantitative comparisons
to OCL, NE-CL, and FS-CL methods in Table 6. The batchsize of example-based
methods is set as 10. As we only align prototypes by finetuning the projection
module, which is much more efficient compared with training the whole network.
Therefore, our method has clear advantages on computation overheads for online
continual learning. Meanwhile, the bi-level optimization quickly converges as
shown in Figure 4.

## References

  * Aljundi et al. (2019) Aljundi, R.; Belilovsky, E.; Tuytelaars, T.; Charlin, L.; Caccia, M.; Lin, M.; and Page-Caccia, L. 2019.  Online Continual Learning with Maximal Interfered Retrieval.  In _NeurIPS_.
  * Fini et al. (2020) Fini, E.; Lathuilière, S.; Sangineto, E.; Nabi, M.; and Ricci, E. 2020.  Online Continual Learning under Extreme Memory Constraints.  In _ECCV_ , 720–735.
  * Gu et al. (2022) Gu, Y.; Yang, X.; Wei, K.; and Deng, C. 2022.  Not Just Selection, but Exploration: Online Class-Incremental Continual Learning via Dual View Consistency.  In _CVPR_ , 7442–7451.
  * Guo, Liu, and Zhao (2022) Guo, Y.; Liu, B.; and Zhao, D. 2022.  Online Continual Learning through Mutual Information Maximization.  In _ICML_.
  * Hou et al. (2019) Hou, S.; Pan, X.; Loy, C. C.; Wang, Z.; and Lin, D. 2019.  Learning a Unified Classifier Incrementally via Rebalancing.  In _CVPR_.
  * Jing and Tian (2021) Jing, L.; and Tian, Y. 2021.  Self-Supervised Visual Feature Learning With Deep Neural Networks: A Survey.  _IEEE TPAMI_ , 43(11): 4037–4058.
  * Kalla and Biswas (2022) Kalla, J.; and Biswas, S. 2022.  S3C: Self-Supervised Stochastic Classifiers for Few-Shot Class-Incremental Learning.  In _ECCV_ , 432–448. Cham.
  * Khosla et al. (2020) Khosla, P.; Teterwak, P.; Wang, C.; Sarna, A.; Tian, Y.; Isola, P.; Maschinot, A.; Liu, C.; and Krishnan, D. 2020.  Supervised Contrastive Learning.  In _NeurIPS_ , 18661–18673.
  * Krizhevsky and Hinton (2009) Krizhevsky, A.; and Hinton, G. 2009.  Learning multiple layers of features from tiny images.  In _Technical Report_.
  * Lee, Hwang, and Shin (2020) Lee, H.; Hwang, S. J.; and Shin, J. 2020.  Self-supervised Label Augmentation via Input Transformations.  In _ICML_ , 5714–5724.
  * Lomonaco and Maltoni (2017) Lomonaco, V.; and Maltoni, D. 2017.  CORe50: a New Dataset and Benchmark for Continuous Object Recognition.  In _CoRL_ , 17–26.
  * Mai et al. (2021) Mai, Z.; Li, R.; Kim, H.; and Sanner, S. 2021.  Supervised Contrastive Replay: Revisiting the Nearest Class Mean Classifier in Online Class-Incremental Continual Learning.  In _CVPR Workshops_ , 3589–3599.
  * Peng et al. (2022) Peng, C.; Zhao, K.; Wang, T.; Li, M.; and Lovell, B. C. 2022.  Few-Shot Class-Incremental Learning from an Open-Set Perspective.  In _ECCV_ , 382–397.
  * Prabhu, Torr, and Dokania (2020) Prabhu, A.; Torr, P. H. S.; and Dokania, P. K. 2020.  GDumb: A Simple Approach that Questions Our Progress in Continual Learning.  In _ECCV_ , 524–540.
  * Shim et al. (2021) Shim, D.; Mai, Z.; Jeong, J.; Sanner, S.; Kim, H.; and Jang, J. 2021.  Online Class-Incremental Continual Learning with Adversarial Shapley Value.  _AAAI_.
  * Vinyals et al. (2016) Vinyals, O.; Blundell, C.; Lillicrap, T.; kavukcuoglu, k.; and Wierstra, D. 2016.  Matching Networks for One Shot Learning.  In Lee, D.; Sugiyama, M.; Luxburg, U.; Guyon, I.; and Garnett, R., eds., _NeurIPS_.
  * Wu et al. (2019) Wu, Y.; Chen, Y.; Wang, L.; Ye, Y.; Liu, Z.; Guo, Y.; and Fu, Y. 2019.  Large Scale Incremental Learning.  In _CVPR_.
  * Zhou et al. (2022a) Zhou, D.-W.; Wang, F.-Y.; Ye, H.-J.; Ma, L.; Pu, S.; and Zhan, D.-C. 2022a.  Forward Compatible Few-Shot Class-Incremental Learning.  In _CVPR_ , 9046–9056.
  * Zhou et al. (2022b) Zhou, D.-W.; Wang, F.-Y.; Ye, H.-J.; and Zhan, D.-C. 2022b.  PyCIL: A Python Toolbox for Class-Incremental Learning.  _SCIENCE CHINA Information Sciences_.
  * Zhu et al. (2021a) Zhu, F.; Cheng, Z.; Zhang, X.-y.; and Liu, C.-l. 2021a.  Class-Incremental Learning via Dual Augmentation.  In _NeurIPS_ , 14306–14318.
  * Zhu et al. (2021b) Zhu, F.; Zhang, X.-Y.; Wang, C.; Yin, F.; and Liu, C.-L. 2021b.  Prototype Augmentation and Self-Supervision for Incremental Learning.  In _CVPR_ , 5871–5880.
  * Zhu et al. (2022) Zhu, K.; Zhai, W.; Cao, Y.; Luo, J.; and Zha, Z.-J. 2022.  Self-Sustaining Representation Expansion for Non-Exemplar Class-Incremental Learning.  In _CVPR_ , 9296–9305.


# A reduced kinetic method for investigating non-local ion heat transport in
ideal multi-species plasmas

N T Mitchell1, D A Chapman2, C J McDevitt3, M P Read2, G Kagan1 1 The Blackett
Laboratory, Imperial College, London SW7 2AZ, UK  2 First Light Fusion Ltd.,
Unit 9/10 Oxford Pioneer Park, Mead Road, Yarnton, Kidlington OX5 1QU, UK  3
Nuclear Engineering Program, Department of Materials Science and Engineering,
University of Florida, Gainesville, Florida 32611, USA
<EMAIL_ADDRESS>
###### Abstract

A reduced kinetic method (RKM) with a first-principle collision operator is
introduced in a 1D2V planar geometry and implemented in a computationally
inexpensive code to investigate non-local ion heat transport in multi-species
plasmas. The RKM successfully reproduces local results for multi-species ion
systems and the important features expected to arise due to non-local effects
on the heat flux are captured. In addition to this, novel features associated
with multi-species, as opposed to single species, case are found. Effects of
non-locality on the heat flux are investigated in mass and charge symmetric
and asymmetric ion mixtures with temperature, pressure, and concentration
gradients. In particular, the enthalpy flux associated with diffusion is found
to be insensitive to sharp pressure and concentration gradients, increasing
its significance in comparison to the conductive heat flux driven by
temperature gradients in non-local scenarios. The RKM code can be used for
investigating other kinetic and non-local effects in a broader plasma physics
context. Due to its relatively low computational cost it can also serve as a
practical non-local ion heat flux closure in hydrodynamic simulations or as a
training tool for machine learning surrogates.

Keywords: non-local transport, reduced kinetic method, multi-species

##  1 Introduction

The collisional transport of energy is a significant issue for the design and
optimization of a wide range of schemes for realizing nuclear fusion in both
magnetized and inertial systems [1]. In the most commonly encountered context
of ICF, where a fuel-containing pellet is irradiated directly by lasers [2] or
indirectly by a laser-produced thermal X-ray bath [3], the influence of heat
flow is prominent in the transport of energy absorbed from the driver to the
ablation front [4] and the hot spot. Here, the dominant carriers of the
absorbed energy are delocalized (or ‘free’) electrons, which transfer their
energy through collisions with ions and other free electrons in the irradiated
material. If the temperature gradient is sufficiently shallow, this energy is
transported deeper into the unheated ablator through conduction, i.e., the
collisional ‘local’ diffusion of heat; the heat flux then obeys Fourier’s law
[5]. However, with high-intensity irradiation the mean free path (MFP) of the
high-energy particles that contribute most to the heat flow routinely exceed
the scale length of the gradient, resulting in transport with quasi-ballistic
features [6]. While ‘non-local’ behaviour may be accurately captured using
full kinetic simulations [7], these methods are impractical for integrated
modeling of the implosion process. Alternatively, and more pragmatically,
numerous simplified models have been constructed which may readily be
implemented as sub-cycled operators in radiation-hydrodynamics codes [8, 9,
10, 11, 12]. At later times, when the fusion fuel is in a state of pre-
ignition, heat flow also acts as a major factor in the power balance that
determines the pathway and probability of success of fuel burn and the onset
of ignition.

The foregoing arguments are usually framed from the perspective of the
electrons, as it is commonly assumed that the ionic heat flux can be neglected
in comparison to the electronic contribution. The assumption relies on
estimating the ratio of the ion to electron heat conductivities by the square
root of the electron-ion mass ratio. However, as recently discussed by Chapman
et al. [13], for fully ionized fusion fuel one finds parity of the
conductivities under modest temperature separations (around $T_{i}/T_{e}\sim
4$ for deuterium), which can be vastly exceeded in burning fusion plasmas.
Even more significant separations may be realized during the ignition and burn
phase as the hallmark of this process is ion temperature runaway until
hydrodynamic decompression occurs and the fuel explodes and rapidly cools.
Non-local effects are also expected for ion heat flux since $\alpha$-particle
deposition will distort the hot tail of the background distribution function
(DF), resulting in an over-abundance of energy-transporting ions relative to a
Maxwell-Boltzmann distribution [14].

Hydrodynamic quantities such as density and temperature are populated
predominantly by thermal particles in the bulk of the velocity distribution,
whereas the heat flux is sourced primarily by suprathermal particles in the
tail of the distribution. In a plasma, faster particles have longer MFPs,
therefore non-local effects are most apparent in the tail of the DF where
particle MFPs may be comparable to length scales on which hydrodynamic
quantities vary. This leads to a regime where the bulk is dominated by
collisions and near-Maxwellian, leaving hydrodynamic quantities undisturbed,
but where the tail of the DF may experience non-local kinetic effects,
perturbing the heat flux. A set of fluid equations describing evolution of
hydrodynamic quantities in such a regime may then remain a valid description
of the plasma, but with a non-local prescription for the heat flux.

In the local limit, ion transport coefficients for a simple plasma of
electrons and a single ion species are provided by Braginskii [15].
Comprehensive formalisms and practical prescriptions for calculations of local
ionic transport processes in an unmagnetized plasma with multiple ion species
are given by Chapman and Cowling [16], Ferziger and Kaper [17], Zhdanov [18],
with compact formulae provided by Kagan and Baalrud [19], and Kagan and Tang
[20, 21].

The local result for the heat flux is generally a linear combination of
gradients of hydrodynamic quantities such as temperature, but extending to
other quantities such as concentration and partial pressure in multi-species
cases. For a given density and temperature, there is a maximum possible heat
flux corresponding to all particles streaming freely at the thermal speed. For
a sufficiently sharp gradient, the local result will exceed the free-streaming
limit, thereby become unphysical. Non-local effects then have a tendency to
inhibit heat flux to remain physical, which is supported by simulations and
Fokker-Planck (FP) codes [22]. This observation motivates a flux-limiter
approach [23, 24], where local transport coefficients are limited in non-local
cases by a fitted function compared to full FP methods. Although
computationally cheap, flux limiters are physically inaccurate and fail to
capture other features of non-local transport such as preheating in cold
regions. More advanced models, such as the Schurtz-Nicolaï-Busquet (SNB)
approach [25], provide a more physically accurate description of non-local
transport for electrons but not ions. For ion transport, this leaves flux
limiters as the usual practical option implemented in hydrodynamics codes.
Chapman et al. [13] indicates that flux limiters may be insufficient for
modelling non-local ion heat transport accurately and therefore may be a major
contributor to discrepancies in modelling ICF systems of interest, in
particular reproducing neutron yields from experiments on uniaxially driven
targets [26, 27].

In this work, we present a reduced kinetic method (RKM) based on approaches by
McDevitt et al. [28, 29, 30] to describe non-local behaviour of the ion DFs,
and therefore ionic heat flux, in a multi-species plasma. The RKM solves only
for the tail of the DF where non-local effects are expected to be most
pronounced. Further noting that fewer particles reside on the tail compared to
the bulk, a linearized Coulomb collision operator is valid even for order-
unity deviations from a Maxwellian distribution in the tail. Transport
quantities such as the heat flux can then be computed from the DFs provided by
the RKM solver. This leaves a novel first-principles model for non-local ion
transport, which as of now has not been widely studied even for a simple,
single ion species plasma. Importantly, this study is conducted in the multi-
species context where ionic heat flux includes additional features not present
in the single ion species or electron cases. In particular, an additional
contribution to the heat flux exists due to relative diffusion of different
ion species which is generally of comparable magnitude to the standard
conductive component. Multi-species scenarios also allow gradients of partial
pressures and concentrations act as drives for transport fluxes along with the
standard drive due to the temperature gradient. The impact of non-locality on
both of these multi-species effects is for the first time investigated. We
consider two components of the heat flux, the reduced heat flux and the
enthalpy flux, and find that non-local deviations of the latter can be
insensitive to sharp pressure and concentration gradients, therefore
increasingly important in comparison to the reduced heat flux in non-local
scenarios. Many hydrodynamics codes treat ion mixtures as a single effective
species, thereby ignore the enthalpy component of the heat flux which may be
particularly significant with sharp concentration gradients and therefore
influence fluid simulations.

The remainder of this paper is organized as follows; in 2, the mathematical
model as derived from the VFP equation and numerical implementation is
presented. In 3, the additional complexity of multi-species heat flux in
comparison to the single-species case is described. In 4, results from the RKM
are presented for a single species and weakly and strongly asymmetric ion
mixtures with gradients in temperature, pressure, and concentration, as well
as consistency checks, convergence studies, and comparisons to SNB and flux-
limited models.

##  2 Reduced kinetic method

Ideal, weakly coupled plasmas may be described by the VFP equation, which is
the Boltzmann equation for a plasma with a Fokker-Planck collision operator as
a consequence of long-range Coulomb collisions with multiple species. The VFP
equation for a species $\alpha$ in laboratory-frame phase-space coordinates
$(t,\boldsymbol{x},\boldsymbol{v})$ for particle position $\boldsymbol{x}$ and
particle velocity $\boldsymbol{v}$ is

$\partial_{t}f_{\alpha}+\boldsymbol{v}\cdot\boldsymbol{\nabla}f_{\alpha}+\boldsymbol{F}_{\alpha}\cdot\boldsymbol{\nabla}_{\boldsymbol{v}}f_{\alpha}=\sum_{\beta}C_{\alpha\beta}\\{f_{\alpha},f_{\beta}\\}\,,$
(1)

where $f_{\alpha}$ is the DF of species $\alpha$ and $\boldsymbol{F}_{\alpha}$
is the force per unit mass acting on species $\alpha$. The FP collision
operator $C_{\alpha\beta}\\{f_{\alpha},f_{\beta}\\}$ represents collisions of
species $\beta$ acting on species $\alpha$, with the sum over all species
$\beta$ including $\alpha$. Ion-electron collisions are almost always
negligible compared to ion-ion collisions and are therefore neglected,
although they may be readily included.

It is convenient to decompose each species’ DF into a Maxwellian component and
its deviation from Maxwellian as $f_{\alpha}=f^{M}_{\alpha}+\delta
f_{\alpha}$. A Maxwellian distribution of a species is given by
$f^{M}_{\alpha}=n_{\alpha}/(\sqrt{\pi}v_{\text{th},\alpha})^{3}\exp(-(|\boldsymbol{v}-\boldsymbol{u}|/v_{\text{th},\alpha})^{2})$,
where $n_{\alpha}$ is the species’ number density, $T_{\alpha}$ is the
species’ temperature in units of energy, and
$v_{\text{th},\alpha}\equiv\sqrt{2T_{\alpha}/m_{\alpha}}$ is the species’
thermal velocity with particle mass $m_{\alpha}$. The centre-of-momentum
velocity, or hydrodynamic velocity, is defined as
$\boldsymbol{u}=\sum_{\alpha}c_{\alpha}\boldsymbol{u}_{\alpha}$, where
$\boldsymbol{u}_{\alpha}=\frac{1}{n_{\alpha}}\int d^{3}\boldsymbol{v}\
\boldsymbol{v}f_{\alpha}$ is the the net flow velocity of a species $\alpha$.
The difference between these different velocities and the hydrodynamic
velocity is typically much smaller than their thermal speeds
$|\boldsymbol{u}_{\alpha}-\boldsymbol{u}|\ll v_{\text{th},\alpha}$ for an ion
mixture near local thermal equilibrium, therefore the Maxwellian drift
velocity for all species may be chosen to be $\boldsymbol{u}$ instead of
$\boldsymbol{u}_{\alpha}$ for convenience. Although different species may have
their own temperatures in the mathematical framework, ion species with
comparable masses and different temperatures is indicative of a system far
from local thermal equilibrium, therefore a common ion temperature
$T_{\alpha}=T$ is taken. Using bilinearity of the FP operator, and noting that
Maxwellians satisfy $C_{\alpha\beta}\\{f^{M}_{\alpha},f^{M}_{\beta}\\}=0$ for
a common temperature and net flow velocity, the collision operator can be
decomposed as

$\displaystyle C_{\alpha\beta}\\{f_{\alpha},f_{\beta}\\}=$
$\displaystyle\,C_{\alpha\beta}\\{\delta
f_{\alpha},f^{M}_{\beta}\\}+C_{\alpha\beta}\\{f^{M}_{\alpha},\delta
f_{\beta}\\}$ (2) $\displaystyle+C_{\alpha\beta}\\{\delta f_{\alpha},\delta
f_{\beta}\\}$

The first term on the RHS of (2) $C_{\alpha\beta}^{T}\\{\delta
f_{\alpha}\\}=C_{\alpha\beta}\\{\delta f_{\alpha},f^{M}_{\beta}\\}$ is the
‘test-particle’ operator which describes collisions of species $\alpha$
against a Maxwellian background of species $\beta$. The second term
$C_{\alpha\beta}^{F}\\{\delta
f_{\beta}\\}=C_{\alpha\beta}\\{f^{M}_{\alpha},\delta f_{\beta}\\}$ is the
‘field-particle’ operator which describes the effect of species $\beta$ on the
Maxwellian part of species $\alpha$. It can be argued that the third, non-
linear term $C_{\alpha\beta}\\{\delta f_{\alpha},\delta f_{\beta}\\}$ is
subleading to the test-particle term $|C_{\alpha\beta}\\{\delta
f_{\alpha},\delta f_{\beta}\\}|\ll|C^{T}_{\alpha\beta}\\{\delta
f_{\alpha}\\}|$ even for order-unity deviations from Maxwellian in the tail of
the distribution, since $\delta f_{\beta}(\boldsymbol{v}^{\prime})$ is smaller
than $f^{M}_{\beta}(\boldsymbol{v}^{\prime})$ globally in
$\boldsymbol{v}^{\prime}$ space, and the second argument of $C_{\alpha\beta}$
is integrated over $\boldsymbol{v}^{\prime}$ [31]. Hence, the full kinetic
equation is reduced to a linear problem with the collision operator

$\displaystyle C_{\alpha\beta}\\{f_{\alpha},f_{\beta}\\}\approx$
$\displaystyle\,C_{\alpha\beta}^{T}\\{\delta
f_{\alpha}\\}+C_{\alpha\beta}^{F}\\{\delta f_{\beta}\\}\,,$ (3)

which is numerically tractable while still capturing relevant kinetic physics.

Figure 1:  Visualization of 1D2V geometry at a boundary. There is an azimuthal
symmetry in the velocity coordinate $\phi$, and translational symmetry
parallel to the boundary.

It is often the case that the length scale of a boundary in a fluid system is
much shorter than its spatial extent in the other two directions parallel to
the boundary, therefore the geometry near the boundary may be considered to be
near one-dimensional planar (1D). This requires that any scalar quantities of
the problem (such as temperature and number densities) only vary in the
direction perpendicular to the boundary. In a one-dimensional spatial planar
geometry, there is an azimuthal symmetry around the single spatial direction,
therefore the velocity space becomes two dimensional (2V), as depicted in
figure 1.

Transport quantities such as heat flux are defined in the centre-of-momentum
frame. A gives the kinetic equation for the deviation of the DF from
Maxwellian $\delta f_{\alpha}=f_{\alpha}-f^{M}_{\alpha}$ in the centre-of-
momentum frame in a 1D2V planar geometry as

$\displaystyle\partial_{t}\delta f_{\alpha}$
$\displaystyle=\sum_{\beta}\bigg{(}C^{T}_{\alpha\beta}\\{\delta
f_{\alpha}\\}+C^{F}_{\alpha\beta}\\{\delta f_{\beta}\\}\bigg{)}-\mu
w\partial_{z}\delta f_{\alpha}$ (4)
$\displaystyle-\bigg{(}\frac{1}{\rho}\partial_{z}p+F_{\alpha}-\sum_{\beta}c_{\beta}F_{\beta}\bigg{)}\bigg{(}\mu\partial_{w}+\frac{1-\mu^{2}}{w}\partial_{\mu}\bigg{)}\delta
f_{\alpha}$ $\displaystyle-\mu
wf^{M}_{\alpha}\bigg{(}\frac{n}{n_{\alpha}}d_{\alpha}+\bigg{(}\frac{m_{\alpha}}{2T}w^{2}-\frac{5}{2}\bigg{)}\partial_{z}\ln
T\bigg{)}\,,$

where $z$ is the single planar spatial coordinate,
$\boldsymbol{w}=\boldsymbol{v}-\boldsymbol{u}$ is the velocity in the centre-
of-momentum frame (peculiar velocity), $w\equiv|\boldsymbol{w}|$ is particle
speed in the centre-of-momentum frame, and
$\mu\equiv\boldsymbol{w}\cdot\boldsymbol{\hat{z}}/w=\cos\vartheta$ is the
cosine of the velocity pitch angle $\vartheta$ in the centre-of-momentum
frame. The diffusive driving force $d_{\alpha}$ is defined in Ferziger and
Kaper [17] as
$d_{\alpha}=(\partial_{z}p_{\alpha}-c_{\alpha}\partial_{z}p)/p-c_{\alpha}\rho/p(F_{\alpha}-\sum_{\beta}c_{\beta}F_{\beta})$.
The concentration $c_{\alpha}=\rho_{\alpha}/\rho$ is the ratio of the species’
mass density $\rho_{\alpha}=n_{\alpha}m_{\alpha}$ to the total mass density of
the mixture $\rho=\sum_{\beta}\rho_{\beta}$. The partial pressure of species
$\alpha$ is given as $p_{\alpha}=n_{\alpha}T$, and $p=\sum_{\alpha}p_{\alpha}$
is the total pressure of the ion mixture. The hydrodynamic velocity has been
neglected by assuming it is much smaller than all thermal speeds
$|\boldsymbol{u}|\ll v_{\text{th},\alpha}$. (4) will then be relaxed to
$\partial_{t}\delta f_{\alpha}=0$ since we solve for transport quantities
being functions only of the instantaneous hydrodynamic profiles.

###  2.1 Test-particle operator

The test-particle collision operator against a Maxwellian background [32] of
species $\beta$ may be written in the form

$C^{T}_{\alpha\beta}\\{f_{\alpha}\\}=\nu_{D}^{\alpha\beta}\hat{\mathcal{L}}_{\mu}f_{\alpha}+\frac{1}{w^{2}}\partial_{w}\bigg{(}\tilde{\nu}_{s}^{\alpha\beta}w^{3}f_{\alpha}+\frac{1}{2}\nu_{\parallel}^{\alpha\beta}w^{4}\partial_{w}f_{\alpha}\bigg{)}\,,$
(5)

where
$\hat{\mathcal{L}}_{\mu}f=\frac{1}{2}\partial_{\mu}[(1-\mu^{2})\partial_{\mu}f]$
is the pitch-angle scattering operator with azimuthal symmetry, and the three
collision frequencies are defined via

$\displaystyle\nu^{\alpha\beta}_{D}$
$\displaystyle=\hat{\nu}_{\alpha\beta}\frac{\text{erf}(\xi_{\beta})-G(\xi_{\beta})}{\xi_{\alpha}^{3}},$
(6) $\displaystyle\tilde{\nu}_{s}^{\alpha\beta}$
$\displaystyle=\frac{m_{\alpha}}{m_{\alpha}+m_{\beta}}\nu^{\alpha\beta}_{s}=2\hat{\nu}_{\alpha\beta}\frac{G(\xi_{\beta})}{\xi_{\alpha}},$
$\displaystyle\nu_{\parallel}^{\alpha\beta}$
$\displaystyle=2\hat{\nu}_{\alpha\beta}\frac{G(\xi_{\beta})}{\xi_{\alpha}^{3}}\,,$

where $\xi_{\alpha}=w/v_{\text{th},\alpha}$, $\text{erf}(x)$ is the error
function [33], $G(x)=[\text{erf}(x)-\frac{2}{\sqrt{\pi}}xe^{-x^{2}}]/2x^{2}$
is the Chandrasekhar function, and the basic collision frequency is

$\hat{\nu}_{\alpha\beta}=\frac{n_{\beta}Z_{\alpha}^{2}Z_{\beta}^{2}e^{4}\ln\Lambda_{\alpha\beta}}{8\pi\sqrt{2m_{\alpha}}\varepsilon_{0}^{2}T^{3/2}_{\alpha}}\,,$
(7)

where we take $\ln\Lambda_{\alpha\beta}=\ln\Lambda$ as the conventional
Coulomb logarithm for simplicity. The latter may take many different forms as
determined from theoretical arguments [34] or fits to high-fidelity plasma
microphysics simulations, however throughout the remainder of this work the
Coulomb logarithm will be treated as a tunable parameter of the RKM framework.
With the foregoing considerations the test-particle collision operator (5) is
then linear in species $\beta$, such that
$C^{T}_{\alpha}=\sum_{\beta}C^{T}_{\alpha\beta}$ can be written in the more
compact form

$C^{T}_{\alpha}\\{\delta
f_{\alpha}\\}=\nu_{D}^{\alpha}\hat{\mathcal{L}}_{\mu}\delta
f_{\alpha}+\frac{1}{w^{2}}\partial_{w}\bigg{(}\tilde{\nu}_{s}^{\alpha}w^{3}\delta
f_{\alpha}+\frac{1}{2}\nu_{\parallel}^{\alpha}w^{4}\partial_{w}\delta
f_{\alpha}\bigg{)}\,,$ (8)

where $\nu^{\alpha}=\sum_{\beta}\nu^{\alpha\beta}$ for each of the three
collision frequencies $\nu_{D}^{\alpha}$, $\tilde{\nu}^{\alpha}_{s}$, and
$\nu_{\parallel}^{\alpha}$. This collision operator appears similar to the
single species form (i.e., if the collision operator (5) had only a single
species present), but where the collision frequencies may have a more
complicated $v$ dependence in asymmetric mass and charge cases due to the
error and Chandrasekhar functions appearing in (6).

###  2.2 Field-particle operator

Rosenbluth, MacDonald, and Judd [35, 36] give the field-particle operator in
1D2V coordinates as

$\displaystyle C^{F}_{\alpha\beta}\\{\delta
f_{\beta}\\}=4\pi\bigg{(}\frac{q_{\alpha}q_{\beta}}{m_{\alpha}}\bigg{)}^{2}\ln\Lambda_{\alpha\beta}\frac{m_{\alpha}f^{M}_{\alpha}}{T_{\alpha}}\Biggl{[}\frac{4\pi
T_{\alpha}}{m_{\beta}}\delta f_{\beta}$ (9)
$\displaystyle-\mathcal{H}_{\beta}+\bigg{(}\frac{m_{\alpha}}{m_{\beta}}-1\bigg{)}w\partial_{w}\mathcal{H}_{\beta}+\frac{m_{\alpha}w^{2}}{2T_{\alpha}}\partial_{w}^{2}\mathcal{G}_{\beta}\Biggl{]},$

where the Rosenbluth potentials $\mathcal{H}_{\beta}$ and
$\mathcal{G}_{\beta}$ and their appropriate derivatives are given in terms of
$\delta f_{\beta}$ as Legendre series over $\mu$ as described in B.

The test-particle operator contains only velocity derivatives on $\delta
f_{\alpha}$, therefore is local in velocity space and does not couple the
deviations of different species. In contrast, the field-particle operator
contains integrals over velocity space and acts on $\delta f_{\beta}$,
therefore is non-local in velocity space and couples the deviations of
different species. This makes the field-particle operator more difficult to
include analytically and numerically than the test-particle operator, yet
retaining the field-particle operator is required for agreement with analytic
results in the local limit such as Braginskii’s heat flux results, as well as
for other important features such as momentum conservation.

Particles in the tail mainly experience collisions with the bulk of the
distribution since it is the most populated region in velocity space. Since
the tail is located far away in velocity space from the bulk, tail particles
do not resolve the non-Maxwellian part of the bulk well and therefore
effectively collide with a Maxwellian background. Such collisions are
represented by the test-particle operator, therefore it is expected that the
test-particle operator dominates over the field-particle operator in the tail
of the distribution, which is the region of phase space which we are primarily
interested in for non-local behaviour. This feature allows for a tractable
numerical implementation of the field-particle operator as discussed in 2.4.

###  2.3 Expansion in Legendre polynomials

Similarly to previous work [37, 28], $\delta f_{\alpha}$ may be expanded in
Legendre polynomials as

$\delta f_{\alpha}(z,w,\mu)=\sum_{l=0}^{\infty}\delta
f_{\alpha}^{(l)}(z,w)P_{l}(\mu)\,,$ (10)

where the Legendre components

$\delta f^{(l)}_{\alpha}(z,w)=\frac{2}{2l+1}\int_{-1}^{1}d\mu\
P_{l}(\mu)\delta f_{\alpha}(z,w,\mu)\,$ (11)

which follows from orthogonality of Legendre polynomials. Legendre polynomials
are a natural basis to expand in over $\mu$ since they are the azimuthally
symmetric versions of spherical harmonics which provide a complete orthogonal
basis over the unit sphere. Integrating (4) over
$\frac{2}{2m+1}\int_{-1}^{1}d\mu\ P_{m}(\mu)$ and making use of the
orthogonality, recursion relations, and eigenvalue property
$\hat{\mathcal{L}}_{\mu}P_{m}=-\frac{1}{2}m(m+1)P_{m}$ of Legendre polynomials
yields the set of PDEs

$\displaystyle\partial_{t}\delta f^{(l)}_{\alpha}$
$\displaystyle=C_{\alpha}^{T,(l)}\\{\delta
f_{\alpha}^{(l)}\\}+\sum_{\beta}C_{\alpha\beta}^{F,(l)}\\{\delta
f_{\beta}^{(l)}\\}$ (12)
$\displaystyle-\biggl{(}w\partial_{z}+\mathcal{F}_{\alpha}\partial_{w}\biggl{)}\biggl{(}\frac{l}{2l-1}\delta
f^{(l-1)}_{\alpha}+\frac{l+1}{2l+3}\delta f^{(l+1)}_{\alpha}\biggl{)}$
$\displaystyle-\frac{1}{w}\mathcal{F}_{\alpha}\biggl{(}\frac{l(l-1)}{2l-1}\delta
f^{(l-1)}_{\alpha}-\frac{(l+2)(l+1)}{2l+3}\delta f^{(l+1)}_{\alpha}\biggl{)}$
$\displaystyle-
wf^{M}_{\alpha}\bigg{(}\frac{n}{n_{\alpha}}d_{\alpha}+\bigg{(}\frac{m_{\alpha}}{2T}w^{2}-\frac{5}{2}\bigg{)}\partial_{z}\ln
T\bigg{)}\delta^{l,1}\,,$

where
$\mathcal{F}_{\alpha}=\frac{1}{\rho}\partial_{z}p+F_{\alpha}-\sum_{\beta}c_{\beta}F_{\beta}$,
the form of $C_{\alpha\beta}^{F,(l)}\\{\delta f_{\beta}^{(l)}\\}$ is given in
B, and $\delta^{l,k}$ is the Kronecker delta.

Truncating the expansion to the first $N_{l}$ Legendre modes, the set of PDEs
may be solved numerically. Fine detail of $\delta f_{\alpha}$ in the pitch
angle $\mu$ is not expected, so the series should converge for a reasonable
number of Legendre modes. This is checked with a convergence study over
$N_{l}$, as well as having been demonstrated in previous work by McDevitt et
al. [29] for a similar problem by comparing the solutions of the system of
PDEs for the Legendre components to the full solution of the 3D equation
without expanding in Legendre polynomials over $\mu$.

Figure 2:  Visualization of solver domain for the set of PDEs for $\delta
f_{\alpha}^{(l)}$.

###  2.4 Numerical implementation

The set of PDEs (12) is solved by discretizing $\delta f_{\alpha}^{(l)}(z,w)$
for all species together to a 4D array $F_{pq}^{s,l}$, where $z=p\Delta z$,
$w=q\Delta w$, $l$ is the Legendre index, and $s$ is the species index.
Derivatives are estimated via the finite differences

$\displaystyle\partial_{z}F_{p,q}^{s,l}=(F_{p+1/2,q}^{s,l}-F_{p-1/2,q}^{s,l})/\Delta
z$ (13)
$\displaystyle\partial_{w}F_{p,q}^{s,l}=(F_{p,q+1/2}^{s,l}-F_{p,q-1/2}^{s,l})/\Delta
w\,,$

where the cell-face values are interpolated via a linear scheme

$\displaystyle F_{p+1/2,q}^{s,l}=$
$\displaystyle\,(1-\delta^{z}_{p+1/2,q})F_{p+1,q}^{s,l}+\delta^{z}_{p+1/2,q}F_{p,q}^{s,l}$
(14) $\displaystyle F_{p,q+1/2}^{s,l}=$
$\displaystyle\,(1-\delta^{w}_{p,q+1/2})F_{p,q+1}^{s,l}+\delta^{w}_{p,q+1/2}F_{p,q}^{s,l}\,.$

The coefficients $\delta^{z}_{p,q}$ and $\delta^{w}_{p,q}$ (which are
functions of grid position) specify the linear interpolation scheme.
$\delta^{z}_{p,q}=1/2$ is chosen for a spatial central difference and a
‘Chang-Cooper’-inspired [38] scheme $\delta^{w}_{p,q}=1/u-1/(e^{u}-1)$, where
$u=2(\Delta w)^{2}q/v_{\text{th},p}^{2}$ is chosen for velocity derivatives.
This choice improves stability and convergence with respect to a simpler
finite difference scheme with fixed $\delta^{w}_{p,q}$. Reshaping the mesh
$F_{pq}^{s,l}$ into a vector $\boldsymbol{F}$, the set of PDEs is discretized
to $\partial_{t}\boldsymbol{F}=M(\boldsymbol{F})+\boldsymbol{P}$, where $M$ is
the linear operator representing the discretized right side of (12) acting on
$\delta f^{(l)}_{\alpha}$, and $\boldsymbol{P}$ is the last inhomogeneous term
reshaped as a vector. The stationary solution may be found via an explicit or
implicit relaxation method. Here, SciPy’s backward-differentiation formulas
implicit integrator [39, 40] is used to relax to the steady-state solution.

As discussed in 2.2, the field-particle term is negligible in the tail
compared to the test-particle operator. This is computationally convenient
since the Jacobian of the field-particle operator is non-sparse in the
velocity index, whereas the Jacobian of the test-particle operator is. For the
solver, the RHS of (12) without the non-sparse integral terms appearing in the
field-particle operator may then be used as an approximate Jacobian, with the
full field-particle operator entering along with the full timestep function.
Although the field-particle operator requires computing
$\Psi_{\beta,(1,2,3,4)}^{(l)}(w)$ integrals as defined in B for each $l$,
species $\beta$, and spatial position, the integrals are cumulative in
velocity, therefore may be computed with reasonable efficiency (see B).
Alternatively, the Jacobian of the field-particle term can be included, and
these two approaches have been found to be in good agreement.

The leading-order Chapman-Enskog (CE) solution [16, 17] is chosen for the the
boundary conditions for the DF at spatial edges (where spatial gradients
become shallow, therefore transport should become local) and at the low tail
cutoff ($w_{\text{min}}=\xi_{\text{min}}v_{\text{th},\alpha}$, with
$\xi_{\text{min}}\sim 1$). The spatial boundary conditions are motivated
physically by arguing that over a sufficient length, hot particles responsible
for non-local effects will have had many collisions and will not penetrate
much further away from the boundary, so the local limit should be restored.
The local solution is exponentially small for the choice of functions of the
bulk profiles, so choosing $f=f^{M}_{\alpha}$ as a spatial boundary condition
is nearly identical. However, the local perturbative solution does provide a
useful boundary condition for the tail cutoff, since $\delta f_{\alpha}$ would
otherwise be zero for $w\leq w_{\text{min}}$, which would remove possibly
important parts of the velocity integrand when calculating fluxes. Since the
bulk particles have shorter MFPs, we still expect the bulk to match the local
solution. At some maximum velocity cutoff above which contributions to the
heat flux integrands are expected to be negligible, $\delta f_{\alpha}=0$ is
chosen due to the fast exponential decay.

##  3 Heat flux in a multi-species plasma

Considering the heat flux in a multi-species plasma has further complications
in comparison to the single species case. A single species may only experience
self-collisions whereas multi-species transport requires also considering
collisions with other species. Furthermore, different species in mixtures
diffuse in the centre-of-momentum frame which gives rise to an additional
source of heat transport of generally comparable magnitude. These additional
complications in the multi-species case are also apparent mathematically since
there is a set of coupled kinetic equations to solve instead of one in the
single species case. Experimental evidence for various types of ICF implosions
indicates that multi-species effects can directly affect the yield [41, 42,
43, 44], further motivating the present study.

For a plasma with $N$ ion species, there are $N$ coupled kinetic equations
which evolve each species’ DF. From this set of kinetic equations, a set of
$5N$ fluid equations may be derived for each species’ mass density, flow
velocity, and temperature. Since the timescales on which different ion species
exchange momentum and energy are much shorter than other (often hydrodynamic)
timescales of interest, the ions typically have a single common temperature
and similar mean flow. There may be cases with highly asymmetric mass mixtures
where such an approximation is not valid, but these are not considered in the
present work. Such extensions could be considered in the future if the need
arises. The separate temperature evolution equations may be combined to form a
single evolution equation for the common ion temperature, where the heat flux
enters only by the species-summed total heat flux
$\boldsymbol{Q}=\sum_{\alpha}\boldsymbol{Q}_{\alpha}$.

The heat flux for a species $\alpha$ is defined in the centre-of-momentum
frame as $\boldsymbol{Q}_{\alpha}=\int d^{3}\boldsymbol{w}\
\frac{1}{2}m_{\alpha}w^{2}\boldsymbol{w}f_{\alpha}$. In a simple plasma
featuring electrons and a single ion species, there is no heat flux associated
with the flow since the only ion species velocity is identical to the
hydrodynamic velocity. On the other hand, in a multi-species plasma, each
species $\alpha$ has an individual diffusive flow defined as
$\boldsymbol{V}_{\alpha}=\boldsymbol{u}_{\alpha}-\boldsymbol{u}$ relative to
the hydrodynamic velocity. This introduces a component to the heat flux in
multi-species cases due to the diffusion of different species in the centre-
of-momentum frame transporting energy, which does not occur in the single-
species case. The heat flux $\boldsymbol{Q}$ may be decomposed into two parts
via $\boldsymbol{Q}=\boldsymbol{q}+\boldsymbol{h}$, the reduced heat flux and
enthalpy flux respectively, defined as

$\displaystyle\boldsymbol{q}=$ $\displaystyle\,\sum_{\alpha}\int
d^{3}\boldsymbol{w}\
\biggl{(}\frac{1}{2}m_{\alpha}w^{2}-\frac{5}{2}T\biggl{)}\boldsymbol{w}f_{\alpha}\,,$
(15) $\displaystyle\boldsymbol{h}=$
$\displaystyle\,\sum_{\alpha}\frac{5}{2}T\int d^{3}\boldsymbol{w}\
\boldsymbol{w}f_{\alpha}\,.$

Although it is not necessary to decompose
$\boldsymbol{Q}=\boldsymbol{q}+\boldsymbol{h}$, it provides more points of
comparison for the RKM with the analytic result in the local limit.

The CE expansion [16, 17] of multi-species Boltzmann equations yields local
results for transport quantities, in particular the heat fluxes in the centre-
of-momentum frame

$\displaystyle\boldsymbol{q}^{\text{local}}=-\lambda^{\prime}\boldsymbol{\nabla}T-nT\sum_{\alpha}D^{(T)}_{\alpha}\boldsymbol{d}_{\alpha}\,,$
(16)

and

$\displaystyle\boldsymbol{h}=\frac{5}{2}\sum_{\alpha}n_{\alpha}T\boldsymbol{V}_{\alpha}\,,$
(17)

in which

$\displaystyle\boldsymbol{V}_{\alpha}^{\text{local}}=-D^{(T)}_{\alpha}\boldsymbol{\nabla}\ln
T-\sum_{\beta}D_{\alpha\beta}\boldsymbol{d}_{\beta}\,,$ (18)

where $\lambda^{\prime}$, $D^{(T)}_{\alpha}$, and $D_{\alpha\beta}$ are the
partial coefficient of thermal conductivity, multi-species thermo-diffusion
coefficients, and multi-species (ordinary) diffusion coefficients
respectively, with the sign convention following Ferziger and Kaper [17].
These transport coefficients are determined from the collision operator
$C_{\alpha\beta}$. For the full Coulomb collision operator with Debye
shielding, these may be provided analytically. Here the weakly coupled limit
results of [19] are used, which provide numerical routines for the multi-
species coefficients.

In the local single ion species case, the only two gradients available to
induce heat flux are the temperature and pressure gradients. By assuming force
balance, the single-species ion heat flux may be expressed purely in terms of
the temperature gradient. In contrast, the multi-species case involves more
gradients on which the heat flux may be dependent, such as concentration
gradients. In addition, for single species cases the force $\boldsymbol{F}$
does not appear explicitly in the kinetic equation for $\delta f$ due to force
balance. For multi-species cases, the force $\boldsymbol{F}_{\alpha}$ may vary
for different species, such as different specific charges
$eZ_{\alpha}/m_{\alpha}$, giving different electrostatic forces on each
species, therefore these forces will appear explicitly in the kinetic equation
for $\delta f_{\alpha}$ and therefore affect the heat flux.

In 1D2V coordinates, the diffusive velocities in the $\boldsymbol{\hat{z}}$
direction are $n_{\alpha}V_{\alpha}=\frac{4\pi}{3}\int_{0}^{\infty}dw\
w^{3}\delta f^{(1)}_{\alpha}$, and the heat fluxes are

$\displaystyle q=$
$\displaystyle\,\frac{2\pi}{3}\,\sum_{\alpha}m_{\alpha}\int_{0}^{\infty}dw\
w^{3}\biggl{(}w^{2}-\frac{5}{2}v_{\text{th}}^{2}\biggl{)}f_{\alpha}^{(1)},$
(19)

and

$\displaystyle h_{\alpha}$
$\displaystyle=\sum_{\alpha}\frac{5}{2}n_{\alpha}TV_{\alpha}\,,$ (20)

which vanishes for a single species since $u_{\alpha}=u$.

##  4 Results

###  4.1 Single species case

Consider a single species with a number density profile $n(z)$ and a
temperature step of the form

$T(z)=\frac{T_{\text{hot}}+T_{\text{cold}}}{2}+\frac{T_{\text{hot}}-T_{\text{cold}}}{2}\text{erf}\left(\dfrac{z}{L}\right)\,,$
(21)

which transitions from $T_{\text{hot}}$ to $T_{\text{cold}}$ with increasing
$z$ over some spatial scale $L$. $n_{0}=n_{\text{total}}(z=-\infty)$,
$T_{0}=T(z=-\infty)$, $m_{0}=\min_{\alpha}(m_{\alpha})$, and
$v_{0}=\sqrt{2T_{0}/m_{0}}$ are introduced as the characteristic number
density, temperature, mass, and velocity to which quantities are appropriately
non-dimensionalized. In particular, heat fluxes are normalized to a
characteristic heat flux $Q_{c}=n_{0}T_{0}v_{0}$, DFs are normalized to
$n_{0}v_{0}^{-3}$, and operators appearing in the kinetic equation are
normalized to a characteristic collision frequency
$\nu_{0}=n_{0}e^{4}\ln\Lambda/4\pi\epsilon_{0}^{2}m_{0}^{2}v_{0}^{3}$. The
Knudsen number is introduced as $N_{K}=\lambda_{0}/L$ as a measure of non-
locality, where
$\lambda_{0}=v_{0}/\nu_{0}=16\pi\varepsilon_{0}^{2}T_{0}^{2}/n_{0}e^{4}\ln\Lambda$
is the MFP of a thermal particle at $z=-\infty$. The RKM solver is run for a
temperature step $T_{\text{hot}}/T_{\text{cold}}=10$ and a variety of Knudsen
numbers from the local limit to non-local cases. The electric field has been
set to zero to focus on the non-local effects due to the temperature profile,
although it can be readily included.

(a) Isobaric case

(b) Isochoric case

Figure 3: Reduced heat flux $q$ (rescaled by the Knudsen number, $N_{K}$) for
a single species isobaric (a) and isochoric (b) temperature step for a variety
of Knudsen numbers from the RKM solver. The local analytic result for the
reduced collision operator and Braginskii’s exact result for the full Coulomb
collision operator are included. The enthalpy flux has not been included since
it is identically zero for a single species, and the temperature and number
density profiles are shown below on the same horizontal scale for reference.

Figure 3 shows the reduced heat flux for an isobaric ($p=nT=\text{constant}$)
and isochoric ($n=\text{constant}$) temperature step. In the local limit the
two cases produce the same result and agree well with the local Braginskii
result. The latter is also the same in the isochoric and isobaric case. In
non-local cases, the isochoric heat flux forms a pedestal in the cold end,
whereas the isobaric heat flux does not. This is expected since the higher
density in the cold end of the isobaric case gives a higher collisionality,
therefore fast particles streaming from the hot end do not penetrate as far.

(a) $f^{(0)}$ in hot end

(b) $f^{(0)}$ in cold end

(c) $f^{(1)}$ in hot end

Figure 4: DF comparison at the hot end $z=-1.5L$ and cold end $z=+1.5L$ for
various Knudsen numbers against local thermal speeds. Note the plots of
$f^{(1)}$ have a symmetric log scale to include negative values. The smaller
subplots show the ratio of each graph to the local Maxwellian with an
appropriate scaling against Knudsen number included. (d) $f^{(1)}$ in cold end
Figure 4: (continued)

Figure 4 shows the first two Legendre components $f^{(0)}=f^{M}+\delta
f^{(0)}$ and $f^{(1)}=\delta f^{(1)}$ of the DF at the hot end and cold end
for the single-species isochoric case. $f^{(0)}$ is the distribution function
averaged over the velocity pitch angle $\theta$. For small Knudsen number
$N_{K}=1/1000$, the RKM result agrees with the local CE solution. For
increasing Knudsen number, the tail of the DF in the hot end is further
depleted and the tail of the DF in the cold end is further enhanced. This is
expected physically at higher Knusden number as fast tail particles in the hot
end may stream over the temperature gradient into the cold end, thereby
reducing the tail population of the DF at the hot end and enhancing it in the
cold end.

Figure 5:  Deviation of the DF $\delta f(z,v,\vartheta)$ at the hot end
$z=-1.5L$ (left) and cold end $z=+1.5L$ (right) in a non-local case with
$N_{K}=0.1$ (top half of each polar plot) and for the local analytic result
(bottom half of each polar plot). Dashed circular contours indicating local
thermal velocities are included. Note the color map uses a symmetric
logarithmic scale with a linear section through zero.

$\delta f(z,w,\vartheta)$ may be reconstructed via the Legendre series to
visualize anisotropy in the DF as shown in figure 5. In the hot end, the tail
of the DF is suppressed in the $\boldsymbol{\hat{z}}\cdot\boldsymbol{w}>0$
portion (noting $\boldsymbol{\hat{z}}\parallel-\boldsymbol{\nabla}T$) as is
expected from faster tail particles with longer MFPs streaming out of the hot
end, thereby depleting its tail population. In the cold end, the tail of the
DF is enhanced in the $\boldsymbol{\hat{z}}\cdot\boldsymbol{w}>0$ portion due
to suprathermal particles streaming into the cold end from the hot end. It may
be noted that this tail enhancement appears particularly anisotropic with the
portion of the tail with $\boldsymbol{\hat{z}}\cdot\boldsymbol{w}<0$ appearing
mostly unchanged in the hot end, since the suprathermal particles from the hot
end must have been travelling in the
$\boldsymbol{\hat{z}}\cdot\boldsymbol{w}>0$ direction originally to enter the
cold end, and they are subject to little collisional scattering in the cold
end since $\nu_{D}$ is small in the tail.

###  4.2 Multi-species cases

For multi-species cases, we define the Knudsen number in a similar manner as
the single-species case, but instead using a maximum MFP
$N_{K}=\lambda_{\text{max}}/L$ defined as
$\lambda_{\text{max}}=\max_{\alpha,z}(\lambda_{0,\alpha}(z))$, with the MFP of
a species $\alpha$ at spatial position $z$ accounting for collisions with all
species via
$\lambda_{0,\alpha}=v_{\text{th},\alpha}/\sum_{\beta}\hat{\nu}_{\alpha\beta}$.
A 50:50 DT mixture and a fully ionized CH2 mixture are chosen as multi-species
examples to investigate both weakly and moderately asymmetric mixtures. Both
of these types of mixtures are relevant to ICF-related systems which often use
DT fuel and plastic layers. A CH2 mixture may be considered to be moderately
asymmetric in mass, since the masses differ by an order of magnitude. Although
the charge ratio is only $Z_{C}/Z_{H}=6$, the charges enter the equations
through the collision frequencies $\nu_{\alpha\beta}\propto
Z_{\alpha}^{2}Z_{\beta}^{2}$, therefore the charge asymmetry is effectively on
an order of magnitude. Three cases are chosen for each mixture: a temperature
gradient at uniform pressure and concentration, a pure pressure gradient at
uniform temperature and concentration, and a pure concentration gradient at
uniform temperature and pressure. These cases with isolated $\nabla T$,
$\nabla p$, and $\nabla c$ are selected to investigate the sensitivity of the
locality of heat fluxes to different driving gradients in temperature,
pressure, and concentration. Similarly to the single-species case, the
electric field has been set to zero to focus on the non-local effects from the
aforementioned gradients. These six cases are shown in figures 6 and 7 for DT
and CH2 respectively.

In a two-species ion mixture ignoring electric fields and electrons, the
diffusive mass flux of the lighter species in the local limit may be written
in the Landau form [45, 21]

$\displaystyle\rho_{l}\boldsymbol{V}_{l}=-\rho
D\bigg{(}\boldsymbol{\nabla}c+k_{p}\boldsymbol{\nabla}\ln
p+k_{T}\boldsymbol{\nabla}\ln T\bigg{)}\propto\boldsymbol{h},$ (22)

where $D$ is the classical diffusion coefficient, and $k_{p}$ and $k_{T}$ are
the baro-diffusion and thermo-diffusion coefficients respectively. It has been
indicated that for a two-species mixture, the enthalpy flux $\boldsymbol{h}$
is directly proportional to the diffusive flux $\boldsymbol{V}_{l}$, therefore
plots which show the enthalpy flux are directly related to the diffusive flux.

(a) DT $\nabla T$

(b) DT $\nabla p$

Figure 6: Components of heat flux for a DT mixture subject to an isolated
temperature gradient (a), isolated pressure gradient (b), and isolated
concentration gradient (c). The shape of the relevant profiles are shown at
the bottom of each subfigure, with undisplayed profiles being fixed constant
at unity or concentrations at $0.5$. Results for $N_{K}\leq 0.01$ are not
visible due to agreement with the local result. (c) DT $\nabla c$ Figure 6:
(continued)

(a) CH2 $\nabla T$

(b) CH2 $\nabla p$

Figure 7: Components of heat flux for a CH2 mixture subject to an isolated
temperature gradient (a), isolated pressure gradient (b), and isolated
concentration gradient (c). The shape of the relevant profiles are shown at
the bottom of each subfigure, with undisplayed pressure or temperature
profiles being fixed constant at unity or concentrations such that
$n_{H}/n_{C}=2$. (c) CH2 $\nabla c$ Figure 7: (continued)

All of the multi-species cases show the same basic qualitiative behaviour as
the single-species cases; increasing Knudsen number inhibits the peak heat
flux and enhances the heat flux in some regions away from the peak. However,
the multi-species cases do demonstrate additional behaviour of interest.
Firstly, non-local effects on the heat flux are demonstrated in the absence of
a temperature gradient, induced instead due to sharp gradients in pressure or
concentration. In many of the cases shown, non-local effects give rise to
large heating pedestals in regions $|z/L|\gtrsim 2$ away from the boundary due
to sharp gradients in $T$, $p$, or $c$. This is in dramatic contrast to the
local heat flux in these regions, which is exponentially small for the chosen
$\sim\tanh(z/L)$ profiles. In most cases this enhancement only appears in one
direction to the boundary, for example the $\nabla p$ cases 6b 7b where the
non-local behaviour is only pronounced for $z/L>0$. In contrast, the DT
$\nabla c$ 6c case has significant non-local behaviour in both directions due
to the near-symmetry, which is not present in any of the other cases.

In the temperature gradient cases 6a 7a, both components of the heat flux are
similarly sensitive to non-local behaviour, and appear similar to the single-
species case. Since the enthalpy flux is proportional to the diffusive flux,
non-local behaviour in the former imply that diffusion can be sensitive to
kinetic effects at similar Knudsen numbers at which heat transport is
sensitive. In contrast, the DT $\nabla p$ case 6b and both $\nabla c$ cases 6c
7c have the reduced heat flux $q$ being significantly affected by non-local
behaviour, whereas the enthalpy flux $h$ is quite insensitive to non-local
behaviour. This leads to cases around the peak flux where $q$ is significantly
suppressed and $h$ is not, so the latter may be increasingly important and
dominant at the peak in non-local cases.

###  4.3 Electron transport

While the main focus of this paper is the ion heat transport, the newly
developed approach can be applied to the electron case as well by treating
electrons as an ion species with $Z_{e}=-1$ and $m_{e}\ll m_{i}$. Since all
other non-local heat flux models available to date are for the electron heat
flux only, we have run this calculation to be able to benchmark our method
against an earlier, electron heat flux approach.

####  4.3.1 Local comparison to Braginskii

In the local limit, Braginskii’s results for the electron heat flux may be
reproduced for a quasi-neutral isobaric simple plasma with an electric field
satisfying local electron force balance

$E=\frac{1}{en_{e}}(-\partial_{z}p_{e}+R_{e}),\quad
R^{\text{local}}_{e}=-B_{e}n_{e}\partial_{z}T_{e}$ (23)

where $R_{e}$ is the electron thermal force, with the local result including
the electron thermal force coefficient $B_{e}$ appearing in (2.9) of [15]. The
terms including flows, friction and viscosity, have been neglected. In the
local limit, the RKM reproduces Braginskii’s results successfully as shown in
figure 8.

Figure 8:  RKM results for the electron heat flux for various ionization
numbers ($Z_{i}=1,2,3,4,100$) matched with Braginskii’s results. Each dashed
black line corresponds to the Z of the nearest colored RKM line. The RKM
solver is run with two species, electron and ion, with $m_{i}/m_{e}=1000$ and
$Z_{e}=-1$. The temperature profile with Knudsen number $N_{K}=1/1000$ has a
temperature drop by a factor of $10$, and the number densities are chosen such
that the mixture is quasi-neutral and isobaric. Similarly to Braginskii, an
electric field ensuring electron force balance using the local result for the
thermal force is included.

####  4.3.2 Non-local comparison to SNB

(a) Full collision operator
Figure 9: Comparison of electron heat flux from the RKM solver, an SNB
approach, and a flux limiter approach for an isochoric simple hydrogen plasma
with a Gaussian temperature profile. (a) uses the full linearized Coulomb
collision operator, and (b) uses the Coulomb collision operator with the
field-particle part fixed at the local result from the Chapman-Enskog
solution. (b) Fixed field-particle part Figure 9: (continued)

The RKM heat flux may be compared to the heat flux as predicted by other non-
local transport models, namely the SNB approach and flux limiters. Figure 9
shows a comparison of the non-local electron heat flux from the RKM solver and
SNB approach for an isochoric simple hydrogen plasma with a Gaussian
temperature profile. For the chosen values $n_{e}=10^{26}\ \text{m}^{-3}$,
$T_{\text{max}}=200\ \text{eV}$, $T_{\text{min}}=100\ \text{eV}$, and
$\ln\Lambda=10$, the characteristic electron mean free path at the hottest
temperature is $\lambda_{0}=6.16\ \mu\text{m}$. The SNB implementation used
here follows the ‘iSNB’ algorithm from Cao et al. [46]. The RKM and SNB
hyperparameters are suitably chosen for convergence.

Braginskii’s result for the electron conductive heat flux $q^{T}_{e}$ in the
local limit chooses an electric field such that electron momentum conservation
is satisfied, which includes the electron friction and thermal forces
determined self-consistently from $\delta f_{e}$. In non-local cases, the
local expressions for electron friction and thermal forces are no longer
valid, therefore the correct electric field yielding zero current is not known
before running the solver. Root-finding methods are applied to find the
electric field in non-local cases such that the electron current, and
therefore enthalpy flux, is zero everywhere.

In addition to SNB, a flux-limited heat flux is included for comparison. The
flux limiter takes the simple harmonic-averaged form

$q_{\text{f.l.}}=\frac{q_{\text{Brag.}}}{1+|q_{\text{Brag.}}/(\alpha
q_{\text{f.s.}})|},$ (24)

where $q_{\text{f.l.}}$ is the flux-limited heat flux, $q_{\text{Brag.}}$ is
the local Braginskii heat flux result,
$q_{\text{f.s.}}=n_{e}T_{\text{min}}v_{\text{th},\text{min}}$ is a
characteristic free-streaming heat flux for the system, and $\alpha$ is the
flux-limiter parameter. Here, a range of values $0.05\leq\alpha\leq 0.15$ are
chosen. The flux-limited heat flux fails to capture enhanced heating in the
colder region in comparison to the SNB and RKM approaches.

The first case uses the full linearized Coulomb collison operator with test-
particle and field-particle operator. The second case uses a modified Coulomb
collision operator where the field-particle operator is fixed at its local
result, i.e. $C\\{\delta f\\}=C^{T}\\{\delta f\\}+C^{F}\\{\delta
f^{\text{local}}\\}$. This modification is proposed for improving
computational performance while retaining some physical accuracy as the
expensive integrals in the field-particle operator are not computed each
timestep, yet agreement in the local limit is retained. As demonstrated, this
approximation yields excellent agreement with SNB. The SNB approach is based
on a simplified collision operator, with this particular implementation based
on a Bhatnagar–Gross–Krook (BGK) type operator $C\\{\delta
f^{(1)}\\}=-\nu_{ei}(w)\delta f^{(1)}$ with $\nu_{ei}(w)\propto w^{-3}$. The
improved agreement between the SNB and RKM approaches when the latter employs
the fixed, rather than first-principle, field-particle operator suggests that
the field-particle operator, which is non-local in velocity space, may play an
important role in non-local behavior that is not captured in simpler collision
models such as BGK that are local in velocity space.

###  4.4 Convergence studies

There are a few hyperparameters of the model which should be chosen for fast
computational performance while converging to the physical solution. The two
hyperparameters of most interest are the number of terms kept in the Legendre
expansion $N_{l}$ and the minimum velocity cutoff
$\xi_{\text{min}}=w_{\text{min}}/v_{\text{th}}$.

The RKM solver is run for an isochoric single-species setup similar to the
isochoric case in 4.1 for a variety of values of $N_{l}$ with
$N_{K}=1/1000,1/20,1/10$ to investigate convergence in different cases of non-
locality. The root-mean-square (RMS) deviation
$\sigma_{q}=\frac{1}{Q_{c}N_{K}|D|}\sqrt{\int_{D}dz\ (q-q_{\text{ref.}})^{2}}$
is introduced to measure convergence of the heat flux, where $D$ is the
spatial domain and $|D|$ is the length of the domain. The reference heat flux
$q_{\text{ref.}}$ is chosen from a solution of the RKM solver with identical
parameters as $q$ apart from the hyperparameter of focus.

Figure 10:  Convergence study for the number of Legendre modes included in the
expansion $\delta f\approx\sum_{l=0}^{N_{l}-1}P_{l}(\mu)\delta f^{(l)}$. The
solver is run identically to the isochoric single species example 3 with
$\xi_{\text{min}}=1$, with the local limit results ($N_{K}=1/1000$) shown in
dashed lines and two significantly non-local cases ($N_{K}=1/20$ and
$N_{K}=1/10$) shown in dotted-dashed lines and dotted lines respectively. In
the local cases, the RKM solutions are difficult to see due to good agreement
with the CE solution. The RMS deviation is also shown. Here, $q_{\text{ref.}}$
uses $N_{l}=6$ results since more Legendre components included in the
expansion should improve the accuracy of the solver.

Figure 10 shows the convergence study for $N_{l}$. Good agreement with the
local analytic solution is practically identical in all cases as expected,
since the local solution for $\delta f^{(l)}$ is only nonzero for $l=1$. In
non-local cases, good convergence is demonstrated even for $N_{l}=4$. Larger
$N_{l}$ increases convergence time, in particular scaling with
$\mathcal{O}(N_{l})$. The default value $N_{l}=6$ is chosen.

For the low velocity cutoff, a value around $\xi_{\text{min}}\approx 1$ should
be taken since a smaller value includes the bulk in the RKM solver domain,
which gives rise to large computational cost and numerical instabilities, and
a larger value does not account for parts of the tail that have significant
contribution to the heat flux. In the region $0.9\leq\xi_{\text{min}}\leq
1.1$, there is a variation of $\sigma_{q}<0.1\%$ for $N_{K}\leq 1/10$,
therefore the heat flux is mostly insensitive to this hyperparameter for a
choice which appropriately defines the tail. This further validates the
approach of solving only for the tail of the distribution beyond some cutoff
since the RKM and local Chapman-Enskog solutions match smoothly at the low
velocity boundary. The default value $\xi_{\text{min}}=1$ is chosen.

###  4.5 Temperature self-consistency

The RKM solver takes input number density and temperature profiles and solves
for the tail of the DF of each species. In extreme non-local cases, a
sufficiently non-Maxwellian tail may result in the density and temperature
moments of the resulting DF to significantly differ from the input profiles,
leaving an inconsistent solution. Since the temperature is a higher moment
than the density, the former is more sensitive to deviation due to non-local
tail behaviour and more likely to be inconsistent. The RKM temperature is
defined as

$T^{\text{RKM}}_{\alpha}(z)=\frac{m_{\alpha}}{3n_{\alpha}}\int
d^{3}\boldsymbol{w}\ w^{2}f_{\alpha}=\frac{4\pi
m_{\alpha}}{3n_{\alpha}}\int_{0}^{\infty}dw\ w^{4}f^{(0)}_{\alpha}.$ (25)

(a) Single species

(b) DT

(c) CH2

Figure 11:  RKM temperature for the single-species isochoric case (a), DT
$\nabla T$ case (b), and CH2 $\nabla T$ case (c). The upper subplots show the
RKM integrated temperature, and the lower subplots shows the ratio of the RKM
temperature to the input temperature profile. For the multi-species cases, the
lighter species ($D$ or $H$) temperature is shown in dashed lines, the heavier
species ($T$ or $C$) temperature is shown in dotted lines, and the total ion
temperature
$T^{\text{RKM}}_{\text{total}}=\sum_{\alpha}n_{\alpha}T^{\text{RKM}}_{\alpha}/\sum_{\alpha}n_{\alpha}$
is shown in solid lines.

Figure 11 shows the RKM temperature for the previously discussed single-
species isochoric case, DT $\nabla T$ case, and CH2 $\nabla T$ case. In all
three cases, the local cases are in good agreement with the local result, with
inconsistency developing with increasing Knudsen number. Even for the extreme
cases of $N_{K}=0.1$, the total temperature is consistent to $\lesssim 5\%$
deviation. In the multi-species cases, it is the temperature of the lighter
species that is more inconsistent with the input profile.

###  4.6 Smallness of field-particle operator

(a) Hot end

(b) Cold end

Figure 12: Comparison of $l=1$ components of the test-particle operator and
field-particle operator at the hot end (a) and cold end (b) from the single-
species isochoric example. The solid lines show results for the $N_{K}=0.1$
results, and the dashed lines show the local analytic results. The lower
subplots show the ratio of field-particle operator to test-particle operator.

Figure 12 shows a comparison of the $l=1$ components of the test-particle
operator and field-particle operator at the hot end ($z/L=-1.5$) and cold end
($z/L=+1.5$) from the single-species isochoric example in 4.1 for the non-
local $N_{K}=0.1$ case against local analytic results. The vertical orange
dashed lines indicates the local low velocity cutoff, chosen to be the local
thermal velocity $w_{\text{min}}=v_{\text{th}}$. The solution for the tail
above the low velocity cutoff $w\geq w_{\text{min}}$ is determined by the RKM
solver. The solution below the low velocity cutoff in the bulk is fixed as the
local Chapman-Enskog solution $\delta f(w<w_{\text{min}})=\delta
f^{\text{local}}(w)$, therefore the test-particle operator is also identical
for the RKM and local results away from the cutoff. In contrast, the RKM and
analytic solutions of the field-particle operator may differ in this region
due to the velocity integrals that sample the solution above the low velocity
cutoff. This deviation is small due to the small population of the tail
compared to the bulk, therefore the local analytic result should still hold as
a good approximation of the bulk.

At the velocity grid point just below the low velocity cutoff there is a spike
in the test-particle operator compared to the local solution. This originates
from the $w^{-2}\partial_{w}(w^{4}\nu_{\parallel}\partial_{w}\delta f)$ term,
suggesting that the numerical scheme as implemented does not constrain
continuity of the second velocity derivative of the DF continuous at the
cutoff. This numerical artifact lies outside the domain of the RKM solver, and
continuity of the second or higher derivatives of the DF are not required,
therefore is not a point of concern.

The red highlight indicates the regions around
$w/v_{\text{th},\text{local}}=\sqrt{5/2}\approx 1.6$ where the local analytic
result $\delta f_{\text{local}}^{(1)}\propto(w^{2}/v_{\text{th}}^{2}-5/2)$
changes sign. The sudden spike of the ratio in this region is therefore not of
concern since both the operators are passing through zero. Apart from these
red regions, the field-particle operator remains small in comparison to the
test-particle operator. The Jacobian of the field-particle operator is then a
small difference between the approximate and true Jacobian, which explains why
the approach of approximating the Jacobian of the collision operator in the
tail without the non-sparse elements of the field-particle operator described
in 2.4 works numerically.

##  5 Discussion

A reduced kinetic method has been presented and applied to investigate non-
local multi-species ion heat transport for the first time. The RKM solves for
the tail of distribution functions according to the Vlasov-Fokker-Planck
equation for a classical, unmagnetized multi-species plasma in a 1D2V
geometry. In the single ion species case, the non-local ion heat flux show
qualitatively similar trends to its well explored electron counterpart. On the
other hand, novel non-local features are recovered for the multi-species ion
heat transport due to additional thermodynamic drives such as concentration
gradients and diffusion of different species in the centre-of-momentum frame.

Results for DFs and heat fluxes in local cases reproduce the analytic results
of the Chapman-Enskog expansion [16, 17, 19] for multi-species ion heat flux
and Braginskii’s [15] results for ion and electron heat flux in a simple
plasma. Non-local electron heat flux results agree well with an SNB
implementation, particularly where the field-particle operator is fixed at the
local solution rather than implemented in a first principle fashion. Expected
qualitative behaviour of non-local effects are demonstrated, in particular
suppression of the peak heat flux and enhancement of heat flux away from
regions of steep gradient. In multi-species cases, non-local effects are
demonstrated from gradients in temperature, concentration, and pressure, with
some cases giving rise to large heat fluxes far away from the boundary where
the local result is exponentially small. The multi-species heat flux has been
decomposed into two components, the reduced heat flux and enthalpy flux, where
the former is the term employed in mainline hydrodynamic codes and the latter
is an inherently multi-species feature due to relative dynamics of the species
in the center-of-momentum frame. As demonstrated here for the first time, in
cases with pressure or concentration gradients, the enthalpy flux is
particularly insensitve to non-local effects, thus becoming increasingly
important in comparison to the reduced heat flux for high Knudsen numbers. The
RKM has been applied here to describe non-local heat flux, but since the
approach solves for the DF of each ion species, it may have future utility in
investigating other kinetic effects such as reactivity reduction in inertial
confinement fusion experiments [29, 47, 24]. Future work relating to the heat
flux may investigate additional effects such as magnetic fields and the
implementation of a RKM prescription to non-local heat flux in a simple 1D
hydrodynamics code.

##  Appendix A Kinetic equation for deviation of distribution function from
Maxwellian

The kinetic equation may be transformed to the centre-of-momentum frame by
transforming the velocity to the peculiar velocity
$\boldsymbol{w}=\boldsymbol{v}-\boldsymbol{u}(\boldsymbol{t},\boldsymbol{x})$,
leaving

$\displaystyle(\partial_{t}-(\partial_{t}\boldsymbol{u})\cdot\boldsymbol{\nabla_{w}}+(\boldsymbol{w}+\boldsymbol{u})\cdot(\boldsymbol{\nabla}-(\boldsymbol{\nabla
u})\cdot\boldsymbol{\nabla_{w}})+\boldsymbol{F}_{\alpha}\cdot\boldsymbol{\nabla_{w}})f_{\alpha}=\sum_{\beta}C_{\alpha\beta}\\{f_{\alpha},f_{\beta}\\}.$
(26)

The DFs are then decomposed via $f_{\alpha}=f^{M}_{\alpha}+\delta f_{\alpha}$,
giving

$\displaystyle(\partial_{t}-(\partial_{t}\boldsymbol{u})\cdot\boldsymbol{\nabla_{w}}+(\boldsymbol{w}+\boldsymbol{u})\cdot(\boldsymbol{\nabla}-(\boldsymbol{\nabla
u})\cdot\boldsymbol{\nabla_{w}})+\boldsymbol{F}_{\alpha}\cdot\boldsymbol{\nabla_{w}})(f^{M}_{\alpha}+\delta
f_{\alpha})=\sum_{\beta}\bigg{(}C^{T}_{\alpha\beta}\\{\delta
f_{\alpha}\\}+C^{F}_{\alpha\beta}\\{\delta f_{\beta}\\}\bigg{)},$ (27)

where the non-linear part of the collision operator has been neglected.

The Maxwellian term may be evaluated identically to Ferziger and Kaper [17]
evaluating $(\mathcal{D}f_{i})^{(0)}$, where the time derivatives acting on
bulk quantities are replaced with the leading-order hydrodynamic (Euler)
equations. Similarly, the
$-(\partial_{t}\boldsymbol{u}+\boldsymbol{u}\cdot\boldsymbol{\nabla
u})\cdot\boldsymbol{\nabla}_{\boldsymbol{w}}\delta f_{\alpha}$ term is
rewritten using the Euler momentum equation
$\partial_{t}\boldsymbol{u}+\boldsymbol{u}\cdot\boldsymbol{\nabla
u}=-\frac{1}{\rho}\boldsymbol{\nabla}p+\sum_{\beta}c_{\beta}\boldsymbol{F}_{\beta}$.
This transforms (27) to

$\displaystyle\biggl{(}\boldsymbol{w}\cdot\bigg{[}\frac{n}{n_{\alpha}}\boldsymbol{d}_{\alpha}+\bigg{(}\frac{m_{\alpha}w^{2}}{2T}-\frac{5}{2}\bigg{)}\boldsymbol{\nabla}\ln
T\bigg{]}+\frac{m_{\alpha}}{T}\bigg{(}\boldsymbol{ww}-\frac{1}{3}w^{2}\mathbb{I}\bigg{)}\colon\boldsymbol{\nabla
u}\biggl{)}f^{M}_{\alpha}$ (28)
$\displaystyle+\biggl{(}\partial_{t}+(\boldsymbol{w}+\boldsymbol{u})\cdot\boldsymbol{\nabla}+(-\boldsymbol{w}\cdot\boldsymbol{\nabla
u}+\frac{1}{\rho}\boldsymbol{\nabla}p+\boldsymbol{F}_{\alpha}-\sum_{\beta}c_{\beta}\boldsymbol{F}_{\beta})\cdot\boldsymbol{\nabla_{w}}\biggl{)}\delta
f_{\alpha}=\sum_{\beta}\bigg{(}C^{T}_{\alpha\beta}\\{\delta
f_{\alpha}\\}+C^{F}_{\alpha\beta}\\{\delta f_{\beta}\\}\bigg{)}.$

Writing 28 in 1D2V coordinates and neglecting the hydrodynamic velocity
assuming $u\ll v_{\text{th}}$ yields the kinetic equation for $\delta
f_{\alpha}$ as

$\displaystyle\partial_{t}\delta f_{\alpha}$
$\displaystyle=\sum_{\beta}\bigg{(}C^{T}_{\alpha\beta}\\{\delta
f_{\alpha}\\}+C^{F}_{\alpha\beta}\\{\delta f_{\beta}\\}\bigg{)}-\biggl{[}\mu
w\partial_{z}+\bigg{(}\frac{1}{\rho}\partial_{z}p+F_{\alpha}-\sum_{\beta}c_{\beta}F_{\beta}\bigg{)}\bigg{(}\mu\partial_{w}+\frac{1-\mu^{2}}{w}\partial_{\mu}\bigg{)}\biggl{]}\delta
f_{\alpha}$ (29) $\displaystyle-\mu
wf^{M}_{\alpha}\bigg{(}\frac{n}{n_{\alpha}}d_{\alpha}+\bigg{(}\frac{m_{\alpha}}{2T}w^{2}-\frac{5}{2}\bigg{)}\partial_{z}\ln
T\bigg{)}\,.$

##  Appendix B Field-particle operator

The Legendre components of the field-particle operator [36] are

$\displaystyle C^{F,(l)}_{\alpha\beta}\\{\delta f^{(l)}_{\beta}\\}$
$\displaystyle=\frac{2}{2l+1}\int_{-1}^{1}d\mu\
P_{l}(\mu)C^{F}_{\alpha\beta}\\{\delta f_{\beta}\\}$ (30)
$\displaystyle=4\pi\bigg{(}\frac{q_{\alpha}q_{\beta}}{m_{\alpha}}\bigg{)}^{2}\ln\Lambda_{\alpha\beta}\frac{m_{\alpha}}{T_{\alpha}}f^{M}_{\alpha}\Biggl{[}\frac{4\pi
T_{\alpha}}{m_{\beta}}\delta
f_{\beta}^{(l)}-\mathcal{H}^{(l)}_{\beta}+\bigg{(}\frac{m_{\alpha}}{m_{\beta}}-1\bigg{)}w\partial_{w}\mathcal{H}^{(l)}_{\beta}+\frac{m_{\alpha}w^{2}}{2T_{\alpha}}\partial_{w}^{2}\mathcal{G}^{(l)}_{\beta}\Biggl{]},$

where

$\displaystyle\mathcal{H}^{(l)}_{\beta}=\frac{4\pi}{2l+1}\biggl{(}\frac{1}{v^{l+1}}\int_{0}^{w}du\
u^{l+2}\delta f_{\beta}^{(l)}(u)+w^{l}\int_{w}^{\infty}du\
\frac{1}{u^{l-1}}\delta
f_{\beta}^{(l)}(u)\biggl{)}=2\pi\biggl{(}\Psi_{\beta,1}^{(l)}(w)+\Psi_{\beta,2}^{(l)}(w)\biggl{)}$
(31)

$\displaystyle
w\partial_{w}\mathcal{H}^{(l)}_{\beta}=\frac{4\pi}{2l+1}\biggl{(}-\frac{l+1}{w^{l+2}}\int_{0}^{w}du\
u^{l+2}\delta f_{\beta}^{(l)}(u)+lw^{l-1}\int_{w}^{\infty}du\
\frac{1}{u^{l-1}}\delta
f_{\beta}^{(l)}(u)\biggl{)}=2\pi\biggl{(}l\Psi_{\beta,1}^{(l)}(w)-(l+1)\Psi_{\beta,2}^{(l)}(w)\biggl{)}$
(32)

$\displaystyle w^{2}\partial_{w}^{2}\mathcal{G}^{(l)}_{\beta}$
$\displaystyle=\frac{4\pi}{2l+1}\biggl{(}\frac{(l+1)(l+2)}{2l+3}\frac{1}{w^{l+3}}\int_{0}^{w}du\
u^{l+4}\delta
f_{\beta}^{(l)}(u)-\frac{l(l-1)}{2l-1}\frac{1}{w^{l+1}}\int_{0}^{w}du\
u^{l+2}\delta f_{\beta}^{(l)}(u)$ (33)
$\displaystyle+\frac{(l+1)(l+2)}{2l+3}w^{l}\int_{w}^{\infty}du\
\frac{1}{u^{l-1}}\delta
f_{\beta}^{(l)}(u)-\frac{l(l-1)}{2l-1}w^{l-2}\int_{w}^{\infty}du\
\frac{1}{u^{l-3}}\delta f_{\beta}^{(l)}(u)\biggl{)}$
$\displaystyle=2\pi\biggl{(}\frac{(l+1)(l+2)}{2l+3}(\Psi_{\beta,4}^{(l)}(w)+w^{2}\Psi_{\beta,1}^{(l)}(w))-\frac{l(l-1)}{2l-1}(\Psi_{\beta,3}^{(l)}(w)+w^{2}\Psi_{\beta,2}^{(l)}(w))\biggl{)}$

where we have introduced the integrals

$\displaystyle\Psi_{\beta,1}^{(l)}(w)=\frac{2}{2l+1}w^{l}\int_{w}^{\infty}du\
u^{1-l}\delta
f_{\beta}^{(l)}(u),\quad\Psi_{\beta,2}^{(l)}(w)=\frac{2}{2l+1}\frac{1}{w^{l+1}}\int_{0}^{w}du\
u^{2+l}\delta f_{\beta}^{(l)}(u),$ (34)
$\displaystyle\Psi_{\beta,3}^{(l)}(w)=\frac{2}{2l+1}w^{l}\int_{w}^{\infty}du\
u^{3-l}\delta
f_{\beta}^{(l)}(u),\quad\Psi_{\beta,4}^{(l)}(w)=\frac{2}{2l+1}\frac{1}{w^{l+1}}\int_{0}^{w}du\
u^{4+l}\delta f_{\beta}^{(l)}(u).$

Substituting into the field-particle operator, we have

$\displaystyle C^{F,(l)}_{\alpha\beta}\\{\delta
f^{(l)}_{\beta}\\}=16\pi^{2}\bigg{(}\frac{q_{\alpha}q_{\beta}}{m_{\alpha}}\bigg{)}^{2}\ln\Lambda_{\alpha\beta}\frac{1}{v_{\text{th},\alpha}^{2}}f^{M}_{\alpha}\Biggl{[}\frac{m_{\alpha}}{m_{\beta}}v_{\text{th},\alpha}^{2}\delta
f_{\beta}^{(l)}+\bigg{(}-1+\bigg{(}\frac{m_{\alpha}}{m_{\beta}}-1\bigg{)}l+\frac{(l+1)(l+2)}{2l+3}\frac{w^{2}}{v_{\text{th},\alpha}^{2}}\bigg{)}\Psi_{\beta,1}^{(l)}(w)$
(35)
$\displaystyle+\bigg{(}-1-\bigg{(}\frac{m_{\alpha}}{m_{\beta}}-1\bigg{)}(l+1)-\frac{l(l-1)}{2l-1}\frac{w^{2}}{v_{\text{th},\alpha}^{2}}\bigg{)}\Psi_{\beta,2}^{(l)}(w)-\frac{l(l-1)}{2l-1}\frac{1}{v_{\text{th},\alpha}^{2}}\Psi_{\beta,3}^{(l)}(w)+\frac{(l+1)(l+2)}{2l+3}\frac{1}{v_{\text{th},\alpha}^{2}}\Psi_{\beta,4}^{(l)}(w)\Biggl{]}.$

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] Callen J D and Kissick M W 1997 Plasma Physics and Controlled Fusion 39 B173 URL https://dx.doi.org/10.1088/0741-3335/39/12B/014
  * [2] Craxton R S, Anderson K S, Boehly T R, Goncharov V N, Harding D R, Knauer J P, McCrory R L, McKenty P W, Meyerhofer D D, Myatt J F, Schmitt A J, Sethian J D, Short R W, Skupsky S, Theobald W, Kruer W L, Tanaka K, Betti R, Collins T J B, Delettrez J A, Hu S X, Marozas J A, Maximov A V, Michel D T, Radha P B, Regan S P, Sangster T C, Seka W, Solodov A A, Soures J M, Stoeckl C and Zuegel J D 2015 Physics of Plasmas 22 110501
  * [3] Meezan N B, Edwards M J, Hurricane O A, Patel P K, Callahan D A, Hsing W W, Town R P J, Albert F, Amendt P A, Hopkins L F B, Bradley D K, Casey D T, Clark D S, Dewald E L, Dittrich T R, Divol L, Döppner T, Field J E, Haan S W, Hall G N, Hammel B A, Hinkel D E, Ho D D, Hohenberger M, Izumi N, Jones O S, Khan S F, Kline J L, Kritcher A L, Landen O L, LePape S, Ma T, MacKinnon A J, MacPhee A G, Masse L, Milovich J L, Nikroo A, Pak A, Park H S, Peterson J L, Robey H F, Ross J S, Salmonson J D, Smalyuk V A, Spears B K, Stadermann M, Suter L J, Thomas C A, Tommasini R, Turnbull D P and Weber C R 2016 Plasma Physics and Controlled Fusion 59 014021 URL https://dx.doi.org/10.1088/0741-3335/59/1/014021
  * [4] Sanmartin J 1993 Electron Heat Conduction in ICF Targets (,) pp 185–196 ISBN 9781003068594
  * [5] Braginskii S I 1958 Soviet Phys. JETP URL https://www.osti.gov/biblio/4344773
  * [6] Rinderknecht H, Amendt P, Wilks S and Collins G 2018 Plasma Physics and Controlled Fusion 60
  * [7] Thomas A, Tzoufras M, Robinson A, Kingham R, Ridgers C, Sherlock M and Bell A 2012 Journal of Computational Physics 231 1051–1079 ISSN 0021-9991 special Issue: Computational Plasma Physics URL https://www.sciencedirect.com/science/article/pii/S0021999111005717
  * [8] Holec M, Nikl J and Weber S 2018 Physics of Plasmas 25 032704 ISSN 1070-664X (Preprint https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.5011818/16154743/032704_1_online.pdf) URL https://doi.org/10.1063/1.5011818
  * [9] Holec M, Nikl J, Vranic M and Weber S 2018 Plasma Physics and Controlled Fusion 60 044019 URL https://dx.doi.org/10.1088/1361-6587/aab05a
  * [10] Chen Z H, Yang X H, Zhang G B, Ma Y Y, Xu H, Luan S X and Zhang J 2023 Physics of Plasmas 30 ISSN 1089-7674 URL http://dx.doi.org/10.1063/5.0146342
  * [11] Brodrick J, Kingham R, Marinak M, Patel M, Chankin A, Omotani J, Umansky M, Del Sorbo D, Dudson B, Parker J, Kerbel G, Sherlock M and Ridgers C 2017 Physics of Plasmas 24 092309
  * [12] Meezan N B, Woods D T, Izumi N, Chen H, Scott H A, Schneider M B, Liedahl D A, Jones O S, Zimmerman G B, Moody J D, Landen O L and Hsing W W 2020 Physics of Plasmas 27 102704 ISSN 1070-664X (Preprint https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/5.0018733/15927783/102704_1_online.pdf) URL https://doi.org/10.1063/5.0018733
  * [13] Chapman D A, Pecover J D, Chaturvedi N, Niasse N, Read M P, Vassilev D H, Chittenden J P, Hawker N and Joiner N 2021 Physics of Plasmas 28 072702
  * [14] Atzeni S, Schiavi A and Bellei C 2007 Physics of Plasmas 14 052702
  * [15] Braginskii S I 1965 Consultants Bureau, New York, 1965, p.205
  * [16] Chapman S, Cowling T G and Park D 1962 American Journal of Physics 30 389–389 ISSN 0002-9505 (Preprint https://pubs.aip.org/aapt/ajp/article-pdf/30/5/389/11912407/389_2_online.pdf) URL https://doi.org/10.1119/1.1942035
  * [17] Ferziger J H and Kaper H G 1972 Mathematical Theory of Transport Processes in Gases (North-Holland)
  * [18] Zhdanov V 2002 Transport Processes in Multicomponent Plasma vol 44 (Taylor and Francis Group)
  * [19] Kagan G and Baalrud S D 2016 arXiv preprint arXiv:1611.09872 URL https://arxiv.org/abs/1611.09872
  * [20] Kagan G and Tang X Z 2012 Physics of Plasmas 19 082709 ISSN 1070-664X (Preprint https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4745869/13888529/082709_1_online.pdf) URL https://doi.org/10.1063/1.4745869
  * [21] Kagan G and Tang X Z 2014 Physics Letters A 378 1531–1535 URL https://doi.org/10.1016%2Fj.physleta.2014.04.005
  * [22] Epperlein E M and Short R W 1991 Physics of Fluids B: Plasma Physics 3 3092–3098 ISSN 0899-8221 (Preprint https://pubs.aip.org/aip/pfb/article-pdf/3/11/3092/12317455/3092_1_online.pdf) URL https://doi.org/10.1063/1.859789
  * [23] Morse R L and Nielson C W 1973 The Physics of Fluids 16 909–920 ISSN 0031-9171 (Preprint https://pubs.aip.org/aip/pfl/article-pdf/16/6/909/12387191/909_1_online.pdf) URL https://doi.org/10.1063/1.1694445
  * [24] Hoffman N M, Zimmerman G B, Molvig K, Rinderknecht H G, Rosenberg M J, Albright B J, Simakov A N, Sio H, Zylstra A B, Gatu Johnson M, Séguin F H, Frenje J A, Li C K, Petrasso R D, Higdon D M, Srinivasan G, Glebov V Y, Stoeckl C, Seka W and Sangster T C 2015 Physics of Plasmas 22 052707 ISSN 1070-664X (Preprint https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4921130/13780497/052707_1_online.pdf) URL https://doi.org/10.1063/1.4921130
  * [25] Schurtz G P, Nicolaï P D and Busquet M 2000 Physics of Plasmas 7 4238–4249 ISSN 1070-664X (Preprint https://pubs.aip.org/aip/pop/article-pdf/7/10/4238/12283913/4238_1_online.pdf) URL https://doi.org/10.1063/1.1289512
  * [26] Derentowicz H, Kaliski S and Ziółkowski Z 1977 Journal of Technical Physics 18 465
  * [27] Derentowicz H, Kaliski S, Wolski J and Ziółkowski Z 1977 Journal of Technical Physics 25 135
  * [28] Tang X Z, Berk H L, Guo Z and McDevitt C J 2014 Physics of Plasmas 21 ISSN 1070-664X 032707 (Preprint https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4868731/15880309/032707_1_online.pdf) URL https://doi.org/10.1063/1.4868731
  * [29] McDevitt C J, Tang X Z, Guo Z and Berk H L 2014 Physics of Plasmas 21 ISSN 1070-664X 032708 (Preprint https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4868732/15880219/032708_1_online.pdf) URL https://doi.org/10.1063/1.4868732
  * [30] Tang X Z, McDevitt C J, Guo Z and Berk H L 2014 Physics of Plasmas 21 ISSN 1070-664X 032706 (Preprint https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4868733/15877954/032706_1_online.pdf) URL https://doi.org/10.1063/1.4868733
  * [31] Landau L D 1936 Physicalische Zeitschrift der Sowjetunion 10 154
  * [32] Helander P 2002 Collisional transport in magnetized plasmas (Cambridge: CUP) ISBN 0521807980
  * [33] Abramowitz M and Stegun I A 1964 Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables ninth dover printing, tenth gpo printing ed (New York: Dover)
  * [34] Spitzer L 1962 Physics of Fully Ionized Gases (New York: Interscience (2nd edition), 1962)
  * [35] Rosenbluth M N, MacDonald W M and Judd D L 1957 Phys. Rev. 107(1) 1–6 URL https://link.aps.org/doi/10.1103/PhysRev.107.1
  * [36] Nishimura S 2015 Physics of Plasmas 22 ISSN 1070-664X 122503 (Preprint https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4936799/13970811/122503_1_online.pdf) URL https://doi.org/10.1063/1.4936799
  * [37] Bell A R, Evans R G and Nicholas D J 1981 Phys. Rev. Lett. 46(4) 243–246 URL https://link.aps.org/doi/10.1103/PhysRevLett.46.243
  * [38] Chang J and Cooper G 1970 Journal of Computational Physics 6 1–16 ISSN 0021-9991 URL https://www.sciencedirect.com/science/article/pii/002199917090001X
  * [39] Virtanen P, Gommers R, Oliphant T E, Haberland M, Reddy T, Cournapeau D, Burovski E, Peterson P, Weckesser W, Bright J, van der Walt S J, Brett M, Wilson J, Millman K J, Mayorov N, Nelson A R J, Jones E, Kern R, Larson E, Carey C J, Polat İ, Feng Y, Moore E W, VanderPlas J, Laxalde D, Perktold J, Cimrman R, Henriksen I, Quintero E A, Harris C R, Archibald A M, Ribeiro A H, Pedregosa F, van Mulbregt P and SciPy 10 Contributors 2020 Nature Methods 17 261–272
  * [40] Byrne G D and Hindmarsh A C 1975 ACM Trans. Math. Softw. 1 71–96 ISSN 0098-3500 URL https://doi.org/10.1145/355626.355636
  * [41] Casey D T, Frenje J A, Gatu Johnson M, Manuel M J E, Rinderknecht H G, Sinenian N, Séguin F H, Li C K, Petrasso R D, Radha P B, Delettrez J A, Glebov V Y, Meyerhofer D D, Sangster T C, McNabb D P, Amendt P A, Boyd R N, Rygg J R, Herrmann H W, Kim Y H and Bacher A D 2012 Phys. Rev. Lett. 108(7) 075002 URL https://link.aps.org/doi/10.1103/PhysRevLett.108.075002
  * [42] Rinderknecht H G, Sio H, Li C K, Zylstra A B, Rosenberg M J, Amendt P, Delettrez J, Bellei C, Frenje J A, Gatu Johnson M, Séguin F H, Petrasso R D, Betti R, Glebov V Y, Meyerhofer D D, Sangster T C, Stoeckl C, Landen O, Smalyuk V A, Wilks S, Greenwood A and Nikroo A 2014 Phys. Rev. Lett. 112(13) 135001 URL https://link.aps.org/doi/10.1103/PhysRevLett.112.135001
  * [43] Rosenberg M J, Rinderknecht H G, Hoffman N M, Amendt P A, Atzeni S, Zylstra A B, Li C K, Séguin F H, Sio H, Johnson M G, Frenje J A, Petrasso R D, Glebov V Y, Stoeckl C, Seka W, Marshall F J, Delettrez J A, Sangster T C, Betti R, Goncharov V N, Meyerhofer D D, Skupsky S, Bellei C, Pino J, Wilks S C, Kagan G, Molvig K and Nikroo A 2014 Phys. Rev. Lett. 112(18) 185001 URL https://link.aps.org/doi/10.1103/PhysRevLett.112.185001
  * [44] Amendt P, Landen O L, Robey H F, Li C K and Petrasso R D 2010 Phys. Rev. Lett. 105(11) 115005 URL https://link.aps.org/doi/10.1103/PhysRevLett.105.115005
  * [45] Landau L D and Lifshitz E M 1987 Fluid Mechanics, Second Edition: Volume 6 (Course of Theoretical Physics) 2nd ed Course of theoretical physics / by L. D. Landau and E. M. Lifshitz, Vol. 6 (Butterworth-Heinemann) ISBN 0750627670 URL http://www.worldcat.org/isbn/0750627670
  * [46] Cao D, Moses G and Delettrez J 2015 Physics of Plasmas 22 082308 ISSN 1070-664X (Preprint https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4928445/15965787/082308_1_online.pdf) URL https://doi.org/10.1063/1.4928445
  * [47] Kagan G, Svyatskiy D, Rinderknecht H G, Rosenberg M J, Zylstra A B, Huang C K and McDevitt C J 2015 Phys. Rev. Lett. 115(10) 105002 URL https://link.aps.org/doi/10.1103/PhysRevLett.115.105002