Numerical modeling of two magnetized counter-propagating weakly collisional plasma flows in arch configuration Open Access

In astrophysics, the interaction of plasma flows with magnetic structures plays a major role. In particular, in the atmospheres of stars and the magnetospheres of planets, magnetic field configurations in the form of an arch, usually filled with plasma, are widespread. Such magnetized plasma arch structures have been observed on the Sun¹ and in the Earth's magnetosphere.² It has been shown that they play a significant role in solar flares,³ as well as in various explosive processes accompanied by the reconnection of magnetic lines.⁴ 

In addition to space observations, magnetized plasma arches are also studied in the laboratory.^5–7 In controlled conditions, it is possible to conduct a more detailed analysis of the processes occurring in the system and give them a physical explanation. However, due to laboratory limitations, studies are conducted on objects much smaller than space ones and at much shorter times. In this case, the correspondence between laboratory modeling and space observations is achieved by comparing dimensionless parameters that allow scaling of physical processes. Nevertheless, not all interaction regimes observed in space can be implemented in the laboratory within the framework of one or two experimental setups, so it is of interest to develop new approaches to laboratory modeling. Recently, our group proposed and implemented an original method for creating a magnetized plasma arch during an arc discharge in a magnetic field created by two coils located at an angle to each other.⁸ The features of the resulting configuration include the supersonic nature of the plasma expansion, magnetic Mach numbers in a wide range $Mm∼0$–10, and the small rate of collisions compared to the gyrofrequency. Consequently, the resulting plasma is non-equilibrium and the interaction is accompanied by the development of instabilities, including those of a kinetic nature.⁹ 

Despite the rich possibilities for performing various diagnostics in a laboratory experiment, many processes occurring in it have too small scales in time and space to be distinguished by modern methods. In this regard, numerical modeling is widely used to plan experiments and explain their results. As an example, we can cite a recent work on the numerical modeling of the interaction process of a plasma flow incident on a magnetic arch¹⁰ in a configuration similar to the experiment on our setup.^11,12 At present, the work is being carried out on the setup to study the collision process of two counterstreaming flows injected from the bases of a magnetic arch along magnetic lines. The present work is aimed at the numerical modeling of this process and the identification of the main interaction modes and the physical processes accompanying them. The main attention is paid to the emerging structures of plasma density, directly observed by its optical glow in the experiment, as well as to the processes of development of various instabilities and, in particular, the Weibel instability.¹³ 

The experiment is carried out in a vacuum chamber with the ambient pressure below $10−6$ Torr. The experimental setup is shown in Fig. 1. The plasma is generated by a specially designed generator as a result of an arc discharge. The cylindrical generator has an outlet opening with a diameter of 2 cm. It is isolated from the chamber walls by a high-voltage insulator. The cathode, which is the plasma source, is made of aluminum. A detailed study of the plasma generator can be found in Ref. 9. In a typical operating mode, it is capable of generating a plasma flow with a velocity of $≈1.5×106$ cm/s and an electron temperature of about 3 eV. The measured average ion charge is $Z=1.7$⁠. Since the energy of the directed ion motion is, thus, about 33 eV, the flow is supersonic with a Mach number $M≈3.5$⁠. By adjusting the discharge capacitor voltage from hundreds of volts to 5 kV, it is possible to change the discharge current value and, as a consequence, the ion concentration in the flow within the range of $1013$– $1016$ cm⁻³. As the voltage changes, the flow velocity, though, remains practically constant. The typical discharge duration is 20  $μ$s, which allows generating a plasma bunch about 30 cm long, which is comparable to the dimensions of a vacuum chamber, so that during the discharge the plasma is able to completely fill it.

The magnetic field is created by two pulsed coils of the same design and placed at a right angle to each other. The current pulse duration in them is about 3 ms. The magnetic field strength in the center of the coils reaches 3.3 T, but in the chamber it is 10–100 times smaller. For instance, for a discharge current of 1 kA, the field in the center of the chamber is 23 mT.

In the experiment, the arc discharge is initiated with some delay relative to the current discharge in the coils. By this time, a stationary magnetic field of an arched configuration is established in the chamber. In the studied configuration, two plasma generators are located inside the coils, and the plasma flows generated by them are injected from the bases of the arch along the magnetic field lines. In preliminary experiments, at low current in the magnetic coils, the chamber was filled with plasma, observed by its optical glow. At high current values in the coils and, consequently, for a strong magnetic field, a glow in the form of an arch is observed as shown in Fig. 2. Thus, one of the goals of this work will be to study the transition between these two regimes.

