Thermal transport plays a critical role in inertial confinement fusion (ICF), leading to fusion target compression and heating. Classically, heat transport in a plasma is described by Spitzer–Härm theory.1 Assuming that the plasma is sufficiently close to thermal equilibrium, the transport coefficients can be derived by linearization of the Vlasov–Fokker–Planck (VFP) equation. However, Spitzer–Härm theory often becomes invalid in the presence of steep temperature gradients in laser-produced plasmas, where nonlocal transport models are required.2–12 Nonlocal thermal transport is also widely encountered in other plasma systems, such as magnetic fusion devices, astrophysical environments, and general laser–plasma interactions.13–24

On the other hand, megagauss magnetic fields can be self-generated in laser-produced plasmas via various mechanisms such as the Biermann battery effect when the temperature and density gradients are noncollinear,25 the Weibel instability due to anisotropic electron velocity distributions,26 and the transport of hot electrons.27 Also, it has been proposed that the application of strong external magnetic fields could improve energy coupling efficiency in ICF, and there have been some experimental implementations of this proposal.28–30 In a magnetized plasma, the transport coefficients are expressed in tensor form and are classically given by the Braginskii theory.31 On the basis of this theory, some more accurate and practical models for arbitrary atomic numbers have been proposed.32,33 Generally, magnetic fields tend to inhibit and divert the heat flux, and nonlocal thermal transport also takes place in magnetized plasmas with steep temperature gradients. In previous studies, the magnetic fields have usually been assumed to be applied in the direction perpendicular to the temperature gradient.10,11,34–37 In real experiments, however, magnetic fields can be found in arbitrary directions.

In this paper, nonlocal electron thermal transport with DC magnetic fields applied in arbitrary directions is investigated theoretically and numerically using a VFP simulation code developed with a spherical harmonics expansion method. In contrast to local heat transport, it is found that the presence of a magnetic field component along the temperature gradient can lead to significant coupling of thermal transport along the directions parallel (longitudinal) and perpendicular (transverse) to the temperature gradient in nonlocal heat transport. The remainder of the paper is structured as follows. In Sec. II, the local theory of thermal transport in a magnetized plasma is reviewed briefly. Section III introduces the basic equations and numerical scheme for the VFP code that is developed here and used in the subsequent investigation. Thermal transport in a magnetized plasma under different configurations is investigated in Sec. IV. Conclusions are presented in Sec. V.

Thermal transport in a fully ionized plasma can be described by the VFP equation$$\frac{\partial f}{\partial t}+\mathbf{v}\cdot \nabla f-\left(\frac{e\mathbf{E}}{m_e}+\mathbf{v}\times \frac{e\mathbf{B}}{m_{e}c}\right)\cdot {\nabla }_{v}f=C_{ee}+C_{ei},$$(1)where f(v, x, t) is the electron distribution function (EDF), m_e and e are the electron mass and charge, and E and B are the electric and magnetic fields. The operators on the right-hand side represent the effects of electron–electron and electron–ion collisions. Under the assumption that the temperature gradient is sufficiently small, with ∇T_e/T_e ≪ 1/λ_e, where λ_e is the electron mean free path, the distribution function can be expanded as a homogeneous zeroth-order term plus a small anisotropic term:31$$f(\mathbf{v},x,t)\simeq f_0(v,x,t)+\frac{\mathbf{v}}{v}\cdot {\mathit{f}}_1(v,x,t),$$(2)where the isotropic term f₀(v, x, t) is initially assumed to be a Maxwellian distribution given by $f_0=n_e{(m_e/2\pi T_e)}^{3/2}\exp (-m_e{v}^2/2T_e)$ and the perturbed term f₁(v, x, t) is related to the thermal transport due to the temperature gradient. Substituting Eq. (2) into Eq. (1), the perturbed term for the steady state can be obtained as$$v\nabla f_0-\frac{e\mathbf{E}}{m_e}{\nabla }_v{f}_0-\frac{e\mathbf{B}}{m_{e}c}\times {\mathit{f}}_1=-{\nu }_{ei}{\mathit{f}}_1,$$(3)where ν_ei is the electron–ion collision frequency. The electron–electron collision term that would otherwise appear on the right-side of Eq. (3) can be neglected in the Lorentz gas approximation in the high-Z limit, which is valid since the characteristic electron–electron collision frequency ν_ee ≈ ν_ei/Z in neutral plasmas.31

