Parametric Study of Magnetosheath Jets in 2D Local Hybrid Simulations

We perform 2D local hybrid simulations of collisionless shocks in order to study the properties of simulated magnetosheath jets as a function of shock properties, namely their Alfvénic Mach number (M A ) and geometry (angle between the upstream magnetic field and the shock normal, θ BN ). In total we perform 15 simulations with inflow speeds of V in = 3.3 C A (Alfvén velocity), 4.5 C A and 5.5 C A and θ BN = 15°, 30°, 45°, 50°, and 65°. Under these conditions, the shock M A varied between 4.28 and 7.42. In order to identify magnetosheath jets in the simulation outputs, we use four different criteria, equivalent to those utilized to identify subsets of magnetosheath jets, called high-speed jets (Plaschke and Hietala and Angelopoulos, Ann. Geophys., 2013, 31, 1877–1889), transient flux enhancements (Němeček et al., Geophys. Res. Lett., 1998, 25, 1273–1276), density plasmoids (Karlsson et al., J. Geophys. Res., 2012, 117, a–n; Karlsson et al., J. Geophys. Res., 2015, 120, 7390–7403) and high-speed plasmoids (Gunell et al., Ann. Geophys., 2014, 32, 991–1009). In our simulations, the density plasmoids were produced only by shocks with M A ≥5.7, while the high-speed plasmoids only formed downstream of shocks with M A ≥6.97. We show that higher M A shocks tend to produce faster jets that tend to have larger surface area, mass, linear momentum and kinetic energy, while these quantities tend to be anticorrelated with θ BN . In general, the increase of θ BN to up to 45° results in increased jet formation rates. In the case of high-speed jets in runs with V in = 3.3 C A and high-speed plasmoids, the jet formation anticorrelates with θ BN . The jet production all but ceases for θ BN = 65° regardless of the shock’s M A . The maximum distances of the magnetosheath jets from the shocks were ≲140 d i (upstream ion intertial lengths), which, estimating 1 d i ∼100–150 km at Earth, corresponds to 2.4–3.3 Earth radii. Thus, none of the simulated jets reached distances equivalent to the average extension of the Earth’s subsolar magnetosheath, which would make them the equivalents of geoeffective jets. Higher M A shocks are probably needed in order to produce such jets.


