<figure><img src="image_1.png" /></figure>

HAL
open science

A Detailed Temperature Map of the Archetypal
Protostellar Shocks in L1157

S. Feng, H. B. Liu, P. Caselli, A. Burkhardt, F. Du, R. Bachiller, C. Codella,
C. Ceccarelli

► To cite this version:

S. Feng, H. B. Liu, P. Caselli, A. Burkhardt, F. Du, et al. A Detailed Temperature Map of the Archetypal Protostellar Shocks in L1157. The Astrophysical journal letters, 2022, 933, 10.3847/2041-8213/ac75d7. insu-03860331

HAL Id: insu-03860331

https://insu.hal.science/insu-03860331

Submitted on 18 Nov 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# A Detailed Temperature Map of the Archetypal Protostellar Shocks in L1157

S. Feng (冯思轶)<fnref n="1" /> ID, H. B. Liu<fnref n="2" /> ID, P. Caselli<fnref n="3" /> ID, A. Burkhardt<fnref n="4" /> ID, F. Du<fnref n="5" />,<fnref n="6" /> ID, R. Bachiller<fnref n="7" /> ID, C. Codella<fnref n="8" />,<fnref n="9" /> ID, and C. Ceccarelli<fnref n="9" /> ID

Received 2022 April 4; revised 2022 June 3; accepted 2022 June 3; published 2022 July 11

## Abstract

We present sensitive NH₃ (1,1)–(7,7) line images from the Karl G. Jansky Very Large Array toward successive shocks, which are associated with the blueshifted outflow lobe driven by the compact protobinary system L1157. Within a projection distance of 0.1 pc, our observations not only trace the quiescent and cold gas in the flattened envelope, but also illustrate the complex physical and chemical processes that take place where the high-velocity jet impinges on its surrounding medium. Specifically, the NH₃ ortho-to-para ratio is enhanced by a factor of 2–2.5 along the jet path, where the velocity offset between the line peak and the blueshifted wing reaches values as high as 10 km s⁻<fnref n="1" />; it also shows a strong spatial correlation with the NH₃ column density, which is enhanced to >10¹<fnref n="6" /> cm⁻<fnref n="2" /> toward the shock cavities. At a linear resolution of 1500 au, our refined temperature map from the seven NH₃ lines shows a gradient from the warm B0 eastern cavity wall (>120 K) to the cool cavity B1 and the earlier shock B2 (<80 K), indicating shock heating.

*Unified Astronomy Thesaurus concepts: Low mass stars (2050); Stellar bow shocks (1586); Astrochemistry (75)*

## 1. Introduction

Protostellar shocks are commonly detected in the earliest stage of star formation. Generated by the impact of supersonic jets from the central protostars on their dense natal parental cloud (e.g., Frank et al. 2014), the shocks create a dense and warm environment in a short timescale by compressing and heating the surrounding gas. Such an environment speeds up chemical processes, which are infeasible in the preshock gas, including, but not limited to, endothermic chemical reactions, ice mantle vaporization, and sputtering (e.g., Viti et al. 2011).

The target region of this work is an archetypal region containing successive shocks (named B0, B1, and B2), which are associated with the blueshifted outflow lobe from the nearby (d ~ 352 pc; Zucker et al. 2019) Class 0 compact protobinary system L1157-mm (≈3L☉; Tobin et al. 2010, 2022). Kinematically, an episodic precessing jet driven by the protobinary system hits the shock cavity wall, generating bright knots (e.g., Gueth et al. 1996, 1998; Podio et al. 2016). Within a projection of 0.1 pc from the protobinary system in the plane of the sky, the kinematic age difference of these shock knots is ~1000 yr, so this region provides us with one of the best space laboratories to study the time-dependent shock chemistry (e.g., Lefloch et al. 2010; Codella et al. 2017).

Multiwavelength line surveys have intensively observed this region, especially B1, over the decades (e.g., Herschel-CHESS, Ceccarelli et al. 2010; IRAM 30 m-ASAI, Lefloch et al. 2018; NOrthern Extended Millimeter Array (NOEMA)-SOLIS,

Ceccarelli et al. 2017) and revealed chemical complexities from diatomic molecules to complex organics (e.g., Tafalla & Bachiller 1995; Bachiller et al. 2001; Benedettini et al. 2007; Arce et al. 2008; Codella et al. 2009, 2010; Lefloch et al. 2012; Benedettini et al. 2013; Busquet et al. 2014; Podio et al. 2014; Fontani et al. 2014a; Codella et al. 2015, 2017; Lefloch et al. 2017, 2018; Codella et al. 2020; Feng et al. 2020; Spezzano et al. 2020). It is clear that sufficient observational data exist to carry out a systematic study of shock chemistry, i.e., examining the origin, excitation, and chemical complexity of different species. However, a crucial input for chemical modeling is the physical structure of this region, which is still missing.

Previous observations attempted to use several inversion lines of NH₃ at centimeter wavelengths to provide constraints on the kinetic temperature of this region, given that these lines have been widely used as an interstellar thermometer<fnref n="10" /> for molecular gas with number density n > 10⁴ cm⁻<fnref n="3" /> (Ho & Townes 1983; Walmsley & Ungerechts 1983; Crapsi et al. 2007; Rosolowsky et al. 2008; Juvela & Ysard 2011; Caselli et al. 2017). However, hindered by the low spatial/velocity resolutions of single-dish point observations (e.g., Bachiller et al. 1993 observed the (1,1)–(4,4) lines; Umemoto et al. 1999 observed the (1,1)–(6,6) lines) and the insufficient sensitivity of interferometric observations (Tafalla & Bachiller 1995 observed only the (1,1)–(3,3) lines), several assumptions were

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Department of Astronomy, Xiamen University, Zengcuo'an West Road, Xiamen, 361005, People's Republic of China; syfeng@xmu.edu.cn</footnote>

<footnote n="2">² Academia Sinica Institute of Astronomy and Astrophysics, No.1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan, Republic of China</footnote>

<footnote n="3">³ Max-Planck-Institut für extraterrestrische Physik, Gießenbachstraße 1, D-85748, Garching bei München, Germany</footnote>

<footnote n="4">⁴ Wellesley College, 106 Central Street, Wellesley, MA 02481, USA</footnote>

<footnote n="5">⁵ Purple Mountain Observatory and Key Laboratory of Radio Astronomy, Chinese Academy of Sciences, Nanjing 210033, People's Republic of China</footnote>

<footnote n="6">⁶ School of Astronomy and Space Science, University of Science and Technology of China, Hefei, Anhui 230026, People's Republic of China</footnote>

<footnote n="7">⁷ Observatorio Astronomico Nacional (OAN, IGN), Calle Alfonso XII, 3, E-28014 Madrid, Spain</footnote>

<footnote n="8">⁸ INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125, Florence, Italy</footnote>

<footnote n="9">⁹ Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="10">¹⁰ The majority of the NH₃ population stays in the metastable levels *J* = *K*, where *J* is the total angular momentum quantum number of NH₃, and *K* is the projection of *J* on the rotational axis. The inversion transitions from different rotational ladders of NH₃ are coupled only collisionally. Having similar frequencies at 1.3 cm wavelength, several inversion lines of NH₃ can be observed simultaneously, so that the calibration uncertainties in the line ratio measurements are minimized.</footnote>
made in their works, leading to large uncertainties in their conclusions.

Interferometric observations have been carried out recently, targeting line ladders of CH<sub>3</sub>CN, H<sub>2</sub>CO, CS, and molecular ions. Although these works provide constraints on the H<sub>2</sub> volume density *n* (e.g., Benedettini et al. 2013; Gómez-Ruiz et al. 2015) and molecular rotation temperature (e.g., Codella et al. 2009; Podio et al. 2014) toward several knots, these attempts are only toward B1.

In this paper, with the newest NH<sub>3</sub> (1,1)–(7,7) observations at high-spatial and high-velocity resolutions, we provide a detailed temperature map over B0–B1–B2 and discuss the complex physical and chemical processes that take place when the high-velocity jet impinges on its surrounding medium.

## 2. Observations and Data Reduction

Using the Karl G. Jansky Very Large Array (JVLA), we have performed observations toward L1157 B0–B1–B2 at the K band in D-array configurations from September to November in 2018. For all observations, a common phase center at 20<sup>h</sup>39<sup>m</sup>10<sup>s</sup>2.7<sup>°</sup>01′10″5 (J2000) and a systemic velocity V<sub>sys</sub> ~ 2.7 km s<sup>-1</sup> were adopted from the Herschel results on the o-NH<sub>3</sub> (1<sub>o</sub>-0<sub>o</sub>) line (Codella et al. 2010). Employing the three-bit sampler, our correlator setup uses 27 independent spectral windows to cover the NH<sub>3</sub> (1,1)–(7,7) transitions (with E<sub>u</sub>/k<sub>B</sub> ranges from 24 to 539 K) and broadband continuum simultaneously.

All epochs of observations used 1331+305 (3C 286), J1642+3948, and J2022+6136 for absolute flux, passband, and complex grain calibrations, respectively.

Following the standard strategy by using the Common Astronomy Software Applications (CASA; McMullin et al. 2007) package release 5.4.0, we manually calibrated the JVLA data. The absolute fluxes of the flux calibrator 3C 286 were referenced from the Perley-Butler 2017 standards (Perley & Butler 2017). For different epochs, a 10%-15% nominal absolute flux calibration accuracy may have to be assumed according to the official documentation.<fnref n="11" /> Nevertheless, the absolute flux calibration errors are factored out when deriving spectral line ratios.

Using the CASA task TCLEAN in CASA 5.4.6, we performed the spectral line imaging (setting the specmode as “cube”) and broadband continuum imaging (setting the specmode as “mfs”) by applying the “multiscale” imaging option with scales values of 0, 5, and 15 times the pixel size. The primary beam (pb) and the maximum recoverable angular scales for the single pointing observations are ~119" and ~61" at 23.7 GHz, respectively. We test the Briggs robust weighting<sup>12</sup> *r* as 2.0, -2.0, and 0.5. The continuum shows >4σ emission only toward the protobinary system, although this system is not resolved from our observations. The continuum intensity peak, measured at different weightings before and after pb correction, is 5σ-10σ, with 1σ rms varying in the range of 0.01-0.04 mJy beam<sup>-1</sup>, which is consistent with Tobin et al. (2022). For all the targeted NH<sub>3</sub> lines, their image properties<sup>13</sup> are listed in Table 1.

## 3. Results

### 3.1. Molecular Spatial Distribution

The hyperfine multiplets of NH<sub>3</sub> (1,1)–(7,7) cover a velocity offset from -30 to +30 km s<sup>-1</sup> with respect to the V<sub>sys</sub> (e.g., Ho & Townes 1983; Mangum & Shirley 2015). We detected all the seven lines toward our target region, with their hyperfine multiplets partially resolvable at a velocity resolution of 0.05–0.2 km s<sup>-1</sup> (Table 1).

Extracting their naturally weighted spectra toward different jet knots from B0 to B2 (Figure A1), we found that the satellite components of each line have a signal-to-noise ratio S/N > 3 in the velocity range of -25 to +25 km s<sup>-1</sup>. Over this velocity range, the integrated intensities of all seven lines trace an extended structure (Figures 1 and B1)—from the rim of the B0 eastern cavity wall to two bow-shock cavities associated with B1 and B2 (denoted as C2 and C1, respectively, in Gueth et al. 1998).

