1. Introduction
The astrosphere is a region of space filled with plasma from stellar wind, surrounded by a boundary called the astropause. This boundary is a tangential discontinuity that separates the stellar wind from the interstellar medium. The most well-known example of an astrosphere is the heliosphere, which surrounds our Sun. The first two-shock model of the heliosphere was proposed by Baranov et al. (Reference Baranov, Krasnobaev and Kulikovskii1970) and has since received widespread recognition and further development (see, for example, Baranov & Malama (Reference Baranov and Malama1993); Malama et al. (Reference Malama, Izmodenov and Chalov2006); Izmodenov et al. (Reference Izmodenov, Malama and Ruderman2009)). For the heliosphere, the astropause is called the heliopause. The shape of the heliopause, according to this model, is a figure of rotation resembling a paraboloid. However, in the presence of magnetic fields or other physical phenomena, the astropause becomes highly asymmetrical (Izmodenov & Alexashov Reference Izmodenov and Alexashov2015; Izmodenov & Alexashov Reference Izmodenov and Alexashov2020). This shape has been observed indirectly in the astrospheres of different stars (see, for example, Decin et al. (Reference Decin, Cox and Royer2012); Kobulnicky et al. (Reference Kobulnicky, Chick and Schurhammer2016); Kobulnicky et al. (Reference Kobulnicky, Schurhammer and Baldwin2017)) and is currently considered the classic shape of an astropause.
However, recently, several scientific groups have challenged this assumption about the shape of the heliopause. In particular, it has been suggested that heliopause may have a two-jet collimated shape (Opher et al. Reference Opher, Drake, Zieger, Zieger and Gombosi2015; Drake et al. Reference Drake, Swisdak and Opher2015). There is a justification for such an exotic shape. The explanation of such a shape can be found in the influence of the azimuthal magnetic field of the Sun, which directs the solar wind towards the north and south poles of the star, creating two astrospheric jets. The first mentions of this effect can be traced back to the work of Yu (Reference Yu1974). Later studies (Golikov et al. Reference Golikov, Izmodenov, Alexashov and Belov2017b; Golikov et al. Reference Golikov, Izmodenov and Alexashov2017a) have shown that the stellar magnetic field can create a two-jet collimated shape of astropause when the interstellar medium is at rest (Figure 1a).
Figure 1. Schematic representation of a tube-like astrosphere in the case of (a) stationary flow and (b) self-oscillations. The figure also shows the astropause (AP) and termination shock (TS).
Subsequent studies of the astrospheric shape with a moving (with respect to the star) interstellar medium were conducted by Korolkov & Izmodenov (Reference Korolkov and Izmodenov2021). The authors showed that the interstellar medium flow pressure first distorts the shape of the tube and then collapses it, transforming it into a classical paraboloid shape. They estimated the Mach number of the interstellar medium at which the collapse occurs. For the Sun, this is approximately 0.3. Stars with strong magnetic fields can resist the interstellar medium flow for a large Mach number, so a two-jet collimated astrosphere is more likely to be seen near these types of stars.
As for the heliopause, there is no clear consensus on its shape at the moment (see Kleimann et al. (Reference Kleimann, Dialynas and Fraternale2022), Section 8). Other physical processes may maintain a two-jet collimated shape even at strong interstellar medium flow. However, the authors of this paper are more inclined towards the classical form, while acknowledging the significant influence of the azimuthal magnetic field of the Sun. This field forms the dense plasma streams in the heliosheath that are observed in our simulations, but the shape of the heliopause remains classical. These jets are directed along the polar regions of the Sun at close distances and then deviate towards the tail region (see Izmodenov & Alexashov (Reference Izmodenov and Alexashov2015), Figure 6; Kornbleuth et al. (Reference Kornbleuth, Opher and Baliukin2021), Figure 1). It should be noted that some studies (Pogorelov Reference Pogorelov2016; Pogorelov et al. Reference Pogorelov, Fichtner and Czechowski2017) question the collimation of the solar wind by the azimuthal magnetic field. These studies suggest that the solar cycle effects cause the density of the solar wind to be highest near the Sun’s equatorial plane.
In this way, the two-jet collimated shape of the astropause is a unique phenomenon. The azimuthal stellar magnetic field can contain astrospheric jets at a large distance from a star, similar to the way plasma is contained in a tokamak using an azimuthal magnetic field. Understanding these processes is essential for advancing our knowledge of plasma physics.
Further studies two-jet collimated shape of the astropause focused on its stability. In Opher et al. (Reference Opher, Drake and Zank2021), it was shown that interstellar atoms could lead to Rayleigh-Taylor instability in a heliospheric jet. In the paper, it was proposed that charge exchange with hydrogen atoms acts as an effective gravity force across the jet, causing instability at the boundary between the dense core and the surrounding rarefied plasma (see dense flow in Figure 1a). Such instability (1) can be the key to the opening of the tail and its Croissant-like shape, (2) can result in mixing of the stellar wind with the interstellar medium, as described by Opher et al. (Reference Opher, Drake, Zieger, Zieger and Gombosi2015), and (3) can also be a source of turbulence in the shock layers of astrospheres or heliospheres.
Recent numerical simulations by Ma et al. (Reference Ma, Opher, Kornbleuth and Toth2024) demonstrate that neutral hydrogen atoms in heliospheric jets can generate a distinct low-velocity solar wind region. The velocity shear at its boundary provides favourable conditions for the development of a Kelvin-Helmholtz (K-H) instability. Quantitative analysis of the simulation data yields a K-H instability growth timescale of 5-7 years.
The present study continues the investigation of the stability of the two-jet collimated astropause. In a sense, it even takes a step back from the work of Opher et al. (Reference Opher, Drake and Zank2021), as we find that astrospheric jets in strong stellar magnetic fields are unstable even without the influence of charge exchange. Instability leads to a deflection of the jet, which then causes the magnetic field to react and a feedback mechanism to arise, trying to bring the jet back to its original position. As a result, undamped oscillations from side to side, or self-oscillations, occur (see Figure 1b). We study the dependence of the self-oscillation period on the problem parameters. We first consider the case when the star is at rest with respect to the interstellar medium and study the instability of a two-jet collimated astropause. In this scenario, the problem has cylindrical symmetry, allowing us to verify the results from a 3D model by using a 2D one. Then in section 3.3 we show that for the stars with strong magnetic fields (approximately 3-4 times stronger than solar), the instability persists even when the motion with respect to the interstellar medium is included. The astropause remains unstable as long as it is collimated into two jets. Once the dynamic pressure of the interstellar medium becomes sufficiently strong, the heliopause transitions to a classical (nearly paraboloidal) shape, and the instability disappears.
The paper structure is as follows: in section 2, the model, boundary conditions, and numerical method are described, and the problem is formulated in dimensionless form; section 3 shows the results of a 2D (cylindrically symmetric) model for the interstellar medium at rest (Subsection 3.1), followed by a 3D model for both at rest and moving interstellar medium (Subsections 3.2 and 3.3, respectively); finally, in the section 4, the results are summarised and discussed.
Figure 2. Streamlines and density isolines for
$\chi = 4$
and various
$M_A$
. OX - rotation axis of the star. The interstellar medium at rest. Cylindrically symmetric 2D calculation with HLLC-type+TVD scheme (4096x2048 cells). Video is available for panel a (“Video-1.avi”, see Data Availability Statement).
2. Model
We consider the plasma within the framework of ideal magnetohydrodynamics (MHD), disregarding any dissipative processes. In our work, we disregard the influence of hydrogen atoms. We aim to demonstrate that instability can arise in the astrospheric jet even without considering the charge exchange effect, which is unexpected. Furthermore, we first assume that the interstellar medium is at rest and non-magnetised. These assumptions allow us to exploit the axial symmetry of the problem and use 2D simulations together with 3D simulations for additional verification. The results of calculations with a moving interstellar medium will be presented in a separate subsection (see Subsection 3.3). The equation of state for the plasma is
$\mathrm{p} = 2\mathrm{n}\mathrm{k}_{\mathrm{B}}\mathrm{T}$
, where p is the pressure, n is the number density of plasma protons,
$\mathrm{k}_{\mathrm{B}}$
is Boltzmann’s constant, and T is the temperature. At the inner boundary, we adopt a Parker spiral solution (unipolar), with:
\begin{equation}B_{\unicode{x03C6}} = B_{\unicode{x03C6}, E} \left( \dfrac{R_E}{R}\right) \sin{\theta}.\end{equation}
The radial component of the magnetic field has been ignored, as its effect is negligible due to its quadratic decrease.
In the case of a moving interstellar medium, we assume that the stellar rotation axis is perpendicular to the direction of motion of the interstellar medium. We define a coordinate system associated with the star: the X-axis aligns with the stellar rotation axis, the Y-axis points in the opposite direction to the interstellar medium flow, and the Z-axis completes a right-handed orthogonal triad.
2.1. Physical parameters
We will consider plasma as a gas with constant heat capacities (
$\gamma = 5/3$
). The solution to the problem formulated above depends on the following five parameters:
\begin{equation}V_0, \ \dot{M}_\star, \ B_{\unicode{x03C6},E} R_E, \ \unicode{x03C1}_{\infty}, \ \mathrm{p}_{\infty}.\end{equation}
Here
$B_{\unicode{x03C6},E}$
is the magnetic field of the stellar wind at the distance
$R_E$
, which can be chosen arbitrarily while the Parker spiral (1) is assumed.
$\dot{M}_\star$
is the rate of mass loss of the star (
$\dot{M}_\star = 4\unicode{x03C0} \unicode{x03C1} V R^2$
).
$\unicode{x03C1}_{\infty}$
and
$\mathrm{p}_{\infty}$
are the density and pressure of the interstellar medium, respectively.
$V_0$
is the terminal velocity of the stellar wind.
When considering the motion of the interstellar medium (see Subsection 3.3), another parameter
$V_\infty$
arises, which represents the velocity of the interstellar medium.
2.2. Dimensionless formulation and boundary conditions
It is much more convenient to formulate the problem in a dimensionless form that has only two dimensionless parameters. Let us choose
$\unicode{x03C1}_{\infty}$
as the characteristic density and
$a_\infty = \sqrt{\gamma \mathrm{p}_\infty/\unicode{x03C1}_\infty}$
as the characteristic velocity. The characteristic distance is determined by
\begin{equation}R_{*} = \sqrt{\dfrac{V_0 \dot{M}_\star}{4\unicode{x03C0} \unicode{x03B3} \mathrm{p}_\infty}}.\end{equation}
The system of MHD equations will remain unchanged, and the boundary conditions (at
$\hat{R}_{\mathrm{in}}$
) will be as follows:
\begin{eqnarray}\hat{V}_{\mathrm{in}} = V_0/a_{\infty} = \chi,\ \, \hat{\unicode{x03C1}}_{\mathrm{in}} =\frac{1}{\chi^2 \hat{R}^2_{\mathrm{in}}},\, \hat{B}_{\mathrm{in}} = \frac{\sqrt{4 \unicode{x03C0}}}{\hat{R}_{\mathrm{in}} M_A },\end{eqnarray}
hats (
$\hat{\cdot}$
) denote dimensionless variables.
As a result, the problem depends on two dimensionless parameters: (1) the Alfvén Mach number of the star (
$M_A = \frac{V_0 \sqrt{4\unicode{x03C0} \unicode{x03C1}_E}}{B_E}$
) and (2) the ratio of the stellar terminal velocity to the speed of sound in the interstellar medium (
$\chi = V_0/a_{\infty}$
).
$\hat{R}_{\mathrm{in}}$
was chosen 0.3 for the 3D model and 0.07 for the 2D model. It could have been chosen arbitrarily, but close enough to the star, while the Parker spiral (1) is still valid.
It is important to note that the parameter
$\chi$
does not affect the stationary solution of the problem, as explained in detail in Korolkov & Izmodenov (Reference Korolkov and Izmodenov2021) (see equation (10) in the paper). However, for a non-stationary problem, it may be important. Specifically, the growth of instability can depend on this parameter, as was the case for the Kelvin–Helmholtz (K-H) instability at the astropause, as seen in figure 2 of Korolkov & Izmodenov (Reference Korolkov and Izmodenov2023). Therefore, we will keep this parameter in the current paper for further study.
In the case of a moving interstellar medium, another dimensionless parameter appears: the Mach number of the incoming flow
$M_\infty = {V_\infty}/{a_\infty}$
. In this case
$\hat{V}_\infty = M_\infty$
.
Below are the values of the dimensionless parameters for the problem in the case of the heliosphere, assuming that the parameters of the solar wind and interstellar medium are as follows:
-
- solar wind (at the Earth orbit):
$n_{p} = 6.5 \textrm{cm}^{-3}, V_0=420 \textrm{km/s}, B_E=37.5 \mu\textrm{G},$
-
- interstellar medium:
$n_p = 0.06 \textrm{cm}^{-3}, T_\infty = 6530 \textrm{K}, V_\infty = 26.4 \textrm{km/s}$
.
Using these values, we can calculate the dimensionless parameters:
$\chi = 31.94,\ \ M_\infty = 1.97,\ \ M_A = 13.34$
(
$T_* = 115$
years).
Let us note immediately that the instability studied in this work will be present for stars with magnetic fields that are 3 or more times stronger than that of the Sun (
$M_A \lesssim 4$
).
2.3. Numerical approach
In the numerical calculations, we used both 2D (cylindrically symmetric) and 3D models.
The 2D model uses a Cartesian grid (in cylindrical coordinates) with a resolution of 4096x2048 cells. This resolution corresponds to approximately 500-550 cells before the termination shock (by radius along the coordinate axis). For the heliosphere, the distance to the shock wave is approximately 90 AU, which means the cell size is approximately 0.16 AU. The equations are solved using the finite volume method with the HLLC-type Riemann solver (see Gurski (Reference Gurski2004)). To increase the spatial order of approximation, a TVD reconstruction with a minmod limiter is used. Calculations are performed on GPUs using CUDA technology.
The 3D model uses a Cartesian grid with refinement, similar to adaptive mesh refinement (AMR). We refine the grid in advance to calculate the stellar spherical wind and the flow in an astrophysical jet in a qualitative manner. Figure 5 is designed to demonstrate grid refinement within the astrospheric jet (dots in the figure represent cell centres). The instability develops precisely in this region. The grid contains 68 cells along the jet radius (corresponding to 136 cells in diameter). We maintained approximately cubical cell shapes to ensure uniform resolution in all directions. Assuming a jet radius of 60 AU for the heliosphere, this corresponds to a cell size of 0.88 AU. Note that the grid refinement is static - the mesh is pre-constructed with higher resolution in critical regions (the jet and supersonic wind) rather than being adapted during simulations. We use the symmetry of the problem relative to the YOZ plane and solve the problem in a half-space. The equations are solved on GPUs using the finite volume method with the HLLD Riemann solver (see Miyoshi & Kusano (Reference Miyoshi and Kusano2005)) and TVD reconstruction (minmod limiter). The grid employed for this paper has about 32 million cells. It should be noted that we also verified the numerical results using HLL and HLLC-type Riemann solvers, as well as first-order accuracy in space. The results with different numerical schemes are very similar.
All calculations (for both 2D and 3D models) were performed using the following methodology. At the initial time step, the stellar wind solution was defined within a small radius (
$\hat{R}_{\mathrm{in}}$
) centred on the star (the origin of the reference system). The remaining computational domain was initialised with interstellar medium parameters. The simulation was advanced until either:
-
• A steady-state solution was achieved (in stable configurations), or
-
• A quasi-periodic regime emerged (in cases where instabilities developed).
To fully resolve the dynamics of the astrospheric jet, including multiple oscillation cycles, the simulation runtime was extended beyond the stabilisation/quasi-periodicity threshold. The total integration time reached 10–12 dimensionless units (equivalent to 1150–1380 years for heliospheric parameters). In section 3.3, calculations were performed up to a dimensionless time of T = 60, which corresponds to 7,000 years for heliospheric parameters. This duration allows disturbances to traverse the jet
$\sim$
90 times at the Alfvén speed. This allowed us to be confident that the solution had fully entered the periodic mode with auto-oscillations.
3. Results and discussion
This section is divided into three parts. Subsection 3.1 presents the results for a model with an interstellar medium at rest. In this case, symmetry in the boundary conditions allows us to assume that the resulting flow is axisymmetric. Therefore, we could search only for an axially symmetric solution, in principle. However, as it will be seen further, our numerical experiments show that axial symmetry holds only until the flow loses stability. After that, asymmetric disturbances appear and grow. An axially symmetrical (2D) model places additional non-physical restrictions on the flow by prohibiting asymmetric disturbances. This restriction may lead to the formation of an unphysical vortex. However, the axisymmetric approach is still attractive because it allows us to perform calculations with a high resolution of numerical grids. This enables us to more accurately determine the onset of instability and explore the flow characteristics up to this moment. So, we use the 2D model until the flow loses stability.
Subsection 3.2 presents the results for a 3D model to investigate asymmetric manifestations of the instability, which are revealed in the form of self-oscillations in the astrospheric jet. Subsection 3.3 presents moving interstellar medium results obtained using a 3D model.
3.1. Interstellar medium at rest: axisymmetric model
Figure 2 shows the streamlines and density isolines for
$\chi = 4$
and different values of
$M_A$
. The smaller the
$M_A$
, the stronger the stellar magnetic field, and the smaller the width of the astrosphere tube. For
$M_A = 8$
(panel a), the width of the tube along the Y-axis is slightly greater than two (dimensionless units), whereas for
$M_A = 3$
(panel d), it is approximately 1.5 in dimensionless units.
For all panels in Figure 2, we can see the Kelvin–Helmholtz (K-H) instability of the astropause. This instability occurs due to the difference in tangential velocity components on both sides of the tangential discontinuity, as described by Landau & Lifshitz (Reference Landau and Lifshitz2013) (section 29; section 84, task 1). Although this instability for a tube-like astrosphere has not been discussed in the literature known to us, we should assume that its properties are similar to those of the instability of classical astrospheres/heliosphere. For these cases, there are numerous works available (for example, Fahr et al. (Reference Fahr, Neutsch, Grzedzielski, Macek and Ratkiewicz-Landowska1986); Baranov et al. (Reference Baranov, Fahr and Ruderman1992); Ruderman & Fahr (Reference Ruderman and Fahr1995); Korolkov & Izmodenov (Reference Korolkov and Izmodenov2023)). In this study, we will not focus on the K-H instability, as its nature is well understood. Note that it does not destroy the astrosphere and grows in space, not in time; in this regard, it is convective. The solution can be described as quasi-stationary or quasi-periodic. Video file “Video-1.avi” (see Data Availability Statement) demonstrates the flow for Panel (a).
Figure 9 from Golikov et al. (Reference Golikov, Izmodenov, Alexashov and Belov2017b) displays results qualitatively similar to those in Figure 2. However, two critical methodological differences explain the absence of Kelvin-Helmholtz (K-H) instabilities in their study: (1) their grid explicitly captured the astropause and (2) they employed a relatively coarse computational grid (60
$\times$
200 cells, corresponding to
$\sim$
8 AU across the astrospheric jet for parameters close to heliosphere) that was insufficient to resolve high-frequency perturbations.
Panel (d) in Figure 2 shows a vortex that is formed near the axis of symmetry at X
$\sim$
1. This is a manifestation of the global instability of the flow, which is clearly seen in the frame of the 3D simulation (see Figure 4 and section 3.2). This instability has asymmetric modes and cannot be correctly described by axisymmetric simulations. Instead, a global-scale vortex appears.
For smaller values of stellar magnetic field (larger values of
$M_A$
), some precursors of the loss of stability are seen in the streamlines, which are close to the axis of symmetry at 1 < X < 2 (see panels (b) and (c) of the Figure). Note also that some features of the instability loss can be seen in the result of calculations with a low spatial resolution of the computational grid. This is confirmed by our numerical tests and also seen in the results of Golikov et al. (Reference Golikov, Izmodenov, Alexashov and Belov2017b) (see Figure 12, C and D).
In our opinion, the above-mentioned global instability is connected with the formation of a high-pressure zone at the axis of symmetry downwind of the termination shock. This zone is clearly visible in Figure 3a, which shows the pressure isolines for
$M_A = 5$
. The high-pressure zone is formed due to the magnetic force acting across the jet, which attracts the flow towards the X-axis. The first mention of this zone can be found in Drake et al. (Reference Drake, Swisdak and Opher2015) (figure 2), and then in Opher et al. (Reference Opher, Drake and Zank2021) (figure 1). Panel (b) in Figure 3 shows the force profile across the channel (at X = 1.5) and the pressure profile. The pressure forms in such a way that its gradient compensates for the magnetic force.
Figure 3. The same calculations as in Figure 2. (a) Streamlines and pressure isolines for
$M_A = 5$
. (b) Pressure (P) and magnetic force (Y-component,
$F_y$
) distributions across the jet along the line X = 1.5 for
$M_A = 5$
. (c) Pressure distributions along the X-axis for different values of
$M_A$
.
Figure 4. (a) Density isolines in the cross-section of the astrospheric jet at two moments in time (5.0 and 5.7). The right panel shows the streamlines in this plane. The simulation was carried out for the parameters
$M_A = 3,\ \chi = 4$
(see “Video-2.avi” in the Data Availability Section). (b) Distributions of magnetic tension
$F_{T_B}$
, magnetic pressure
$F_{p_B}$
, and total magnetic force
$\mathbf{F}_{\mathrm{mag}}$
along cut-line 1. (c) Distributions of plasma density
$\unicode{x03C1}$
, x-axial velocity component
$V_x$
, velocity magnitude
$|\mathbf{V}|$
, magnetic field magnitude along cut-line 1. (d) Distributions of plasma density
$\unicode{x03C1}$
, magnetic field magnitude
$|\mathbf{B}|$
, pressure p, velocity magnitude
$|\mathbf{V}|$
along cut-line 2.
An increase in pressure near the X-axis leads to the formation of a pressure gradient along the X-axis. The gradient acts as an effective force hindering flow motion. Pressure profiles along the X-axis for different Alfven Mach numbers are shown in Figure 3c. The lower the
$M_A$
, the higher the pressure gradient behind the termination shock (TS). This gradient causes the plasma to flow around the high-pressure area. As a consequence, instability arises. This is somewhat similar to how instability can occur behind a streamlined body or vortex (Kármán vortex street).
To confirm this hypothesis, we developed a simplified toy model that simulated the flow of plasma/gas through a flat channel with a transverse force field. We found that this force creates a region of high pressure on the channel’s axis, which causes instability similar to that of an astrophysical jet. While the results of this numerical experiment are not presented in this paper, they can be found in the paper by Korolkov & Izmodenov (Reference Korolkov and Izmodenov2024).
As seen from the toy model or further research (see Subsection 3.2), this instability is asymmetric and cannot be described using a cylindrically symmetric model. Therefore, a 3D model is needed for further research.
It should be emphasised that the high-frequency oscillations observed in the profiles in Figure 3b,c can be attributed to the K-H instability at the astropause. The instability generates disturbances that propagate against the flow since the flow is subsonic.
3.2. Interstellar medium at rest: 3D model
This section describes the results of 3D modelling for the problem with the interstellar medium at rest.
Figure 4a (and “Video-2.avi”, see Data Availability Statement) shows snapshots of the density isolines and the streamlines in the XOY plane (longitudinal section of the astrospheric jet). For this calculation,
$M_A = 3$
and
$\chi = 4$
. Snapshots are presented at two different times (
$\approx 5$
and
$\approx 5.7$
in dimensionless units) to show the jet’s self-oscillations. We can see that the jet fluctuates from one side of the jet to the other. When converted to physical units for a heliosphere (
$T_* = 115$
years), these time markers correspond to 575 and 655 years of simulation time.
Figure 4b displays the transverse force distribution within the astrospheric jet (cut-line 1 in Figure 4a). The figure shows: magnetic tension
$F_{T_B} = (B\cdot\nabla)B/(4\unicode{x03C0})$
, magnetic pressure
$F_{p_B}=-(\nabla B^2)/(8\unicode{x03C0})$
, and total magnetic force
$\mathbf{F}_{\mathrm{mag}} = ([\nabla \times \mathbf{B}] \times \mathbf{B})/(4\unicode{x03C0})$
. Moreover,
$\mathbf{F}_{\mathrm{mag}} = F_{T_B} + F_{p_B}$
. The magnetic force exhibits a restoring tendency, acting to return the plasma flow to its initial symmetric configuration after deviations from the symmetry axis. Notably, magnetic tension and pressure demonstrate comparable magnitudes.
Figure 4c presents the parameter distribution along cut-line 1 in Figure 4a. It shows plasma density
$\unicode{x03C1}$
, x-axial velocity component
$V_x$
, velocity magnitude
$|\mathbf{V}|$
, and magnetic field magnitude
$|\mathbf{B}|$
. The plots clearly reveal parameter asymmetries induced by plasma flow deflection. Based on the small difference between
$V_x$
and
$|\mathbf{V}|$
, it can be concluded that the stream flows predominantly along the astrospheric jet.
Figure 4d presents the parameter distribution along cut-line 2 in Figure 4a. It shows plasma density
$\unicode{x03C1}$
, magnetic field magnitude
$|\mathbf{B}|$
, pressure p, velocity magnitude
$|\mathbf{V}|$
. The magnetic field oscillations enable estimation of the instability wavelength at 0.7 dimensionless units. For heliospheric parameters, this corresponds to approximately 90 AU. These observations confirm that the investigated self-oscillations occur in the low-frequency regime. The observed self-oscillation period is 1.6 in dimensionless units or
$\sim$
26 years for parameters close to heliospheric (see “Video-2.avi” and Figure 6).
Figure 5 shows the same self-oscillations in the YOZ plane (for X = 1, the simulation is the same as in Figure 4a, but the time points are different). Thin black curves show the magnetic field lines. A colour map indicates the density. In this plane, the instability is seen as the oscillation of the point of density maximum. The magnetic field lines are also concentrated around this point. The oscillations are diagonal and not directed along the coordinate axis. The choice of the oscillation direction is physically equiprobable but numerically related to the computational grid. The cell centres are marked with dots in Figure 5 to show the mesh resolution across the jet. Based on this figure, it is possible to calculate, for example, that the diameter of the astrospheric jet corresponds to 136 grid cells.
Figure 5. The same calculations as in Figure 4, but on the YOZ plane for X = 1 (the time points are different). The streamlines represent the magnetic field lines, and the density isolines are marked with colour. Dots show the mesh resolution across the jet, and arrows schematically indicate the direction and magnitude of the magnetic force (
$F_{\mathrm{mag}}$
).
The arrows in the left panel illustrate schematically the direction of the magnetic force (
$\mathbf{F}_{\mathrm{mag}}$
). This force also has a component along the jet (X-component), but it is much smaller than the components in the YOZ plane. The main direction of this force is towards the centre of the jet. The force is stronger on the side where the jet is deflected. This provides a feedback mechanism that helps the jet return to its original position.
We emphasise that our simulations do not reveal significant velocity shear across the jet (Figure 4c). This key difference from Ma et al. (Reference Ma, Opher, Kornbleuth and Toth2024) results stems from our exclusion of neutral hydrogen atoms in the model, which eliminates the physical conditions necessary for Kelvin-Helmholtz (K-H) instability development. The observed instability instead manifests as large-scale self-oscillations of the astrospheric jet, exhibiting periodic deflections in two directions (most clearly visible in the supplementary video materials “Video-2.avi”). When scaled to heliospheric parameters, these oscillations exhibit a period of approximately 26 years. Notably, the visual characteristics of this instability differ fundamentally from classical K-H instability. Careful examination of the oscillation dynamics in the video materials shows no recognisable K-H features (e.g. characteristic vortex roll-up or wave-like patterns at the shear interface).
Thus, the jet undergoes periodic self-oscillations. We calculated the period of the global oscillations of the jet for different values of parameters (
$M_A$
and
$\chi$
). To do this, we expanded the density at a specific point inside the jet as a function of time into a Fourier series. One dominant harmonic (with its multiples) exists, which corresponds to the oscillation period of the jet. Similarly, we can determine the oscillation period from a video recording of the oscillations by counting the number of deviations over a certain period. Both methods yield the same result. Figure 6 illustrates the relationship between the oscillation period (T) and the
$\chi$
. The dots indicate the calculation results. Each curve represents a different value of
$M_A$
. We found an inverse dependence of the oscillation period on
$\chi$
. The curves in the figure represent
$T = 4.178/\chi$
,
$6.4/\chi$
, and
$7.2/\chi$
for
$M_A$
= 2, 3, and 4, respectively. When
$M_A$
is greater than 4, no oscillations are observed.
The inverse relationship between the oscillation period and
$\chi$
is because the Alfvén velocity is directly proportional to
$\chi$
, so an increase in
$\chi$
leads to an increase in the speed of sound. This results in faster jet oscillations as the feedback mechanism becomes more rapid.
We did not observe any high-frequency harmonics in the Fourier series expansion of the density. This may be due to two possible reasons: (1) the insufficient grid resolution used in 3D simulations, which may have cut off high-frequency harmonics or (2) the presence of a dominant harmonic can suppress secondary oscillation modes, a phenomenon commonly observed in dynamic systems with forced oscillations. This effect has been demonstrated in heliospheric studies, where Korolkov & Izmodenov (Reference Korolkov and Izmodenov2023) showed that periodic solar wind disturbances can suppress Kelvin-Helmholtz instabilities.
Figure 6. Dependence of the self-oscillation period on the
$\chi$
for different values of
$M_A$
with the interstellar medium at rest. The dots indicate the calculation results, while the lines show their approximations. The curves in the figure represent
$T = 4.178/\chi$
,
$6.4/\chi$
, and
$7.2/\chi$
for
$M_A$
= 2, 3, and 4, respectively.
Figure 7. The results for a moving interstellar medium with the Mach number
$M_\infty$
equal to 0.6 (a - d2) and 2 (e, f);
$M_A$
= 3 in all cases. The interstellar medium flows from right to left. The coordinate system remains the same as before: the X-axis is the stellar rotation axis, and the Y-axis is directed against the flow of the interstellar medium. Panel (a) shows density isolines and streamlines on the Z = 0 plane. Panel (b) is a closer view (to the star) of the same result. Panel (c) shows density isolines and streamlines on the X = 0 plane. Panels (d2) and (d1) show the density dependence on time and its Fourier series at the point (X = 1.5, Y,Z = 0), respectively. Panels (e) and (f) show the density isolines and streamlines on the Z = 0 and X = 0 planes, respectively.
3.3. Moving interstellar medium: 3D model
This section demonstrates that the jet instability persists not only in the case of an interstellar medium at rest but also when the interstellar medium flows around the star. A series of calculations were carried out for different velocities (the Mach numbers,
$M_\infty$
) of the interstellar medium and
$M_A$
= 3,
$ \chi$
= 4. For ease of understanding, we have rotated all the figures in this section so that the interstellar medium flows from right to left. The coordinate system remains the same as before: the X-axis is the stellar rotation axis, and the Y-axis is directed against the flow of the interstellar medium. All results are quasi-stationary, obtained after long-time evolution. Figure 7d2 shows calculations carried out up to a dimensionless time of T=60 (where for parameters close to heliospheric, this corresponds to approximately 7000 years).
Figure 7 shows the results for
$M_\infty$
= 0.6 (Panels a-d) and
$M_\infty$
= 2 (Panels e, f). In the first case (
$M_\infty$
= 0.6), the astropause retains its tube-like shape. In the second case (
$M_\infty$
= 2), bifurcation of the global structure of the flow occurs, and the astropause becomes open in the direction of the tail.
Panel (a) in Figure 7 shows density isolines and streamlines at the Z = 0 plane. Panel (b) gives a zoomed view of Panel (a). It is seen that the astropause is slightly tilted due to the dynamic pressure of the interstellar medium. Similar to section 3.2, jet self-oscillations are seen in Panel (b), particularly in the dense plasma flow behind the termination shock near (X = 1, Y = 0). Panel (c) shows density isolines and streamlines at the X = 0 plane. The flow behind the astropause has an asymmetric and unsteady nature due to the self-oscillations of the jet. A reverse flow is seen behind the tube. This effect has been discussed in Korolkov & Izmodenov (Reference Korolkov and Izmodenov2021) (for example, figure 3, Panels a3, b3).
Panel (d1) in Figure 7 shows a Fourier series expansion of the density at the point (X = 1.5, Y,Z = 0). The X-axis represents the frequency of oscillation, which is equal to one divided by the period of the oscillation. The Y-axis shows the values of coefficients in the Fourier expansion of the density. There is one main period of oscillation that is equal to 1.4. For the interstellar medium at rest, this period is 1.6 (see Figure 6). Therefore, as the velocity of the interstellar medium increases, the period decreases slightly. Panel (d2) shows the density as a function of time at the same point. The magnitude of density varies by about
$\pm14\%$
.
Panels (e) and (f) in Figure 7 illustrate the density isolines and streamlines in the XOY and YOZ planes, respectively. These panels clearly show the global flow pattern, including the open tail of the astropause. Despite this, the influence of the stellar magnetic field remains strong, causing the astropause to be significantly flattened along the Z-axis. Nevertheless, there are no self-oscillations in this case. This is the main conclusion from the results in this case.
4. Conclusion
In this paper, we investigated the global instability of astrospheric jets for stars having strong magnetic fields (
$\gtrsim 3$
times stronger than the Sun). We demonstrate that instability can occur even when the charge exchange process is not considered. The main results of this work are summarised as follows:
-
• For stars with
$M_A \lesssim 4$
, there are self-oscillations of the astrophysical jet (see Figure 4 and “Video-2.avi” in Data Availability Statement). This happens if the interstellar medium is at rest or moves slowly enough so that the tube-like shape of the astropause is still preserved. -
• If the interstellar medium is at rest, the self-oscillation period of a jet decreases as the stellar magnetic field increases. Additionally, the period is inversely related to the
$\chi$
number (
$\sim 1/\chi$
). This means that if the stellar winds have the same dynamic pressure (
$\unicode{x03C1}_1V_1^2=\unicode{x03C1}_2V_2^2$
), but
$\unicode{x03C1}_1 \neq \unicode{x03C1}_2$
, then the self-oscillating periods of their jets are related by:
$T_1 = T_2\dfrac{\chi_1}{\chi_2}$
(see section 2.2). -
• If the interstellar medium has a relative motion with respect to the star and the astrosphere has a tube-like shape, then the period of self-oscillation decreases slightly compared to that of the interstellar medium at rest. The period dependence on
$\chi$
and
$M_A$
has not been established yet, but similar dependencies are expected. -
• The mechanism of jet instability has not been fully understood, but the authors attribute it to the magnetic force acting transversely to the jet. A flow-hindering pressure gradient develops behind the termination shock on the jet’s axis (see Figure 3a,c). Flows with such pressure gradients can be unstable, as shown in a separate study by Korolkov & Izmodenov (Reference Korolkov and Izmodenov2024). On the other hand, the feedback mechanism is clear. It arises from the tension in magnetic field lines that tend to pull the rejected jet back to its original position (see Figure 5), leading to its self-oscillations.
It should be noted that the self-oscillations of the jet have an analogy with very well-studied self-oscillations of a fountain jet Karlikov & Trushina (Reference Karlikov and Trushina1998); Karlikov et al. (Reference Karlikov, Tolokonnikov and Trushina2009). In this context, the pressure gradient along the jet’s axis can be seen as an analogue of the destabilising effect of gravity. At the same time, the reaction of the magnetic field acts as a feedback mechanism similar to the propagation of perturbations in a fluid to the base of a fountain jet.
Acknowledgement
The research is carried out using the equipment of the shared research facilities of HPC Resources of the Higher School of Economics (Kostenetskiy et al. Reference Kostenetskiy, Chulkevich and Kozyrev2021).
Funding statement
SK thanks the Theoretical Physics and Mathematics Advancement Foundation ‘BASIS’ grant 22-8-3-42-1.
Competing interests
None.
Data availability statement
The videos released with this paper have been posted on Zenodo under a Creative Commons Attribution license at DOI: https://doi.org/10.5281/zenodo.14203184.
The other data underlying this article will be shared on reasonable request to the corresponding author.
Open access