INTRODUCTION
Magnetosheath jets (see Plaschke et al., 2018, and the references therein) is an umbrella term for different kinds of phenomena discovered in the Earth's magnetosheath (Lucek et al., 2005) that produce localized increases of the solar wind (SW) dynamic pressure (P dyn ). They were first reported by Němeček et al. (1998) as increases of ion flux (product of ion density and velocity) located downstream of the quasi-parallel bow-shock of Earth. The authors named these structures transient flux enhancements (TFEs) and defined them as regions, where the ion flux exceeds the ambient value, averaged over a 20 min time interval, by at least 50%.
Later, different authors used different quantities and thresholds in order to identify and study magnetosheath jets. Here we mention four of them. Karlsson et al. (2012) and Karlsson et al. (2015) used downstream ion density in order to identify density plasmoids. These were defined as structures in which the ion density exceeds the ambient value, obtained by averaging during a 500 s interval, by at least 50%. Karlsson et al. (2012) classified plasmoids into four groups depending on their magnetic field and velocity signatures. Plasmoids were classified as diamagnetic (paramagnetic), if the magnetic field magnitude inside them diminished (increased) and embedded (fast) if the velocity inside them remained the same (increased) compared to the ambient values. It was proposed by Karlsson et al. (2015) that diamagnetic plasmoids may arrive from the SW and then cross into the magnetosheath, while the paramagnetic plasmoids are produced at the bow-shock. One possible cause of paramagnetic density plasmoids discussed by Karlsson et al. (2015) are the short-large amplitude magnetic structures (SLAMS, e.g., Schwartz and Burgess, 1991) transmitted into the magnetosheath. On the other hand, in their local hybrid simulations of collisionless shocks, Preisser et al. (2020a) found that plasmoids may form due to magnetic reconnection just behind the quasiparallel shocks. Plaschke et al. (2013) used the magnetosheath dynamic pressure calculated with the x (radial, along Sun-Earth line) component of the plasma velocity, P dyn,x , in order to study another subset of magnetosheath jets, called high-speed jets. These were defined as regions where three conditions are satisfied: 1. the P dyn,x exceeds the value of 25% of upstream SW dynamic pressure (P dyn,x ≥ 0.25 · P dyn,SW ) during the entire jet interval, 2. the maximum P dyn,x during the same interval must exceed half of the P dyn,SW and 3. the plasma velocity component along Sun-Earth line inside high-speed jets is negative (antisunward). Gunell et al. (2014) studied another subset of magnetosheath jets, called high-speed plasmoids. These were defined as entities that propagate antisunward in the magnetosheath with velocities whose absolute values exceed those in the ambient plasma by a factor of 2 or more.
A possible formation mechanism of high-speed jets, TFEs and high-speed plasmoids was proposed by Hietala et al. (2009). According to them, the SW is processed differently in different locations on the shock's surface due to shock rippling. It was suggested by Hietala et al. (2009) that the SW is decelerated less and is focused downstream of the shock ripples, which gives rise to localized increases of magnetosheath velocity and density and therefore of P dyn . Hietala and Plaschke (2013) suggested that~97% of magnetosheath jets form due to shock rippling. Later global hybrid simulation results by Hao et al. (2016) were consistent with the jet formation mechanism proposed by Hietala et al. (2009).
In the past, large statistical studies of magnetosheath jets have been performed (Plaschke et al., , 2016Archer and Horbury, 2013;Raptis et al., 2020a,b;Liu et al., 2020) in which physical properties of magnetosheath jets, such as their occurrence and sizes were studied. In the latest of such studies, Plaschke et al. (2020) calculated the average jet scale sizes to be 0.12 R E (Earth radius) and 0.15 R E in directions perpendicular and parallel to their direction of propagation. This is about an order of magnitude less than some previous estimations.
Some past investigations have clearly shown the correlation between the IMF orientation and the jet occurrence rate. The statistical analysis of high-speed jets by Plaschke et al. (2013) has shown that these events occur predominantly downstream of quasi-parallel bow-shock, for low IMF cone angles (below~30°). Raptis et al. (2020b) have shown that magnetosheath jets can sometimes also be found downstream of quasi-perpendicular bow-shock. These may have very different origin than jets downstream of quasi-parallel bow shock, as has been discussed by Blanco-Cano et al. (2020) and  (these authors talk about the jets in the quasi-parallel and quasiperpendicular magnetosheath). Raptis et al. (2020b) showed that on average, the jets in the quasi-parallel magnetosheath occur for higher solar wind speeds than jets in the quasiperpendicular magnetosheath.
To our best knowledge, no studies have yet shown any correlation between jet occurrence rate or jet properties, such as their size, mass, speed, etc. and the upstream solar wind and IMF properties (e.g., β, Mach number, etc.) Only  have suggested that the jet occurrence is at most weakly dependent on other solar wind conditions or solar wind variability. Directly relating the upstream solar wind properties and the θ BN angle to the jet properties is a task that should in principle be possible, due to existence of large amount of data gathered by several space missions. This is mostly due to the fact that observationally it is not always possible to determine upstream conditions and shock properties that led to production of individual jets. In order to fill this gap, we perform a parametric study of magnetosheath jet properties as a function of the shock properties. We use the results of 15 local hybrid simulations in which we vary the shock's Alfvénic Mach number, M A (defined as the upstream solar wind speed in the shock's rest frame expressed in units of C A ), and their geometry, determined by the angle θ BN between the upstream magnetic field and the shock normal vector. In the following sections we describe the numeric setup of the simulations and then show how jet properties vary as a function of θ BN , M A and their distance from the shock, Δx.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org February 2022 | Volume 9 | Article 793195 2 NUMERICAL SETUP AND JET IDENTIFICATION CRITERIA We perform local hybrid (kinetic ions, fluid electrons) simulations with the 2D HYPSI code (Burgess and Scholer, 2015;Gingell et al., 2017;Trotta and Burgess, 2019;Preisser et al., 2020a,b). The simulations were performed on a grid of N x × N y = 1,000 × 800 cells with dimensions Δ x = Δ y = 0.5 d i (d i = c/Ω is the upstream ion inertial length) in both directions. The SW was injected from the left-hand, open simulation border along the xaxis with inflow velocities of 3.3 C A , 4.5 C A , and 5.5 C A (initial Alfvén speed). Thus, the x-axis in our simulations is equivalent to the x-axis in the Geocentric Solar Magnetospheric coordinate system. The upstream electron and ion βs (ratio between thermal and magnetic pressures) were 0.5 and the adiabatic coefficient γ was 5/3. Initially, there were 100 particles per cell, to keep the statistical noise at a reasonable level. Periodic boundary conditions were applied along the y-direction, while the right hand side of the simulation acts as a reflecting wall. A shock is therefore initiated with the injection method (Quest, 1985), and propagates along the negative x-direction The upstream magnetic field lied in the xy plane and its orientation was varied so that shocks with θ BN = 15°, 30°, 45°, and 50°were produced. Since in our simulations the shock normals point in the negative x direction, the θ BN is equal to the angle between the upstream magnetic field and the x-axis. We note here that we opted not to simulate exactly parallel shocks with θ BN = 0°. This is because it is very unlikely for a real shock to be exactly parallel. Given that perturbations make θ BN vary, the angle θ BN = 15°is representative for shock geometries spanning, roughly between~5 and~25°. θ BN = 0°would a special case in terms of upstream waves since they produce almost transverse waves, with little compression in the upstream. The case with θ BN = 15°reproduces most effects of almost parallel geometry, and is a better choice than exactly parallel shock. Simulation time was measured in units of the inverse of proton gyrofrequency (Ω −1 ). The simulation timestep was of 0.005 Ω −1 and the outputs were produced every 2.5 Ω −1 . In total, we performed fourteen simulations with initial conditions and the resulting shock Mach numbers (equal to the upstream solar wind velocity in the shock rest frame expressed in units of C A ) that are summarized in Table 1. The physical quantities in this paper, such as velocity and dynamic pressure, are calculated in the shock-rest frame. The running times of the simulations range between and 250-325 Ω −1 .
The search for jets was performed in a rectangular domain with a size of δx × δy = 200 d i ×400 d i with its upstream edge located 12 d i downstream of the shock. By taking into account that the typical d i upstream of the Earth's bow-shock is of the order of 150 km, these dimensions roughly correspond to 4.7 × 9.4 R E . The first condition was chosen because in the outputs from earlier simulation times, the extension of the region downstream of the shocks in the x-direction was only slightly larger than 200 d i . The second condition was used because the simulated shocks are rippled and choosing a domain to close to them meant that occasionally parts of the shock could enter the domain and interfere with the jet identification process. Downstream plasma moments, used in Eqs 3-5, were calculated by averaging over this domain. This is different from observations where averaging in time is used. The initial conditions are used as a proxy for upstream plasma moments in Eqs 1, 2. For each run, we selected 21 consecutive outputs with well developed shocks. In order to identify the jets we use four different criteria, equivalent of those used by Němeček et al. (1998) for TFEs, Plaschke et al. (2013) for high-speed jets, Karlsson et al. (2012) for density plasmoids and (Gunell et al., 2014) for high-speed plasmoids. Some of these works use timeaveraged downstream plasma quantities in order to define thresholds for magnetosheath jets. For example, Němeček et al. (1998) and Karlsson et al. (2012) avereage the ion flux density and magnetosheath plasma density during 20-min and 500 s time intervals, respectively. However for our simulations, such time averaging would require simulation outputs with very high time cadence. Hence, we opt for spatial averaging of the required downstream plasma properties in the described rectangular domain. As for the upstream quantities, we use initial plasma and magnetic field values.
In the case of high-speed jets, the threshold is composed of two conditions: where the subscripts d and u refer to downstream and upstream regions, respectively, the subscript 0 means initial values at t = 0, max refers to the maximum value inside high-speed jets and x refers to velocity component along the x-axis. P dyn,x,d , n and v stand for the dynamic pressure, ion number density and ion velocity. For TFEs the threshold writes as for density plasmoids it is while for high-speed plasmoids it is expressed as The triangular brackets denote quantities that were spatially averaged in the domain described in the previous paragraph. Additionally, in order to be selected, the jets had to have surface areas larger than a threshold, determined separately for each type of jets. This was done so that the selected events are physically meaningful structures and not just background noise. After careful inspection, the appropriate minimum areas for highspeed jets, TFEs, density plasmoids and high-speed plasmoids were decided to be 25 d 2 i (100 cells), 7.5 d 2 i (30 cells), 2 d 2 i (8 cells) and 2 d 2 i (8 cells), respectively. These areas have been chosen so that the selected magnetosheath jets exhibit P dyn , ion flux, density and velocity signatures that stand out from the rest of the downstream region as opposed to being part of downstream fluctuations that populate the whole downstream region.