Comparing the total fluxes of the (1,1)–(4,4) lines from our JVLA observations with those given by the single-dish observations (Bachiller et al. 1993; Umemoto et al. 1999)<fnref n="14" /> we found 20% and 15% of the flux missing toward B1 and B2, respectively. Note that the angular resolution of Effelsberg observations (Bachiller et al. 1993) is coarser than that of JVLA by a factor of 10, so effects such as pointing errors and the attenuation at the edge of the Gaussian-shaped Effelsberg beam, may bring uncertainty in the comparison. Although we cannot precisely recover the flux per pixel from the Effelsberg single-point data, such comparison indicates that the images in this work cover the NH<sub>3</sub> extended emissions toward B1 and B2 at a similar level.

Limited by the largest recoverable angular scale with one pointing in these observations, only the central region (i.e., within the gray circle in Figure 1) of the image has high fidelity. In order to apply the kinetic temperature map derived from NH<sub>3</sub> lines to the chemical property measurements (i.e., column density, abundance) of the other molecular lines at a spatial resolution of ~1500 au, we labeled seven clumpy substructures (B0e, B1a, B1b, B1c, B1f, B1i, and B2a), which were identified from previous observations (e.g., by using line emissions from CH<sub>3</sub>CN, HCN, CH<sub>3</sub>OH, SiO, CS, H<sub>2</sub>CO, HCO<sup>+</sup>, HNCO, SO, and SO<sub>2</sub> in Benedettini et al. 2007; Codella et al. 2009; Gómez-Ruiz et al. 2013; Burkhardt et al. 2016; Feng et al. 2020). Moreover, the large pb of JVLA allows us to investigate the complete B0 cavity wall and the B2 shock front together with the phase center B1. In the maps with pb correction (Figure B1), we found two clumpy structures showing stronger emissions than the rest of B0 and B2. Both structures are along the jet path derived by Podio et al. (2016) and were previously denoted B0a and B2b in Burkhardt et al. (2016).

Each substructure may be traced by more than one molecule, and the absolute coordinate of these intensity peaks (listed in Table A1) are different from tracer to tracer<fnref n="15" /> by 1"-2". Nevertheless, these nine clumpy substructures disentangled three chemical layers (e.g., Fontani et al. 2014b; Codella et al. 2017).

<footnote n="11">¹¹ https://science.nrao.edu/facilities/vla/docs/manuals/oss/performance/fdscale</footnote>

<footnote n="12">¹² The naturally weighted ($r = 2.0$) image has high sensitivity but low spatial resolution, the uniformly weighted ($r = -2.0$) image has low sensitivity but high spatial resolution, and the image with $r = 0.5$ is a good trade-off between resolution and sensitivity.</footnote>

<footnote n="13">¹³ Line image properties have negligible changes before and after continuum subtraction.</footnote>

<footnote n="14">¹⁴ The line emission peaks observed with Nobeyama 45 m (Umemoto et al. 1999) are in general lower than those observed with Effelsberg 100 m (Bachiller et al. 1993) by a factor of 1.5–2, probably due to the beam dilution.</footnote>

<footnote n="15">¹⁵ B2a corresponds to the emission peak of CO toward B2 in the present work, which is ~5" offset to the SO and SO<sub>2</sub> emission peak reported in Feng et al. (2020) and is ~2" offset to the coordinate reported in Burkhardt et al. (2016).</footnote>
<table><caption>**Table 1**
Spectroscopic Parameter of the Targeted NH₃ Lines and Their Image Properties from the JVLA Observations</caption><thead><tr><th rowspan="3">Mol.</th><th rowspan="3">Freq.<br>(GHz)</th><th rowspan="3">Transition</th><th rowspan="3">Sμ<sup>2a</sup><br>(D<sup>2</sup>)</th><th rowspan="3">E<sub>u</sub>/k<sub>B</sub><sup>b</sup><br>(K)</th><th rowspan="3">No. hfs<sup>c</sup></th><th rowspan="3">ΔV<sup>d</sup><br>(km s<sup>-1</sup>)</th><th colspan="3">Δθ (P.A.)<sup>e</sup></th><th colspan="3">1σ rms<sup>f</sup></th></tr><tr><th colspan="3">(" × ",°)</th><th colspan="3">(mJy beam<sup>-1</sup> ch<sup>-1</sup>)</th></tr><tr><th>r = 2.0</th><th>r = 0.5</th><th>r = -2.0</th><th>r = 2.0</th><th>r = 0.5</th><th>r = -2.0</th></tr></thead><tbody><tr><td>p-NH<sub>3</sub></td><td>23.694</td><td>1(1)0a-1(1)0s</td><td>6.6</td><td>24</td><td>18<sup>g</sup></td><td>0.049</td><td>4"63 × 3"77 (17°)</td><td>4"24 × 3"19 (26°)</td><td>3"71 × 2"63 (31°)</td><td>6.3</td><td>5.1</td><td>6.1</td></tr><tr><td>p-NH<sub>3</sub></td><td>23.723</td><td>2(2)0a-2(2)0s</td><td>14.7</td><td>65</td><td>21<sup>g</sup></td><td>0.197</td><td>4"90 × 3"81 (19°)</td><td>4"60 × 3"36 (23°)</td><td>4"34 × 2"93 (29°)</td><td>2.2</td><td>2.5</td><td>2.3</td></tr><tr><td>o-NH<sub>3</sub></td><td>23.870</td><td>3(3)0a-3(3)0s</td><td>46.4</td><td>124</td><td>26<sup>h</sup></td><td>0.196</td><td>4"20 × 3"60 (−6°)</td><td>3"67 × 3"01 (2°)</td><td>3"20 × 2"46 (10°)</td><td>2.8</td><td>2.9</td><td>2.8</td></tr><tr><td>p-NH<sub>3</sub></td><td>24.139</td><td>4(4)0a-4(4)0s</td><td>31.8</td><td>201</td><td>7<sup>i</sup></td><td>0.194</td><td>4"12 × 3"53 (−1°)</td><td>3"62 × 2"97 (3°)</td><td>3"15 × 2"43 (10°)</td><td>1.3</td><td>1.3</td><td>2.2</td></tr><tr><td>p-NH<sub>3</sub></td><td>24.533</td><td>5(5)0a-5(5)0s</td><td>40.5</td><td>296</td><td>7<sup>i</sup></td><td>0.191</td><td>4"25 × 3"65 (0°)</td><td>3"74 × 3"05 (10°)</td><td>3"31 × 2"48 (18°)</td><td>1.2</td><td>1.3</td><td>2.2</td></tr><tr><td>o-NH<sub>3</sub></td><td>25.056</td><td>6(6)0a-6(6)0s</td><td>98.5</td><td>409</td><td>7<sup>i</sup></td><td>0.187</td><td>4"50 × 3"65 (20°)</td><td>4"09 × 3"10 (26°)</td><td>3"62 × 2"51 (28°)</td><td>1.3</td><td>1.3</td><td>2.3</td></tr><tr><td>p-NH<sub>3</sub></td><td>25.715</td><td>7(7)0a-7(7)0s</td><td>56.3</td><td>539</td><td>7<sup>i</sup></td><td>0.182</td><td>4"09 × 3"46 (2°)</td><td>3"61 × 2"90 (8°)</td><td>3"19 × 2"34 (16°)</td><td>1.1</td><td>1.2</td><td>2.2</td></tr></tbody></table>

**Notes.**

a Sum of all hyperfine multiplets.

b Higher upper energy level of the transition.

c Number of hyperfine multiplets used for deriving the optical depth of the main component.

d Channel width.

e Synthesized beam.

f Without pb correction (with uniform noise level over the entire map), without angular resolution smoothing, measured with the referred beam size ("beam") per channel ("ch").

g Number of hyperfine multiplets, adopted from Mangum & Shirley (2015).

h Number of hyperfine multiplets, adopted from PySpecKit (Ginsburg & Mirocha 2011).