In classical local transport models, it is usually assumed that the magnetic field B is perpendicular to the temperature gradient, because the magnetic field does not change the transport along the magnetic field for the first order of perturbations.32 Under a transverse magnetic field, the solution to Eq. (3) is given straightforwardly by$$\begin{array}{l}\hfill {\mathit{f}}_1=\frac{1}{{\nu }_{ei}^2+{\omega }^2}(-{\nu }_{ei}\mathbf{g}-\mathit{\omega }\times \mathbf{g}),\end{array}$$(4)$$\begin{array}{l}\hfill \mathbf{g}=v\nabla f_0-\frac{e\mathbf{E}}{m_e}{\nabla }_v{f}_0,\end{array}$$(5)where ω = eB/m_ec is the gyrofrequency of electrons in the magnetic field. Equation (4) shows that there are two terms contributing to the perturbation of the distribution function, the first of which is due to the temperature gradient and the second to the presence of the magnetic field.34 It should be noted that f₁ is rotated by the Lorentz force, and so the heat flux is no longer parallel to the temperature gradient. In addition, the “effective” electron collision frequency is reduced by the magnetic field via the factor $1/(1+{\omega }^2/{\nu }_{ei}^2)$.36

To calculate the electric field and heat flux, Eq. (4) is substituted into the Ampère–Maxwell law:$$\mathbf{J}=-\frac{4\pi e}{3}\int {\mathit{f}}_1{v}^{3}dv=\frac{c}{4\pi }\nabla \times \mathbf{B}.$$(6)The heat flux is then given by$$\begin{array}{ll}\hfill \mathbf{Q}& =\frac{2\pi m_e}{3}\int {\mathit{f}}_1{v}^{5}dv=\frac{2\pi m_e}{3}\int \frac{v^5}{{\nu }_{ei}^2+{\omega }^2}(-{\nu }_{ei}\mathbf{g}-\mathit{\omega }\times \mathbf{g})dv.\hfill \end{array}$$(7)In the notation adopted by Braginskii,31 the above integrations yield the following transport relations:$$\begin{array}{ll}\hfill e{n}_e\mathbf{E}=& -\nabla P_e+\mathbf{J}\times \mathbf{B}/c+{\alpha }_{\perp }\mathbf{J}/n_e-{\alpha }_{\wedge }\mathbf{b}\times \mathbf{J}/n_e\hfill \\ \hfill & -{\beta }_{\perp }\nabla T_e-{\beta }_{\wedge }\mathbf{b}\times \nabla T_e,\hfill \end{array}$$(8)$$\begin{array}{l}\hfill \mathbf{Q}=-\left({\beta }_{\perp }+\frac{5}{2}\right)\frac{T_e}{e}\mathbf{J}-{\beta }_{\wedge }\frac{T_e}{e}\mathbf{b}\times \mathbf{J}-{\kappa }_{\perp }\nabla T_e-{\kappa }_{\wedge }\mathbf{b}\times \nabla T_e,\end{array}$$(9)where P_e is the scalar intrinsic pressure, n_e is the electron number density, and b is the unit vector in the magnetic field direction. The κ coefficients are for the thermal heat flux from temperature gradients, whereas the β coefficients give the heat flux driven by electric currents. These quantities are expressed in terms of components parallel to the driving terms, which are perpendicular to the magnetic field and indicated by subscript ⊥, and components in a crossed direction, which are perpendicular to both the magnetic field and the driving term and indicated by subscript ∧. The driving terms are the current density and temperature gradient. The transport coefficients read as follows:35$$\begin{array}{cc}\hfill & \hfill {\kappa }_{\perp }=\frac{2T_e{n}_e}{m_e{v}_T}\left|\mathrm{\Omega }\right|\frac{{\varphi }_{11}{\varphi }_7^2-{\varphi }_7{\varphi }_9^2+{\mathrm{\Omega }}^2({\varphi }_{11}{\varphi }_{10}^2+{\varphi }_7{\varphi }_{12}^2-2{\varphi }_9{\varphi }_{10}{\varphi }_{12})}{{\varphi }_7^2+{\mathrm{\Omega }}^2{\varphi }_{10}^2},\hfill \\ \hfill & \hfill {\kappa }_{\wedge }=\frac{2T_e{n}_e}{m_e{v}_T}\frac{{\varphi }_{14}{\varphi }_7^2+{\varphi }_{10}{\varphi }_9^2+2{\varphi }_7{\varphi }_9{\varphi }_{12}+{\mathrm{\Omega }}^2({\varphi }_{14}{\varphi }_{10}^2-{\varphi }_{10}{\varphi }_{12}^2)}{{\varphi }_7^2+{\mathrm{\Omega }}^2{\varphi }_{10}^2},\hfill \\ \hfill & \hfill {\beta }_{\perp }=\left|\mathrm{\Omega }\right|\frac{{\varphi }_9{\varphi }_7+{\mathrm{\Omega }}^2{\varphi }_{12}{\varphi }_{10}}{{\varphi }_7^2+{\mathrm{\Omega }}^2{\varphi }_{10}^2}-\frac{5}{2},{\beta }_{\wedge }=\frac{{\varphi }_{12}{\varphi }_7-{\varphi }_9{\varphi }_{10}}{{\varphi }_7^2+{\mathrm{\Omega }}^2{\varphi }_{10}^2},\hfill \\ \hfill & \hfill {\alpha }_{\perp }=\frac{1}{2}{m}_e{n}_e{v}_T\frac{{\varphi }_7}{{\varphi }_7^2+{\mathrm{\Omega }}^2{\varphi }_{10}^2},\hfill \\ \hfill & \hfill {\alpha }_{\wedge }=\frac{1}{2}{m}_e{n}_e{v}_T\left|\mathrm{\Omega }\right|\left(1-\frac{{\varphi }_{10}}{{\varphi }_7^2+{\mathrm{\Omega }}^2{\varphi }_{10}^2}\right),\hfill \end{array}$$(10)where Ω = ωτ_ei is the Hall coefficient, τ_ei is the electron–ion collision time, u = v/v_te is the normalized velocity with $v_{te}=\sqrt{2T_e/m_e}$, and the dimensionless functions ϕ_n are$${\varphi }_n(\mathrm{\Omega })=\frac{4}{3\sqrt{\pi }}\int \frac{u^n{e}^{-u^2}du}{1+{\mathrm{\Omega }}^2{u}^6}.$$(11)