RESULTS
In Figure 1 we show an example of the output from the model with θ BN = 15°and V in = 3.3 C A at t = 200 Ω −1 . The colors represent the magnetic field magnitude. Shock is located at x~300 d i and can be distinguished by red trace (high B-field magnitude). We can see that the upstream and the downstream regions are highly turbulent.
In Figure 2 we show examples of each type of event studied here. These events were selected for illustration purposes from different runs. In general, some of the jets were found by more than one criteria, but there are also many events that have been identified by one criteria only. We did not look into exactly how often different jet definitions overlap, but the survey performed by Archer and Horbury (2013) suggests that this is something that commonly occurs. The colorscales on panels a), b), c) and d) represent the P dyn,x , F = ρ d v d , number density and velocity along x-axis, respectively. The black contours on these panels outline a high-speed jet, a TFE, a density plasmoid and a high-speed plasmoid. Below we examine the behaviour of the statistical properties of magnetosheath jets as a function of M A , θ BN and the distance from the shock, Δx.

High-Speed Jets
In Figure 3 we show statistical properties of high-speed jets as a function of θ BN and the inflow velocity of the runs. On panels a-f we plot the average number of jets, their average area, mass, linear momentum, kinetic energy and center of mass velocity. The latter is normalized to the initial upstream C A . The blue, the red and the FIGURE 1 | An example of the output from the model with θ BN = 15°and V in = 3.3 C A at t = 200 Ω −1 . The colors represent the magnetic field magnitude, B tot in units of initial upstream magnetic field strengths, B 0 . Shock is located at x~300 d i and can be distinguished by red trace (high B-field magnitude). The upstream plasma is propagating from the left along the x direction and is specularly reflected at the right boundary. Periodic boundary contidions are applied at the bottom and top boundaries. We can see in Figure 3A) that the average number of jets increases with V in and with θ BN for values of θ BN ≤ 45°. On the panels b-f we observe that the average jet properties behave differently as a function of θ BN for different V in . Area, mass, momentum and kinetic energy of jets produced in runs with V in = 3.3 C A , 4.5 C A , and 5.5 C A peak at θ BN = 50°, 30°, and 15°, respectively. The exception is the kinetic energy of jets in models with V in ≤ 4.5 C A that increases monotonically for θ BN ≤ 50°. In the case of the runs with V in = 5.5 C A , these properties exhibit another, smaller peak at θ BN = 45°. The center of mass velocity of the jets shows very weak dependence on θ BN for ≤50°.

