Criss-cross Nebula: A Case Study of Shock Regions with Resolved Microstructures at Scales of ∼1000 au

Using integral field spectroscopy (IFS) from MaNGA, we study the resolved microstructures in a shocked region in the Criss-cross Nebula (CCN), with an unprecedentedly high resolution of ≲1000 au. We measure surface brightness maps for 34 emission lines, which can be broadly divided into three categories: (1) the [O iii] λ5007-like group including seven high-ionization lines and two [O ii] auroral lines that uniformly present a remarkable lane structure, (2) the Hα λ6563-like group, including 23 low-ionization or recombination lines that present a clump-like structure, and (3) [O ii] λ3726 and [O ii] λ3729 showing high densities at both the [O iii] λ5007 lane and the Hα clump. We use these measurements to constrain resolved shock models implemented in MAPPINGS V. We find our data can be reasonably well fitted by a model that includes a plane-parallel shock with a velocity of 133 ± 5 km s −1, plus an isotropic two-dimensional Gaussian component, which is likely another clump of gas ionized by photons from the shocked region, and a constant background. We compare the electron density and temperature profiles as predicted by our model with those calculated using observed emission-line ratios. We find different line ratios to provide inconsistent temperature maps, and the discrepancies can be attributed to observational effects caused by limited spatial resolution and projection of the shock geometry, as well as contamination of the additional Gaussian component. Implications on shock properties and perspectives on future IFS-based studies of the CCN are discussed.