Figure 1 shows the transport coefficients as functions of the applied magnetic field strength ωτ_ei, which indicates both the suppression and deflection of the heat flux by the magnetic field. As the magnetic field increases, the transport coefficients κ_⊥ and β_⊥ along the temperature gradient decay rapidly, while the crossed transport coefficients κ_∧ and β_∧ perpendicular to the gradient and the magnetic field first increase before ωτ_ei ∼ 0.1 and then decrease.

Since the diffusion approximation is used, the above theory is limited to local thermal transport. To study nonlocal thermal transport under magnetic fields along arbitrary directions, we have developed a VFP code with a longitudinal spatial coordinate and three velocity components (1D3V) for electrons, while the ions are fixed as a neutral background. Generally, the VFP equation can be written as$$\frac{\partial f}{\partial t}+v_x\frac{\partial f}{\partial x}-\left(\frac{e\mathbf{E}}{m_e}+\mathbf{v}\times \frac{e\mathbf{B}}{m_{e}c}\right)\cdot {\nabla }_{v}f=C_{ee}+C_{ei}.$$(12)

Even though the VFP equation gives a complete description of electron dynamics, direct solution of this equation is difficult. An efficient way is to reduce the computational cost by expanding the distribution function in spherical harmonics $Y_{\ell }^m$, with the expansion being truncated at a certain order ℓ:38–40$$f(\mathbf{v},x,t)=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{f}_{\ell }^m(v,x,t)Y_{\ell }^m.$$(13)