Transient Flux Enhancements
We show the properties of TFEs in the Figure 4. As in the previous case, the average number of TFEs increases with V in for θ BN ≤ 45°and V in ≥4.5. The number of TFEs in runs with V in = 3.3 C A diminishes monotonically with θ BN . The area, mass, momentum and kinetic energy of TFEs in runs with V in = 5.5 C A peak at θ BN = 15°with a second peak at 45°. In runs with V in = 4.5 C A these quantities present a peak at θ BN = 30°. Finally, in runs with V in = 3.3 C A these quantities exhibit largest values at θ BN = 15°and 50°except for the mass which peaks at θ BN = 45°. The center of mass velocity of TFEs positively correlates with V in and with θ BN (panel f). Figure 5 shows statistical properties of density plasmoids. We can see that no plasmoids were produced in runs with V in = 3.3 C A . For θ BN ≤ 45°their average number increases with θ BN . Plasmoid properties (area, mass, momentum, kinetic energy) show complex behaviour with θ BN . In the case of plasmoids produced in runs with V in = 5.5 C A , they peak at θ BN = 45°a nd exhibit a smaller peak for θ BN = 15°. Their center of mass velocity peaks at θ BN = 15°and exhibits a second peak at θ BN = 45°. Values of the properties of the jets produced in runs with V in = 4.5 C A exhibit approximately monotonic growth at θ BN ≤ 50°and then decline for larger angles. Figure 6 shows properties of high-speed plasmoids. These were produced only in runs with V in = 5.5 C A and only for θ BN = 15°, 30°and 45°. We can see that their number, area, mass, momentum, and kinetic energy decline with θ BN , while their center of mass velocity peaks at θ BN = 30°.

High-Speed Plasmoids
Finally we comment on the average number of different types of events, which we consider as a proxy of jet formation rates, by looking at panels 3a-6a. We can see that the most numerous events, on average, were TFEs, followed by high-speed jets, density plasmoids and high-speed plasmoids. Although density plasmoids in the run with V in = 5.5 C A were on average more numerous than the high-speed jets in the same runs, their numbers were smaller in runs with V in = 4.5 C A and they were not produced in runs with V in = 3.3 C A .

