Magnetohydrodynamic Simulations of Spicular Jet Propagation Applied to Lower Solar Atmosphere Model

We report a series of numerical experiments for the propagation of a momentum pulse representing a chromospheric jet, simulated using an idealized magnetohydrodynamic model. The jet in a stratified lower solar atmosphere is subjected to a varied initial driver (amplitude, period) and magnetic field conditions to examine the parameter influence over jet morphology and kinematics. The simulated jet captured key observed spicule characteristics including maximum heights, field-aligned mass motions/trajectories, and cross-sectional width deformations. Next, the jet features also show a prominent bright, bulb-like apex, similar to reported observed chromospheric jets, formed due to the higher density of plasma and/or waves. Furthermore, the simulations highlight the presence of not yet observed internal crisscross/knots substructures generated by shock waves reflected within the jet structure. Therefore we suggest verifying these predicted fine-scale structures in highly localized lower solar atmospheric jets, e.g., in spicules or fibrils by high-resolution observations, offered by the Daniel K. Inoyue Solar Telescope or otherwise.


Introduction
Understanding the physical mechanisms that maintain the observed thermal profile at the solar interface region between the photosphere and the low corona remains an intriguing problem in solar physics. This is where the plasma pressuredominated lower solar atmosphere (photosphere and chromosphere) becomes magnetic pressure-dominated (low corona). Observations suggest a ubiquitous thin magnetic flux tube (MFT) features in the region that topologically links the photosphere with the corona and transport mass, energy, and momentum across this transitional zone balance the radiative losses. Moreover, these MFT structures act as conduits for magnetohydrodynamic (MHD) wave modes to propagate from the lower solar atmosphere to the corona and transfer energy between these regions.
Depending on the observed wavelengths, physical characteristics (e.g., length, lifetime, velocity), regions (quiet Sun, active region, or coronal hole) and/or locations (at limb or on disk), the MFT structures are further categorized as different features (e.g., spicules, mottles, dynamic fibrils, macrospicules, x-ray jets, EUV jets, coronal jets) that permeate the solar atmosphere (see reviews by Beckers 1968Beckers , 1972Tsiropoula et al. 2012). These MFT structures are key to understanding the wave and pulse propagation in the solar atmosphere with their observed dynamics (transverse, rotational, width perturbations) reflecting the presence of confined MHD wave modes (sausage, kink, torsional Alfvén, fluting).
In recent years, the field-aligned/longitudinal propagation of dominant MHD waves/pulses in these MFT structures become a key focus of research to understand chromospheric heating (Narain & Ulmschneider 1990;Zaqarashvili & Erdélyi 2009;Jess et al. 2015). Understanding wave/pulse propagation in a dynamic, gravitationally strongly stratified, and inhomogeneous environment can provide vital clues about the dissipation mechanisms and/or formation of shock-like behavior associated with these thin MFT structures. Observations suggest upward flow velocities in the range of 15-110 km s −1 for typical lifetimes of around 50-400 s for off-limb spicular structures. The visible apex of these observed features attains a maximum height of ∼7 Mm and descends back on the surface following a parabolic trajectory (Pereira et al. 2012(Pereira et al. , 2016. Several theoretical and numerical approaches to gain an insight into the morphology and dynamical properties of jets in the solar atmosphere led to a plethora of models (early review by Sterling 2000) that included physical mechanisms such as pressure/velocity pulse originating in the middle-upper chromosphere leading to the formation of shocks (Hollweg 1982;Shibata et al. 1982;Suematsu et al. 1982;Sterling & Mariska 1990;Heggland et al. 2007;Kuźma et al. 2017), Alfvén waves in an MFT waveguide that nonlinearly couples and generates shocks Hollweg 1992;Kudoh & Shibata 1999;Matsumoto & Shibata 2010), magnetic reconnection driven jets (Yokoyama & Shibata 1995, 1996Archontis et al. 2005;Isobe et al. 2008;Nishizuka et al. 2008;Sterling et al. 2010;González-Avilés et al. 2017, and/or Joule heating due to ion-neutral collision damping (Haerendel 1992;James et al. 2003). More realistic 3D simulations with complex physics mimicking the lower solar atmosphere were employed by  using the Bifrost code by Gudiksen et al. 2011). By combining a multitude of physical mechanisms that involve radiative losses, thermal conduction along the magnetic field, and ionneutral and nonlocal thermal equilibrium effects, the potential role of spicular jets in heating at the solar transition zone was highlighted.
Despite the numerous observational, theoretical, and numerical studies used to understand the morphology and formation of jets in the lower solar atmosphere, qualitative investigation on parameters that influence the observed dynamical characteristics of these features remains somewhat missing. Here, we aim to address the fundamental problem regarding the propagation of a momentum pulse/wave originating at the top of the photosphere and propagating through the chromosphere and the transition region (TR) into the lower corona in an idealized/stratified solar atmosphere using a simplistic numerical model. The parameter space comprises of driver times, amplitudes, and magnetic field strength, which are examined to assess their role in determining the height, width, internal (sub)structures, and dominant MHD wave mode.
The paper is organized as follows. Section 2 describes the numerical methods used to simulate the chromospheric jet, Section 3 discusses the parameter scan on the jets simulations regarding their morphology and dynamics, along with the resulting substructures and possible energy/mass transfer mechanisms, and Section 4summarizes the key findings of our investigation and their implications, together with possible scope for future studies.

Magnetohydrodynamic Simulations
Our computational model for the numerical simulations is based on a 2D setting, using the grid-adaptive MPI-AMRVAC version 2.0 software (Tóth 1996;Keppens et al. 2012;Porth et al. 2014;Xia et al. 2018). The governing MHD equations are solved in the following form (including gravity as an additional source term), Above is the full system of MHD equations employed in our model with variables denoted as: time (t), mass density (ρ), plasma velocity v, and magnetic field (B). The total pressure (p tot ) and energy density (e) terms are given as, where γ = 5/3 denotes the ratio of specific heats. Solar gravitational acceleration is given by, where G is the gravitational constant and M e denotes the solar mass.

Equilibrium
In our model, the TR layer is located at 2 Mm, where the temperature connects the 8000 K, upper chromosphere to a 1.5 MK corona, as shown in Figure 1. The temperature in this region as a function of height is mathematically represented as, where y corresponds to the vertical direction, T 0 (0) = T ch = 8 × 10 3 K (T c = 1.8 × 10 6 K) is the chromospheric (coronal) temperature, and y 2 tr = Mm (w 0.02 tr = Mm) is the TR height (width). At t = 0 the total pressure varies only in the vertical direction. The magnetohydrostatic equilibrium is achieved with, where j is the current density. Since the magnetic field is uniform in the vertical direction, the Lorentz force (j × B) is zero. The temperature and pressure are related through the ideal gas law, where M is the mean atomic weight, T is the temperature, and  is the universal gas constant. Therefore, the relation for ρ and p can be written in the following form, ρ 0 and p 0 are the initial mass density and pressure equilibrium, respectively, at the lower boundary of the computational domain, while the pressure scale heights are given by, Using the temperature profile given by Equation (8) and taking ρ 0 = 2.34 × 10 −4 kg m −3 , we create pressure and mass density as a function of height by means of Equations (11) and (12).
To simulate chromospheric jets in a lower solar atmospheric environment, we set the domain size to 50 × 30 Mm with a level one resolution of 32 × 24 (giving a physical resolution of ∼1.56 × 1.25 Mm). In addition, unidirectional grid stretching is applied to the vertical direction where, from the bottom/ origin, the grid cells vary by a constant factor of 1.1 from cell to cell. The combination of coarse level one resolution and grid stretching allowed us to investigate the entire lifetime of the simulated jet(s) with strong signals traveling upwards becoming increasingly diffusive. To accurately capture the sharpness of the TR, we applied seven levels of AMR giving a spatial resolution of ∼12 × 10 km. The advantage of using AMR is that extra grid cells are allocated dynamically where smallscale features occur. Therefore one may cost-effectively resolve steep gradients such as the TR or the boundary of a jet as it quickly evolves in the computational domain. Due to the choice of the discretization scheme for the boundary conditions, we employed ghost cells of two grid layers encasing the physical domain in each direction. The numerical scheme applied for spatial discretization is the Harten-Lax-van Leer solver (Harten et al. 1983) and a third-order Čada limiter (Čada and Torrilhon 2009) for the time discretization. We used a CFL number of 0.8 and the generalized Lagrange multiplier (GLM)-MHD method to maintain ∇ · B = 0 (Dedner et al. 2002). For both the left and right boundaries, we utilized periodic boundary conditions. In ghost cells at the lower boundary of the computational domain, we set ρ, e, and B to their initial values. In ghost cells for the upper boundary, the values for ρ and e are determined by the gravitational stratification and B was extrapolated assuming zero as the normal gradient. For both the upper and lower boundaries, we took an antisymmetric boundary for the velocity components.

Driver
The jet structure in the series of our simulations is initiated using a momentum pulse in the ghost cells of the computational domain for a specified time (ranging from 50-300 s). The jet is launched symmetrically by a driver at the center of the computational domain, which varies both spatially and temporally as shown in Figure 1. In the x direction, the vertical component of jet velocity has a Gaussian profile with the full width at half maximum (FWHM) of 187.5 km (w j ). As the jet evolves in time, this function will attain a switch-off phase where it extinguishes with a hyperbolic tangent, where v j is the velocity of the jet, A is the amplitude of the driver, t is time, P duration of the driver time, and x 0 is central location of the jet initiation. The width of the driver is determined by Δx, which is the FWHM based on the jet width (w j ) (see Figure 1),

Parameter Scan
Using the numerical setup described above, we conducted multiple simulations to determine the role of three key parameters (magnetic field strength, driver period, and amplitude) on jet characteristics in the lower solar atmosphere. All possible combinations with respective magnitudes B = 20, 40, 60, and 80 G, P = 50, 200, and 300 s, and A = 20, 40, 60, and 80 km s −1 were investigated to qualitatively estimate their impact on jet morphology (maximum height and width) and kinematics (trajectory and transverse displacement). In our simulations, the choice of initial velocities are based on observed speeds of spicules (∼15-120 km s −1 ) and the momentum needed to lift material from the low solar atmosphere to spicular heights (∼4-8 Mm). The driver periods coincide with 5 minute photospheric oscillations due to the pmode (Leighton et al. 1962), which has already been shown as a potential driver for chromospheric jets (De Pontieu et al. 2004). The magnetic field estimates are consistent with stretched computational grid at t = 0 s with seven levels of AMR. The panel at the top right shows the initial temperature (red) and density (blue) stratification for the lower solar atmosphere (0-10 Mm). In the lower panels (left), an example of the Gaussian function, with A = 60 km s −1 and w j = 187.5 km (marked by red dots) at t = 0 s, is shown, which is used as the driver for jets in our simulations. The driver velocities for amplitude (A = 60 km s −1 ) is shown for periods in the range P = 50-300 s in the lower right panel.
For accurate identification of parameter influence on jet morphology and kinematics, we developed a jet tracking software ( Figure 2). The software tracks the jet apex (yellow triangle) and the cross-sectional widths (blue circular dots) in both spatial and temporal domains, using a tracer parameter that tracts the advection of selected fluid in the computational box (Porth et al. 2014). We assign a high arbitrary value (∼100) to the jet material, while the ambient environment is set to zero. By using a threshold of 15% of the jet material, the simulation is converted into a binary image with one identifying cell associated with the numerical jet and zero to the ambient medium. The cross-sectional widths are estimated as a function of the distance between the opposite jet cells at edges, taken from a horizontal slice across the simulated jet at a specified height(s). Furthermore, the visible apex of the jet is selected as the highest index of a jet cell in the simulation.
The visible apex of the simulated jet indicates a parabolic motion over time, as shown in Figure 3, with each panel highlighting the dynamics based on the initial velocity of the driver. The estimated jet trajectories reveal the nonballistic path (s) of the apex, which is consistent with reported observations (Hansteen et al. 2006;De Pontieu et al. 2007;Rouppe van der Voort et al. 2007) for chromospheric jets (e.g., mottles, fibrils, spicules). In our numerical setup, the minimum amplitude for the driver is estimated to be >20 km s −1 for the nearphotospheric plasma to attain spicular heights. Also, for most cases, the maximum apex for the simulated jets was found to be proportional to the driver amplitude.
The influence of key parameters on the jet apex (panels to the left) and mean jet cross-sectional widths (panels to the right) are displayed in Figure 4. Here, the magnetic field strength (aligned with jet injection) has a minimal impact on the field-aligned flow; however, it must be noted that any deviation between magnetic field direction and flow can severely influence the apex height attained by the jet. Interestingly, the driver periods influence the maximum height of the jet for a shorter duration (50 s); though, any increase over 200 s has no/minimal impact over the jet height(s). In our study, we found that the main parameter that strongly affects the maximum height of the jet apex is the driver amplitude. The apex height is determined by the magnitude of the driver amplitude with the jet apex following a parabolic trajectory. Based on these indications, we modeled the jet height with a power-law function given as h Cv . 1 6 j n max

( ) =
By taking an average value for each velocity data point and then fitting an optimal curve using least squares, we estimated values of C = 10 −2.21 and n = 1.72, as shown in the lower (left) panel in Figure 4. The estimates from the power-law approximation suggest a nonlinear relation between apex height and pulse strength, in contrast to those reported by, e.g., Singh et al. (2019). It must, however, be noted that Singh et al. (2019) used an instant pressure pulse close to the TR as a driver in their simulation, which might not be directly comparable with our investigation.
The panels (right) in Figure 4 show effects of different parameter combinations on cross-sectional widths of the simulated jet structure. The widths were estimated at each height separated by 1 Mm, with the mean width taken as an average of values at all heights for the lifetime of the jet. Our simulations suggest a strong influence of magnetic field strength in determining the width(s) of the jet structure with estimated cross-sectional widths inversely related to the magnetic field magnitudes. This might be due to the higher internal pressure of the chromospheric jet, as compared to the ambient environment, during the rise of the feature through the stratified solar atmosphere. This results in expansion of the jet structure and subsequent bending/distortion of the surrounding field. However, in the case of higher magnetic field magnitudes, the jet experiences strong tension forces that result in greater collimation of the jet structure. These results further indicate the possibility of lower cross-sectional widths of jet features located near the strong magnetic field environment than those in the quiet Sun regions.
Another important aspect that influences the jet width is the lifetime of the driver of the jet. The simulations suggest a linear relationship between the driver periods and the jet crosssectional widths. This could be due to the longer supply of plasma to the jet by a sustained driver resulting in increased widths. However, the amplitude of the driver had a minimal effect on the width for the jets aligned with radial magnetic fields, though this might not be the case with the jet direction misaligned with the background magnetic fields.
Based on the results of the parameter scan, we define a "standard" jet from our simulations, which has characteristics more akin to the classical spicules. The result is shown in Figure 5, with P = 300 s, B = 60 G, and A = 60 km s −1 where the top, middle, and bottom panels highlight variations in the density, temperature, and numerical Schlieren (indicative of normalized density gradient magnitude). The numerical Schlieren is defined as, where c 0 = 5, c 1 = 0.05, and c 2 = − 0.001. It must be noted that even though the numerical Schlieren is a physically defined term, its use in our analysis is purely from a qualitative perspective due to their ability to visualize significant variations of a chosen parameter. An important facet of our simulation is related to the field-aligned jet motions ( Figure 5) that show distinct kinematic and morphological characteristics during the rise and fall phase. In the rise phase (left panels of Figure 5), the simulated jet show complex internal substructures and sausage-like deformation of crosssectional widths in density and Schlieren data. However, during the fall phase (right panels of Figure 5), the complex/internal beam substructure cease due to driver switch-off. At this stage, the boundary deformation is no longer symmetric, though the jet axis has noticeable transverse displacement akin to an ideal MHD kink wave mode. Modification in sausage-like wave properties associated with cross-sectional width estimates during the rise and fall phases of simulated jets were recently reported by e.g., (Dover et al. 2020).
The temperature profile of the numerical jet suggests periodic distortions in the TR layer over the jet lifetime. The TR deforms as the jet penetrates this layer while generating waves, also known as TR quakes (Scullion et al. 2011). During its entire lifetime, the jet structure remains isothermal/cool without much variations for any combinations of the selected parameters. Figures 6 and 7, showcases for the simulated jet with a density structure highlighting variations in the parameter space in each row for similar time-instance as shown in Figure 5. The top and middle panels in Figure 6 exhibits the effect of strong (B = 80 G) and weak (B = 20 G) magnetic field strengths, respectively. For the strong magnetic field, the jet was found to be more collimated with higher density, mostly concentrated near the apex. In this case, there is an increase in the number of crisscrosses and no sausage-like jet boundary deformations. However, in the case of the weak magnetic field, the jet is more diffusive as it balloons out in higher solar atmospheric layers. Here, the simulations do suggest a clear complex beam structure and sausage-like deformation.  velocity magnitudes on jet morphology. For higher velocity (A = 80 km s −1 ), the simulated jet had an internal substructure similar to the standard jet, but with more crisscross/knot features at matching time steps. The jet shows kinking behavior from t = 223 s onward, possibly due to the initial high velocity in the slender structure, which is more susceptible to kink instability. However, at t = 288 s the internal substructures disappear with the standard jet, indicating a close association between the driver time and knot lifetimes. For A = 20 km s −1 (top panel, Figure 7), the simulations suggest that for a lower velocity the jet heights are reduced along with the absence of any significant substructures.
The middle and bottom panels of Figure 7 highlight the role of the driver period for two characteristic cases, one with magnitudes P = 50 s another with 200 s. In the case of the shorter driver period (P = 50 s), the jet at t = 21 s shows a similar evolutionary trend as for those with longer duration (P = 200 and 300 s), with a lifetime of around t = 223 s. For long duration drivers (P = 200 s), the jet boundaries are smooth and show no kink-like deformation in the fall phase. A key difference between the two driver periods is the amount of mass injection in the jet, which has direct implications for the cross-sectional widths. At the end of the lifetime of the jet (t = 223 s), the jet has no complex internal substructures, suggesting a stronger role of initial velocities in generating these substructures.

Jet Beam Substructures
Our simulations for chromospheric jet structures revealed a peculiar aspect with complex internal substructures. Multiple cases with varied driver amplitudes/periods and magnetic field strength ( Figures 5 and 6) highlighted crisscross/knot patterns within the jet boundary, along with the cross-sectional distortion of the feature. Similar behavior appears to be ubiquitous and was reported in previous studies associated with astrophysical (van Putten 1996;de Gouveia Dal Pino 2005;Hada et al. 2013;Cohen et al. 2014;2017AnA606A103H et al. 2017) and laboratory jets (Menon & Skews 2010;Edgington-Mitchell et al. 2014;Ono et al. 2014).
Simulation results indicate the presence of these patterns at the sites of higher densities/pressure within the jet structure. We refer to these sites as knots in our investigation and propose two possible mechanisms for their formation: 1. Knots are a manifestation of nonequilibrium pressure forces prominent within and around the jet structure. The internal pressure results in expansion of the jet crosssectional area, which is then counterbalanced by the tension forces of the jet. Eventually, the magnetic pressure (force) dominates the internal plasma pressure and the jet boundary decreases while resulting in the compression of the jet structure. These processes over height and time appear as cross-sectional deformations of the jet, along with the sites of high density/pressure as knots within the jet structure. 2. The knot features could also be the sites within the jet structure where internal shock waves (Norman et al. 1982) become reflected from the jet boundary. If the jet has supersonic velocities (simulated jet Mach numbers are around ∼1-3), this can give rise to a myriad of internal substructures due to high velocities.
Continuing in the framework of (2), a schematic overview of these beam (sub)structures due to supersonic flow are given in Figure 8, highlighting the osculating jet boundary and cross patterns occurring inside the jet. The jet boundary undulates as the gas periodically overexpands and collapses inwards, as it tries to reach an equilibrium with the ambient atmospheric pressure. The jet structure repeatedly overshoots its equilibrium primarily due to the effects of the boundary, communicated to the interior by sound waves, which are traveling slower than the supersonic flow of the jet. This leads to the formation of low/high regions of pressure in the jet being out of phase with the cross-sectional width (Figure 8).
The mechanism for the formation of these low/high pressure points within the jet structure could be explained in terms of the perturbed pressure balance between the jet and its surrounding Figure 5. Panels show snapshots of temporal evolution of the "standard" numerical jet with uniform radial magnetic field in stratified atmosphere with a driver period (P = 300 s), magnetic field (B = 60 G), and amplitude (A = 60 km s −1 ). The density (top), temperature (middle), and Schlieren (bottom) images highlight complex internal substructures (knots) of the simulated jet with bright apex, possibly due to high-density concentrations. The temperature (middle) panel indicates an isothermal nature of the jet structure with colder plasma components from the photospheric layer. The 9 s animation shows the dynamics of the jet. The animation runs from t = 0.0 s to t = 605.4 s. (An animation of this figure is available.) Figure 6. Evolution of the simulated jet with different combination of magnetic field and driver amplitude parameters, in a stratified solar atmosphere with a uniform radial magnetic field. Parameters are varied for the "standard" jet configuration to identify the effects of a particular parameter on jet morphology and kinematics. Panels area. Since the jet internal pressure is greater than the ambient atmospheric pressure, the dynamics of the jet become similar to an underexpanded jet (Norman et al. 1982;Edgington-Mitchell et al. 2014). As a result, crisscross pattern appear due to a series of internal shock waves and subsequent expanding structures, identified as knots in our simulations. These expanding fans or knots (blue dashed lines in Figure 8) are outward flows formed at the sites of enhanced pressure differences between the jet (internal shocks) and ambient atmosphere.
The Mach lines of these expansion waves reflect inwards at the jet boundary while forming compression waves/fans due to the pressure continuity (red lines in Figure 8). These compression waves are reflected at a near-constant angle at the jet boundary and, due to the curvature at the edges, the Mach lines of the reflected compression waves tend to converge in a conical shock wave before arriving at the jet axis. Depending on the angle between the incident shock and jet axis, the shock could either be reflected back or form a Mach disk for small and large angles, respectively (see black lines, Figure 8). As the plasma flow passes through this shock region, it will increase the local pressure at the site. This increase in pressure, together with the outward expanding fan, allows repetition of the process.
This explains the formation of a stationary crisscross pattern and the knots along with cross-sectional deformations, density enhancements, and cavities in the jet structure. Also, this provides a vital clue about the disappearance of these complex substructures during the drivers switch-off phase. In our simulations, the standard jet ( Figure 5) clearly showcases these knot-like structures and, with the driver switch-off at t = 288 s, these internal substructures disappear. Similar behavior is evident from Figures 6 and 7; though, the presence of knots appears to be sensitive to particular parameters like a strong magnetic field, low driving speed, and shorter driving times.

Bright apex of the Jet
A peculiar feature of our simulations (see, e.g., Figures 5(b)-(c)) is the presence of a bright, bulb-like apex with higher density/emission as compared to the rest of the jet structure. This feature is also prominent in multiple simulated jets with varied driver amplitude, period, and/or magnetic field strength magnitudes (see, Figures 6(b), (c), (k), (l), 7(f), (g)). Schlieren images of the jet apex ( Figure 9) reveal regions with strong pressure imbalances within the (sub)structure formed due to the internal shockwave patterns. These regions, with alternate high and low pressures, suggest possible sites for superimposition of internal shock waves with leading high density/pressure concentration due to upward flows of plasma and/or waves.
Observational evidence for these bright bulb-like substructures at the jet apex was recently reported in a few studies  observations for spicule structures and identified bright leading substructures as "heating fronts." The authors proposed the formation of these regions as a consequence of upward propagating magnetic-hydrodynamic waves and/or the dissipation of electric currents due to ambipolar diffusion between ions and neutrals. Srivastava et al. (2018) reported ubiquitous observations of chromospheric jets with a leading bright apex, trailed by a dark region in running-difference images. They labeled these jets as "tadpoles" and suggested the role of pseudo shocks in the formation of these observed bright (sub) structures.

Conclusion
In this paper, numerical simulations for the propagation of momentum pulse, confined in a chromospheric jet structure, are performed and investigated. The simulated jet is localized, mimicking the environment of a stratified lower solar atmosphere. The obtained computational jet revealed complex internal substructures with overall field-aligned/longitudinal plasma motions akin to observed spicule dynamics. The jet, for a varied driver (amplitude, period) and magnetic field configurations showed distinct changes in morphological (cross-sectional width, maximum length) and kinematic (transverse displacement, longitudinal trajectory) profiles. The key findings of our investigation are: 1. The visible apex of the jet structure followed a nonballistic parabolic trajectory with heights consistent with observed and reported spicule motions. A parameter scan suggests a strong influence of the driver amplitude in determining the longitudinal dynamics of a jet structure. 2. Our simulations suggest that jets emanating from regions with higher background magnetic field strengths are more collimated/dense as compared to those originating in conditions typical of the quiet Sun. 3. Complex internal substructures are formed during the rise phase of the jet structure. So-called crisscross/knot-like patterns are formed possibly due to the reflection of internal shock waves at the jet boundaries. 4. Temperature maps from our simulations reveal periodic distortions of the transition region due to the penetration of the jet structure at this solar atmospheric layer.
However, the jet remained isothermal for its entire lifetime for all combinations of the initial parameters. 5. The simulated jets in our numerical experiment show a bright, bulb-like apex, similar to the observed cases of chromospheric jets. These could be the regions with high density/pressure due to upflow of plasma and/or waves.
Our study provides an initial overview of the role of a range of driving parameters and background magnetic field conditions on the observed dynamics of chromospheric structures (spicules). As a next step, we plan to examine the role of misalignment between jet propagation and background magnetic field topology to further identify its impact on jet morphology and dynamics. Here, for the first time in solar jets, the possible presence of internal substructures in these jets are highlighted, which can provide important insights into the generation mechanisms of chromospheric jets (pulse-, reconnection-, convection-, solar global oscillation-, or vortexdriven, etc.). Furthermore, the (simulated) substructures found in these highly localized chromospheric jets may give a new window and impetus to investigate the jet behavior at the highest possible resolution available by using data, e.g., from the Daniel K. Inouye Solar Telescope (DKIST) (Rast et al. 2019;Rimmele et al. 2020) and, in the future, the upcoming European Solar Telescope (