i Number of hyperfine multiplets, adopted from TopModel in "Splatalogue" database (https://www.cv.nrao.edu/php/splat).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The velocity structure from the cavity wall B0 to shocks B1 and B2 (the same for both cases with or without pb correction). The three left columns show the velocity structure maps of the (2,2) and (3,3) lines. The background in each panel shows the velocity map $V_p$, where the line intensity peaks (the first column), $V_b$, where the blueshifted wing extends to (the second column), and the peak-to-bluest velocity offset $V_p - V_b$ (the third column). All maps are naturally weighted. The black contours are the same continuum emissions as those in Figure 1. The purple contour indicates the $5\sigma$ emission from the (1,1) integrated intensity toward the flattened envelope around the protobinary system, which is perpendicular to the outflow direction. The gray dashed contours and the labeled clumpy substructures are the same as shown in Figure 1, and the synthesized beam is plotted in the bottom-left corner. The pb-corrected line emissions with $< 5\sigma$ rms are blank in each panel. The black curve guides the path of the precessing jet from L1157-mm as modeled by Podio et al. (2016). The last column shows the (2,2) and (3,3) line profiles extracted from clumpy substructures along the jet path. The gray dashed vertical line indicates the system velocity of the cloud $V_{sys} = 2.7$ km s$^{-1}$, and the red dashed vertical line indicates the velocity where the main component of each line blueward $V_{lsr} = -10$ km s$^{-1}$. The line intensity Jy beam$^{-1}$ is converted to the brightness temperature $T_B$, and the horizontal cyan line indicates the baseline ($T_B = 0$ K). A narrow feature appears toward the protobinary system mm on the spectrum, probably tracing the quiescent and cold gas in the flattened envelope.</figcaption></figure>

deceleration effect of the underlying jet or wind after traveling large distances (Zhang et al. 2000). In each shock, the pattern of blueshifted wings becoming broader with the distance from the

outflow source indicates the “Hubble-law”, which is modelled as a consequence of the bow shock (e.g., Smith et al. 1997; Downes & Ray 1999).
When comparing the peak-to-bluest velocity offset ($V_p - V_b$) in the line of sight, we note that a high offset of 10 km s⁻¹ is present along the curved jet path.

## 4. Analysis and Discussion

To determine the temperature and density profile over a region in the ISM, the rotational diagram (RD) method was previously widely used with two assumptions: (1) All lines are in the local thermal equilibrium (LTE) condition; (2) the lines are optically thin or the optical depth is known to allow correction. However, our test of the RD method shows that significant blueshifted emissions on the line profiles toward this shocked region lead to large parameter uncertainties (see Appendix C for details).

Instead, with a proper assumption of the source geometry, the large velocity gradient (LVG) escape probability approximation is preferred to constrain the gas properties toward a particular location, when the collision rates of a particular species are known (Sobolev 1960; Goldreich & Scoville 1976).

### 4.1. LVG Approximation

Using the non-LTE statistical equilibrium radiative transfer code RADEX (van der Tak et al. 2007) and a related solver (Fujun Du's myRadex),<fnref n="18" /> we apply the MultiNest Algorithm (Feroz & Hobson 2008; Feroz et al. 2009, 2019) and derive the probability density function (PDF) of the kinetic temperature $T_{\text{kin}}$, the H₂ volume density $n$, and the column density $N_T$ for p-NH₃ and o-NH₃ toward all pixels.

Radiative and nonreactive collisional transitions in the gas phase do not change the molecular spin orientations. Therefore, we treat *o*/p-NH₃ as distinct species and assume that the transitions between them are forbidden (interconversion processes between *o*- and *p*-transitions are ignored). In this work, the collision rates are relative to the sum of all the hyperfine components and between NH₃ (1,1)–(6,6) and *p*-H₂, which are provided by the Leiden Atomic and Molecular Database (LAMDA, Schöier et al. 2005).<fnref n="19" />

All the pb-corrected images are smoothed to the same angular resolution (i.e., 4′′90 × 3′′81) and the same pixel size (i.e., 0′′5). The intensities of the lines are integrated over the velocity range of −25 to +25 km s⁻¹ by assuming a single velocity component. The geometry is set as “LVG” (“Slab” was tested and no appreciable differences were shown). The FWHM line widths of all lines are adopted as 4 km s⁻¹ based on the median from all spectrum fittings (see Table A2). In the case where multivelocity components are smeared within the VLA’s synthesized beam, an FWHM line width of 9 km s⁻¹ is adopted from Bachiller et al. (1993) and Lefloch et al. (2012) for the test. Two different values have been considered for the beam-filling factor: (i) unity for all lines; and (ii) 0.5 or 0.1 for transitions higher than (3,3). The best-fit parameters are searched over relatively large ranges for the H₂ number density $n$ (10³–10⁷ cm⁻³), $T_{\text{kin}}$ (5–10⁴ K), and $N_T$ (10¹⁴–10¹⁷ cm⁻²).

An example of the PDF from the LVG approximation is shown toward B1a in Figure D1, and examples of parameter maps are given in Figure D2. With different line-width and beam-filling factor combinations, the best fits toward all nine

clumps are listed in Table D1. Although degeneracy among the above three parameters is inevitable, we found that the $T_{\text{kin}}$ and the $N_T$ for *o*- and *p*-NH₃ are well constrained toward these representative positions, with uncertainties less than 30%.

From Figure D1 we also note that the *n* of each clump cannot be constrained from the LVG fittings, and thus, varying its value in the aforementioned range does not change the best-fit result of $T_{\text{kin}}$ and the $N_T$. In this warm shocked environment, the critical densities of all seven NH₃ lines are in the range of 5 × 10³–3 × 10⁴ cm⁻³ (estimated by using the Einstein coefficients and the collision rates provided by LAMDA in the temperature range of 50–300 K, which are consistent with the effective critical density given by, e.g., Shirley 2015). Previous CO and CS observations from Lefloch et al. (2012), Benedettini et al. (2013), and Gómez-Ruiz et al. (2015) indicate $n > 10^4$ cm⁻³ at an angular resolution in the range of 3′′–20′′, which should validate the LTE assumption for the NH₃ lines. Because this number density cannot be confirmed from our LVG fittings, the RD method would need to be applied with caution.

Our tests show that a smaller beam-filling factor for higher transitions result in better consistency of the $T_{\text{kin}}$ derived from *o*/p-NH₃ separately (Table D1). The collisional coefficients, the line width, and the integrated intensity of each line adopted in the fittings change the best-fit absolute value of $T_{\text{kin}}$ and $N_T$ toward individual pixels, though they do not change the contrast between pixels for both parameter maps (see Appendix D for details).

### 4.2. The Gas Kinetic Temperature Map

Assuming a line width of 4 km s⁻¹, a unity beam-filling factor for (1,1)–(3,3), and 0.1 for (4,4)–(7,7), Figure 3 shows the maps of gas kinetic temperature, the relative abundance ratio between *o*- and *p*-NH₃ and their total column density derived from the above-mentioned LVG best fit.

In a 0.1 pc scale field of view, good agreement on $T_{\text{kin}}$ derived from *o*/p-NH₃ separately (relative ratio ~1, shown in yellow in the first column of Figure 3, with <10% uncertainty) indicates the appropriateness of these assumptions.

Throughout the entire southern outflow shocks, the mean $T_{\text{kin}}$ map reveals an intrinsic gradient (second column of Figure 3): warm components (>120 K) appear toward the spots where the jet impinges on the eastern cavity walls (B0a–B0e–B1a), cool components (<80 K) toward the cavity B1b–B1c and the older shock B2.

Although the temperature gradient from B0 to B2 has been roughly indicated by previous observations (e.g., Umemoto et al. 1999; Tafalla & Bachiller 1995) and was proposed by the model (e.g., Lefloch et al. 2012; Podio et al. 2014 suggested a slightly larger temperature gradient from B0 with > 200 K to B1 with ~60 K and B2 with ~20 K), it has never been confirmed by a detailed map until this work. At a linear resolution of ~1500 au, our refined results indicate that the older components have experienced more postshock cooling.

### 4.3. The Column Density Map and the Ortho-to-para Ratio Map of NH₃

Given that no continuum at millimeter and centimeter wavelengths is detected toward B0–B1–B2 from the existing data and that the molecular abundance with respect to H₂

<footnote n="18">¹⁸ https://github.com/fjdu/myRadex</footnote>

<footnote n="19">¹⁹ No data are provided for the collision rate between NH₃ (7,7) and *p*-H₂ in LAMDA.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Kinetic temperature maps and column density maps for *o*- and *p*-NH₃, derived from LVG MultiNest best fitting, by assuming a single component (with an FWHM line width of 4 km s⁻¹) and two filling factors (1 for (1,1)–(3,3) lines and 0.1 for (4,4)–(7,7) lines) toward each position. Panels from left to right in the upper row show the relative ratio and the mean of *T*<sub>kin</sub> from *o*/*p*-NH₃, as well as the total column density and the OPR of NH₃. Panels in the lower row show the uncertainty with respect to the value of the above parameters. The pixels where the pb-corrected NH₃ (5,5) and (6,6) lines < 5σ emission and CO(1–0) < 8σ emission are blank. The gray dashed contours and the labeled clumpy substructures are the same as shown in Figure 1. The black curve guides the path of the precessing jet from L1157-5 mm as modeled by Podio et al. (2016). These maps are derived by smoothing all lines to the same pixel size and angular resolution, i.e., 4.99′′90 × 3.77′′81, shown as the synthesized beam in the bottom left of each panel.</figcaption></figure>

measured from different observations has a large uncertainty especially in the shocked region, the column density and ortho-to-para ratio (OPR) of NH₃ are more reliable physical indicators than its abundance with respect to H₂ in this work.

A U-shape structure appears on both the *o*- and *p*-NH₃ column density map toward B1, connecting the spots B1a–B1c–B1b, with the total value reaching >1 × 10¹⁶ cm⁻² to the eastern wall. Two U-shape structures also appear toward B2a and B2b, with the column density peak twisting toward the western wall (third column of Figure 3). These U-shapes are significant on the integrated intensity maps, especially for the lines observed with a larger pb.

Measuring the Spearman’s rank correlation<fnref n="20" /> coefficient ρ (Cohen 1988), the total column density shows moderate (0.4–0.5) spatial correlation with the velocity offset in the line of sight (*V*<sub>p</sub> − *V*<sub>b</sub>) toward B0 and B2 locally, as well as strong (0.63) correlation toward B1 locally (Figure D3). However, the correlation is weak (0.18) if treating the successively shocks as one entity. This could imply that the sputtering of ammonia off

<footnote n="20">²⁰ The Spearman rank correlation coefficient ρ is a nonparametric measure of statistical dependence between two variables. This coefficient can assess how well a monotonic function (no matter whether it is linear or not) can describe the relationship between two variables. The coefficient ρ is in the range from −1 (decreasing monotonic relation) to 1 (increasing monotonic relation), with zero indicating no correlation.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Left column: the $V_p$ map of the main hyperfine component of the (1,1) and (2,2) lines toward the flattened envelope around the marked protobinary system, which should elongate in the direction shown as the gray two-headed arrow (Chiang et al. 2010), perpendicular to the outflow. In each panel, the black contours indicate the integrated intensity (with pb correction), starting from $5\sigma$ and increasing by $3\sigma$. The gray dashed contours, starting from $4\sigma$ ($\sigma = 0.27$ Jy beam$^{-1}$) and increasing with the step of $5\sigma$, show the CO (1-0) emission (Gueth et al. 1996), and the synthesized beam is plotted in the bottom-left corner. The pb-corrected (1,1) line emissions with $<5\sigma$ rms are blank in each panel. Right column: beam-averaged line profiles extracted from the protobinary system and two positions perpendicular to the outflow direction in the flattened structure (labeled as $\alpha$ and $\beta$ on the $V_p$ map). The gray dashed vertical line indicates the $V_{sys} = 2.7$ km s$^{-1}$ of the cloud. The line intensity Jy beam$^{-1}$ is converted to the brightness temperature $T_B$, and the horizontal cyan line indicates the baseline ($T_B = 0$ K).</figcaption></figure>

the surface of ice grains is not very sensitive to shock at a velocity up to ~10 km s$^{-1}$, possibly indicating that NH$_3$ is efficiently returning in the gas phase already at the lowest values of $V_p - V_b$ (which may then be related to the shock velocity threshold for ice sputtering).

The OPR seems more enhanced along the curved jet path (B0a–B0e–B1a–B1c–B2b) compared to the rest by a factor of 2–2.5 (the fourth column of Figure 3). Interestingly, these locations where the OPR reaches a maximum (as high as 4–5 or 2.8–3.5 when adopting a line width of 4 km s$^{-1}$ or 9 km s$^{-1}$, respectively) are spatially coincident with the largest velocity offset in the line of sight ($V_p - V_b \sim 10$ km s$^{-1}$) in Figure 2, showing a strong Spearman’s rank correlation ($\rho > 0.5$) toward the entire mapping region (Figure D3). Although the locations of the maximum OPR and the total column density peak are not identical, the Spearman’s rank correlation between these parameters is strong ($\rho > 0.5$) toward B0, B1, and B2 locally. Such spatial correlation has also been reported toward other protostellar shocked regions (e.g., shocks associated with Orion-KL, Goddi et al. 2011).

Considering that the uncertainties from observation, data reduction, and LVG estimations ($<15\%$) are uniform for the entire map, we believe the enhancement of column density toward the cavity, the enhancement of OPR toward the jet path, as well as its spatial correlation with the velocity offset are intrinsic dynamic effects.

The OPR gradient revealed by our observations at high angular resolution improves the single-dish result from Umemoto et al. (1999).<fnref n="21" /> Enlightened by the explanation from the comprehensive gas-grain chemical network in Harju et al. (2017), surface and gas-phase processes may have different contributions along the cavity walls and within the cavity. According to Faure et al. (2013), OPR < 1 can be reproduced if H$_2$ is mainly in para form. Therefore, a high OPR (2–5) given from our study may be the result of o-H$_2$, which is not included in our LVG fitting but whose abundance may indeed increase as expected in shocked gas (see discussion in Appendix D).

<footnote n="21">²¹ Umemoto et al. (1999) assumed all lines are optically thin and reported the OPR toward B1 to be ~1.7$^{+0.2}_{-0.3}$, while toward B2 to be ~1.7 ± 0.2.</footnote>
### 4.4. The Flattened Envelope

Previous studies have revealed that the flattened structure surrounding L1157-mm is a mixture of outflow components with the inner envelope on the scale of thousands of astronomical units (e.g., Looney et al. 2007; Chiang et al. 2010; Tobin et al. 2010; Lee et al. 2012). Although the protobinary system mm is at the edge of our pb where the noise level is high from the pb-corrected images (Figures A1–B1), line emissions especially those from (1,1) and (2,2) have $S/N > 3$.

To study the velocity gradient across the flattened envelope, we present the $V_p$ maps of both lines within the region where the (1,1) line shows $>5\sigma$ emission after pb correction (Figure 4). Comparing the beam-averaged line profiles from the protobinary system mm, a northeast position $\alpha$, and a southwest position $\beta$ perpendicular to the outflow, we find that a $1 \text{ km s}^{-1}$ width narrow feature appears toward the mm on the spectra of the (2,2) line. Such feature appears on the (3,3) line as well (Figure 2), probably tracing the quiescent and cold gas in the flattened envelope. The blueshifted wing may trace the warmer material partially belonging to the outflow.

A velocity gradient of the $N_2H^+$ (1–0) line is reported by Chiang et al. (2010) and Tobin et al. (2010), which is normal to the outflow elongation and indicates the presence of rotation. In contrast, our $NH_3$ lines present a velocity gradient from the northeast (~2 km s$^{-1}$) to the southwest (~3 km s$^{-1}$). Nevertheless, this gradient is the same as what was shown on the JVLA $NH_3$ (1,1) map by Tobin et al. (2010), covering ±30" east and west offset to the protobinary system mm. Another consistency with Tobin et al. (2010) is the redshifted (~2.7 km s$^{-1}$) gas of $NH_3$, which seems to be curved down toward B0. Therefore, instead of an effect of sidelobe contamination (e.g., in dark blue on the (2,2) map of Figure 4), we believe the northeast to southwest velocity gradient of the $NH_3$ lines indicates the interaction between the outflow and the envelope.

### 5. Conclusions

Our newest JVLA observations provide high-sensitivity images of $NH_3$ (1,1)–(7,7) lines toward the archetypal protostellar shock, which is associated with the chemically rich blueshifted outflow in L1157. In the 0.1 pc scale field of view, we draw detailed kinetic temperature and $NH_3$ column density maps of this successively shocked region for the first time at a linear resolution of 1500 au. Our conclusions are as follows:

1. The emissions of all seven lines highlight the curved precessing jet path, from the eastern cavity wall of B0 and B1 to the cavity B2. In the line of sight, the map of the peak-to-bluest velocity offset reaches as high as 10 km s$^{-1}$ on this path.

2. Treating *o*- and *p*-NH$_3$ as distinct species, the kinetic temperature maps from LVG fittings show good agreement. The high-precision (<10% uncertainty) temperature map reveals an intrinsic gradient from the warm B0 eastern cavity wall (>120 K) to the cool cavity B1 and the earlier shock B2 (<80 K).

3. Both *o*- and *p*-NH$_3$ show a column density enhancement toward three U-shape cavities in B1 and B2, reaching as high as $>1 \times 10^{16} \text{ cm}^{-2}$. The OPR is enhanced by a factor of 2–2.5 along the curved jet path compared to the

rest, showing a strong spatial correlation with the peak-to-bluest velocity offset and with the total column density. All of these chemical gradients may be linked to the shocks.

4. A flattened envelope surrounding the protobinary system appears at the edge of our pb. We find a $1 \text{ km s}^{-1}$ width narrow feature on the line spectrum, probably tracing the quiescent and cold gas in the flattened envelope. The velocity map for the line intensity peak velocity shows a gradient from the northeast (~2 km s$^{-1}$) to the southwest (~3 km s$^{-1}$) of this structure, with the redshifted (~2.7 km s$^{-1}$) gas extended toward B0, indicating the interaction between the outflow and the envelope.

We gratefully acknowledge the useful and complimentary comments of the anonymous referee.

H.B.L. is supported by the Ministry of Science and Technology (MoST) of Taiwan (grant Nos. 108-2112-M-001-002-MY3).

F.D. is supported by the National Natural Science Foundation of China (NSFC) through grant Nos. 11873094 and 12041305.

C.Co. acknowledges the PRIN-INAF 2016 “The Cradle of Life—GENESIS-SKA (General Conditions in Early Planetary Systems for the rise of life with SKA).”

C.C. acknowledges receiving funds from the European Union’s Horizon 2020 research and innovation programme from the European Research Council (ERC) for the project “The Dawn of Organic Chemistry” (DOC), grant agreement No. 741002, and from the Marie Skłodowska-Curie for the project “Astro-Chemical Origins” (ACO), grant agreement No. 811312.

This research made use of NASA’s Astrophysics Data System.

*Software:* astropy (Astropy Collaboration et al. 2013), Numpy (van der Walt et al. 2011), CASA (v5.4.0; v5.6.1; McMullin et al. 2007), bettermoments (v1.0.0; Teague & Foreman-Mackey 2018), HfS (Estalella 2017).

## Appendix A
Line Profile

Along the jet paths, the protobinary system L1157-mm and nine representative clumpy substructures from B0 to B1 and B2 are denoted, with their coordinates listed in Table A1. The beam-averaged line profiles of $NH_3$ (1,1)–(7,7) are extracted from these nine positions, shown in Figure A1.

<table><caption>**Table A1**
Positions Corresponding to the Protostellar Object and the Clumpy Substructures</caption><thead><tr><th></th><th>R.A. (J2000)</th><th>Decl. (J2000)</th></tr></thead><tbody><tr><td>mm</td><td>20<sup>h</sup>39<sup>m</sup>06<sup>s</sup>170</td><td>68°02′15.07s<sup>′′</sup></td></tr><tr><td>B0a</td><td>20<sup>h</sup>39<sup>m</sup>10<sup>s</sup>026</td><td>68°01′41.30s<sup>′′</sup></td></tr><tr><td>B0e</td><td>20<sup>h</sup>39<sup>m</sup>10<sup>s</sup>365</td><td>68°01′23.80s<sup>′′</sup></td></tr><tr><td>B1a</td><td>20<sup>h</sup>39<sup>m</sup>10<sup>s</sup>301</td><td>68°01′13.51s<sup>′′</sup></td></tr><tr><td>B1b</td><td>20<sup>h</sup>39<sup>m</sup>08<sup>s</sup>734</td><td>68°01′10.54s<sup>′′</sup></td></tr><tr><td>B1c</td><td>20<sup>h</sup>39<sup>m</sup>09<sup>s</sup>677</td><td>68°01′06.45s<sup>′′</sup></td></tr><tr><td>B1f</td><td>20<sup>h</sup>39<sup>m</sup>09<sup>s</sup>444</td><td>68°01′16.70s<sup>′′</sup></td></tr><tr><td>B1i</td><td>20<sup>h</sup>39<sup>m</sup>10<sup>s</sup>000</td><td>68°01′01.11s<sup>′′</sup></td></tr><tr><td>B2a</td><td>20<sup>h</sup>39<sup>m</sup>12<sup>s</sup>755</td><td>68°00′38.44s<sup>′′</sup></td></tr><tr><td>B2b</td><td>20<sup>h</sup>39<sup>m</sup>14<sup>s</sup>217</td><td>68°00′24.73s<sup>′′</sup></td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Figure A1.** Beam-averaged line profiles (the line intensity Jy beam⁻¹ is converted to the brightness temperature T_b), with the center toward nine clumpy substructures (Table A1) in the plane of the sky. The HfS fit is shown in red overlying each spectrum. All lines are extracted from images smoothed to the same pixel size (0.′′5), and their velocity resolutions are listed in Table 1. In each panel, two red vertical lines at -25 and + 25 km s⁻¹ indicate the velocity range within which we integrate the intensity; the gray dashed vertical line indicates the V_sys = 2.7 km s⁻¹ of the cloud. The horizontal cyan line indicates the baseline (T_B = 0 K).</figcaption></figure>

Applying the Monte Carlo fitting tool HfS developed by Estalella (2017), we fit the hyperfine multiplets of each line with a Gaussian shape model using the “hfs_fit” procedure and assuming one velocity component. The best-fit results of the line parameters are listed in Table A2 and shown as a red curve overlying the observation spectrum in Figure A1.

Testing the other mainstream fitting packages, such as CLASS/GILDAS (Pety 2005) and PySpecKit (Ginsburg & Mirocha 2011), we found similar fitting results.

Given that the following assumptions were made in these mainstream fitting packages, some caveats need to be noted:

(1) The mainstream fitting packages take the Gaussian profile as a model, but our target lines show significant blueshifted emissions in this successive shocked region.

Therefore, the centroid velocity of the main component is not strictly identical to the velocity of the peak intensity, and the FWHM line width from fitting is slightly narrower than the observed line wing, especially for the (3,3)–(6,6) lines (Figure A1).

(2) The mainstream fitting packages provide the optical depth of the main component τ_m. This optical depth is close to that of the entire line τ, but only when the satellites have <30% contribution to the line-integrated intensity, i.e., for the (2,2)–(7,7) lines. As to the (1,1) line where satellites may also be optically thick, τ_m may be significantly smaller than τ.

(3) As a compromise of the fit to hyperfine multiplets, one velocity component is assumed in the fitting. This may bring in an overestimation of the FWHM line width.
<table><caption>Table A2
The Best-fit Parameters for Spectra in Figure A1</caption><thead><tr><th>Position</th><th>Best-fit parameter<sup>a</sup></th><th>(1,1)</th><th>(2,2)</th><th>(3,3)</th><th>(4,4)</th><th>(5,5)</th><th>(6,6)</th><th>(7,7)</th></tr></thead><tbody><tr><td rowspan="5">B0a</td><td>δV(km s<sup>-1</sup>)<sup>b</sup></td><td>3.59 ± 0.56</td><td>4.47 ± 0.50</td><td>4.34 ± 0.21</td><td>4.64 ± 1.75</td><td>4.80 ± 0.65</td><td>9.60 ± 2.79</td><td><sup>h</sup></td></tr><tr><td>V<sub>lsr</sub>(km s<sup>-1</sup>)<sup>c</sup></td><td>0.56 ± 0.20</td><td>1.02 ± 0.19</td><td>1.19 ± 0.08</td><td>-0.61 ± 0.89</td><td>0.38 ± 0.44</td><td>-0.93 ± 1.18</td><td><sup>h</sup></td></tr><tr><td>A(1 - e<sup>-τm</sup>)(K)<sup>d,e</sup></td><td>1.04 ± 0.11</td><td>1.07 ± 0.10</td><td>3.60 ± 0.13</td><td>0.46 ± 0.11</td><td>0.37 ± 0.05</td><td>0.31 ± 0.07</td><td><sup>h</sup></td></tr><tr><td>τ<sup>*</sup><sub>m</sub> = 1 - e<sup>-τm</sup></td><td>0.26 ± 0.43</td><td>(4.0 ± 4.0)E-4</td><td>(2.0 ± 2.0)E-4</td><td>1.00 ± 0.07<sup>g</sup></td><td>0.72 ± 0.46<sup>g</sup></td><td>(4.0 ± 4.0)E-4</td><td><sup>h</sup></td></tr><tr><td>I(K km s<sup>-1</sup>)<sup>f</sup></td><td>9.54<sup>+0.36</sup><sub>-3.64</sub></td><td>7.00<sup>+0.33</sup><sub>-1.40</sub></td><td>23.44<sup>+0.49</sup><sub>-2.02</sub></td><td>5.28<sup>+0.27</sup><sub>-1.09</sub></td><td>2.94<sup>+0.07</sup><sub>-0.26</sub></td><td>3.66<sup>+0.12</sup><sub>-0.78</sub></td><td>1.17<sup>+0.17</sup><sub>-0.08</sub></td></tr><tr><td rowspan="5">B0e</td><td>δV(km s<sup>-1</sup>)<sup>b</sup></td><td>1.57 ± 0.53</td><td>4.94 ± 1.67</td><td>4.11 ± 0.25</td><td>4.31 ± 0.45</td><td>5.29 ± 4.98</td><td>3.82 ± 0.63</td><td><sup>h</sup></td></tr><tr><td>V<sub>lsr</sub>(km s<sup>-1</sup>)<sup>c</sup></td><td>0.77 ± 0.01</td><td>0.88 ± 0.24</td><td>1.49 ± 0.09</td><td>0.66 ± 0.33</td><td>-0.53 ± 0.98</td><td>0.20 ± 0.34</td><td><sup>h</sup></td></tr><tr><td>A(1 - e<sup>-τm</sup>)(K)<sup>d,e</sup></td><td>0.63 ± 0.14</td><td>0.75 ± 0.08</td><td>3.15 ± 0.14</td><td>0.49 ± 0.05</td><td>0.22 ± 0.07</td><td>0.29 ± 0.03</td><td><sup>h</sup></td></tr><tr><td>τ<sup>*</sup><sub>m</sub> = 1 - e<sup>-τm</sup></td><td>0.80 ± 0.16</td><td>0.45 ± 1.49<sup>g</sup></td><td>(5.0 ± 5.0)E-4</td><td>0.90 ± 0.14<sup>g</sup></td><td>0.99 ± 0.44<sup>g</sup></td><td>1.00 ± 0.02<sup>g</sup></td><td><sup>h</sup></td></tr><tr><td>I(K km s<sup>-1</sup>)<sup>f</sup></td><td>4.94<sup>+0.43</sup><sub>-1.52</sub></td><td>6.46<sup>+0.42</sup><sub>-1.65</sub></td><td>19.88<sup>+0.44</sup><sub>-1.66</sub></td><td>3.86<sup>+0.16</sup><sub>-0.42</sub></td><td>2.23<sup>+0.13</sup><sub>-0.19</sub></td><td>2.37<sup>+0.07</sup><sub>-0.09</sub></td><td>0.57<sup>+0.16</sup><sub>-0.07</sub></td></tr><tr><td rowspan="5">B1a</td><td>δV(km s<sup>-1</sup>)<sup>b</sup></td><td>5.46 ± 0.39</td><td>5.57 ± 0.49</td><td>5.34 ± 0.17</td><td>5.36 ± 0.33</td><td>6.87 ± 0.65</td><td>6.40 ± 0.39</td><td>1.78 ± 4.53</td></tr><tr><td>V<sub>lsr</sub>(km s<sup>-1</sup>)<sup>c</sup></td><td>-0.22 ± 0.18</td><td>0.13 ± 0.14</td><td>0.58 ± 0.07</td><td>-0.21 ± 0.22</td><td>-0.17 ± 0.27</td><td>-0.21 ± 0.23</td><td>-2.37 ± 0.08</td></tr><tr><td>A(1 - e<sup>-τm</sup>)(K)<sup>d,e</sup></td><td>1.00 ± 0.07</td><td>1.28 ± 0.06</td><td>4.67 ± 0.12</td><td>0.71 ± 0.04</td><td>0.50 ± 0.04</td><td>0.57 ± 0.03</td><td>0.18 ± 0.04</td></tr><tr><td>τ<sup>*</sup><sub>m</sub> = 1 - e<sup>-τm</sup></td><td>0.27 ± 0.28</td><td>0.50 ± 0.33</td><td>(3.1 ± 3.1)E-3</td><td>0.83 ± 0.10</td><td>0.01 ± 0.01</td><td>0.53 ± 0.22</td><td>1.00 ± 0.01</td></tr><tr><td>I(K km s<sup>-1</sup>)<sup>f</sup></td><td>12.46<sup>+0.14</sup><sub>-4.57</sub></td><td>10.68<sup>+0.42</sup><sub>-2.16</sub></td><td>32.76<sup>+0.14</sup><sub>-2.35</sub></td><td>5.35<sup>+0.11</sup><sub>-0.17</sub></td><td>4.01<sup>+0.06</sup><sub>-0.35</sub></td><td>4.58<sup>+0.08</sup><sub>-0.12</sub></td><td>1.53<sup>+0.02</sup><sub>-0.18</sub></td></tr><tr><td rowspan="5">B1b</td><td>δV(km s<sup>-1</sup>)<sup>b</sup></td><td>2.72 ± 4.35</td><td>5.30 ± 0.44</td><td>4.83 ± 0.28</td><td>5.72 ± 0.64</td><td>3.62 ± 0.73</td><td>4.46 ± 0.39</td><td>3.56 ± 1.26</td></tr><tr><td>V<sub>lsr</sub>(km s<sup>-1</sup>)<sup>c</sup></td><td>1.41 ± 0.01</td><td>-0.41 ± 0.13</td><td>-0.37 ± 0.07</td><td>-1.07 ± 0.28</td><td>-0.97 ± 0.50</td><td>-1.24 ± 0.30</td><td>-1.68 ± 0.62</td></tr><tr><td>A(1 - e<sup>-τm</sup>)(K)<sup>d,e</sup></td><td>1.11 ± 0.14</td><td>1.24 ± 0.06</td><td>3.50 ± 0.08</td><td>0.60 ± 0.06</td><td>0.28 ± 0.04</td><td>0.35 ± 0.03</td><td>0.14 ± 0.04</td></tr><tr><td>τ<sup>*</sup><sub>m</sub> = 1 - e<sup>-τm</sup></td><td>0.84 ± 0.13</td><td>0.55 ± 0.28</td><td>0.76 ± 0.13</td><td>(1.8 ± 1.8)E-3</td><td>1.00 ± 0.03<sup>g</sup></td><td>0.95 ± 0.08<sup>g</sup></td><td>0.95 ± 0.72<sup>g</sup></td></tr><tr><td>I(K km s<sup>-1</sup>)<sup>f</sup></td><td>15.15<sup>+0.04</sup><sub>-5.97</sub></td><td>9.00<sup>+0.41</sup><sub>-1.24</sub></td><td>23.90<sup>+0.49</sup><sub>-0.93</sub></td><td>3.42<sup>+0.20</sup><sub>-0.25</sub></td><td>2.59<sup>+0.22</sup><sub>-0.57</sub></td><td>2.59<sup>+0.05</sup><sub>-0.21</sub></td><td>0.84<sup>+0.15</sup><sub>-0.11</sub></td></tr><tr><td rowspan="5">B1c</td><td>δV(km s<sup>-1</sup>)<sup>b</sup></td><td>5.30 ± 0.25</td><td>5.37 ± 0.22</td><td>4.97 ± 0.13</td><td>4.64 ± 0.28</td><td>4.15 ± 0.31</td><td>5.75 ± 0.71</td><td><sup>h</sup></td></tr><tr><td>V<sub>lsr</sub>(km s<sup>-1</sup>)<sup>c</sup></td><td>0.00 ± 0.12</td><td>0.00 ± 0.09</td><td>0.11 ± 0.05</td><td>-0.60 ± 0.19</td><td>-0.53 ± 0.24</td><td>-0.79 ± 0.29</td><td><sup>h</sup></td></tr><tr><td>A(1 - e<sup>-τm</sup>)(K)<sup>d,e</sup></td><td>1.51 ± 0.07</td><td>1.65 ± 0.06</td><td>4.67 ± 0.10</td><td>0.82 ± 0.05</td><td>0.44 ± 0.03</td><td>0.50 ± 0.05</td><td><sup>h</sup></td></tr><tr><td>τ<sup>*</sup><sub>m</sub> = 1 - e<sup>-τm</sup></td><td>0.32 ± 0.17</td><td>0.03 ± 0.03</td><td>0.02 ± 0.02</td><td>0.76 ± 0.13</td><td>0.96 ± 0.05<sup>g</sup></td><td>(1.0 ± 1.0)E-4</td><td><sup>h</sup></td></tr><tr><td>I(K km s<sup>-1</sup>)<sup>f</sup></td><td>19.08<sup>+0.09</sup><sub>-7.36</sub></td><td>11.20<sup>+0.70</sup><sub>-1.55</sub></td><td>29.59<sup>+0.59</sup><sub>-1.61</sub></td><td>5.38<sup>+0.03</sup><sub>-0.37</sub></td><td>3.22<sup>+0.14</sup><sub>-0.44</sub></td><td>3.20<sup>+0.01</sup><sub>-0.04</sub></td><td>1.10<sup>+0.09</sup><sub>-0.18</sub></td></tr><tr><td rowspan="5">B1f</td><td>δV(km s<sup>-1</sup>)<sup>b</sup></td><td>3.88 ± 0.36</td><td>5.32 ± 0.53</td><td>5.00 ± 0.52</td><td>4.02 ± 0.33</td><td>6.37 ± 0.84</td><td>4.15 ± 0.32</td><td><sup>h</sup></td></tr><tr><td>V<sub>lsr</sub>(km s<sup>-1</sup>)<sup>c</sup></td><td>-0.66 ± 0.01</td><td>-0.27 ± 0.16</td><td>0.24 ± 0.09</td><td>-0.37 ± 0.27</td><td>-0.11 ± 0.35</td><td>-0.38 ± 0.28</td><td><sup>h</sup></td></tr><tr><td>A(1 - e<sup>-τm</sup>)(K)<sup>d,e</sup></td><td>0.76 ± 0.08</td><td>0.92 ± 0.05</td><td>3.26 ± 0.09</td><td>0.50 ± 0.04</td><td>0.37 ± 0.04</td><td>0.37 ± 0.03</td><td><sup>h</sup></td></tr><tr><td>τ<sup>*</sup><sub>m</sub> = 1 - e<sup>-τm</sup></td><td>0.58 ± 0.28</td><td>0.65 ± 0.27<sup>g</sup></td><td>0.44 ± 0.43</td><td>0.97 ± 0.04<sup>g</sup></td><td>0.09 ± 0.09</td><td>0.99 ± 0.02<sup>g</sup></td><td><sup>h</sup></td></tr><tr><td>I(K km s<sup>-1</sup>)<sup>f</sup></td><td>7.44<sup>+0.03</sup><sub>-2.69</sub></td><td>7.44<sup>+0.25</sup><sub>-1.36</sub></td><td>22.54<sup>+0.35</sup><sub>-1.74</sub></td><td>3.22<sup>+0.35</sup><sub>-0.14</sub></td><td>2.59<sup>+0.08</sup><sub>-0.12</sub></td><td>2.95<sup>+0.03</sup><sub>-0.18</sub></td><td>0.51<sup>+0.03</sup><sub>-0.19</sub></td></tr><tr><td rowspan="5">B1i</td><td>δV(km s<sup>-1</sup>)<sup>b</sup></td><td>1.52 ± 5.84</td><td>5.18 ± 0.34</td><td>4.18 ± 0.38</td><td>3.64 ± 0.34</td><td>4.26 ± 0.80</td><td>3.96 ± 0.82</td><td><sup>h</sup></td></tr><tr><td>V<sub>lsr</sub>(km s<sup>-1</sup>)<sup>c</sup></td><td>1.86 ± 0.01</td><td>-0.26 ± 0.15</td><td>-0.21 ± 0.07</td><td>-0.84 ± 0.26</td><td>0.11 ± 0.34</td><td>-0.45 ± 0.34</td><td><sup>h</sup></td></tr><tr><td>A(1 - e<sup>-τm</sup>)(K)<sup>d,e</sup></td><td>0.56 ± 0.07</td><td>0.93 ± 0.05</td><td>2.44 ± 0.07</td><td>0.53 ± 0.05</td><td>0.33 ± 0.05</td><td>0.32 ± 0.05</td><td><sup>h</sup></td></tr><tr><td>τ<sup>*</sup><sub>m</sub> = 1 - e<sup>-τm</sup></td><td>0.99 ± 0.02<sup>g</sup></td><td>(4.0 ± 4.0)E-4</td><td>0.39 ± 0.39</td><td>0.95 ± 0.07<sup>g</sup></td><td>(3.0 ± 3.0)E-4</td><td>(5.4 ± 5.4)E-3</td><td><sup>h</sup></td></tr><tr><td>I(K km s<sup>-1</sup>)<sup>f</sup></td><td>9.61<sup>+0.04</sup><sub>-3.68</sub></td><td>5.73<sup>+0.24</sup><sub>-0.68</sub></td><td>13.96<sup>+0.37</sup><sub>-0.93</sub></td><td>3.73<sup>+0.35</sup><sub>-0.71</sub></td><td>1.56<sup>+0.11</sup><sub>-0.14</sub></td><td>1.38<sup>+0.05</sup><sub>-0.01</sub></td><td>0.61<sup>+0.07</sup><sub>-0.04</sub></td></tr><tr><td rowspan="5">B2a</td><td>δV(km s<sup>-1</sup>)<sup>b</sup></td><td>2.70 ± 0.10</td><td>3.21 ± 0.11</td><td>3.15 ± 0.05</td><td>3.68 ± 0.26</td><td>3.82 ± 0.37</td><td>4.38 ± 0.33</td><td>3.43 ± 3.09</td></tr><tr><td>V<sub>lsr</sub>(km s<sup>-1</sup>)<sup>c</sup></td><td>1.99 ± 0.03</td><td>1.81 ± 0.04</td><td>1.75 ± 0.02</td><td>1.63 ± 0.10</td><td>1.47 ± 0.15</td><td>1.52 ± 0.12</td><td>1.09 ± 0.15</td></tr><tr><td>A(1 - e<sup>-τm</sup>)(K)<sup>d,e</sup></td><td>4.62 ± 0.13</td><td>3.76 ± 0.11</td><td>10.29 ± 0.15</td><td>1.67 ± 0.10</td><td>0.90 ± 0.07</td><td>1.05 ± 0.06</td><td>0.34 ± 0.12</td></tr><tr><td>τ<sup>*</sup><sub>m</sub> = 1 - e<sup>-τm</sup></td><td>0.31 ± 0.11</td><td>(1.0 ± 1.0)E-3</td><td>(2.8 ± 2.8)E-3</td><td>(1.2 ± 1.2)E-3</td><td>(7.0 ± 7.0)E-4</td><td>(2.0 ± 2.0)E-4</td><td>0.97 ± 0.86<sup>g</sup></td></tr><tr><td>I(K km s<sup>-1</sup>)<sup>f</sup></td><td>30.82<sup>+0.18</sup><sub>-11.73</sub></td><td>15.80<sup>+1.12</sup><sub>-2.93</sub></td><td>41.63<sup>+1.37</sup><sub>-3.51</sub></td><td>7.14<sup>+0.06</sup><sub>-0.94</sub></td><td>3.64<sup>+0.16</sup><sub>-0.03</sub></td><td>5.54<sup>+0.02</sup><sub>-0.58</sub></td><td>1.62<sup>+0.04</sup><sub>-0.02</sub></td></tr></tbody></table>
<table><caption>Table A2
(Continued)</caption><thead><tr><th>Position</th><th>Best-fit parameter<sup>a</sup></th><th>(1,1)</th><th>(2,2)</th><th>(3,3)</th><th>(4,4)</th><th>(5,5)</th><th>(6,6)</th><th>(7,7)</th></tr></thead><tbody><tr><td rowspan="5">B2b</td><td>δV(km s<sup>-1</sup>)<sup>b</sup></td><td>3.86 ± 0.25</td><td>3.90 ± 0.22</td><td>3.81 ± 0.10</td><td>4.24 ± 0.54</td><td>3.95 ± 0.68</td><td>4.31 ± 0.69</td><td>2.74 ± 1.20</td></tr><tr><td>V<sub>lsr</sub>(km s<sup>-1</sup>)<sup>c</sup></td><td>1.81 ± 0.16</td><td>1.65 ± 0.08</td><td>1.60 ± 0.04</td><td>1.33 ± 0.21</td><td>1.48 ± 0.26</td><td>1.41 ± 0.28</td><td>1.44 ± 0.54</td></tr><tr><td>A(1 − e<sup>-τ<sub>m</sub></sup>)(K)<sup>d,e</sup></td><td>2.57 ± 0.16</td><td>2.68 ± 0.12</td><td>7.88 ± 0.18</td><td>1.23 ± 0.13</td><td>0.78 ± 0.11</td><td>0.76 ± 0.10</td><td>0.25 ± 0.08</td></tr><tr><td>τ<sub>m</sub><sup>*</sup> = 1 − e<sup>-τ<sub>m</sub></sup></td><td>0.11 ± 0.12</td><td>(1.3 ± 1.3)E − 3</td><td>(3.0 ± 3.0)E − 4</td><td>(1.0 ± 1.0)E − 4</td><td>(3.3 ± 3.3)E − 3</td><td>0.02 ± 0.02</td><td>0.97 ± 0.93<sup>g</sup></td></tr><tr><td>I(K km s<sup>-1</sup>)<sup>f</sup></td><td>22.30<sup>+0.10</sup><sub>-8.67</sub></td><td>14.57<sup>+0.99</sup><sub>-3.07</sub></td><td>39.51<sup>+0.65</sup><sub>-2.39</sub></td><td>5.88<sup>+0.59</sup><sub>-0.11</sub></td><td>4.41<sup>+0.17</sup><sub>-0.58</sub></td><td>3.60<sup>+0.12</sup><sub>-0.07</sub></td><td>1.88<sup>+0.1</sup><sub>-0.72</sub></td></tr></tbody></table>

**Notes.**

a Given by using the “hfs_fit” procedure in the HfS package, fit by assuming a Gaussian line profile; all beam-averaged lines are extracted from images with the same pixel size.

b Intrinsic FWHM line width by taking into account the opacity.

c Centroid velocity of the main component.

d Peak intensity of the main component (for lines broader than the channel width).

e Optical depth of the main component τ<sub>m</sub>, provided by the HfS as τ<sub>m</sub><sup>*</sup>, τ<sub>m</sub> = −ln(1 − τ<sub>m</sub><sup>*</sup>).

f Integrated intensity over the velocity range from −25 to +25km s<sup>-1</sup>, the uncertainty is given by the difference between this range and the velocity range of −10~+5 km s<sup>-1</sup> (the main component of the hyperfine multiplets) and −27~+33 km s<sup>-1</sup> (line intensity down to zero).

g Fitting results need to be taken cautiously for τ<sub>m</sub> > 1 and with an uncertainty larger than the value.

h Line emission is less than 3σ.
## Appendix B
### Images with Primary Beam Correction

Given that the beam responses are different for each antenna in a synthesis array, the flux density at the beam edge is biased low. Dividing the image by the pb, we apply the pb correction to each

line datacube and provide the integrated intensity maps over the velocity range of −25 to +25 km s⁻¹ in Figure B1. Compared with Figure 1, line emission toward the edge of the pb (e.g., mm, B2b) are recovered. However, noise levels are increasing toward these positions (see example spectrum in Figure 2).

<figure><img src="image_9.png" /></figure>

**Figure B1.** The outline of the southern (blueshifted) outflow lobe from the Class 0 compact protobinary system L1157-mm (images with different weightings are given in different rows). Color maps show the intensity maps of NH₃ (1,1)–(7,7) integrated over the velocity range of −25 to +25 km s⁻¹ (with pb correction). The black contours, starting from 4σ (σ = 0.03 mJy beam⁻¹) and increasing with a step of 4σ, show the continuum emission with the same weighting as the line emissions (with pb correction). The rest of the contours and all the labels are the same as in Figure 1. The solid gray circle indicates the largest recoverable scale and the dashed gray circle indicates the primary beam.
## Appendix C
Line Optical Depth Correction

Even though the mainstream fitting packages bring in a centroid velocity map a bit different from the peak velocity map (Section A), they can constrain the lower limit of the optical depth for a particular line $\tau$ for each pixel. Therein, to correct the molecular column density, which is underestimated, we can multiply the line-integrated intensity from observations by a factor of $\tau/(1 - e^{-\tau})$.

Applying the “hfs cube mp” procedure in the HfS fitting tool, we focus on deriving the optical depth $\tau_m$ map and its uncertainty map for each line. Because the parameter directly

given by the fitting program is $\tau_m^* = 1 - e^{-\tau_m}$, the error propagation algorithm enhances the uncertainty of $\tau_m$ due to the fitting caveats listed in Section A. Therefore, we set a threshold for the pixels with a high fidelity, i.e., when the uncertainty with respect to the value is less than 1 dex.

As shown in Figure C1, the optical depths of all seven lines are negligible ($\tau_m < 1$) toward B2, while those toward the B1 cavity wall (B1a, B1f, B1i, B1b) are larger than 3, except for the (3,3) line. The $\tau_m$ of (4,4)–(6,6) may be underestimated toward some pixels, where the uncertainties are higher than the value by a factor of more than 5.

<figure><img src="image_10.png" /><figcaption>**Figure C1.** Optical depth $\tau_m$ maps of the (1,1)–(6,6) lines. The gray dashed contours and the labeled clumpy substructures are the same as shown in Figure 1. The integrated intensity of the line within the red solid contour shows >4$\sigma$ emissions. The pixels where lines show <4$\sigma$ emissions or uncertainty 1 dex higher than the value of $\tau_m$ are blank.</figcaption></figure>
## Appendix D
LVG Fitting Example and Caveats

The optical depth of each line is significant, uncertain, and not uniform throughout the entire region. Even correcting each line-integrated intensity with the optical depth $\tau_m$, the rotation diagram (RD) method is not an optimized approach to deriving the temperature and density structure of our source. Instead, the LVG approximation is suitable under such a circumstance, though it is computationally expensive.

Figure D1 gives a PDF example of the $T_{\text{kin}}$, $N_T$, and $n$ from LVG MultiNest fitting toward B1a by treating the *o*- and *p*-NH$_3$ as distinct species. Table D1 lists the mean, standard deviation, and the best fit of $T_{\text{kin}}$ and $N_T$ for *o*/p-NH$_3$ toward nine clumps. Figure D2 gives a combination of $T_{\text{kin}}$ and $N_T$ maps for *o*/p-NH$_3$ by testing the cases assuming the beam-filling factor to be unity for all lines and the line width to be 4 km s$^{-1}$ or 9 km s$^{-1}$. Figure D3 gives possible spatial correlations between the variables derived from LVG fittings.

One major caveat in the above LVG fitting is that LAMDA only provides the collision rates between the *p*/-*o*-NH$_3$ and *p*-H$_2$ (Danby et al. 1988). Although we note that some *o*-H$_2$ cross sections were computed at a few selected energies (e.g., Offer & Flower 1989; Rist et al. 1993; Bouhafs et al. 2017), no collision rates between *p*/-*o*-NH$_3$ and *o*-H$_2$ are available in LAMDA. Unlike the case in cold dense cores, the abundance of *o*-H$_2$ can become significant in warm and shocked regions (e.g., Flower et al. 2006; Neufeld et al. 2019).

A second caveat is the beam-filling factor, which we take as unity for all lines in the fitting, given that the line emissions show more extended spatial distributions toward all the clumpy substructures than the synthesized beam. However, considering that shock discontinuities cannot be spatially resolved by observations, substructures with high column density and high temperature might be beam smeared. This is likely the case after checking the peak brightness temperatures of all

transitions toward the representative positions (Figure A1) and their optical depth maps (Figure C1). For testing, we set it as unity for (1,1)–(3,3) and 0.5, or 0.1, or 0.01 for (4,4)–(7,7), and then compare the testing results with the results by assuming the beam-filling factor of unity for all lines. When the filling factor for higher transition levels is adopted in the range of 0.1–1, a smaller filling factor decreases the kinetic temperature significantly, by 30% for *p*-NH$_3$ but by <10% for *o*-NH$_3$. In contrast, the column density of *o*/p-NH$_3$ only shows imperceptible changes (<5%) when a different beam-filling factor is adopted.

A third caveat is the FWMH line width we adopted, i.e., the median from all spectrum fitting as 4 km s$^{-1}$. This seems to be consistent with the FWHM line width of (1,1)–(3,3) at an angular resolution of 40", which was reported to be 8–9 km s$^{-1}$ by Bachiller et al. (1993). However, the peak-to-bluest velocity offset varies pixel by pixel (Figure 2), being broader along the jet path than within the cavity by a factor of 2. As a compromise to the NH$_3$ hyperfine multiplets fitting, we assume one velocity component. In the case of multivelocity components being smeared within the VLA's synthesized beam, an FWHM line width of 9 km s$^{-1}$ is adopted for the test. We find that such a line width leads to a 15% and 20% increase in the best-fit $T_{\text{kin}}$ for *p*-NH$_3$ and *o*-NH$_3$, as well as a 15% and 45% decrease in their $N_T$, respectively (Table D1). Nevertheless, these increase and decrease factors are globally consistent throughout the entire map pixels, so the contrast (local gradient) of each map does not change.

A minor fourth caveat is the velocity range over which we integrate the line intensity. In Table A2, we list the uncertainty of the line-integrated intensity $I$, which is given by comparing the integration (i) where the line emission has $S/N > 3$ (–25 to +25 km s$^{-1}$); (ii) where the line emission is above zero (–27 to +33 km s$^{-1}$), and (iii) includes only the main component of the hyperfine multiplets (–10 to +5 km s$^{-1}$). Due to strong

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure D1. The PDF of *o*/p-NH$_3$ parameters toward B1a, derived by using RADEX and the MultiNest algorithm. A line width of 4 km s$^{-1}$ is adopted for all lines, and the beam-filling factor is assumed to be unity for the (1,1)–(3,3) lines while to be 0.1 for the (4,4)–(7,7) lines. The fitting parameters, including the H$_2$ volume density (*n*), *o*/p-NH$_3$ column density *N*(*o*/p-NH$_3$) and the gas kinetic temperature (*T*$_{kin}$), are listed.</figcaption></figure>
<table><caption>Table D1
Mean and Standard Deviation of the PDF and Best-fit Results Using the LVG Method</caption><thead><tr><th rowspan="2">Location</th><th rowspan="2">Species<sup>a</sup></th><th rowspan="2">ff<sup>b</sup></th><th colspan="2"><i>N</i><sub>T,n</sub><sup>c</sup> (10<sup>15</sup> cm<sup>-2</sup>)</th><th colspan="2"><i>T</i><sub>kin,n</sub><sup>c</sup> (K)</th><th colspan="2"><i>N</i><sub>T,b</sub><sup>d</sup> (10<sup>15</sup> cm<sup>-2</sup>)</th><th colspan="2"><i>T</i><sub>kin,b</sub><sup>d</sup> (K)</th><th rowspan="2"><i>n</i><sup>e</sup> (10<sup>5</sup> cm<sup>-3</sup>)</th></tr><tr><th>mean ± std</th><th>best fit</th><th>mean ± std</th><th>best fit</th><th>mean ± std</th><th>best fit</th><th>mean ± std</th><th>best fit</th></tr></thead><tbody><tr><td rowspan="4">B0a</td><td><i>p</i>-NH<sub>3</sub></td><td>I</td><td>1.0 ± 0.0</td><td>1.1</td><td>436 ± 158</td><td>327</td><td>0.8 ± 0.1</td><td>0.9</td><td>525 ± 174</td><td>802</td><td>0.1–8.1</td></tr><tr><td></td><td>II</td><td>1.0 ± 0.0</td><td>1.0</td><td>305 ± 124</td><td>204</td><td>0.9 ± 0.1</td><td>0.9</td><td>369 ± 145</td><td>271</td><td>0.03–33.7</td></tr><tr><td><i>o</i>-NH<sub>3</sub></td><td>I</td><td>4.2 ± 0.8</td><td>4.4</td><td>102 ± 8</td><td>100</td><td>2.6 ± 0.4</td><td>2.6</td><td>122 ± 11</td><td>121</td><td>15.2–73.4</td></tr><tr><td></td><td>II</td><td>4.2 ± 0.8</td><td>4.4</td><td>102 ± 8</td><td>100</td><td>2.6 ± 0.4</td><td>2.6</td><td>122 ± 11</td><td>121</td><td>0.5–5.0</td></tr><tr><td rowspan="4">B0e</td><td><i>p</i>-NH<sub>3</sub></td><td>I</td><td>0.8 ± 0.0</td><td>0.9</td><td>278 ± 91</td><td>233</td><td>0.8 ± 0.1</td><td>0.8</td><td>386 ± 132</td><td>292</td><td>4.8–61.4</td></tr><tr><td></td><td>II</td><td>1.0 ± 0.0</td><td>0.9</td><td>294 ± 78</td><td>196</td><td>0.8 ± 0.1</td><td>0.8</td><td>439 ± 181</td><td>198</td><td>0.03–0.1</td></tr><tr><td><i>o</i>-NH<sub>3</sub></td><td>I</td><td>4.0 ± 0.7</td><td>4.2</td><td>103 ± 9</td><td>102</td><td>2.5 ± 0.4</td><td>2.6</td><td>124 ± 11</td><td>122</td><td>0.3–4.9</td></tr><tr><td></td><td>I</td><td>4.1 ± 0.7</td><td>4.2</td><td>103 ± 9</td><td>102</td><td>2.5 ± 0.4</td><td>2.6</td><td>124 ± 11</td><td>123</td><td>0.8–92</td></tr><tr><td rowspan="4">B1f</td><td><i>p</i>-NH<sub>3</sub></td><td>I</td><td>0.8 ± 0.1</td><td>0.8</td><td>175 ± 41</td><td>154</td><td>0.6 ± 0.1</td><td>0.7</td><td>206 ± 57</td><td>178</td><td>5.8–46.0</td></tr><tr><td></td><td>II</td><td>0.8 ± 0.1</td><td>0.8</td><td>122 ± 18</td><td>117</td><td>0.6 ± 0.1</td><td>0.6</td><td>139 ± 24</td><td>130</td><td>6.4–29</td></tr><tr><td><i>o</i>-NH<sub>3</sub></td><td>I</td><td>3.7 ± 0.7</td><td>3.9</td><td>106 ± 9</td><td>104</td><td>2.3 ± 0.3</td><td>2.4</td><td>127 ± 11</td><td>126</td><td>0.1–0.2</td></tr><tr><td></td><td>II</td><td>3.7 ± 0.7</td><td>3.9</td><td>106 ± 9</td><td>105</td><td>2.3 ± 0.3</td><td>2.4</td><td>128 ± 12</td><td>126</td><td>11–84</td></tr><tr><td rowspan="4">B1a</td><td><i>p</i>-NH<sub>3</sub></td><td>I</td><td>1.6 ± 0.1</td><td>1.6</td><td>275 ± 98</td><td>211</td><td>1.3 ± 0.0</td><td>1.3</td><td>400 ± 144</td><td>319</td><td>1.1–32.2</td></tr><tr><td></td><td>II</td><td>1.6 ± 0.2</td><td>1.6</td><td>182 ± 61</td><td>146</td><td>1.3 ± 0.1</td><td>1.3</td><td>243 ± 90</td><td>184</td><td>5.2–6.4</td></tr><tr><td><i>o</i>-NH<sub>3</sub></td><td>I</td><td>7.7 ± 1.5</td><td>8.1</td><td>98 ± 8</td><td>97</td><td>4.3 ± 0.7</td><td>4.4</td><td>119 ± 11</td><td>117</td><td>10.3–14.1</td></tr><tr><td></td><td>II</td><td>7.8 ± 1.5</td><td>8.1</td><td>98 ± 8</td><td>97</td><td>4.3 ± 0.7</td><td>4.4</td><td>119 ± 11</td><td>117</td><td>1.7–2.8</td></tr><tr><td rowspan="4">B1b</td><td><i>p</i>-NH<sub>3</sub></td><td>I</td><td>1.5 ± 0.1</td><td>1.6</td><td>103 ± 11</td><td>100</td><td>1.3 ± 0.0</td><td>1.3</td><td>118 ± 13</td><td>115</td><td>1.9–2.9</td></tr><tr><td></td><td>II</td><td>1.6 ± 0.2</td><td>1.7</td><td>83 ± 8</td><td>82</td><td>1.3 ± 0.1</td><td>1.3</td><td>96 ± 11</td><td>96</td><td>3–98</td></tr><tr><td><i>o</i>-NH<sub>3</sub></td><td>I</td><td>4.7 ± 0.9</td><td>4.9</td><td>94 ± 7</td><td>93</td><td>2.8 ± 0.4</td><td>2.9</td><td>112 ± 9</td><td>110</td><td>1.4–4.7</td></tr><tr><td></td><td>II</td><td>4.7 ± 0.9</td><td>4.9</td><td>94 ± 7</td><td>93</td><td>2.8 ± 0.4</td><td>2.9</td><td>112 ± 9</td><td>110</td><td>0.01–13</td></tr><tr><td rowspan="4">B1c</td><td><i>p</i>-NH<sub>3</sub></td><td>I</td><td>1.7 ± 0.2</td><td>1.7</td><td>116 ± 16</td><td>112</td><td>1.4 ± 0.1</td><td>1.4</td><td>136 ± 20</td><td>131</td><td>63.1–70.4</td></tr><tr><td></td><td>II</td><td>1.8 ± 0.2</td><td>1.9</td><td>87 ± 9</td><td>85</td><td>1.4 ± 0.1</td><td>1.4</td><td>102 ± 12</td><td>100</td><td>0.4–2.1</td></tr><tr><td><i>o</i>-NH<sub>3</sub></td><td>I</td><td>5.4 ± 1.0</td><td>5.6</td><td>86 ± 6</td><td>86</td><td>3.1 ± 0.5</td><td>3.1</td><td>102 ± 8</td><td>102</td><td>0.2–27.6</td></tr><tr><td></td><td>II</td><td>5.4 ± 1.0</td><td>5.6</td><td>86 ± 6</td><td>86</td><td>3.1 ± 0.5</td><td>3.2</td><td>102 ± 8</td><td>101</td><td>0.1–0.2</td></tr><tr><td rowspan="4">B1i</td><td><i>p</i>-NH<sub>3</sub></td><td>I</td><td>0.6 ± 0.1</td><td>0.7</td><td>329 ± 117</td><td>253</td><td>0.6 ± 0.1</td><td>0.6</td><td>433 ± 152</td><td>319</td><td>1.8–2.4</td></tr><tr><td></td><td>II</td><td>0.7 ± 0.1</td><td>0.7</td><td>476 ± 181</td><td>92</td><td>0.6 ± 0.1</td><td>0.6</td><td>523 ± 189</td><td>92</td><td>0.01–11.2</td></tr><tr><td><i>o</i>-NH<sub>3</sub></td><td>I</td><td>1.5 ± 0.2</td><td>1.5</td><td>98 ± 7</td><td>98</td><td>1.1 ± 0.1</td><td>1.1</td><td>109 ± 7</td><td>108</td><td>38.2–42.7</td></tr><tr><td></td><td>II</td><td>1.5 ± 0.2</td><td>1.5</td><td>98 ± 7</td><td>98</td><td>1.1 ± 0.1</td><td>1.1</td><td>109 ± 7</td><td>108</td><td>9–34</td></tr><tr><td rowspan="4">B2a</td><td><i>p</i>-NH<sub>3</sub></td><td>I</td><td>3.2 ± 0.4</td><td>3.3</td><td>86 ± 8</td><td>84</td><td>2.3 ± 0.2</td><td>2.3</td><td>101 ± 10</td><td>99</td><td>1.7–23.2</td></tr><tr><td></td><td>II</td><td>3.7 ± 0.4</td><td>3.8</td><td>68 ± 6</td><td>67</td><td>2.5 ± 0.2</td><td>2.5</td><td>81 ± 7</td><td>80</td><td>0.01–0.02</td></tr><tr><td><i>o</i>-NH<sub>3</sub></td><td>I</td><td>11.7 ± 2.3</td><td>12.1</td><td>83 ± 6</td><td>83</td><td>6.0 ± 1.0</td><td>6.2</td><td>100 ± 8</td><td>99</td><td>1.3–13.6</td></tr><tr><td></td><td>II</td><td>11.7 ± 2.3</td><td>12.1</td><td>84 ± 6</td><td>83</td><td>6.0 ± 1.0</td><td>6.2</td><td>99 ± 7</td><td>98</td><td>0.2–10</td></tr><tr><td rowspan="4">B2b</td><td><i>p</i>-NH<sub>3</sub></td><td>I</td><td>2.4 ± 0.2</td><td>2.4</td><td>95 ± 11</td><td>93</td><td>1.8 ± 0.1</td><td>1.8</td><td>111 ± 13</td><td>109</td><td>6.7–11.2</td></tr><tr><td></td><td>II</td><td>2.7 ± 0.3</td><td>2.7</td><td>71 ± 6</td><td>71</td><td>1.9 ± 0.2</td><td>2.0</td><td>84 ± 8</td><td>83</td><td>0.01–6</td></tr><tr><td><i>o</i>-NH<sub>3</sub></td><td>I</td><td>11.9 ± 2.4</td><td>12.4</td><td>75 ± 4</td><td>75</td><td>6.1 ± 1.1</td><td>6.4</td><td>89 ± 6</td><td>88</td><td>1.1–51.0</td></tr><tr><td></td><td>II</td><td>12.6 ± 2.4</td><td>12.4</td><td>75 ± 4</td><td>75</td><td>6.1 ± 1.1</td><td>6.4</td><td>89 ± 6</td><td>88</td><td>0.5–3</td></tr></tbody></table>

Notes.

a Estimated by excluding the (7,7) line because no collision rate is given in LAMDA.

b Case I corresponds to a unity filling factor for all lines at any pixel; Case II corresponds to a unity filling factor for (1,1)–(3,3) and as 0.1 for (4,4)–(7,7) at any pixel.

c Assuming a narrow (“n”) FWHM line width of 4 km s⁻¹ (the median from the fittings to the observed spectrum in Table 1).

d Assuming a broad (“b”) FWHM line width of 9 km s⁻¹ adopted from Bachiller et al. (1993) and Lefloch et al. (2012).

e The PDF is not obvious, so the range of the best fits based on different line width assumptions is listed.
<figure><img src="image_13.png" /><figcaption>**Figure D2.** The best fit of the kinetic temperature maps and column density maps for *o*- and *p*-NH₃, derived from LVG MultiNest, by assuming a single component and unity filling factor for all lines toward each position. The results assuming the line width as 4 km s⁻¹ and 9 km s⁻¹ are shown in the upper and lower rows, respectively. The pixels where the NH₃ (5,5) and (6,6) lines show <5σ emission and CO(1–0) show < 8σ emission are blank for *p*- and *o*-NH₃, respectively. The gray dashed contours and the labeled clumpy substructures are the same as shown in Figure 1. These maps are derived by smoothing all lines to the same pixel size and angular resolution, i.e., 4''90 × 3''81, shown as the synthesized beam in the bottom left of each panel.</figcaption></figure>
<figure><img src="image_14.png" /></figure>

**Figure D3.** Possible correlation between variables. Values are extracted from pixels after smoothing the parameter maps to the same angular resolution and plotted with a bivariate Gaussian kernel density estimate as contours. B0, B1, and B2 are plotted, with the Spearman's rank correlation coefficient ρ given in red, blue, and yellowish green, respectively, while the correlation coefficient of the entire target is given in black. The pixels where the NH₃ (5,5) and (6,6) lines show <5σ emission are blank.

<figure><img src="image_15.png" /></figure>

satellite emission, the uncertainty for the (1,1) line is >20%. Nevertheless, it does not change our best-fit result significantly because the uncertainty for the (2,2) line is <30% and for the rest lines is <10% (less than the systematic uncertainty).

## ORCID iDs

S. Feng (冯思轶) https://orcid.org/0000-0002-4707-8409

H. B. Liu https://orcid.org/0000-0003-2300-2626

P. Caselli https://orcid.org/0000-0003-1481-7911

A. Burkhardt https://orcid.org/0000-0003-0799-0927

F. Du https://orcid.org/0000-0002-7489-0179

R. Bachiller https://orcid.org/0000-0002-5331-5386

C. Codella https://orcid.org/0000-0003-1514-3074

C. Ceccarelli https://orcid.org/0000-0001-9664-6292

## References

Arce, H. G., Santiago-García, J., Jørgensen, J. K., Tafalla, M., & Bachiller, R. 2008, ApJL, 681, L21

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bachiller, R., Martin-Pintado, J., & Fuente, A. 1993, ApJL, 417, L45

Bachiller, R., Pérez Gutiérrez, M., Kumar, M. S. N., & Tafalla, M. 2001, A&A, 372, 899

Benedettini, M., Viti, S., Codella, C., et al. 2007, MNRAS, 381, 1127

Benedettini, M., Viti, S., Codella, C., et al. 2013, MNRAS, 436, 179

Bouhafs, N., Rist, C., Daniel, F., et al. 2017, MNRAS, 470, 2204

Burkhardt, A. M., Dollhopf, N. M., Corby, J. F., et al. 2016, ApJ, 827, 21

Busquet, G., Lefloch, B., Benedettini, M., et al. 2014, A&A, 561, A120

Caselli, P., Bizzocchi, L., Keto, E., et al. 2017, A&A, 603, L1

Ceccarelli, C., Bacmann, A., Boogert, A., et al. 2010, A&A, 521, L22

Ceccarelli, C., Caselli, P., Fontani, F., et al. 2017, ApJ, 850, 176

Chiang, H.-F., Looney, L. W., Tobin, J. J., & Hartmann, L. 2010, ApJ, 709, 470

Codella, C., Benedettini, M., Beltrán, M. T., et al. 2009, A&A, 507, L25

Codella, C., Ceccarelli, C., Bianchi, E., et al. 2020, A&A, 635, A17

Codella, C., Ceccarelli, C., Caselli, P., et al. 2017, A&A, 605, L3

Codella, C., Fontani, F., Ceccarelli, C., et al. 2015, MNRAS, 449, L11

Codella, C., Lefloch, B., Ceccarelli, C., et al. 2010, A&A, 518, L112

Cohen, J. 1988, Statistical Power Analysis for the Behavioral Sciences (London: Routledge)

Crapsi, A., Caselli, P., Walmsley, M. C., & Tafalla, M. 2007, A&A, 470, 221

Danby, G., Flower, D. R., Valiron, P., Schilke, P., & Walmsley, C. M. 1988, MNRAS, 235, 229

Downes, T. P., & Ray, T. P. 1999, A&A, 345, 977

Estalella, R. 2017, PASP, 129, 025003

Faure, A., Hily-Blant, P., Le Gal, R., Rist, C., & Pineau des Forêts, G. 2013, ApJL, 770, L2

Feng, S., Codella, C., Ceccarelli, C., et al. 2020, ApJ, 896, 37
Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, OJAp, 2, 10

Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2006, A&A, 449, 621

Fontani, F., Codella, C., Ceccarelli, C., et al. 2014a, ApJL, 788, L43

Fontani, F., Sakai, T., Furuya, K., et al. 2014b, MNRAS, 440, 448

Frank, A., Ray, T. P., Cabrit, S., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 451

Ginsburg, A., & Mirocha, J. 2011, PySpecKit: Python Spectroscopic Toolkit, Astrophysics Source Code Library, ascl:1109.001

Goddi, C., Greenhill, L. J., Humphreys, E. M. L., Chandler, C. J., & Matthews, L. D. 2011, ApJL, 739, L13

Goldreich, P., & Scoville, N. 1976, ApJ, 205, 144

Gómez-Ruiz, A. I., Codella, C., Lefloch, B., et al. 2015, MNRAS, 446, 3346

Gómez-Ruiz, A. I., Hirano, N., Leurini, S., & Liu, S.-Y. 2013, A&A, 558, A94

Gueth, F., Guilloteau, S., & Bachiller, R. 1996, A&A, 307, 891

Gueth, F., Guilloteau, S., & Bachiller, R. 1998, A&A, 333, 287

Harju, J., Daniel, F., Sipilä, O., et al. 2017, A&A, 600, A61

Ho, P. T. P., & Townes, C. H. 1983, ARA&A, 21, 239

Juvela, M., & Ysard, N. 2011, ApJ, 739, 63

Lee, K., Looney, L., Johnstone, D., & Tobin, J. 2012, ApJ, 761, 171

Lefloch, B., Bachiller, R., Ceccarelli, C., et al. 2018, MNRAS, 477, 4792

Lefloch, B., Cabrit, S., Busquet, G., et al. 2012, ApJL, 757, L25

Lefloch, B., Cabrit, S., Codella, C., et al. 2010, A&A, 518, L113

Lefloch, B., Ceccarelli, C., Codella, C., et al. 2017, MNRAS, 469, L73

Looney, L. W., Tobin, J. J., & Kwon, W. 2007, ApJL, 670, L131

Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Neufeld, D. A., DeWitt, C., Lesaffre, P., et al. 2019, ApJL, 878, L18

Offer, A., & Flower, D. R. 1989, JPhB, 22, L439

Perley, R. A., & Butler, B. J. 2017, ApJS, 230, 7

Pety, J. 2005, in SF2A-2005: Semaine de l'Astrophysique Francaise, ed. F. Casoli (Les Ulis: EdP-Sciences), 721

Podio, L., Codella, C., Gueth, F., et al. 2016, A&A, 593, L4

Podio, L., Lefloch, B., Ceccarelli, C., Codella, C., & Bachiller, R. 2014, A&A, 565, A64

Rist, C., Alexander, M. H., & Valiron, P. 1993, JChPh, 98, 4662

Rosolowsky, E. W., Pineda, J. E., Foster, J. B., et al. 2008, ApJS, 175, 509

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369

Shirley, Y. L. 2015, PASP, 127, 299

Smith, M. D., Suttner, G., & Yorke, H. W. 1997, A&A, 323, 223

Sobolev, V. V. 1960, Moving Envelopes of Stars (Cambridge, MA: Harvard Univ. Press)

Spezzano, S., Codella, C., Podio, L., et al. 2020, A&A, 640, A74

Tafalla, M., & Bachiller, R. 1995, ApJL, 443, L37

Teague, R., & Foreman-Mackey, D. 2018, RNAAS, 2, 173

Tobin, J. J., Cox, E. G., & Looney, L. W. 2022, ApJ, 928, 61

Tobin, J. J., Hartmann, L., Looney, L. W., & Chiang, H.-F. 2010, ApJ, 712, 1010

Umemoto, T., Mikami, H., Yamamoto, S., & Hirano, N. 1999, ApJL, 525, L105

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Viti, S., Jimenez-Serra, I., Yates, J. A., et al. 2011, ApJL, 740, L3

Walmsley, C. M., & Ungerechts, H. 1983, A&A, 122, 164

Zhang, Q., Ho, P. T. P., & Wright, M. C. H. 2000, AJ, 119, 1345

Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2019, ApJ, 879, 125