Properties of Magnetosheath Jets as a Function of Their Distance From the Shock
In this section we discuss the evolution of jets. Ideally we would study the evolution of average jet properties as a function of time. However, there are several reasons for why this is not always possible. For one, we would have to define the exact moment of formation for every jet and it is not clear how to do this. Also, our time between the outputs is 2.5 Ω −1 , which limits our ability to track the jet formation. Thus, we choose to study the average jet  Figure 3. The blue, the red and the black traces represent properties of jets produced in runs with V in = 3.3 C A , 4.5 C A and 5.5 C A respectively. The units are given in brackets. The subindex u refers to the upstream region, while the subindex 0 means the initial value.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org February 2022 | Volume 9 | Article 793195 properties as a function of their distance from the shock. Here we make two suppositions: that all the jets are formed at the shock and not further downstream and that their average distance is proportional to their age. Since this distance is different for jets with different velocities, in Figures 7-10 we group the jets according to the V in of the corresponding runs, so that jet velocities in each Figure do not vary a lot.

High-Speed Jets
Figure 7 exhibits scatter plots of jet properties as a function of their distance from the shock, Δx. The red triangles, blue squares, pink asterisks and black squares represent jets produced in runs with θ BN = 15°, 30°, 45°, and 50°, respectively. Figure 6A-C show results from the runs with V in = 3.3 C A , 4.5 C A , and 5.5 C A , respectively. All the shown quantities are normalized to the highest value on each panel. Panels i-v show plots of normalized area, mass, momentum, kinetic energy and center of mass velocity. Figures 8-10 have the same format.
The first thing we see on all panels in Figure 7 is that the number of high-speed jets rapidly decreases as a function of Δx. On velocity panels v we see that jets at largest distances from the shock have velocities that are somewhere in the middle between the largest and the smallest velocity near the shock. This is due to the fact that the velocity of the jets diminishes with the distance because the jets exchange their momentum with the surrounding plasma and thus decelerate. The minimum jet velocity is either maintained or increased slightly with Δx. This happens because the jets that start with small velocities decelerate further downstream and no longer fulfill the established threshold needed for their identification. Also, on all panels v we can see that average jet velocity positively correlates with θ BN , as has been seen in previous Figures.
In the Figure 7A we can observe that there are six jets produced in runs with θ BN = 50°that stand out as they exhibit values of normalized area, mass, momentum and kinetic energy above 0.5. These jets are observed at distances of up to 100 d i . The  Figure 3. The red and the black traces represent properties of jets produced in runs with 4.5 C A and 5.5 C A respectively. The units are given in brackets. The subindex u refers to the upstream region, while the subindex 0 means the initial value.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org February 2022 | Volume 9 | Article 793195 rest of the jets all occupy roughly the domain with the normalized values of ≤0.5. In the Figure 7B there are no outliers, while in the Figure 7C there are three jets, produced in runs with θ BN = 15°, 30°, and 45°, that stand out from the rest, however this time they only appear at small Δx. It can be seen on all panels that the maximum Δx was reached by the jets produced in runs with θ BN = 45°and 50°. The maximum distance is~130-140 d i . Among the jets from the run with V in = 5.5 C A there is only one jet at distances larger than 100 d i .
For a bulk of the jets, maximum values of their area, mass, momentum and energy tend to diminish monotonically with Δx, the only exception being the six outliers in Figure 7A.

Transient Flux Enhancements
Scatter plots of TFE properties are shown in Figures 8A-C. We see that these events are more numerous than high-speed jets. The maximum distance at which they were found arẽ 120 d i for those from runs with V in = 3.3 C A and~140 d i for the other runs. For TFEs in runs with V in = 4.5 C A and 5.5 C A , the maximum values of their properties (area, mass, momentum, kinetic energy and velocity) diminish with Δx. For jets form runs with V in = 3.3 C A , this is true for the majority of events that start out with normalized values of ≲0.5, while for those events that exhibit higher initial values, monotonic decrease is not obvious.