This approach is acceptable because the collisions ensure that the distribution function tends to isotropy. On substitution of the above expansion into Eq. (12), the equation for the harmonic components can be written as$$\frac{\partial f_{\ell }^m}{\partial t}{+\mathcal{A}}_{\ell }^m{+\mathcal{E}}_{\ell }^m{+\mathcal{B}}_{\ell }^m{=\mathcal{C}}_{\ell }^m.$$(14)On the left-hand side of Eq. (14), ${\mathcal{A}}_{\ell }^m$ represents the contribution of spatial advection, ${\mathcal{E}}_{\ell }^m$ that of electric fields, and ${\mathcal{B}}_{\ell }^m$ that of magnetic fields, while the term on the right-hand side is due to collisions (see Ref. 39 for more details of the terms on both sides). The electric field is calculated via the quasineutral approximation or by Ampère’s law, and we ignore the evolution of the magnetic fields. The timescale of electron transport is much shorter than that of ion motion, and so the ions are treated as a cold background.

The numerical implementations for solving the equations are identical to the methods in Ref. 40. The central difference approach is used for the derivatives in real space, and the Vlasov operator on the left-hand side is integrated by the Runge–Kutta method. Following the Chang–Cooper scheme,41 the Fokker–Planck collision operator is differenced using an implicit and energy-conserving scheme according to Ref. 42.

To simplify the study, we normalize the variables as follows:$$\begin{array}{cccccc}\hfill v& \to \frac{v}{v_n},\hfill & \hfill t& \to \frac{t}{{\tau }_n},\hfill & \hfill x& \to \frac{x}{{\lambda }_n},\hfill \\ \hfill f& \to \frac{v_n^{3}f}{n_n},\hfill & \hfill \mathbf{E}& \to \frac{e{\tau }_n^2\mathbf{E}}{m_e{\lambda }_n},\hfill & \hfill \mathbf{B}& \to \frac{e\mathbf{B}{\tau }_n}{m_e},\hfill \end{array}$$(15)where the temperature is normalized to a reference value T_n and the density to n_n. Correspondingly, $v_n=\sqrt{2k_B{T}_n/m_e}$ is the thermal velocity, ${\tau }_n=m_e^2{v}_n^3/(4\pi Z{e}^4{n}_n\ln \mathrm{\Lambda })$ is the electron–ion collision time, and λ_n = v_nτ_n. For simplicity, the Coulomb logarithm ln Λ is assumed to be constant everywhere in spatial space.

To validate the code, we have carried out benchmark studies by comparing the simulation results with those of the classical local transport theory given in Sec. II. The initial temperature profile is$$T_e(x)=\left\{\begin{array}{cc}{T}_0,\hfill & x>2x_{\max }/3,\hfill \\ T_0+2\mathrm{\Delta }T,\hfill & x<x_{\max }/3,\hfill \\ T_0+\mathrm{\Delta }T-\mathrm{\Delta }T\cos \left(\frac{3x}{x_{\max }}\pi \right),\hfill & \mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}.\hfill \end{array}\right.$$(16)