INTRODUCTION
Shocks may be produced by multiple processes that are key in galaxy formation and evolution, such as stellar winds, supernova feedback, AGN feedback, and gas accretion/outflow at different scales.There have been many studies of shocks and shock-induced structures in both the Milky Way and external galaxies, e.g.shockinduced filaments in the Orion-Eridanus superbubble (Pon et al. 2014), the shock-excited Herbig-Haro objects (Dopita & Sutherland 2017), and shock ionization in passive galaxies (Yan 2018), star forming galaxies (Shopbell & Bland-Hawthorn 1998;Veilleux et al. 2005;Sharp & Bland-Hawthorn 2010;Westmoquette et al. 2011;Ho et al. 2014), galaxy mergers (Farage et al. 2010;Soto et al. 2012;Soto & Martin 2012;Medling et al. 2021), AGN host galaxies (Dopita et al. 1997(Dopita et al. , 2015;;Terao et al. 2016;Molina et al. 2018Molina et al. , 2021;;Holden et al. 2023) and post starburst galaxies (Alatalo et al. 2016).These studies have deepened our understanding of shocks in different aspects.However, previous observations of shocked regions are mostly limited to either narrow-band imaging of emission lines with very low spectral resolution and no kinematics information, or spectroscopy with limited spatial coverage and resolution.It is crucial to have resolved spectroscopy with high spatial and spectral resolutions in order to fully understand not only the shocks themselves but also their roles in the dynamic heating of interstellar medium (ISM) and destruction of interstellar clouds.
There have also been many studies of nearby shocked regions, such as Crab Nebula at ∼ 2 kpc (Chevalier & Gull 1975;Clark & Stephenson 1977;Davidson & Fesen 1985;Nomoto 1987;Fesen et al. 1997;Hester 2008;Martin et al. 2021) and IC443 at ∼ 2 kpc (Rho & Petre 1998;Jones et al. 1998;Chevalier 1999;Carpenter et al. 1995;Olbert et al. 2001;Troja et al. 2008;Alarie & Drissen 2019;Deng et al. 2022), which usually present regular shell-like structures and have high shock velocities and electron densities.In this work we study the Criss Cross Nebula (CCN), which has different structures and physical conditions, compared with other well studied nearby shock regions.Located at (α = 62.55 deg, δ = -4.98 deg) or (l = 196.97 deg, b = -37.83 deg) and with an angular size of 6 ′ × 3 ′ , CCN is dubbed by its unique morphology and diverse microstructures and is believed to be produced by one or multiple shocks (see Figure 1).This nebula was firstly studied by Zanin & Weinberger (1997, hereafter Z97) and Zanin & Weinberger (1999, hereafter Z99), who suggested that the CCN is ionized by a shock with a velocity of about 40 km/s and electron density in the range 1−100 cm −3 based on narrow-band images and long-slit spectra.A following study by Temporin et al. (2007, hereafter T07) found a photoionized arc-like structure of Hα emission near the CCN, which was called "canopy" by the authors considering its shape and the relative position to the CCN.Though still uncertain, the distance to CCN should be within ∼ 400 pc and could be even more nearby, at ∼ 150 pc as argued in T07.CCN is projected at the same area of the Orion-Eridanus superbubble and is very close to Arc A, a major structure within the bubble.The small distance, cross-like morphology, low shock velocity, and complex large-scale environment make the CCN a unique object for studies of shocks.
In this work we analyze the integral field spectroscopy (IFS) data for a small region of the CCN (∼ 20 ′′ × 20 ′′ ), which was observed as a filler target by the Mapping Nearby Galaxies at Apache Point Observatory (MaNGA; Bundy et al. 2015) survey.At a distance of 150-400 pc, the spatial resolution ∼ 2.5 ′′ of MaNGA corresponds to a physical scale of ∼ 375 AU to 1000 AU.Although covering only a small part of the nebula, the IFS data from MaNGA allows for the first time a case study of resolved microstructures in shock regions with unprecedentedly-high resolution.We will measure all the emission lines in the MaNGA FoV and use the measurements to constrain both physical parameters and geometry parameters of the underlying shock.
The paper is structured as follows.In section 2 we describe our data and measure the emission lines.In section 3 we use the emission line measurements to constrain shock models.We discuss in section 4 and summarize our work in section 5.

MaNGA Data
A small part of the CCN has been observed as a filler target by the MaNGA survey, which is one of the three major experiments of the fourth generation Sloan Digital Sky Survey (SDSS-IV; Blanton et al. 2017).During its six-year operation from July 2014 through August 2020, MaNGA obtained integral-field spectroscopy (IFS) for ≳ 10 4 nearby galaxies selected from the NASA Sloan Atlas (NSA; Blanton et al. 2011).The IFS data were obtained using 17 science integral-field units (IFUs) with a field of view (FoV) ranging from 12 ′′ to 32 ′′ (Drory et al. 2015).The MaNGA spectra cover a wavelength coverage from 3600 Å to 10300 Å with a spectral resolution of R ∼ 2000, taken with the two dual-channel BOSS spectrographs (Smee et al. 2013) on the Sloan 2.5-meter telescope at the Apache Point Observatory (Gunn et al. 2006).Details of MaNGA sample selection, flux calibration, survey execution and observing strategy can be found in Wake et al. (2017), Law et al. (2015), and Yan et al. (2016a,b).MaNGA raw data were reduced using the MaNGA data reduction pipeline (DRP; Law et al. 2016) which produces a 3D datacube for each target with an effective spatial resolution that can be described by a Gaussian with a full width at half maximum (FWHM) of ∼ 2.5 ′′ .All the MaNGA data, including raw data and DRP products, are released as part of the final data release of SDSS-IV (DR171 ; Abdurro 'uf et al. 2022).
We use the datacube produced by the MaNGA DRP for this study.The upper-right panel in Figure 1 shows the MaNGA FoV, a hexagon covering approximately 20 ′′ ×20 ′′ , thus only a small part of the CCN as indicated by the cyan hexagon in the upper-left and upper-middle panels.As an example, the lower panel of the same figure shows the spectrum of the central pixel (indicated by the black cross in the upper-right panel), with the short black lines marking some of the emission lines.
The MaNGA DRP creates a sky model for each plate based on data from the sky fibers, in a way that is designed and optimized for nearby galaxies.For the CCN, however, the "background" measured by the sky fibers is contaminated by the emission from diffuse ionized gas surrounding the CCN.We apply our spectral fitting method described in Section 2.3 to the sky model constructed by DRP, finding some emission lines to be oversubtracted.We do not attempt to rebuild the sky model.Rather, we measure the flux of these over-subtracted lines in the sky model, and we take into account their effect in the following analysis.A total of 4 emission lines to what we measured from the DRP datacube.We notice that the over-subtracted flux is typically one order of magnitude smaller than the flux measured from the DRP spectrum, and so the flux correction is generally a tiny effect on our results.

Ancillary data
We have obtained the [Oiii] λλ4959, 5007 narrowband image for the whole CCN, using the China Near Earth Object Survey Telescope at Xuyi Observational Station operated by Purple Mountain Observatory (PMO).The observation was taken on December 30, 2021, with a total exposure time of 1.5 hours.In addition, kindly provided by Bringfried Stecklum, the narrowband images of [Sii] λλ6717, 6731 and Hα λ6563 used in T07 are also available for our study.The upper-left panel in Figure 1 shows the RGB image of the whole CCN, constructed by combining the narrowband images of [Oiii], Hα and [Sii].
We observed the 12 CO(1-0) line for the CCN using the PMO 13.7-m millimeter telescope (Zuo et al. 2004) located in Delingha, China, over the period from April 26, 2022 through May 6, 2022.A total of 136 scans were taken, and each scan was reduced to produce a datacube with a pixel size of 30 ′′ by applying the standard reduction pipeline of this telescope.After fitting and subtracting the baseline in the spectrum and masking any heavily affected spectra, all unmasked spectra were stacked to create a final spectrum, which has a rms of T CO = 0.039 K in 0.15 km/s.This spectrum provides no detection of 12 CO(1-0) with a 3σ upper limit of 0.14 K km/s when assuming a line width of 10 km/s.If we adopt a conversion factor of α CO = 4.35 (M ⊙ /pc 2 )/(K km/s), this upper limit corresponds to an H 2 surface density lower than 0.6 M ⊙ /pc 2 .

Emission line measurements
The MaNGA spectra present numerous emission lines but almost no continuum over the whole wavelength coverage, as can be seen from the example spectrum in Figure 1.Therefore, for each spaxel we do not attempt to obtain an overall continuum component, and we directly apply our emission line fitting code to measure all the detected lines.These include a total of 34 emission lines, as listed in Table 1.Although the continuum is close to zero in all the spectra, the baseline beside each emission line always deviates to some degree from a zero horizontal line due to the residual in sky subtraction and/or flux calibration in the MaNGA DRP.We determine a local baseline with a linear fitting of the spectrum over both sides of each emission line, excluding obvious emission features and outliers by sigma clipping.
Each emission line is modeled with a single Gaussian.Most of the lines are well separated from each other, and each line is fitted independently.In some special cases as listed in Table 1, two or more neighboring lines are fitted simultaneously.In these cases some Gaussian function parameters may be tied for reasonable considerations, which are also listed and explained in the table.For a given line or a group of lines, the best-fit model is determined by minimizing the χ 2 , defined as where n is the number of wavelength pixels over the emission line fitting window centered at the restframe wavelength of the line, flux i and model i are the observed and model fluxes of the i-th pixel, and err i is the error of the observed flux of the i-th pixel.The width of the fitting window ranges from 20 Å to 100 Å, mainly depending on the number of lines included in the fitting.
As an example, Figure 2 displays the result of the fitting for all the 34 lines for the central spaxel in the MaNGA datacube.As can be seen, all the lines are well-fitted.
In fact, the majority of the lines in the whole MaNGA datacube are well-fitted, except that for some weak lines such as [Ariii] λ7135, the emission line signal in some spaxels is overwhelmed by noise.In such cases, our fitting code gives best-fit models with small fluxes and relatively large uncertainties, which can be identified as non-detections given a signal-to-noise (S/N) threshold.
Figure 3 displays our measurements of the surface brightness map for all the 34 lines.In each panel we show the spaxels with S/N > 3 only.According to the morphology of the maps, we find all the lines can be broadly divided into three groups.In the figure, the three groups of maps are shown separately with a blank gap inserted in between every two neighboring groups, while within each group the lines are ordered by increasing wavelength.As can be seen, the maps in the first group present a remarkable lane of high brightness across the upper-left part of the MaNGA FoV.This group consists of seven high-ionization lines and two [Oii] auroral lines: In what fol-lows this group is referred to as "Group I", or "[Oiii] λ5007-like group" to reflect the strongest fluxes of the [Oiii] λ5007 line in this group.We note that, for [Ariii] λ7135 the lane feature is less pronounced than other lines of the same group, but we still include it in Group I considering that the emission in the lane region is the strongest and the only significant feature of this line.
The second group is referred to as "Group II" or "Hα λ6563-like group", which is made up of 23 low-ionization or recombination lines: Hη λ3836, Hζ λ3889, Hϵ λ3970, The main feature of this group is a clump-like emission region located in the lower-right part of the MaNGA FoV.Most of the maps in this group show a continuous area of clump, while some weak lines show more clumpy and less concentrated structures.We include both cases in this group, as long as the peak and the majority of the high brightness area are located in the lower-right part of the MaNGA FoV.The third group includes only two lines, [Oii] λ3726 and [Oii] λ3729, which show clump-like features at both the position of the [Oiii] λ5007 lane and the position of the Hα clump, with a bridge-like structure connecting the two clumps.
To more clearly see the similarities of the maps within each group and the differences between different groups, in every panel we repeatedly show the brightness map of both [Oiii] λ5007 (the solid black contours) and Hα λ6563 (the dashed contours).In the first group the lanelike feature shows up in all the maps with highly similar shapes and locations, while the clump-like feature seen in the Hα map also similarly shows up in most of the maps in Group II.When comparing the two groups, it is noticeable that the lane-like feature in Group I has little spatial overlap with the clump-like feature in Group II.The spatial offset between the two groups of lines implies that the MaNGA region of the CCN is likely to be shocked by a blast wave moving from lower-right to upper-left.On the other hand, however, the fact that the Hα bright region is clump-shaped instead of laneshaped cannot be simply explained by a single shock.
Figure 4 displays the velocity map for eight emission lines which have reliable velocity measurements, typically with uncertainties less than 10 km/s.Overall, all the lines show relatively positive velocities in the lowerright region and negative velocities in the upper-left region.When comparing the velocity maps more closely, one can easily see both the high similarities within each group and the apparent differences between different b Generally each emission line is modeled by a single Gaussian, denoted as G1(a1, λ1, σ1) with three free parameters: amplitude a1, center wavelength λ1, and the line width σ1 resulting from kinematical and instrumental broadening.In some cases, neighbouring lines are fitted simultaneously by a linear combination of multiple Gaussians with some parameters tied during fitting, for the reasons provided in column "Notes".c Classification group according to the surface brightness morphology of the emission line as described in detail in subsection 2.3.We use notation ⟨1⟩ for the [Oiii] λ5007 like group, ⟨2⟩ for the Hα λ6563 like group, and ⟨3⟩ for the group of [Oii]   [NeIII] λ3869 [OIII] λ4363 [OIII] λ4959 [OIII] λ5007 [ArIII] λ7135 [OII] λ7320 groups.For instance, the velocity maps of [Oiii] λ5007 and [Oiii] λ4959 show a lane-like structure in the same position as their surface brightness maps, a unique feature which is not seen in the velocity maps of other groups.These velocity maps should in principle provide useful constraints to the physical processes behind the CCN.However, due to the limited spectral resolution of MaNGA, it is hard to decompose the different velocity components from the spectra.Therefore, we will use only the surface brightness measurements when constraining the shock models below in section 3.

The shock model
In this section, we attempt to fit the MaNGA data by applying the shock model implemented in the MAPPINGS V code (Sutherland & Dopita 2017;Sutherland et al. 2018).The basic ideal of MAPPINGS V is to calculate the  To compare with the resolved MaNGA data, we use shock profiles calculated by the model, i.e. the physical properties and emissivity as functions of the separation r from the shock front, instead of the integrated spectrum which is the default output of the model.In this case, the shock front locates at r = 0, the region with r > 0 has already been shocked (thus referred to as "shocked region"), and the region with r < 0 is expected to be shocked in future (thus referred to as "precursor").The precursor may have been ionized by the photons escaping from the shocked region.Since MAPPINGS V adopts parallel-slab geometry, in the direction perpendicular to the shock direction the physical properties and emissivity are homogenous in the model.Therefore, the model is essentially a one-dimensional model, and we need to simplify the geometry of the observed region in consideration in order for appropriate comparisons between model and data.Our geometry model is illustrated in Figure 5.We assume the region is a "slab" of ionized gas with a certain thickness, oriented with an inclination angle of i (the angle between the "slab" and the line of sight).As a result, the gas slab has a line-of- sight depth of D los , and so the slab thickness is given by D = D los cos(i).
For a given emission line, the observed surface brightness profile, F obs (r a ) where r a is the angular separation from the shock front, can then be modelled by In the equation, the symbol "⋆" represents the convolution operation.Here, f (r a ) is the intrinsic emissivity profile given by the model.For the calculation of f (r a ), we adopt a distance of 150 pc to CCN in order to convert r a to physical separations with the consideration of projection effect (Uncertainties in CCN distance and the effects on our modelling will be discussed in subsection 4.1.).Due to the thickness and inclination of the gas region, the observed surface brightness at given r a is effectively an integral (along line of sight) of intrinsic emissivity at different physical separations.To take into account this effect, we first obtain an average surface brightness at r a by convolving f (r a ) with the averaging filter A(K t ) with a kernel size of K t = D los sin(i), which is then multiplied by the line-of-sight depth D los to give the total surface brightness, as can be seen from the first term at the right hand of Equation 2. This term is further convolved by G 2.5 , a Gaussian filter with FWHM = 2.5 ′′ , which models the point spread function (PSF) of the MaNGA data.Finally, the last term of the equation represents an additional component, F b , to account for the background emission which could be either a constant background/foreground in the line of sight (signal or noise, or both), or produced by ionizing mechanisms other than the shock.As we will show in subsection 3.2, the model with a constant background cannot explain our data.Thus, we have an extended model in subsection 3.3.
In addition, we also need to determine the location of the shock front and the projected direction of the shock.For simplicity, we fix the shock direction as the diagonal line going from the lower-right corner to the upper-left corner in the MaNGA FoV.This direction is perpendicular to the lane-like structure in the [OIII] surface brightness map, at which we have also seen sharp jumps in some line ratios (see Figure 6; the shock direction is indicated by the arrow in each panel).The shock should have moved from the Hα clump region to the [Oiii] lane, so that the latter should be closer to the shock front.The location of the shock front cannot be determined by simple considerations, and this is modeled by a free parameter, ∆S, the angular distance of the shock front to the central spaxel of the MaNGA FoV.
The physical parameters of the shock considered in this work include the pre-shock hydrogen density n pre , the shock velocity v sh , and the transverse magnetic field B. Considering the fact that CCN is projected in and also possibly located in the Orion-Eridanus superbubble, we adopt the metal abundance of the Orion nebula provided in MAPPINGS V. We note that changes in abundance can result in complex and non-linear modifications of shock properties.For instance, under certain parameters, a shock with a high metallicity can produce a higher [Oiii] λ5007/Hβ λ4861 ratio when compared to its low metallicity counterpart.However, we find that, if the metallicity and abundance are not fixed, the model will have too many degrees of freedom and will be very hard to be constrained well with the available data.Therefore, we opt to fix the metallicity and abundance in this work and leave the effect of potential variations in metallicity or abundance to be studied in future when more data becomes available.

Model 1: shock with a constant background
We first consider the model with a constant background component.We have attempted to fit the shock and the background components simultaneously using all available line and line ratio measurements from the full MaNGA FoV, but finding ourselves unable to obtain a reasonably well constrained solution.We thus opt for a two-step approach, in which we first constrain the shock parameters, including both the physical parameters (n pre , v sh , B) and the geometric parameters (K t , D, ∆S), by applying the model to a selected set of emission lines and line ratios from a "clean region" selected out of the MaNGA FoV, and then we constrain the background component F b for all the emission lines from the full MaNGA FoV but fixing the parameters of the shock obtained from the first step.
The measurements used in the first step include the surface brightness profile of Hβ λ4861 representing the Hα λ6563-like group, the profile of [Oiii] λ5007 representing the [Oiii] λ5007-like group, and the profiles of [Oii] λ3729 and [Oii] λ3726 from the third group.We include both lines from the [Oii] doublet because they are comparably strong, unlike the other two groups in which one or two lines are much stronger than the other lines.For the Hα λ6563-like group, we use Hβ rather than Hα because the former is closer to [Oiii] λ5007 in wavelength, and so the effect of dust extinction can be mitigated (although the dust extinction in CCN should be very weak; see subsection 4.2 for discussion).In fact, we have repeated the fitting procedure using Hα to replace Hβ in the first step, and we find very similar results (see Table 2).
The "clean region" used in the first step is indicated by a black quadrangle in each panel of Figure 6.This region is selected to have simple shock profiles in the line ratios along the projected direction of the shock (as indicated by the arrow).The exact boundary of this "clean region" is manually specified to maximize the number of spaxels included in the modelling, and at the same time to minimize the contamination of more complex features.We notice that the apparent difference between regions inside and outside the "clean region" is a hint that the whole region cannot be simply modeled with a single shock.We will come back and discuss this point in subsection 4.4.
It is time consuming to calculate the shock properties for a given set of model parameters, and so the calculation would be very expensive when the geometry parameters are fitted simultaneously with the shock parameters.In practice, we use a "zoom-in" grid fitting approach.We start with a coarse grid of the shock model  parameters, and find the best-fit geometry parameters for each node of the grid by minimizing the reduced χ 2 .
The best-fit model in this grid is then given by the node at which the reduced χ 2 is the smallest over the whole grid, χ 2 min .A weight is then given to each of the nodes by W = exp(− 1 2 (χ 2 − χ 2 min )), and the 1 σ uncertainty of each parameter is calculated from the weighted average of all the nodes.If the uncertainties are significantly smaller than the step sizes of the grid, we construct a finer grid around the best-fit node and repeat the above process until the uncertainties are comparable with the step sizes of the grid.
The best-fit results from the first step are shown in Figure 7, where we plot the observed profiles from the clean region and the best-fit shock profiles.The bestfit model parameters and uncertainties are listed in Table 2.The model parameters constrained by using the Hα line instead of Hβ are also listed in the table.Generally, all the profiles are well-fitted in the whole range of r a probed, except that the model slightly overpredicts (underpredicts) the ratio of [Oii] λ3729/Hβ λ4861 at the smallest (largest) r a (see the third panel of the figure) and the peak of [Oiii] λ5007/Hβ λ4861 is a bit shallower in model.It is encouraging that the models with Hβ and Hα give rise to very similar constraints for all the model parameters.Accordingly, the shock has a velocity of v sh ∼ 110 km/s, a pre-shock hydrogen density of n pre ∼ 1.3 cm −3 , and a transverse magnetic field which is very close to zero (though with relatively large uncertainties).The shock velocity is larger than the value found by Z97, which is 40 km/s.This discrepancy may not be surprising, considering that the targeted position of the spectroscopic study in Z97 is outside the MaNGA FoV.The geometry parameters of our best-fit model show that the direction of shock is almost perpendicular to our line of sight, with an inclination angle of 22 deg (with Hβ) or 19 deg (with Hα).The shock front is about 7 arcsec from the central spaxel, thus close to the upper-left edge of the MaNGA FoV.
Figure 8 shows the best-fit results for some lines as obtained from the second step, in which we fix the shock model parameters and further constrain the background component for individual lines.For a given line, the figure displays the maps of the observed surface brightness, the surface brightness predicted by the model, the residual defined as the difference between the observation and model prediction (Res), and the relative residual (RR) defined as the ratio of the residual relative to the geometric mean of the surface brightness measurement and its error for each spaxel, i.e.RR = Res/ √ S 2 + err 2 , where S is the surface brightness of the emission line in the given spaxel, and err is the 1σ uncertainty of the surface brightness.For spaxels with high S/N, the relative residual is effectively equal to fractional error, while for low S/N spaxels, the relative residual is reduced to the ratio between residual and measurement error, thus telling whether the residual can be explained by data uncertainties.
As can be seen, our model successfully captures some main features in the data, such as the [Oiii] lane-like structure and the overall difference in surface brightness between the upper-left and lower-right regions.However, the model substantially underpredicts or overpredicts the surface brightness for certain regions and for all the lines considered, which can be easily identified from the relative residual maps.

Model 2: shock with a 2D Gaussian component
The above result strongly suggests that the F b component in our model cannot be simply a constant background.We have considered a different model in which we add another shock, that is, the observed region is jointly ionized by two distinct shocks, but found such kind of model can not fit the data well (see subsection 4.5).
In this section we present our second model (Model 2 hereafter) in which we include a two-dimensional (2D) Gaussian component in addition to the shock component and the constant background.In this case, Equation 2 is rewritten as (3) where the new term F G (x, y) is an isotropic 2D Gaussian function (σ x = σ y ).Note that the observed surface brightness and the intrinsic emissivity are also written as two-dimensional functions of spaxel position (x, y).This is because the total surface brightness in both data and model can no longer be described by a one-dimensional profile after the Gaussian component is included.Different emission lines share the same Gaussian center (x, y), but have their own width σ x and amplitude A. The fitting procedure is divided into two successive steps as for the first model (Model 1 hereafter).Different from Model 1, for the first step of Model 2 we use more emission lines and all available spaxels from the MaNGA FoV to constrain the shock parameters, considering the larger number of model parameters and the more complex model components.Specifically, we use the surface brightness profiles of This model can fit the data much better than Model 1, although the physical origin of the Gaussian component still remains unclear.The best-fit parameters and uncertainties of the shock component are listed in Table 2. Compared with Model 1, the shock in Model 2 has a slightly higher pre-shock hydrogen density (1.4 cm −3 versus 1.3 cm −3 ), slightly larger velocity (133 km/s versus ∼110 km/s), and a much stronger magnetic field (1.6 µG versus ≲ 10 −4 µG).For the geometric parameters, both K t and D are slightly smaller than those from Model 1, while ∆S increases from ∼7 arcsec to 17.2 arcsec.The inclination angles from the two models are quite similar, at around 19 deg.
The fitting results from the second step are shown in Figure 9, for the same set of emission lines as in the previous figure.As can be seen, the RR is significantly reduced in most spaxels and for all the lines when compared to the fitting results of Model 1. From the model-predicted maps, one can easily find both the differences between different groups of lines and the similarities within the same group: the three lines from the [Oiii] λ5007-like group are uniformly dominated by the lane-like structure, the Gaussian component dominates all the lines from the Hα λ6563-like group, and both components are presented in the map of [Oii] λ3729 with comparable amplitudes.We have examined the ratios of different lines for the Gaussian component alone, finding them to be comparable to the expected line ratios of a clump of gas ionized by the photons from the shocked region (subsection 4.5).Although the fits are substantially improved with respect to Model 1, one can still see patterns of significant residuals.For instance, the maps tend to present positive values of RR in the left corner for the [Oiii] λ5007-like lines as well as the upper-right edge for the Hα λ6563-like lines.At the border of the MaNGA FoV, these regions could be adjacent to other regions of the CCN that may be ionized in a different way.In addition, a lane-like residual pattern can be seen in/around the original [Oiii] lane structure, for all the emission lines but with different orientations for different line groups.This pattern is likely to be induced by another shock along the upper-right to lower-left direction.We will come back and discuss on this point in subsection 4.4.We evaluate the correlations between model parameters by calculating the weighted Pearson correlation coefficient r pearson based on the grid generated during the fitting process.The result is shown in Figure 10 for the six free parameters describing the shock and the geometry.The shock parameters (v sh , n pre , B) show weak correlations between each other, but they all show strong correlations with one or more geometric parameters.This analysis suggests that, though not very strong, our model suffers from degeneracies caused by the projection effect and the uncertainties in the geometry of the gas region.We would like to point out that the degeneracies between model parameters are taken into account in the 1σ uncertainties of our best-fit model parameters, which are actually the weighted standard deviation of the marginalized distributions.

Shock profiles of electron temperature and density
The top-left panel of Figure 11 displays the profiles of both electron temperature T e (the red curve) and density n e (the blue curve) in the shock component as predicted by our best-fit model (Model 2).Plotted in the other panels are the emissivity profiles (solid lines) and the corresponding ion fraction profiles (dashed lines) of hydrogen (H), carbon (C), nitrogen (N), oxygen (O) and sulfur (S).As can be seen, the MaNGA FoV (indicated by the gray band) happens to capture the jump in density and temperature (and subsequently, the line emissivities) at r ∼ 3500 AU.
An alternative and more commonly-adopted way to measure temperatures and densities of ionized gas is to use a set of selected emission line ratios that are sensitive to n e or T e as diagnostics, based on a given atomic database.It is interesting to compare the T e and n e profiles predicted by our model with those calculated directly from emission line ratios.We use PyNeb (Luridiana et al. 2015) to calculate T e and n e for each spaxel within the MaNGA FoV based on our measurements of emission line ratios.For consistency, we use CHIANTI 8 (Del Zanna et al. 2015), the same atomic database as used in MAPPINGS V.In fact, we have explored alternative databases available in PyNeb, finding them to provide similar n e and T e .We use the line ratio of [Sii] 6731/6717 as the diagnostic of n e , and use five different line ratios as diagnostics of T e : [Nii] 5755/6548, [Oiii] 4363/(4959+5007), [Oii] (7320+7330)/(3726+3729), [Sii] (4076+4069)/(6731+6716), and [Ci] 8729/9852.For a given pair of T e -sensitive and n e -sensitive diagnostics, we obtain the measurements of T e and n e simultaneously for each spaxel in the MaNGA FoV.Therefore, five combinations of T e -sensitive line ratios and the n esensitive line ratio give rise to a set of five measurements of the T e and n e maps, which are shown in the top and middle rows in Figure 12.
Each combination utilizes different T e -sensitive line ratios while keeping the n e -sensitive line ratio fixed as [Sii] 6731/6717.The alternative n e -sensitive line ratio [Oii] 3726/3729 is not utilized due to its significant uncertainty, as illustrated in the second panel of Figure 7.This uncertainty obscures the true density structure and generates a misleading clumpy artificial density distribution.The primary cause of this significant uncertainty stems from that the [Oii] λλ 3726,3729 doublets is not well resolved in the MaNGA spectra, as can be seen in the first two panels of Figure 2. Therefore, using [Oii] 3726/3729 would risk over-interpreting the data.However, it would be worthwhile to investigate the consistency between the density measured by [Oii] 3726/3729 and [Sii] 6731/6717 when higher quality data becomes available.It is important to note that, although n e is mainly constrained by [Sii] 6731/6717, the different T esensitive lines lead to slight differences in the derived n e .This is why the n e maps in the first row of Figure 12 differ quantitatively, although they qualitatively show good agreement.The bottom row in the same figure shows the corresponding χ 2 maps.As can be seen, the χ 2 values are close to zero in most spaxels of a given map, indicating that all the lines relevant for measuring T e and n e are well-fitted.We note that, we have added back the over-subtracted sky flux when calculating [SII] λλ 6717,6731, although this correction has no significant effect on our results.From Figure 12, all the density maps present a prominent jump from the upper-left low density region (∼ 1 cm −3 ) to the lower-right high density region (∼ 100 cm −3 ), and the boundary of the two regions coincides with the location of the lane-like structure in the surface brightness maps of the [Oiii] λ5007-like group.This result qualitatively agrees with the predicted density profile in Figure 11.Unlike n e , the measurements of T e are different, depending on what temperature diagnostics are used.In this case, only [Oii] (7320+7330)/(3726+3729) captures the overall trend of the predicted temperature profile in Figure 11.However, we should note that it is not suitable to directly compare the profiles shown in Figure 11 with the measurements in Figure 12, as the former are the intrinsic shock profiles without the Gaussian component and observational effects (e.g.limited spatial resolution and projection effect) in the real data.
In Figure 13 we compare the observed maps of the diagnosing line ratios (panels in the top row) with the predicted maps of the same line ratios from Model 2 in which we have included both the Gaussian component and the observational effects.Overall, our model can well reproduce the distribution of all the line ratios (the second row), although with subtle structures remaining in the maps of absolute residual (the third row) and relative residual (the bottom row).This result demonstrates that the discrepancies in T e as measured from different line ratios are caused by observational effects as well as the contamination of non-shock components such as the Gaussian component in our case.It appears that different line ratios can be affected by these effects in different ways, thus leading to inconsistent measurements.We have examined the effect of the Gaussian component and the observational effects separately, and we find both effects to have contributed significantly to the total discrepancies.Our analysis thus strongly warns that electronic density and temperature measurements based on observed emission line ratios can be seriously biased due to observational effects and non-shock components, and should be taken with caution.If assumed to be a part of the Orion-Eridanus superbubble, CCN should be at a distance ranging from 150 pc to 400 pc, corresponding to the close and far ends of the superbubble (Pon et al. 2016;Kounkel et al. 2017;Zari et al. 2017).Using the three-dimensional extinction map constructed by Vergely et al. (2022) 2 , we find the differential extinction along the line of sight to CCN to present two peaks, located at 150 pc and 475 pc, respectively.The two peaks coincide with the two edges of the Orion-Eridanus superbubble, and the CCN is likely to be associated with one of the peaks.In addition, we have attempted to constrain the distance to CCN with our data, by including the distance as a free parameter in Model 2, D CCN .We obtain a best-fit distance of D CCN = 170 ± 94 pc.This value is still consistent with the values from the first two considerations, given the large uncertainties.This result implies that the distance cannot be well constrained given the currently available data.Fortunately, we found that both the B and v sh have a weak correlation with the distance of CCN.On the other hand, although n pre has a stronger correlation with D CCN , the marginalized uncertainty of n pre only increases from 1.4 ± 0.2 cm −3 to 1.4 ± 0.3 cm −3 .The geometry parameters are indeed correlated with the distance, but our test shows that the data still prefer a small inclination angle, with models of i < 30 deg showing a larger Bayesian evidence than those of i ≥ 30 deg.

The dust extinction in CCN
The dust extinction in CCN is also uncertain.An upper limit of the dust extinction along the line of sight towards CCN can be estimated by integrating the differential extinction from Vergely et al. (2022) over the distance range from 0 up to 475 pc.This gives rise to A 5500 Å = 0.14 mag, corresponding to a color excess of E(B − V) < 0.05 mag assuming the Milky Way extinction curve of Cardelli et al. (1989) with R V = 3.1.This is smaller than the average value estimated from Hα λ6563/Hβ λ4861 in the MaNGA FoV, E(B − V) = 0.12 mag, using the electronic temperature for each spaxel as constrained by our best-fit model and assuming the Case B recombination coefficient.We find that, however, different Balmer line ratios lead to inconsistent E(B − V) values, implying that the extinction curve of Cardelli et al. (1989) may not be applicable for CCN, or that the Case B assumption is not appropriate in shocked regions (e.g.see Figure 4 of Dopita & Sutherland 2017).In fact, our best-fit shock model does predict higher Hα λ6563/Hβ λ4861 ratios than that of Case B at fixed temperature.We have also attempted to include the dust extinction in our model.We find no significant improvement in goodness-of-fit, and the null hypothesis of E(B − V) = 0 can not be rejected at a 3σ level.In summary, though without direct measurements of the dust content, constraints from different aspects consistently suggest that the dust extinction in CCN is rather limited.Our modelling presented in the previous section should have not been affected significantly by dust extinction.

4.3.
A single shock interacting with a cold gas cloud, or multiple shocks?
The MaNGA data and our best-fit model suggest that the CCN region observed by MaNGA is ionized by a shock which moves across the MaNGA FoV from the lower-right corner to the upper-left corner.This direction is opposite to the shock direction indicated by the protruding double-arc structure of the whole CCN, which goes from left to right as suggested by Z97 and T07.A natural explanation on this discrepancy is that CCN is a shocked cloud of cold gas (atomic gas, or molecular gas, or both).In this case, the overall curvature in the emission line surface brightness maps traces the shape of the gas cloud, but not the shock front.The protruding arcs to the right are actually the right edge of the cloud, which has collided head-on with the shock wave coming from the right.The shock has already swept into but not through the cloud, and consequently the shock front is found on the left side of the arc and still within the cloud.In a recent theoretical study by Kupilas et al. (2022), the simulation of a shocked Hi cloud produces a double-arc structure which is similar to CCN (see their Figure 4).Accordingly, one would expect to detect Hi on the left side of the double arc.We have examined the Hi data of the Galactic All Sky Survey (GASS; McClure-Griffiths et al. 2009;Kalberla et al. 2010), finding no Hi enhancement at the position of CCN.Considering the low spatial resolution of GASS (14.4 ′ ), however, the Hi cloud may suffer from smearing effect and so the scenario of a shocked Hi cloud cannot be firmly ruled out with the current data.As abovementioned, we have also used the Purple Mountain Observatory (PMO) 13.7-m millimeter telescope to observe the 12 CO(1-0) line in CCN (subsection 2.2).This observation provides an upper limit of the H 2 surface density, which is 0.6 M ⊙ /pc 2 if we assume a conversion factor of α CO = 4.35 (M ⊙ /pc 2 )/(K km/s) and a line width of 10 km/s.This result largely rules out the scenario of a shocked molecular cloud, consistent with the simulation in Kupilas et al. (2022) where the CCN-like double arc can not be formed in a shocked H 2 cloud (see their Figure 9).
The double-arc structure of the CCN could also be produced by multiple shocks, rather than the interaction between a cold gas cloud and a single shock as discussed above.In fact, the idea of multiple shocks has its own supporting evidences, including significant differences in the line ratio distributions in/outside the "clean region" in Figure 6, residuals of the best-fit model in Figure 9, the double-arc structure of the whole CCN, and the vari- Both IFU observations covering the whole CCN and high-resolution Hi observations would be needed if one were to discriminate the scenario of gas cloud with a single shock and the picture of multiple shocks.

Source of the shock
No matter which scenario is at work, as Z97 and T07 suggested, the source of the shock is most likely to be a supernova, which may be the one responsible for the Orion-Eridanus superbubble.In this case, however, the shock direction should go from left to right, thus in disfavor of the Hi cloud scenario.In addition, we failed to find any candidates of source supernova after having searched known supernovae and pulsars around CCN. Thus, the source of CCN is still an open question.

The Gaussian component
The physical origin of the Gaussian component is not immediately clear from our results.This component is possibly caused by another resolved shock.We have attempted to fit the residual of the shock component in Model 2 with another shock moving in the same direction, i.e. from lower-right to upper-left in the MaNGA FoV.We find that, however, some emission lines such as [Sii] λ6717 and [Nii] λ6583 cannot be fitted by this model.We have also tried with a shock direction from upper-right to lower-left, which fails to fit the lines in Group II.Another possible case is that the second shock propagates along the line of sight.We thus use an integrated shock model (the default output of the MAPPINGS V) to fit the line ratios in the Gaussian component.We find that a shock with n pre = 0.1 cm −3 , v sh = 450 km/s, and B = 0.1 µG can explain all the line ratios except [Sii]λ6717/Hα and [Sii]λ4069/Hα, as can be seen from Figure 14.Given the high velocity of this shock, we would expect the velocity field at the location of the Gaussian component to differ significantly from the other regions.However, this is not observed in the MaNGA FoV (see Figure 4).Therefore, another shock is unlikely to be the origin of the Gaussian component.Alternatively, the Gaussian component could be a clump of gas near the shocked region.In this scenario, the gas in this clump can be ionized by photons from shocked region.The most likely candidate that can be modeled with MAPPINGS V is the precursor of the shock in Model 2. Although this gas clump may have a different geometry than the precursor, most of the line ratios of the Gaussian component agree with the line ratio profiles of the shock precursor as predicted by Model 2 with a precusor-to-shock front distance of 2-4 pc, as shown in Figure 15.Two of the line ratios, [Ci] λ9852/Hα λ6563 and [Oiii] 4363/Hβ λ4861 show discrepancies between the data and the model, which could be attributed to a variety of reasons, e.g. the homogenous gas assumption adopted in MAPPINGS V, differences in the geometry of this gas clump compared to the precursor component in MAPPINGS V, differences in the metal abundance of the gas compared to our assumptions, or the contribution of photons from the shocked region outside the MaNGA FoV.However, investigation of these possibilities is beyond the scope of the current paper.Considering that the majority of the line ratios can be explained with no additional free parameters, we suggest that the Gaussian component is very likely a clump of gas ionized by photons from the shocked region, thus physically related to the shocked region.Further studies in future are needed to pin down this scenario.

SUMMARY
In this work, we have analyzed the integral field spectroscopy data for a small part of the Criss Cross Nebula (CCN), which was observed as a filler target by the MaNGA survey (Figure 1).At a distance of ∼ 150−475 pc, the spatial resolution of MaNGA (∼ 2.5 ′′ ) corresponds to a physical scale less than 1000 AU.Within the MaNGA FoV, a hexagon covering ∼ 20 ′′ × 20 ′′ , we have obtained the surface brightness map for 34 emission lines in the optical band, by fitting Gaussian profiles to the observed spectrum of each spaxel (Table 1 and Figure 2).We find the emission lines can be broadly divided into three groups according to their surface brightness maps (Figure 3): (1) the [Oiii] λ5007-like group including seven high-ionization lines and two [Oii] auroral lines which uniformly present a remarkable lane of high brightness in the upper-left part of the MaNGA FoV, (2) the Hα λ6563-like group including 23 lowionization or recombination lines which mostly present a clump of high brightness in the lower-right part of the MaNGA FoV, and (3) the third group including only two lines, [Oii] λ3726 and [Oii] λ3729, which show clumplike structures at positions of both the [Oiii] λ5007 lane and the Hα λ6563 clump.The velocity maps of different lines in a given group show high similarities, but vary from group to group (Figure 4).The three categories of microstructures indicate both regularity and complexity in the CCN.
We assume the lane structure in the [Oiii] λ5007-like group is produced by a shock which is moving across the CCN from the lower-right to upper-left direction.We use the measurements of emission line brightness and surface brightness ratios between different lines to constrain both physical parameters and geometry parameters of the shock model.We have examined two models: a simple model with a single shock plus a constant background (Model 1; Equation 2), and a more complex model with an additional Gaussian component (Model 2; Equation 3).We show that the MaNGA data can be reasonably well-fitted with Model 2 (Figure 8; Figure 9).The physical counterpart of the Gaussian component is likely another clump of gas ionized by photons from the shocked region (Figure 15), though more data is needed to further confirm this scenario.Our best-fit model parameters indicate that the shock has a pre-shock hydrogen density of n pre = 1.4 ± 0.2 cm −3 , a velocity of 133 ± 5 km/s, and a transverse magnetic field of B = 1.6 ± 0.3 µG, with the shock direction orientated to the line of sight by an inclination angle of 19 ± 14 deg (Table 2; Figure 5).
Finally, we have compared the electronic density and temperature profiles as predicted by Model 2 (Figure 11) with those calculated directly using observed emission line ratios as commonly-done in the literature (Figure 12).We find different line ratios to give rise to similar density maps but significantly inconsistent temperature maps.The discrepancies can be attributed to observational effects caused by the limited spatial resolution of the data and the projection of the shock geometry, as well as the contamination of the additional Gaussian component.After including all these effects, our model can well reproduce all the emission line ratios that are relevant for calculating density and temperature (Figure 13).Thus, electronic densities and temperatures measured from observed line ratios should be taken with caution (see also Cameron et al. 2022 for a similar caveat).
The CCN is an ideal lab for studying microstructures in shock regions and testing existing shock models.Considering its nearby distance, relatively high surface brightness and unique morphology, CCN will be a good template for studying similar nebulae in greater distance, e.g.those being/to be observed by the Local Volume Mapper of the SDSS V project (Konidaris et al. 2020) and the AMAZE project (Yan et al. 2020).For the CCN itself, as mentioned, future studies would benefit from high-resolution Hi and IFS observations, both covering the whole CCN.In fact, we have recently done both observations, using MeerKAT and SITELLE on CFHT, and we expect to come back and present more scientific results soon.

Figure 1 .
Figure 1.Upper: Criss-Cross Nebula (CCN) as seen at different scales.The left two images are constructed using narrowband images of [Sii] λλ6717, 6731 and Hα λ6563 from Temporin et al. (2007) and [Oiii] λλ4959, 5007 obtained in this work.The rightmost image is based on MaNGA data of the same lines.Lower: the MaNGA spectrum of the central spaxel in the MaNGA field of view (marked as a black cross in the upper-right image).The short black lines below the spectrum indicate the position of 34 emission lines that are measured and studied in this work.Some of the lines are more closely shown in the insets.areconfirmed in the sky model including [Oii] λλ3726, 3729 and [Sii] λλ6717, 6731, whose redshifts are close to what we measured from the DRP datacube.We notice that the over-subtracted flux is typically one order of magnitude smaller than the flux measured from the DRP spectrum, and so the flux correction is generally a tiny effect on our results.
λ3726 and [Oii] λ3729.d These two lines are only marginally resolved under MaNGA resolution.e It is hard to find a clean left and right window for these lines separately, and [Sii] λ4076 is sometimes too weak.f This configuration aims to alleviate the confusion effect of [Caii] λ7324 on [Oii] λλ7320, 7330.gThis configuration aims for more robust results when [Cii] λ9826 is weak.

Figure 2 .
Figure 2. Fitting results for the 34 emission lines of one example spaxel.The grey histogram is the observed spectrum with baseline subtracted, and the shaded grey region is the 1 σ error.The solid red line is the best-fit model, with the line center marked by the dashed vertical line.

]Figure 3 .
Figure3.The surface brightness map of all the 34 emission lines measured within the MaNGA FoV in the unit of 10 38 erg/s/kpc 2 .The solid and dashed gray contours represent the 75%, 90% (bolded), 95%, and 99% quantiles in the surface brightness map of [Oiii] λ5007 and Hα λ6563 lines, respectively, which are repeated in every panel for reference.

Figure 4 .
Figure 4.The velocity map of eight emission lines in the unit of km/s as measured from the MaNGA data.The bolded contours show 90% quantiles of [Oiii] λ5007 (solid line) and Hα λ6563 (dashed line) surface brightness.flowsolution of one-dimensional radiative shock models with parallel-slab geometry.The cooling function is established based on the CHIANTI 8 atomic database(Del Zanna et al. 2015), and the pre-ionization effect is addressed via an iterative approach.Other than the usually used Rankine-Hugoniot flow equations, in MAPPINGS V a unique choosing of integral form of the conservation laws is employed to mitigate round-off errors.The reader is referred toSutherland & Dopita (2017)  for more details about MAPPINGS V.To compare with the resolved MaNGA data, we use shock profiles calculated by the model, i.e. the physical properties and emissivity as functions of the separation r from the shock front, instead of the integrated spectrum which is the default output of the model.In this case, the shock front locates at r = 0, the region with r > 0 has already been shocked (thus referred to as "shocked region"), and the region with r < 0 is expected to be shocked in future (thus referred to as "precursor").The precursor may have been ionized by the photons escaping from the shocked region.Since MAPPINGS V adopts parallel-slab geometry, in the direction perpendicular to the shock direction the physical properties and emissivity are homogenous in the model.Therefore, the model is essentially a one-dimensional model, and we need to simplify the geometry of the observed region in consideration in order for appropriate comparisons between model and data.Our geometry model is illustrated in Figure5.We assume the region is a "slab" of ionized gas with a certain thickness, oriented with an inclination angle of i (the angle between the "slab" and the line of sight).As a result, the gas slab has a line-of-

Figure 5 .
Figure 5.A sketch of the shock model adopted in our work.The shock is assumed to have a parallel-slab geometry with a depth of D los along the line of sight, and is orientated with an inclination angle of i.The vertical purple line marks the shock front, and the orange arrow in left side indicates the shock direction.

Figure 8 .Figure 9 .
Figure 8. Surface density maps of 12 emission lines (panels from top to bottom).For each line, panels from left to right show the observed map, the predicted map by Model 1, and the absolute and relative residual maps.All the maps are color-coded by surface brightness in units of 10 38 erg/s/kpc 2 .The gray contours show 90% quantiles of [Oiii] λ5007 (solid line) and Hα λ6563 (dashed line) surface brightness.

Figure 11 .Figure 12 .
Figure 11.In the top left panel, we present the intrinsic shock profiles of the electronic density (ne) in blue and the temperature (Te) in red, as predicted by Model 2. The remaining panels display the ion fractions of hydrogen (H), carbon (C), nitrogen (N), oxygen (O), and sulfur (S), as well as the emissivity of the corresponding emission lines.These panels are arranged from left to right and top to bottom.The shaded region in each panel represents the coverage of the MaNGA FoV.

Figure 13 .
Figure13.Shown in the top panels are the observed surface density maps of six emission line ratios that are used to measure electronic density and temperature maps in the previous figure (Figure12).The predicted maps of the same lines by Model 2 are shown in the panels in the second row.The third and last rows show the absolute and relative residual maps.

Figure 15 .
Figure 15.The line ratio profile of the precursor of Model 2 (solid line) is compared to the fitted line ratio of the Gaussian component (represented by the horizontal dashed line).The colorful shaded region represents the 1 σ uncertainty of the fitted line ratio of the Gaussian component.Additionally, the grey shaded region indicates the coverage of the MaNGA FoV.

Table 1 .
All the lines fitted in MaNGA data

Table 2 .
The best-fit parameters of shock models with 1 σ uncertainty The residual of the best-fit integrated shock model obtained by comparing it with the line ratio of the Gaussian component.The error bar represents the 1 σ uncertainty of the line ratio of the Gaussian component.ation of shock velocity in different regions of CCN (e.g.133 km/s in the MaNGA FoV studied in our work versus 40 km/s in the region studied by Z97).It's likely that there are two or even more shocks with different orientations in the CCN region.