This system is usually treated numerically by means of the particle-in-cell method¹⁵ or some kind of grid-based methods.¹⁶ However, performing full-scale modeling with realistic parameters using a fully kinetic electromagnetic code turns out to be technically challenging even in two-dimensional geometry due to relatively low velocities of plasmas: a typical plasma frequency of electrons is of the order of $1012$ s⁻¹, the corresponding time step should be less than 1 ps, i.e., modeling processes with a duration of about 10  $μ$s requires tens of millions of iterations, which is thousand times greater than the typical capabilities of modern programs.

The two-dimensional hybrid simulation was carried out by the AKA code.²¹ This code was previously successfully tested on a model problem for analyzing the magnetic reconnection process in the perturbed Harris sheet on electron spatial scales,²² in the problem of long-term dynamics of the Weibel instability,²⁰ as well as in modeling experiments on laser ablation of plasma in an external magnetic field^23–25 and in self-generated magnetic fields.²⁶ 

The hybrid simulations were close to the conditions of the real experiment. The simulation domain size corresponded to a region of $50×50$ cm, the initial ion concentration in the plasma flows was $1016$ cm⁻³, and the plasma was assumed to be fully ionized and consist of singly ionized aluminum. Plasma flows were generated at the boundary of the simulation domain. The initial plasma velocity was varied in different simulations near the value of $106$ cm/s, the temperature of electrons and ions was assumed to be initially zero. The total simulation time was about 80  $μ$s. The external magnetic field was considered to be given, its configuration was obtained by means of a separate electromagnetic calculation for the known geometry of real coils. Its value in the simulation domain was about 0.25 T, which corresponds to the estimate of the magnetic field in the center of the chamber in the experiment. The boundaries were supposed to be open for particles and fields.

For a qualitative examination of electron kinetic effects, a fully kinetic electromagnetic simulation was also carried out by the particle-in-cell method.¹⁵ This required scaling the system to reduce the computation time. Two main approximations were made: first, the flow velocity was increased to bring it closer to the speed of light, while remaining nonrelativistic—in the simulation it was equal to $109$ cm/s; and second, the ion–electron mass ratio was reduced to 16, which allowed the ion and electron scales to be brought closer to each other. This reduced mass ratio is a common approach in kinetic modeling that allows capturing electron kinetic effects while maintaining a computationally feasible separation between the ion and electron scales. The remaining quantities were chosen so as to preserve the main dimensionless parameters without significant changes: the ratio of the magnetic Mach number is on the order of unity, the ratios of the ion gyroradius and ion inertial length to the arch width are on the order of unity, and the ratio of the electron gyroradius to the arch width is less than unity. We would like to emphasize that maintaining a magnetic Mach number of order unity for higher flow velocities requires setting much higher magnetic fields, however, thanks to scaling theory, this can be reduced to laboratory values, if only the plasma can be considered collisionless as it is in our case. The detailed comparison of dimensional and dimensionless quantities in hybrid and kinetic simulations are given in Table I.

To facilitate comparison between the hybrid and kinetic simulations, the results presented in the following sections are given in the normalized form. In the hybrid simulations, distances are normalized to the ion inertial length $d0=cε0mi/qi2N0$⁠, time is normalized to inverse plasma cyclotron frequency $ωci−1=mi/qiB0$⁠, velocities are normalized to Alfvén velocity $V0=d0ωci=B0/μ0miN0$⁠, magnetic fields are normalized to $B0$⁠, and electric fields are normalized to $E0=V0B0$⁠. Ion charge is supposed to be equal to $qi=e$⁠. The values $N0$ and $B0$ are given in Table I. Specific values of normalized quantities are given in the corresponding figure captions. In the kinetic simulations due to the focus on electron-scale physics, time is normalized to the inverse electron plasma frequency $ωe0−1=ε0me/e2N0$⁠, distances are normalized to the electron inertial length, $de0=c/ωe0$⁠, and velocities are normalized to the speed of light $c=de0ωe0$⁠. We also would like to underline that quantitative comparison of hybrid and fully kinetic simulations is not possible due to computational restrictions in the kinetic case, so the comparison is merely qualitative.

