Suppression of electron thermal conduction by whistler turbulence in a sustained thermal gradient
The dynamics of weakly magnetized collisionless plasmas in the presence of an imposed temperature gradient along an ambient magnetic field is explored with particle-in-cell simulations and modeling. Two thermal reservoirs at different temperatures drive an electron heat flux that destabilizes off-angle whistler-type modes. The whistlers grow to large amplitude, , and resonantly scatter the electrons, significantly reducing the heat flux. A surprise is that the resulting steady state heat flux is largely independent of the thermal gradient. The rate of thermal conduction is instead controlled by the finite propagation speed of the whistlers, which act as mobile scattering centers that convect the thermal energy of the hot reservoir. The results are relevant to thermal transport in high astrophysical plasmas such as hot accretion flows and the intracluster medium of galaxy clusters.
Introduction. Thermal conduction is integral to a wide variety of phenomena occurring in space and astrophysical plasmas. Ascertaining the rate of thermal conduction in such systems is therefore of fundamental importance. In particular, the microphysics of weakly collisional, weakly magnetized plasmas, which is not fully understood, may play a pivotal role in determining the transport properties of the global system Rincon et al. (2014). A magnetic field makes thermal conduction anisotropic and when the plasma is greater than order unity, the system is susceptible to a host of microscale kinetic instabilities which tend to suppress thermal fluxes via particle scattering Ramani and Laval (1978); Levinson and Eichler (1992); Gary and Li (2000); Kunz et al. (2014); Komarov et al. (2016); Riquelme et al. (2016). Such instabilities are expected to operate in rarefied plasma environments such as the Intracluster Medium (ICM) of galaxy clusters Schekochihin et al. (2010); Kunz et al. (2011); Roberg-Clark et al. (2016) as well as hot accretion flows Quataert (1997); Sharma et al. (2007) and the solar wind Hellinger and Trávníček (2013); Hollweg (1974); Gary et al. (1994).
The impact of instabilities on transport is tied to their nonlinear evolution, which in turn is influenced by the input of free energy from the surrounding astrophysical environment that can be included through appropriate boundary conditions in numerical models. Examples include those used in shearing-box Kunz et al. (2016); Riquelme et al. (2016) and compressing-box Sironi and Narayan (2015); Sironi (2015) simulations of instabilities driven by pressure anisotropies and their impact on the transport of heat and momentum in accretion flows. Here we focus on the dynamics of a system in contact with two thermal reservoirs at different temperatures, driving a heat flux parallel to an ambient magnetic field. We use particle-in-cell (PIC) simulations to model the resulting heat flux instability and steady-state suppressed thermal conduction with a sustained thermal gradient. This model is in contrast to previous work in which the heat flux instability was studied as an initial value problem Roberg-Clark et al. (2016).
Numerical Scheme. We carry out two-dimensional (2D) simulations using the PIC code Zeiler et al. (2002) to model thermal conduction along an imposed temperature gradient in a magnetized, collisionless plasma with open boundaries. calculates particle trajectories using the relativistic Newton-Lorentz equations and the electromagnetic fields are advanced using Maxwell’s equations. The ends of the simulation domain act as thermal reservoirs at two different temperatures separated by a distance , forming a temperature gradient and driving a heat flux. An initially uniform magnetic field threads the plasma along the gradient and is free to evolve in time. The initial particle distribution function is chosen to model the free-streaming of particles from each thermal reservoir and has the form
where is the initial density, is the Heaviside step function, is the thermal speed, and the parallel and perpendicular directions are with respect to . The cold particles are given a parallel drift speed to ensure zero net current () in the initial state while the error function makes the density of hot and cold particles equal. also has nonzero pressure anisotropy () and a heat flux . is not unstable in a 1D system since only off-angle modes resonate with particles near the large phase space discontinuity in at (Fig. S1a).
When particles exit the open boundaries they are re-injected with velocities pulled from (at ) or (). The drift velocity is then recalculated at each time step to ensure that the current of re-injected particles cancels the current of outgoing particles at the cold reservoir. The electromagnetic field components at the thermal reservoir boundaries are where . Periodic boundary conditions are used for both particles and fields in the direction. Ions in the simulation are not evolved in time and act as a charge-neutralizing background. The subscript denotes an electron quantity.
Simulation Parameters. We have performed six simulations in which , and are varied independently so as to change and . The baseline simulation has , , , , and , where is the gyroradius, is the cyclotron frequency, and is the plasma frequency. The parameters for each simulation are listed in Table 1. Each simulation uses 560 particles per cell, has a transverse length of , and is run to . The largest simulation () has a spatial domain of by grid cells.
Whistler Turbulence. Initializing the simulations with leads to an impulse of transient fluctuations in the out-of-plane magnetic field that propagate towards the hot thermal reservoir (evidence for this is shown later). These fluctuations are driven by the initial pressure anisotropy and quickly lead to a sharp drop in the anisotropy to the marginally stable level for firehose-type modes (not shown). The fluctuations rapidly damp and become dynamically unimportant in the simulations and are not discussed further.
The reinjection and mixing of hot and cold particles results in a continuous source of heat flux in the simulation domain. The heat flux drives off-angle (), slowly propagating (, elliptically polarized whistler modes that reach large amplitude, (fig. 1a), and strongly scatter electrons, isotropizing the electron distribution function (see the supplementary material). The heat flux drops well below its initial value . Some reflection of waves occurs at the cold plate boundary but the heat flux is insensitive to the length of the simulation domain, confirming that such reflection does not impact the integrated results.
Strong scattering by the whistlers causes inherently 2D structures to develop in quantities such as the temperature (figure 1b) and heat flux (figure 1c). In figure 1b the trajectories of four electron macro-particles from the simulation, tracked starting from an initial position for a period of in steady state, are overlaid over , which does not vary appreciably during the time of the orbits. Some particles reverse their parallel velocity several times as a result of scattering in the strong magnetic fluctuations. Because the system is 2D the particle out-of-plane canonical momentum, , is a conserved quantity. Since and kinetic energy is mostly conserved in the magnetic fluctuations, the electrons are confined to relatively narrow channels in .
Suppression of Thermal Conduction. Suppression of the heat flux develops over a time of hundreds of resulting in a steady state in which a continuous temperature profile has formed between the hot and cold reservoirs (fig. 2a) and the heat flux has leveled off to a nearly constant value (fig. 2b). Fig. 2c shows the time profiles of average heat flux for six simulations.
The expectation for a system subject to Coulomb scattering (or another scattering process) is that the heat flux is diffusive, . We find instead that the final heat flux is insensitive to the ambient gradient. The black lines in fig. 2c correspond to simulations with a fixed but differing box lengths or hot to cold temperature jumps. For all of these runs the heat flux settles at around . Thus, the heat flux rather than the gradient controls the dynamics. As long as is significantly greater than , the hot plate controls the final heat flux. However, the two simulations with differing have noticeably different asymptotic heat fluxes that follow the scaling (fig 2c inset). To explain this result we turn to the physics of scattering by large-amplitude whistler waves.
Scattering by Whistlers. The physics of resonant interactions of particles with elliptically polarized whistlers is well-documented in the literature (see e.g. Roberg-Clark et al. (2016) and references therein). In the frame of a single off-angle whistler wave, total kinetic energy is conserved and particles which satisfy the various resonance criteria , are trapped Karimabadi et al. (1992). For , resonant particles experience small oscillations in the plane. For large-amplitude whistlers () resonances can overlap, leading to irreversible diffusive behavior along circular, constant energy curves in the whistler wave frame Karimabadi et al. (1992). In the presence of multiple whistlers with differing parallel phase speeds some diffusion may also occur perpendicular to circles of constant energy Karimabadi et al. (1992). Resonance overlap is an effective mechanism for heat flux suppression since it causes large deflections in the particle pitch angle , quenching the parallel heat flux Roberg-Clark et al. (2016).
To demonstrate that this is the physics at play in our simulations, in fig. 3a we show a resonance diagram in for four trapped particles with differing energy in the simulation with at steady state. Particle energy is mostly conserved and the primary diffusion is in pitch angle Karimabadi et al. (1992). All the particles display significant deflection so the bulk of particles undergo trapping by the whistlers. Also of note is that the nearly-circular contours in velocity space are effectively centered about , indicating that the whistler phase speed is small compared to the thermal speed .
To quantify the rate of scattering by the whistlers we calculate the quantity by averaging over individual trajectories of roughly particles for the simulation (fig. 3b). The diffusion rate is half the linear slope of at late time, where is the scattering time. We find . We plot versus time for particles in fig. 4 to illustrate the particle motion. Some particles are diverted back towards the initial particle location at once scattering becomes significant while others maintain their initial direction of propagation. The linear trend of mean-squared displacement in 3b is evidence for diffusive behavior. Pitch angle scattering in a spectrum of whistler turbulence was also reported by Keenan and Medvedev (2016).
Steady State Heat Flux. The results of fig. 2a have demonstrated that the asymptotic rate of thermal conduction in the presence of large-amplitude whistler waves is largely independent of the temperature gradient and instead follows a scaling . A simple explanation for this result, consistent with a comment in Levinson and Eichler (1992), is that whistlers act as particle scattering centers that propagate at their phase speed and control the net flow of high-energy particles carrying the bulk of the heat flux. The resulting heat flux is simply the product of the phase speed and the thermal energy of the hot plasma, .
The whistler wave phase speed is determined via the cold plasma dispersion relation, . Taking (as in Roberg-Clark et al. (2016)) for whistlers at high , we find
In figure 5a we show a spacetime diagram ( versus ) of the out-of-plane at a single value of . After a transient associated with the anisotropy-driven waves of the initial distribution that was discussed earlier, the whistlers propagate at a nearly uniform speed in the direction of (). To confirm that the unstable modes have , we show the power spectrum for the runs with at and in fig. 5b. The spectra are nearly isotropic in the 2D Fourier space (not shown) so in the spectra shown the energy has been summed over . We find a spectral index of for the modes near although we note that the more important point is to establish that the spectrum peaks near even as varies. A more complete exploration of the spectrum requires simulations with a third spatial dimension. In addition we find that for each of the six simulations in Table 1, , where was measured in the middle of the simulation domain. These results strongly support the scaling
where is a coefficient of order unity. Equation (3) reveals the crucial role of the background magnetic field in facilitating thermal transport since it controls the propagation of whistlers. In the case of a very small magnetic field the whistlers barely propagate and the thermal conduction is virtually shut off. However, no whistler growth was found in a simulation with (not shown), indicating that heat flux suppression by whistlers requires a finite ambient magnetic field. Recent PIC simulations with an imposed thermal gradient suggest that pressure anisotropy driven modes are at play when there is no initial ambient magnetic field Schoeffler et al. (2017). Those results are consistent with the transient growth of fluctuations seen in our simulations in the case of . These reach finite amplitude but then rapidly decay on time scales short compared with the development of the heat-flux instability.
Discussion. A caveat of our model is that the imposed thermal gradient is much larger than that measured in environments such as the ICM Levinson and Eichler (1992). However, the present simulations suggest that the transport is insensitive to the imposed temperature gradient (although the sign of the parallel heat flux is determined by the sign of through the whistler phase speed). The point is that heat flux instability is directly driven by the collisionless heat flux, which depends only on the temperature difference across a domain, rather than the ambient gradient. It seems likely, therefore, that the current results apply to cases in which the temperature gradient is far weaker. A full treatment of the ICM also requires the inclusion of weak collisions not present in our kinetic model.
A question is how the microphysics of whistler scattering will affect heating and thermal conduction in the intracluster medium. The scaling of heat flux in (3) with implies a suppression factor of roughly below the free-streaming thermal conduction. The functional dependence is a noticeable departure from the Spitzer conductivity Spitzer (1962) proportional to often used in hydrodynamic or MHD models of the ICM (e.g. Ruszkowski and Oh (2010), Yang and Reynolds (2016)). Our results may therefore significantly alter the equilibria associated with clusters of galaxies, which result from a balance between thermal conduction and radiative cooling.
Our results show promising similarities with the observations of thermal conduction in the solar wind by Bale et al. Bale et al. (2013) in which the heat flux takes on a constant value, independent of collisionality and the ambient temperature gradient, in the weak collisionality regime where the collisional mean-free-path exceeds the temperature scale length. However, much of their data is in a regime of much lower than in the present simulations. The exploration of the transition from high to low with analysis and simulations is underway so that more detailed comparisons with solar wind observations can be made.
- F. Rincon, A. A. Schekochihin, and S. C. Cowley, MNRAS 447, L45 (2014).
- A. Ramani and G. Laval, Physics of Fluids 21, 980 (1978).
- A. Levinson and D. Eichler, The Astrophysical Journal 387, 212 (1992).
- S. P. Gary and H. Li, The Astrophysical Journal 529, 1131 (2000).
- M. W. Kunz, A. A. Schekochihin, and J. M. Stone, Physical Review Letters 112, 1 (2014).
- S. V. Komarov, E. M. Churazov, M. W. Kunz, and A. A. Schekochihin, MNRAS 477, 467 (2016).
- M. A. Riquelme, E. Quataert, and D. Verscharen, The Astrophysical Journal 824 (2016).
- A. A. Schekochihin, S. C. Cowley, F. Rincon, and M. S. Rosin, MNRAS 405, 291 (2010).
- M. W. Kunz, A. A. Schekochihin, S. C. Cowley, J. J. Binney, and J. S. Sanders, MNRAS 410, 2446 (2011).
- G. T. Roberg-Clark, J. F. Drake, C. S. Reynolds, and M. Swisdak, The Astrophysical Journal 830, L9 (2016).
- E. Quataert, The Astrophysical Journal, Volume 500, Issue 2, pp. 978-991. 500, 978 (1997).
- P. Sharma, E. Quataert, G. W. Hammett, and J. M. Stone, The Astrophysical Journal 667, 714 (2007).
- P. Hellinger and P. M. Trávníček, Journal of Geophysical Research: Space Physics 118, 5421 (2013).
- J. V. Hollweg, Journal of Geophysical Research 79, 3845 (1974).
- S. P. Gary, E. E. Scime, J. L. Phillips, and W. C. Feldman, Journal of Geophysical Research 99, 23391 (1994).
- M. W. Kunz, J. M. Stone, and E. Quataert, Physical Review Letters 117, 1 (2016).
- L. Sironi and R. Narayan, The Astrophysical Journal 800, 88 (2015).
- L. Sironi, The Astrophysical Journal 800 (2015).
- A. Zeiler, D. Biskamp, J. F. Drake, B. N. Rogers, M. A. Shay, and M. Scholer, Journal of Geophysical Research 107, 1230 (2002).
- H. Karimabadi, D. Krauss-Varban, and T. Terasawa, Journal of Geophysical Research 97, 13853 (1992).
- B. D. Keenan and M. V. Medvedev, Journal of Plasma Physics 82, 905820207 (2016).
- K. M. Schoeffler, N. F. Loureiro, and L. O. Silva, arXiv preprints (2017), arXiv:1707.06069v1 .
- L. Spitzer, Physics of Fully Ionized Gases, 2nd ed. (New York: Interscience, 1962).
- M. Ruszkowski and S. P. Oh, The Astrophysical Journal 713, 1332 (2010).
- H.-Y. K. Yang and C. S. Reynolds, The Astrophysical Journal 818, 181 (2016).
- S. D. Bale, M. Pulupa, C. Salem, C. H. K. Chen, and E. Quataert, The Astrophysical Journal 769, L22 (2013).