Density Plasmoids
In Figure 9 we show the properties of density plasmoids. These were not produced in runs with V in = 3.3 C A . Something that is very obvious is that the density plasmoids reached maximum distances of~70 d i and 100 d i in runs with V in = 4.5 C A and 5.5 C A , respectively. Thus, not only are these events more difficult to produce, they merge with the surrounding plasma much faster than the high-speed jets and TFEs. In the case of the V in = 4.5 C A , it was the events in runs with θ BN = 50°that reached the highest normalized values of area, mass, momentum and kinetic energy, while in runs with V in = 5.5 C A , it was the plasmoids from runs with θ BN = 45°.

High-Speed Plasmoids
The properties of high-speed plasmoids are shown in Figure 10. These were only produced in runs with V in = 5.5 C A and θ BN = 15°, 30°, and 45°. High-speed plasmoids survived to maximum distances of~50 d i , less than any other type of events. Thus, they are difficult to form and they quickly merge into the ambient plasma.

DISCUSSION
In this work we perform a parametric study of magnetosheath jets that formed in fifteen 2D local hybrid simulations of collisionless shocks carried out with the HYPSI numerical code.
Two parameters were varied in the simulations: the inflow velocity, V in , and the angle θ BN . Both variations result in variations of the shock's M A . In order to identify jets in the simulations, we applied four different criteria, equivalent to those used for identifying four different subsets of magnetosheath jets: high-speed jets , transient flux enhancements (Němeček et al., 1998), density plasmoids (Karlsson et al., 2012) and high-speed plasmoids (Gunell et al., 2014). Some past investigations (e.g., Plaschke et al., 2013;Raptis et al., 2020b) have clearly shown that the upstream IMF cone angle and θ BN impact the jet occurrence rate. Higher rates occur during smaller IMF cone angles (≲30°) and θ BN . However clear correlation between the jet occurence and the SW parameters, such as its upstream β, M A and the IMF cone angle, has not yet been established (e.g., Plaschke et al., 2018). It has been suggested by Plaschke et al. (2013) that the jet occurrence is only weakly dependent on other solar wind conditions or solar wind variability, if at all. Raptis et al. (2020b) have shown that a minority of the magnetosheath jets can also be found downstream of the quasi-perpendicular bow-shock and that they are on average related to lower solar wind velocities than the jets observed downstream of the quasi-parallel bow-shock. It has been shown by Blanco-Cano et al. (2020) and  that the jets downstream of the quasi-perpendicular bowshock may have very different origin. As such, they are not produced in our simulations. There has not yet been any investigation that would to relate the upstream solar wind properties and the θ BN angle to the jet properties, such as their size, mass, speed, momentum, kinetic energy etc. as we have done here. Although local simulations allow us to study magnetosheath jet properties as a function of selected shock parameters, they do have some limitations. Compared to the global simulations and real bow-shocks, our shocks and foreshocks are much simpler. In our simulations, the upstream region is quite limited in size, which results in a relatively small number of upstream ULF waves (between 5 and 20, depending on the model, see . These waves also exist for shorter time intervals, which may lead to less wave steepening and smaller number of magnetic structures such as shocklets and SLAMS . This could result in smaller number of magnetosheath jets. Also, in real bow-shocks and in global simulations, the upstream ULF waves are a mix of waves with different properties produced upstream of the sections of the bow-shock that have different nominal geometries and Mach numbers. Another thing is that in our simulations the upstream solar wind always impinges perpendicularly to the shock surface, which, in the case of the bow-shocks is only true close to the subsolar point. All these facts may influence the properties of the shock ripples and consequently of magnetosheath jets. Additionally, draping of the IMF lines downstream may cause the regions deep in the magnetosheath to be magnetically connected to the regions of the FIGURE 7 | Scatter plots of properties of high-speed jets as a function of their distance from the shock. Panels (A-C) show results for runs with V in = 3.3 C A , 4.5 C A and 5.5 C A . Panels (i-v) show normalized area, mass, linear momentum, kinetic energy and velocity. All the quantities are normalized to the highest value on each panel. Red triangles, blue squares, purple asterisks and black dots represent jets from runs with θ BN = 15°, 30°, 45°, and 50°, respectively.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org February 2022 | Volume 9 | Article 793195 bow-shock that are much farther away and have different geometries than is the case in our simulations. The SW flow in the magnetosheath also progressively deviates from the Sun-Earth line as one approaches the magnetopause, which may affect the trajectories of the magnetosheath jets. Our finding may thus be more representative of the regions just behind bow-shocks and less so of the regions closer to the magnetopause. For each run, we analyzed 21 consecutive outputs. We study jet properties, such as their area, mass, linear momentum, kinetic energy and center of mass velocity, first as a function of the shock's M A and θ BN and then as a function of their distance from the shock, Δx. The latter is chosen as a proxy of the jet's age in order to study their evolution.
We can see in Figures 3A-6A that in our simulations the number of TFEs was the largest and the number of high-speed plasmoids the smallest among all fours types of events. Although the number of identified events to some extent depends on the minimum surface area chosen for each subset of magnetosheath jets, a visual inspection showed that this does not significantly change the relative occurrences of different types of magnetosheath jets.
We can see in Figures 5F, 6 that density plasmoids were produced at shocks with M A ≥5.7, while high-speed plasmoids were only produced at shocks with M A ≥6.97. It thus seems that these events require higher M A shocks for their formation than high-speed jets and TFEs.
The number of jets in general increases with M A for θ BN values ≤45°. The exceptions are the TFEs in runs with C A = 3.3 C A and high-speed plasmoids, whose number decreases with θ BN . Considering the fact that increasing θ BN results in higher M A , we thus conclude that, in general, higher M A favors the formation rates of magnetosheath jets. However, for some critical θ BN,c > 50°, the jet production all but ceases regardless of the shock's M A .
On the other hand, the magnetosheath jet area, mass, linear momentum and kinetic energy tend to diminish with increasing θ BN , although they tend to be larger in runs with larger V in . Thus, the M A of the shocks plays an important role in jet properties-larger M A mean higher values of the mentioned jet properties, but their anticorrelation with θ BN is strong also. The velocity of magnetosheath jets tends to be positively correlated with M A , except in the case of high-speed jets in runs with θ BN ≤ 45°, where there is almost no dependence of the jet center of mass velocity on θ BN . Positive correlation of jet velocity on θ BN is strong for TFEs, while in the case of both type of plasmoids, the dependence of the velocity on θ BN is complex.
If the jets are produced by the mechanisms proposed by Hietala et al. (2009) and Karlsson et al. (2012), which has been claimed by, for example Plaschke et al., 2013), it is reasonable to expect the jet properties to depend on the properties of shock ripples, which in turn have been shown by  to depend on the properties of upstream ULF waves. These authors observed a complex dependence of the upstream ULF wave amplitude and wavelength on M A and θ BN and showed that it correlates strongly with increasing θ BN and weakly with M A (see their Figure 6).
The behaviour of the amplitude of the upstream ULF waves was explained in terms of the theory of generation of upstream waves proposed by Barnes (1970) and numerical simulation results obtained by Burgess (1987). Barnes (1970) showed that the amplitude of upstream ULF waves is proportional to the product of the beam velocity and density: Burgess (1987) on the other hand showed that 1) the fieldaligned beam velocity increases and 2) its density decreases strongly with increasing θ BN for 40°< θ BN < 60°.  suggested that the interplay between V b and N b is probably the reason why the amplitudes of the upstream waves exhibit a complicated dependence on the shock's M A and θ BN .
It was concluded by these authors that the properties of the shock ripples are influenced by the properties of the upstream ULF waves and thus they themselves exhibit a complex dependence on M A and θ BN .  divided the simulated shocks into weakly, intermediately and strongly rippled. Among the strongly rippled shocks were those with.
The shocks classified as intermediately rippled were those from models with.
The rest were classified as weakly rippled. By comparing this with Figures 3-6, there does not seem to be any obvious correlation between the degree to which the shocks are rippled FIGURE 9 | Scatter plots of properties of density plasmoids as a function of their distance from the shock. The format is the same as in Figure 7 except that runs with V in = 3.3 C A are missing since no density plasmoids were formed in them.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org February 2022 | Volume 9 | Article 793195 and the simulated jet properties. The detailed analysis of correlation beteween the shock ripples and magnetosheath jets will be the topic of future investigations.
Here we make one observation about the upstream ULF wave wavelength (λ) and jet properties.  showed that the λ of upstream ULF waves increases strongly with θ BN , however the values of different jet properties (area, mass, momentum, energy and velocity) in most cases show the opposite trend, as they tend to diminish with θ BN . The exceptions are, as can be seen in Figure 3A)-6A), density plasmoids with V in = 4.5 C A .  also showed that the ULF wave λ positively correlates with M A , similar to average values of jet properties for most jets (Figures 3-6).
The properties of different types of jets also evolve with a distance from the shock, Δx. In general, the maximum values of the jet properties diminish more or less monotonically with Δx. The exception to the rule are a small number of simulated highspeed jets and TFEs produced in runs with V in = 3.3 C A (see Figures 3A, 4A).
In our simulations, high-speed jet and TFEs reached maximum distances of~130-140 d i . Density plasmoids were identified up to~100 d i from the shock, while maximum distances for high-speed plasmoids were~50 d i . In local hybrid simulations performed by Preisser et al. (2020a), the density plasmoids survived to large distances from the shock, because inside them the plasma is magnetically isolated from the ambient magnetosheath plasma. This could mean that in our simulations these events also exist at large distances from the shock, however they no longer fulfill the corresponding identification criteria and are thus not detected anymore.
We can compare these distances to those at the bow-shock of Earth. At 1 A.U. typical values for ion inertal length in the pristine SW are between 100 d i and 150 d i . Thus, 140 d i corresponds to 1.4 ·10 4 km-2.1 ·10 4 km, which can be expressed as 2.4 R E -3.3 R E . This is less than 4 R E , which is a typical width of the magnetosheath along the Sun-Earth line. This means that the simulated jets do not reach sufficient distances from the shock so that they would be equivalent to the so called geoeffective jets, which hit the magnetopause. Perhaps simulations with higher M A which would produce faster jets, and/or with longer simulation times are needed in order to produce such magnetosheath jets.