By varying the values of x_max and ΔT, different temperature gradients can be obtained. The magnetic field is uniform, and so J = 0 and the heat flux driven by the current is zero. The plasma is initialized under the assumptions that the density is spatially uniform, the electron distribution is Maxwellian, and Z = 16. A continuous boundary condition is applied on two boundaries along the x direction. The simulation is stopped when the heat wave front reaches the boundary.

In the case of an unmagnetized plasma, according to the local transport theory, the electrons with velocities near v ≃ 3.7v_te make the greatest contribution to the heat flux. Correspondingly, f₁ ≈ 533(λ_e/L_T)f_M at v ≃ 3.7v_te, where L_T = T_e/∇T_e is the scale length. Therefore, the classical electron thermal transport theory only holds when λ_e/L_T ≪ 10⁻³. We have carried out a numerical simulation in the local transport regime. Figure 2 presents the heat transport along the temperature gradient with ΔT = 1.5T₀, where the degree of nonlocality at x = x_max/2 is about λ_e/L_T ∼ 5 × 10⁻⁵, and the magnetic field B = B_xe_x + B_ze_z, with eB_xτ_n/m_e = ω_xτ_n = 0.1 and eB_zτ_n/m_e = ω_zτ_n = 0.1. It is confirmed that the simulation results obtained with our code agree well with those of the local transport theory. It is worth mentioning that the Braginskii heat fluxes are calculated using the instantaneous density and temperature profiles obtained from the VFP simulations.

Using the VFP code, we investigate the coupling of thermal transports along the directions parallel and perpendicular to the temperature gradient in the presence of a DC magnetic field along different directions.

First, we consider a magnetic field applied along the transverse (z) direction, which will lead to thermal transport in the y direction besides the transport along the temperature gradient in the x direction. When the condition λ_e/L_T ≪ 10⁻³ is not satisfied, the electron distribution function (EDF) will no longer be Maxwellian. Furthermore, the perturbation f₁ may be greater than f₀ in some velocity region, which will cause breakdown of the local theory of thermal transport. As a result, the heat flux calculated by the classical theoretical models will be significantly overestimated in the case of a large temperature gradient. Figure 3 shows the heat flux along the directions parallel and perpendicular to a large temperature gradient. At a moderately large temperature gradient λ_e/L_T ≈ 0.01, Fig. 3(a) clearly shows that nonlocal transport appears in an unmagnetized plasma. Our VFP simulation demonstrates that the peak heat flux Q_x with B_z = 0 is much smaller than the saturated heat flux predicted by Spitzer–Härm theory. Moreover, the heat flux obtained from our VFP simulation is distributed over a wider area, which implies the existence of a preheating effect due to the nonlocal thermal transport.

When a transverse magnetic field is applied along the z direction with a normalized field strength ω_zτ_n = 0.05, which is about 16 T for a plasma with Z = 16, n_e = 10²¹ cm⁻³, and T_e = 2.5 keV. It is found that the heat flux Q_x along the temperature gradient can be significantly reduced in comparison with the unmagnetized case. In the case ω_zτ_n = 0.05, the heat flux Q_x obtained from the VFP simulation almost coincides with that estimated by the classical Braginskii model, indicating that thermal transport with a large temperature gradient may become local again under a strong transverse magnetic field. However, Fig. 3(b) shows that a transverse heat flux Q_y along the y direction will be present in the case of a transverse magnetic field B_z. As illustrated in Fig. 3(b), this transverse heat flux Q_y is well predicted by the Braginskii theory for a moderately large temperature gradient λ_e/L_T ≈ 0.01. This is because the magnetic field will be able to localize the heat flux if the Larmor radius is much shorter than the scale length of the temperature gradient.36 For the magnetic field ωτ_n = 0.05 and temperature gradient λ_e/L_T ≃ 0.01 used in Fig. 3(b), we have r_L/L_T ≈ 0.2.