The fully kinetic simulations were performed using the Smilei code,²⁷ which implements a fully electrodynamic, relativistic particle-in-cell method. Unlike the hybrid approach, this code solves the full system of Maxwell's equations, including the displacement current. Previously, we have performed a comparison of the Smilei code against hybrid AKA code for the problem of Weibel instability at ion time scales.²⁰ It showed the qualitative agreement between the two approaches, with similar trends observed in the behavior of the system, despite the exact values of physical quantities were different. Specifically, the maximal magnitude of the self-generated magnetic field and the spectrum of the instability were found to be similar, as well as the residual electron pressure anisotropy level. The code allows one to inject particles during the calculation, but there is no way to specify an external constant field. Therefore, the external magnetic field was simulated by a system of superheavy charged particle currents. The size of the simulation region in both coordinates was 2.4 ion inertial lengths = 9.6 electron inertial lengths. The temperature of the injected particles was zero. The boundary conditions were open for both particles and fields.

A series of calculations were performed using a hybrid numerical modeling method, in which the number of flows and their speed were varied. Note that for the selected parameters of plasma flows, the pressure of an individual flow and the magnetic field pressure in the center of the interaction region are comparable. Accordingly, the global nature of the interaction depends on whether the flow pressure exceeds the magnetic pressure or not. Two typical regimes were observed in the modeling: subcritical at lower flow pressure and overcritical at higher flow pressure.

Figure 3 shows the calculation result for one flow with a velocity of $106$ cm/s. It is evident that the flow as a whole follows the magnetic field lines. At the same time, the plasma concentration during its expansion decreases, since the magnetic field in the center of the calculation region is lower than at the plasma ejection point. This plasma expansion is accompanied by a compression of the magnetic field at both boundaries of the plasma tube being formed. However, at its apex, at late times, a partial plasma ejection is observed.

The transverse plasma escape is explained by the $E×B$ drift, which occurs as a result of an inductive effect: as the plasma tends to push magnetic fields aside due to the diamagnetic effect, the time-dependent magnetic field induces a transversal electric field $Ez$⁠. This field is clearly visible in Fig. 4 at the time of 5  $μ$s.

In the same figure at later times, it is seen that a strong alternating field arises at the arch bottom boundary, polarized perpendicular to the plane and propagating along the magnetic field. This wave has a wavelength close to the ion inertial length $di=cε0mi/e2ni$⁠, and therefore, it can be seen as a surface low-frequency mode excited at frequency close to ion cyclotron frequency $ωci=eB/mi$⁠. Indeed, in the limit of zero frequency for Alfvén waves, we have $ω=kVA$⁠, where $ω$ is a wave frequency, k is a wave number, and $VA=B/μ0mini$ is Alfvén velocity. For $k=1/di$⁠, we get $ω=ωci$⁠. This ion cyclotron wave is expected to grow in the plasma with the anisotropic ion distribution function. Particularly, it should be unstable if $T∥>T⊥$ (⁠ $T∥$ and $T⊥$ being dispersions of particle velocities along and perpendicular to local magnetic field).²⁸ Here, due to short times and low loss rate of particles through arc ends, we have for ions $T∥>T⊥$ as they are injected along magnetic lines with zero transverse momentum so it is, indeed, unstable and can grow.

Figure 5 shows the simulation results for two counter-propagating flows with the same interaction parameters. Note that for the first 15–20  $μ$s the dynamics with two flows is almost the same as the dynamics with one flow. A plasma-filled tube is formed, at the boundaries of which the magnetic field is compressed by the plasma pressure. However, at later times, the expansion of the outer side of the tube occurs much faster than in the case of one flow, which is due to the doubled pressure of the plasma flows.

As can be seen from Fig. 6, the ions in such a configuration are demagnetized, since their gyroradius is comparable to the size of the tube, and the electrons are effectively magnetized—the plasma beta, although large in the center of the tube, is close to unity at its boundary, where the field is compressed.

At later times, the tube is stretched and forms a configuration with oppositely directed magnetic fields. The zero-magnetic-field line is formed along the diagonal of the simulation box. Magnetic reconnection, however, is slow at these parameters. The tearing instability is not observed, so the formation of plasmoids and the separation of the top of the tube from the base do not occur, the configuration slowly evolves but remains stable even 70  $μ$s after the onset of interaction. Such stability is associated with a relatively low pressure of the plasma flows.