CONCLUSION
In this section we briefly summarize the results of this investigation in which we study the properties of magnetosheath jets as a function of upstream solar wind and IMF conditions in fourteen local hybrid simulations of collisionless shocks.
• In the case of shocks produced with V in = 4.5 C A and 5.5 C A , the magnetosheath jet formation rate positively correlates with θ BN (and thus with M A ) for θ BN ≤ 45°. For larger θ BN it diminishes strongly. The exceptions are the high-speed plasmoids and high-speed jets in runs with V in = 3.3 C A , whose formation rate decreases with increasing θ BN . • The average jet properties, such as their area, mass, linear momentum kinetic energy and velocity tend to increase with shock's M A and diminish with increasing θ BN . However the latter dependence is often not monotonic. • There seems to exist a critical angle θ BN between 50°and 65°w hen jets are no longer produced regardless of the shock's M A . This may imply that magnetosheath jets detected downstream of the quasi-perpendicular shocks could be of different origin than that proposed by Hietala et al. (2009). This result is in agreement with past observational works, such as Plaschke et al. (2013), Raptis et al. (2020b) and others who show that the magnetosheath jets mostly occur for low θ BN angles. • The number of identified jets in our simulations is different for different subsets. It is the largest for TFEs, followed by high-speed jets, density plasmoids and high-speed plasmoids. • Density plasmoids were not produced in runs with V in = 3.3 C A , while high-speed plasmoids were not produced in runs FIGURE 10 | Scatter plots of properties of high-speed plasmoids as a function of their distance from the shock. The format is the same as in Figure 7. The format is the same as in Figure 7 except that runs with V in = 3.3 C A and 4.5 C A are missing since no high-speed plasmoids were formed in them.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org February 2022 | Volume 9 | Article 793195 with V in = 3.3 C A and 4.5 C A . It is therefore expected that these subsets of magnetosheath jets will be more difficult to detect, especially downstream of low M A shocks. • Smaller number of density plasmoids could also be partially due to the fact that due to relatively small simulation domain, the size of the upstream region is limited, so that the ULF waves undergo less steepening, which in turn may result in fewer SLAMS formed in the foreshock. • The number of the magnetosheath jets diminishes with the distance from the shocks, Δx. This is because the velocity of the jets diminishes with the distance form the shock, so they no longer comply with the established criteria. • The events that were found at largest Δx from the shocks are the high-speed jets and the TFEs. No jet was found beyond 140 d i . By taking typical values of upstream d i at 1 A.U. to be in the range of 100-150 km, 140 d i correspond to 2.4 R E -3.3 R E . This is much less than the typical width of the magnetosheath along the Sun-Earth line, which is about 4 R E .
• It thus follows from our simulations that of the four types of magnetosheath jets studied here, only high-speed jets and TFEs could eventualy hit the magnetopause and become geoeffective. • Runs with M A > 8 are needed in order to produce geoeffective magnetosheath jets.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
A.-A-analysis and interpretation of simulation results. PK-analysis and interpretation of simulation results and supevision (principal adivser). LP-analysis and interpretation of simulation results. DT-analysis and interpretation of simulation results. DB-author of the hypsi code, analysis and interpretation of simulation results.