For comparison, Figs. 3(c) and 3(d) show the heat flux with a large temperature gradient λ_e/L_T ≈ 0.05. In this case, nonlocal transport becomes significant for both the unmagnetized and magnetized plasmas. More importantly, the nonlocal thermal transport in this case occurs not only along the temperature gradient (Q_x) but also along the transverse direction (Q_y) with B_z ≠ 0. The thermal transport along the y direction exhibits a stronger nonlocality than that along the x direction. Usually, the so-called heat flux limiter $\left|Q_{\text{VFP}}/m_e{n}_e{v}_{te}^3\right|$ can be used to estimate the degree of nonlinearity of the transport, and a flux limiter of 0.15 is often used in MHD fluid codes.43,44 The limiter is around 0.7 for the longitudinal heat flux in Fig. 3(c), while for the Righi–Leduc heat flux in Fig. 3(d) it is about 0.18. That is to say, the VFP simulations show that the heat flux limiter for the transverse (Righi–Leduc) heat flux is much smaller than that for the standard component of the longitudinal heat flux.

We now consider the case in which the DC magnetic field has nonzero components along both the transverse (z) and longitudinal (x) directions: B = B_xe_x + B_ze_z. The thermal flux Q_z along the z direction is now no longer zero, since the magnetic field B_x along the x direction can divert the heat flux Q_y along the y direction to the z direction. According to the local heat transport theory, the transport parallel to the magnetic field will not be affected by the applied magnetic field. Taking the scalar product of both sides of Eq. (3) with B gives B · [v∇f₀ − (eE/m_e)∇_vf₀ + ν_eif₁] = 0. Therefore, without loss of generality, it can be assumed that the magnetic field is perpendicular to the heat flux, i.e., f₁ · ω = 0. In the nonlocal transport regime, however, the isotropic part of the electron distribution can be far from Maxwellian, and its evolution is associated with all the components of f₁. Furthermore, the first-order expansion is no longer sufficiently accurate in the nonlocal transport regime. Actually, high-order spherical harmonic terms play an important role in the transport process, and the first-order harmonic f₁ is strongly coupled with those high-order terms, as shown in the following.

Figure 4 compares the heat fluxes along different directions with a large temperature gradient λ_e/L_T ≈ 0.05 under a magnetic field B = B_xe_x + B_ze_z. We truncate the spherical harmonic expansion at order ℓ = 8 here. Figure 4(a) shows that in comparison with the case of a purely transverse magnetic field B = B_ze_z, the heat flux Q_x along the temperature gradient can be slightly enhanced by the x component of the magnetic field B = B_xe_x + B_ze_z, and the Braginskii theory significantly overestimates Q_x. By contrast, Fig. 4(b) shows that the heat flux Q_y along the cross direction will be significantly reduced by the B_x component of the magnetic field. This is because B_x will deflect the heat flux Q_y partially into the heat flux Q_z. As another consequence, a nonzero Q_z is induced, as shown in Fig. 4(c). As can be seen from Figs. 4(b) and 4(c), the Braginskii theory significantly overestimates Q_y as well as Q_z.