Figure 7 shows the simulation result for a plasma flow velocity doubled and, accordingly, four times higher pressures of these flows. It is evident that in this case, a quasi-stationary tube is not formed. Almost immediately, the process of separation of a part of the flows with the formation of plasmoids begins. At later times, a configuration with an elongated tube is formed, but its shape is deformed, and the magnetic field pressure inside is higher. In contrast to the subcritical regime, here, the magnetic reconnection process is more intense which is associated with the turbulent nature of the interaction in this regime. Clearly visible plasmoids with a closed magnetic field configuration break away from the top of the tube (see the upper left corner at time instant 30  $μ$s), which also manifests the highly nonequilibrium and intensely turbulent nature of the observed process. So, we can assume that the magnetic reconnection in this regime is driven mainly by ion dynamics and the Hall term in the generalized Ohm law (8). It is also supported by the fact that the characteristic scale of the field heterogeneities inside the elongated tube is of the order of ion inertial length. Inside the plasma, we see an intense mixing of magnetic field lines with a quasi-turbulent state arising. The hybrid code, though, does not allow us to investigate the turbulization as it does not resolve electron spatial scales.

As in the case of one flow, the formation of a plasma tube in the case of two flows is accompanied by resonant excitation of surface waves at ion cyclotron frequency. This can be seen in Fig. 8. Note that the excitation of these waves is observed for both relatively slow flows and faster ones. There the in-plane component of the electric field is also shown. It is also alternating at tube boundaries which is expected as the waves at this range of frequencies should be elliptically polarized. We also note that the in-plane component exceeds the z-component which means that this is a mode connected to Alfvén wave at lower frequencies. In other words, this is an ordinary elliptically polarized wave.

A feature of the two-flow interaction compared to the case of one flow is a much more pronounced density filamentation along the flows (cf. Figures 3 and 5 at 15  $μ$s). We attribute it to Weibel instability. In the simulation performed, electrons are described hydrodynamically, but taking into account the second-order moments — the components of the pressure tensor. In this case, the Weibel instability manifests itself in the presence of anisotropy of this tensor.

Figure 9 shows the distribution of the anisotropy of the pressure tensor obtained during the simulation. It is evident that a highly anisotropic plasma is formed when the flows collide. At later times, however, the anisotropy in this region of space drops to unity, which may be further evidence that Weibel instability had developed and reached saturation by this time. Initially, the instability leads to filamentation of currents, which then give rise to filamentation of density.

The claim that the observed filamentation is a manifestation of Weibel instability can be further supported by the following estimates. The Weibel instability arises from the anisotropy of electron temperature, leading to countercurrent electric currents being attracted to create periodic filaments. An external magnetic field suppresses this instability by preventing electron shifts across field lines, particularly when the plasma pressure to magnetic pressure ratio is less than a critical value.²⁹ Specifically, in the limiting case of zero transverse temperature, the instability is suppressed if the longitudinal plasma pressure is less than twice the magnetic pressure. Based on Fig. 5 at $t=4$ $μ$s, we can estimate the magnetic pressure at the point of collision of flows to be less than $0.1B02$⁠. The ram pressure of the flows in this case are of the order of $B02/2$ which gives an electron pressure in the moment of collision to be of the order of $B02$ (implying that the thermal energy of electrons in collided flows is equal to the sum of their ram pressures just before the collision). It means that the conditions for the development of Weibel instability are met.

Weibel instability, however, in these conditions cannot be correctly described in a hydrodynamic approximation for electrons, since their distribution function, at least in the early stages of the development of instability, is far from equilibrium. For a qualitative consideration of the kinetically determined features of the developed Weibel instability, fully kinetic modeling was carried out by the particle-in-cell method in the model conditions described in Sec. III.

Kinetic simulation results are shown in Fig. 10. It can be seen that the structure of the tube as a whole coincides with the results of hybrid modeling. At the same time, filamentation of electron density is much stronger, and the filament widths are generally smaller due to the better resolution of electron scales, which is compromised by the smaller size of the overall system compared to ion scales and the smaller ion-to-electron mass ratio. We also note that the filaments are longer and more regular compared to the hybrid simulations, which we attribute to the more pronounced anisotropy of the distribution function in the case of the kinetic plasma description.

It should be noted that in our setup Weibel instability is developing in the presence of an external magnetic field. In this case, the growth of modes with wavelengths lower than the gyroradius is inhibited so it is expected that the mode with a wavenumber $k=1/rce$ will be the most pronounced. This is consistent with the filament width observed in the simulation.

From the anisotropy distribution in Fig. 11, it is evident that a plasma strongly anisotropic along the field is formed in the tube. Electrons and ions circulate along the tube, reflecting from its ends. As a result of this movement, their flow velocity transforms into a thermal velocity longitudinal to the magnetic field, which remains high for a long time relative to the transverse velocity due to the absence of collisions.