In Fig. 5, we plot the time-averaged ratio $\left|Q_{\text{VFP}}/Q_{\text{Brag}}\right|$ of the heat flux obtained from the VFP simulation to that estimated by the Braginskii theory under different magnetic fields at x_max/2. When the magnetic field is perpendicular to the temperature gradient (ω_xτ_n = 0), the transverse magnetic field B_z inhibits the nonlocal effect, and the heat flux in the x direction is almost local when ω_zτ_n ≥ 0.1. The heat flux Q_y in the y direction from the VFP simulation is significantly smaller than that estimated by the Braginskii theory, and the radio $\left|Q_{\text{VFP}}/Q_{\text{Brag}}\right|$ is always smaller than 0.4 for ω_zτ_n ≤ 0.125. This indicates that the transverse magnetic field B_z can localize the heat flux Q_x along the temperature gradient with ω_xτ_n ≥ 0.1 (the corresponding r_L/L_T ≤ 0.05). However, this transverse magnetic field will induce a nonzero heat flux Q_y in the y direction, and the classical Braginskii theory will overestimate the value of Q_y. The heat fluxes Q_y and Q_z in the two transverse directions can be localized by an increase in the magnetic field component B_x along the temperature gradient, as shown in Figs. 5(b) and 5(c), since electron transverse motion will be suppressed by the longitudinal magnetic field B_x, thereby reducing the effective electron mean free path in the transverse directions.34 However, Fig. 5(a) shows that the Braginskii theory overestimates the heat flux Q_x more with increasing B_x, and the heat transport along the temperature gradient exhibits a stronger nonlocality. With the increase of the longitudinal magnetic field B_x, the electrons will be guided along the x direction and the transverse electron motion will be suppressed. Consequently, the heat flux Q_x is enhanced, while the heat fluxes Q_y and Q_z in the transverse directions are localized. Therefore, with increasing B_x, the heat transport along the x direction exhibits a stronger nonlocality with a smaller Q_x,VFP/Q_x,Brag, as shown in Fig. 5(a). In the NIF, the typical L_T is about 5 mm,45,46 and so magnetic fields of the order of several tesla are high enough to affect the transport process. The experiment also shows that an external field along the hohlraum results in a significant increase in plasma temperature.47

As shown in Fig. 4(c), a nonzero Q_z will be induced owing to the rotation of Q_y by the B_x component of the magnetic field. Figure 5(c) further shows that the Q_z obtained from the VFP simulations is always much lower than that estimated by the classical Braginskii theory $(\left|Q_{z,\mathrm{V}\mathrm{F}\mathrm{P}}/Q_{z,\mathrm{B}\mathrm{r}\mathrm{a}\mathrm{g}}\right|\le 0.1)$. That is to say, the nonlocality of the heat flux component Q_z is more significant, and hence Q_z should be treated more carefully.

In our simulations, the initial electron distribution function is assumed to be Maxwellian, and therefore a certain response time is required to generate the heat flux. In the local transport theory, the Braginskii transport is a quasistatic state. When the temperature gradient is steep, however, we find that the heat flux cannot reach a quasistatic value before the temperature profile changes significantly. For an initial steep temperature gradient λ_e/L_T ≈ 0.05, Fig. 6 shows that the growth rates of the heat flux components along different directions are significantly different. The heat flux Q_x in the x direction increases most rapidly at the beginning, reaching a peak at t ∼ 5τ_n. It then decays rapidly during 5τ_n ≤ t ≤ 20τ_n and slowly after t ∼ 30τ_n. In comparison, the growth rates of the heat flux components in the y and z directions are slower. The heat flux Q_x is due mainly to the temperature gradient along the x direction, while the heat flux components Q_y and Q_z in the directions perpendicular to the temperature gradient come from deflection by the magnetic field, which naturally lags behind the growth of Q_x. The result is that the nonlocal heat flux depends not only on the current temperature profile, but also on the history of the temperature profile. This could be one reason why some nonlocal heat flux models fail, since they only consider the T_e profile at the current time.

In the local transport theory, the diffusion approximation is usually adopted, with only the first-order spherical harmonic being considered. However, the situation is much more complicated in the nonlocal transport regime, particularly in magnetized plasmas, where transverse thermal transport is induced as well. Figure 7 shows the ratios of the higher-order spherical harmonic terms $f_1^0$, $f_2^0$, $f_3^0$, and $f_4^0$ to the zeroth-order (isotropic) term $f_0^0$ of the EDF in velocity space at the points x = x_max/3, x_max/2, and 2x_max/3 in the nonlocal regime. For a relatively large temperature gradient, it can be seen that the first-order term $f_1^0$ can be greater than the isotropic term $f_0^0$ in the high-energy region around v ≈ 3.7v_n that makes the greatest contribution to the heat flux. Also, the high-order terms appear to be more important in the spatial region around x = 2x_max/3 where the relatively cold plasma is preheated by the electrons at the tail of the electron distribution function. By contrast, the high-order terms are not so important at the top of the temperature profile (x ≤ x_max/3). In particular, the high-order spherical harmonic terms can become comparable to or even larger than the zeroth-order term $f_0^0$ in the high-energy velocity region. Therefore, linear perturbation theory is no longer valid, and the contributions of those high-order terms must be considered precisely in the nonlocal transport model.