Note that here we also observe a slow plasma outflow from the top of the tube and it is mainly associated with the large gyroradius of the ions. In Fig. 12, it is seen that it is of the order of the transverse size of the tube, and accordingly, along with the electric drift, a fairly strong gradient drift is observed.

Let us summarize the main conclusions obtained in the modeling. In the experiment, two interaction regimes should be observed: subcritical, corresponding to a magnetic Mach number $Mm<1$⁠, at which the plasma flow pressure is less than the magnetic pressure, and overcritical, corresponding to $Mm>1$⁠, at which the plasma flow pressure is higher than the magnetic pressure. In the subcritical regime, the magnetic arch is filled with plasma with the formation of a plasma arch, and in the overcritical mode, the plasma breaks through the magnetic lines and fills the entire volume of the chamber. Note that these regimes have already been observed in preliminary experiments with a discharge in magnetic fields of different magnitudes.

In the subcritical regime, due to the diamagnetic effect, the magnetic field at the boundary of the arch is compressed, which is also accompanied by the formation of denser plasma regions. In the experiment, such thickenings can be observed as a brighter glow of the edges of the arch as seen at Fig. 2. Along the boundaries of the arch, excitation of a wave near the ion cyclotron frequency is observed, which can likely be recorded in the experiment by a MHz-range receiver. The top of the plasma arch slowly moves outward due to the $E×B$ drift, which can be recorded in the experiment by a high-speed camera with a submicrosecond shutter, filming the plasma with a variable delay after the start of the discharge. At low values of the magnetic Mach number, however, the speed of this drift will be small due to the relative weakness of the diamagnetic effect.

The plasma arch demonstrates the most interesting dynamics in the near-critical mode. In this case, the arch expansion rate is comparable to the initial flow velocity and the Alfvén velocity. During the discharge, a significant expansion of the plasma is observed with the formation of a region with oppositely directed magnetic fields. At magnetic Mach numbers slightly greater than unity, intense turbulence is observed in this region, accompanied, among other things, by the formation of plasmoids breaking away from the arch. It can be expected that in this case, intense high-frequency radiation will be observed at frequencies near the electron cyclotron resonance, which can be recorded by a GHz-range receiver. In addition, it can be expected that high-energy electrons will be generated, which can also be recorded by an electron spectrometer.

In the plasma tube, the development of Weibel instability with the formation of plasma filaments is also expected. These filaments are also visible in experiment at Fig. 2 by their enhanced optical glow.

A two-dimensional numerical simulation of the process of propagation of one and two counterstreaming flows in an external magnetic field with an arched configuration similar to the experimental one was carried out. Full-scale simulation was carried out using a hybrid method, in which ions are described kinetically, and electrons are described hydrodynamically taking into account the evolution of their pressure tensor. Electron kinetic effects were also investigated in a fully kinetic simulations using the particle-in-cell method in a scaled geometry with a reduced electron-to-ion mass ratio, increased flow velocities and magnetic field strength, and reduced scales of the interaction region. It was shown that as a result of plasma propagation, an arched plasma tube is formed, in which the magnetic field is compressed at the edges. The plasma propagating along the tube is polarized, resulting in an $E×B$ drift outward from the tube. At low flow velocities, this leads to a slow expansion of the tube and its quasi-stationary evolution. As a result, a region with oppositely directed magnetic field lines is formed, but magnetic reconnection in such a configuration is slow. At higher flow velocities, however, the tube evolution is much faster, stationarity is not achieved, and magnetic reconnection is more intense. This highly nonequlibrium process is accompanied by formation of plasmoids detached from the base of the arch.

In both cases, generation of elliptically polarized surface waves along tube boundaries is observed at a frequency close to the ion cyclotron frequency. Inside the tube, density filamentation is observed, caused by the development of Weibel instability, due to the occurrence of electron pressure anisotropy during the interpenetration of plasma flows. Fully kinetic modeling at scaled parameters confirmed the main conclusions of the hybrid modeling.

This research was funded by Russian Science Foundation (Grant No. 23-12-00317). The simulations were performed on resources provided by the Joint Supercomputer Center of the Russian Academy of Sciences.

The authors have no conflicts to disclose.

Artem V. Korzhimanov: Conceptualization (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Sergey A. Koryagin: Writing – review & editing (equal). Andrey D. Sladkov: Investigation (equal); Software (equal); Visualization (equal). Mikhail E. Viktorov: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal).

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