Considering higher-order expansion terms, the components of $f_1^0$ and $f_1^1$ will be coupled together through $f_2^m$ (for m = 0, 1, 2). Therefore, the magnetic field component B_x along the temperature gradient will also affect the transport process. Figure 8 shows the time evolution of the heat flux at the point x = x_max/2 when the magnetic field is not completely perpendicular to the temperature gradient. As can be seen, the truncated order of the spherical harmonic expansion used in the VFP simulation has a relatively weak effect on the heat flux Q_x in the x direction, but a more significant effect on the heat fluxes Q_y and Q_z. With the first-order expansion term only, all the components Q_x, Q_y, and Q_z of the heat flux are overestimated, while with the first two expansion terms, all the components are underestimated. More interestingly, the higher-order terms have a more significant effect on the heat flux components Q_y and Q_z that are perpendicular to the temperature gradient. Generally, we find that the VFP simulation results converge when the expansion is taken up to order ℓ = 8, as has been done in the simulations presented above.

The effect of high-order terms on the distribution of the heat flux is shown in Fig. 9. It is found that the strong preheating effect at the heat front of Q_y and Q_z can be exactly preserved only when the high-order spherical harmonic terms are retained, while the first-order spherical harmonic expansion is already enough to treat the preheating effect at the heat front of Q_x.

The nonlocal thermal transport process in magnetized plasmas has been studied theoretically and numerically with the VFP model, in which the magnetic field has nonzero components both perpendicular to and along the temperature gradient. For this purpose, a VFP code with one-dimensional spatial coordinate and three-dimensional velocity components has been developed. Generally, the magnetic field component perpendicular to the temperature gradient tends to restrain the heat flux along this gradient and induces a transverse heat flux. When the magnetic field has two components along the transverse and longitudinal directions, it is found that the heat flux along the third direction appears via a coupling between the longitudinal magnetic field and the transverse heat flux. Nonlocal heat transport is found in both the longitudinal and transverse directions, provided the temperature gradients are sufficiently large. For nonlocal transport under a magnetic field along an arbitrary direction, the magnetic field will reduce the nonlocality of the heat transport in the direction perpendicular to the magnetic field, i.e., the difference between the heat fluxes predicted by the Braginskii theory and the VFP simulations tends to decrease with increasing magnetic field strength. In real experiments, transverse magnetic fields are induced at laser ablation fronts. We note that external magnetic fields are now being introduced on purpose to control plasma temperature and density profiles.48–50 Generally, a transverse heat flux may be beneficial to the formation of a more uniform ablation front with a higher temperature.51 Furthermore, a stronger nonlocal effect in the transverse direction can reduce heat transport and produce plasmas with different geometrical features, and so it is important to implement an appropriate thermal transport model in MHD simulations.

More importantly, the nonlocal heat flux depends not only on the current but also on the preceding temperature profiles. In the process of heat flux evolution, the response time of the heat flux component along the temperature gradient is much shorter than that of the component perpendicular to the gradient. When the temperature gradient is steep, the contribution of higher-order terms in the spherical harmonic expansion of the electron distribution function becomes important even for weakly magnetized plasmas, especially for the thermal transport in the direction perpendicular to the temperature gradient.

Acknowledgment. This work is supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA25050100), the National Natural Science Foundation of China (Grant Nos. 12135009 and 11975154), and the Science Challenge (Project No. TZ2018005).