Searching for further evidence for cloud-cloud collisions in L1188

In order to search for further observational evidence of cloud-cloud collisions in one of the promising candidates, L1188, we carried out observations of multiple molecular lines toward the intersection region of the two nearly orthogonal filamentary molecular clouds in L1188. Based on these observations, we find two parallel filamentary structures, both of which have at least two velocity components being connected with broad bridging features. We also found a spatially complementary distribution between the two molecular clouds, as well as enhanced $^{13}$CO emission and $^{12}$CO self-absorption toward their abutting regions. At the most blueshifted velocities, we unveil a 1~pc-long arc ubiquitously showing $^{12}$CO line wings. We discover two 22 GHz water masers, which are the first maser detections in L1188. An analysis of line ratios at a linear resolution of 0.2 pc suggests that L1188 is characterised by kinetic temperatures of 13--23~K and H$_{2}$ number densities of 10$^{3}$--10$^{3.6}$ cm$^{-3}$. On the basis of previous theoretical predictions and simulations, we suggest that these observational features can be naturally explained by the scenario of a cloud-cloud collision in L1188, although an additional contribution of stellar feedback from low-mass young stellar objects cannot be ruled out.


Introduction
Cloud-cloud collisions have been proposed to play an important role in molecular cloud evolution, star formation, and, on a grander scale, in galaxy evolution (e.g., McKee & Ostriker 2007). Theoretically, many studies suggest that the formation of giant molecular clouds could be due to agglomeration via successive inelastic cloud-cloud collisions (e.g., Kwan 1979;Roberts & Stewart 1987). On the other hand, cloud-cloud collisions are thought to be important for star formation. Previous studies suggested that a cloud-cloud collision could induce star-forming seeds in the shock-compressed interface (e.g., Habe & Ohta 1992;Anathpindika 2010). Such a process would favor the formation of massive molecular filaments (Inoue et al. 2018) from which massive stars can fragment (Inoue & Fukui 2013). Numerical simulations suggest that such a mechanism greatly enhances star formation rates and efficiencies (Wu et al. 2017). Furthermore, cloud-cloud collisions may make significant contributions to global star formation rates in spiral galaxies because such collisions are expected to occur frequently (once every 8-10 Myr in spiral arms, Dobbs et al. 2015). In addition, Li (2017) argued that kinematic dissipation by cloud-cloud collisions might lead The reduced datacubes are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http:// cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/632/A115 to the formation of giant massive clumpy structures in highredshift gas-rich galaxies. Identifying cloud-cloud collisions is indispensable in order to verify these theories, while it is still challenging.
From an observational point of view, a few nearby molecular clouds are believed to exhibit features of cloud-cloud collisions (Duarte-Cabral et al. 2011;Nakamura et al. 2012). These candidates are mainly identified by blueshifted and redshifted velocity fields. Toward the rather distant infrared dark cloud G035.39−00.33 (∼2.9 kpc), parsec-scale SiO emission (e.g., Jiménez-Serra et al. 2010) and the spatio-kinematic offset between different tracers (e.g., Henshaw et al. 2013;Bisbas et al. 2018) are also proposed to result from cloud-cloud collisions. Numerical simulations have shown that broad bridge features in position-velocity diagrams can serve as an isolated signature of cloud-cloud collisions (Haworth et al. 2015a,b). These criteria have already led to many cloud-cloud collision candidates, such as NGC 3603 (Fukui et al. 2014), RCW 38 (Fukui et al. 2016), M 20 (Torii et al. 2011(Torii et al. , 2017, L1188 (Gong et al. 2017), RCW 79 (Ohama et al. 2018), NGC 6618/M17 (Nishimura et al. 2018), GM24 (Fukui et al. 2018a), NGC 6334/NGC 6357 (Fukui et al. 2018b), and the references therein. However, there are still some uncertainties as to whether these candidates are true collisions or not; because cloud morphologies and internal motions are usually complex and a competitive process, the feedback A&A 632, A115 (2019) by massive stars, complicates the picture. Therefore, additional observational features of cloud-cloud collisions is beneficial for robust identifications.
L1188, located at a distance of ∼800 pc (Gong et al. 2017;Szegedi-Elek et al. 2019), was first mapped in 13 CO (1-0) with the Nagoya-4 m telescope (Abraham et al. 1995). A subsequent study suggests that L1188 is one of the most promising cloudcloud collision candidates with a total molecular gas mass of 3.9 × 10 3 M and without massive stars (M > 8 M ) inside (Gong et al. 2017). L1188 is thus much less affected by the energetic feedback from massive stars when compared with the aforementioned candidates associated with massive stars. Furthermore, previous CO observations demonstrate that L1188 has a relatively simple structure, which consists of two nearly orthogonal filamentary molecular clouds, namely L1188a and L1188b (Gong et al. 2017). These properties allow the assignation of observational features to different phenomena. In order to search for further evidence of cloud-cloud collisions, we performed a multiple molecular line study toward the intersection region in L1188, covering a size of 20 × 20 (4.7 pc × 4.7 pc, see Fig. 1).

Observations and data reduction
We performed mapping observations of five different molecules with four different telescopes. Table 1 summarizes these observations.

PMO-13.7 m observations
The intersection region was observed in HCO + (1-0), SiO (2-1) and CH 3 OH (8 0 -7 1 A + ) with the Purple Mountain Observatory 13.7 m telescope (PMO-13.7 m) from 2018 May 1 to 15 (project code: 18A003). We employed the 3 × 3 beam sideband separation Superconducting Spectroscopic Array Receiver (SSAR) and fast Fourier transform spectrometers (FFTSs) for our observations (Shan et al. 2012). Each of the 18 FFTS modules provides 16 384 channels, covering an instantaneous bandwidth of 1 GHz. This results in a channel spacing of 61 kHz, that is, 0.20 km s −1 at 86 GHz. The region was observed in the on-the-fly mode with a scanning speed of 50 per second and a dump time of 0.3 s. These observations took around 40 hour in total.
The antenna temperature, T * A , was obtained with the standard chopper wheel method (Ulich & Haas 1976). Throughout this work, the antenna temperature is converted into the main beam brightness temperature scale, T mb , with the relation, where η mb is the main beam efficiency. The main beam efficiency is 57% at 86 GHz according to the telescope's status report 1 . The flux calibration error is estimated to be roughly 10%. During the observations, system temperatures were 145-301 K on a T * A scale. The half-power beam width (HPBW) is about 60 at the observed frequencies. The pointing was found to have an uncertainty of 5 , which is less than 1/10 of the HPBW.
The GILDAS 2 software was used to reduce the data (Pety 2005). Raw data were gridded into regular separations of 30 between two adjacent pixels. The achieved sensitivities are ∼0.08 K on a T mb scale for SiO (2-1) at a channel width of 0.20 km s −1 .
CO (1-0) and 13 CO (1-0) data, also observed with the PMO-13.7 m telescope as part of the Milky Way Imaging Scroll Painting (MWISP 3 ; Su et al. 2019) project, have already been presented in Gong et al. (2017) and are also used in our analysis.

IRAM-30 m observations
CO (2-1) observations (project code: 138-16, 023-18) were performed with the Institut de Radioastronomie Millimétrique 30 m (IRAM-30 m) telescope and the 9 dual-polarization pixel HEterodyne Receiver Array (HERA, Schuster et al. 2004) in 2017 March and 2018 August. The selected region was observed in the on-the-fly mode, scanning along right ascension and declination directions. The dump time is typically 0.5 s with a scanning speed of 6 per second. The calibration was done every 6-8 min. FFTSs were used as backends to record signals from both linear polarizations. These FFTSs provide a channel spacing of 195 kHz, corresponding to a channel width of 0.25 km s −1 at 230 GHz. The focus was adjusted via observations of the radio galaxy 2021+614. Pointing, which was regularly checked on nearby strong continuum sources, was found to have a typical rms uncertainty of ∼5 . The spectral data are expressed in units of main beam brightness temperature by adopting a forward efficiency of 92% and a main beam efficiency of 57%. The HBPW is about 11 at 230 GHz. The typical noise is about 0.5 K at a channel width of 0.25 km s −1 . These IRAM-30 m data were also analyzed with the GILDAS software.

JCMT-15 m observations
We carried out 12 CO (3-2) and 13 CO (3-2) observations (project code: M17AP002) toward the central 20 × 20 region in L1188 with the James Clerk Maxwell Telescope (JCMT) from 2017 May 31 to June 28. The 16-element Heterodyne Array Receiver Program (HARP) was used as front-end, while the Auto-Correlation Spectral Imaging System (ACSIS) was employed as back-end (Buckle et al. 2009). ACSIS provides 7172 channels,  each of which is 30.5 kHz-wide. This is equivalent to a width of 0.026 km s −1 at 345 GHz. We used the on-the-fly mode to scan the observed region in orthogonal directions. At 345 GHz, the HPBW is about 14. 5 and the main beam efficiency is 61%.
Pointing was checked using JCMT standard sources, resulting in a typical rms uncertainty of ∼3 . The observed data were reduced with the ORAC Data Reduction pipeline (ORAC-DR, Jenness et al. 2015) of the Starlink 4 software (Currie et al. 2014). For the following analysis, all spectra have been binned to a channel width of 0.27 km s −1 to improve signal-to-noise ratios (S/N). The typical noise levels are about 1 and 0.7 K at a channel width of 0.27 km s −1 for 12 CO (3-2) and 13 CO (3-2), respectively. Velocities are given with respect to the local standard of rest (LSR) throughout this work.

Effelsberg-100 m observations
We carried out H 2 O (6 1,6 -5 2,3 ) measurements (project code: 22-16, 59-17) in a position-switching mode with the doublebeam and dual-polarization K-band receiver of the 100-m telescope at Effelsberg/Germany 5 during 2019 April. The two beams are at the same declination but separated by 192 in right ascension. We performed a raster map toward the region showing 12 CO wing emission (see results given below) with a grid spacing of 15 , leading to a coverage of 360 × 75 . The FFTS provides a bandwidth of 300 MHz, consisting of 65 536 channels with a channel spacing of 4.6 kHz which is equivalent to 0.06 km s −1 at 22 GHz. The full width at half maximum (FWHM) beamsizes are about 40 at 22 GHz. NGC 7027 was used for pointing, focus, and flux calibration, and the flux calibration accuracy is estimated to be ±20%. Nearby pointing results were obtained toward the radio galaxy 2021+614, and the rms pointing uncertainty was found to be 5 . Typical system temperatures on a T * A scale are 49-80 K. The main beam efficiency was 59.7%. The conversion factor from flux density (S ν ) to main beam temperature (T mb ) is T mb /S ν = 1.75 K Jy −1 at 22 GHz. These Effelsberg-100 m data were also analyzed with the GILDAS software.

Molecular distribution and kinematics
Our new observations have led to the detection of extended emission in 12 CO (2-1), 12 CO (3-2), HCO + (1-0), and compact filamentary structures in 13 CO (3-2). The integrated intensity maps are shown in Fig. 2. The spatial distributions in 12 CO (2-1) and 12 CO (3-2) are similar to 12 CO (1-0), but both have a factor of ∼4 higher angular resolutions compared with previous 12 CO (1-0) observations (e.g., Gong et al. 2017). Previously, only five positions had been observed in HCO + (1-0) (Verebélyi et al. 2013), so the distribution of HCO + (1-0) emission is revealed for the first time by our measurements. We designate the three clumpy regions, which are indicated in Fig. 2c, as C1, C2, and C3. The 13 CO (3-2) data have rather low S/N, so their resolution has been smoothed to 55 to achieve a higher dynamic range hereafter. It turns out that the YSO candidates appear to cluster toward region C2 where strong and clumpy HCO + (1-0) and 13 CO (3-2) emission is present. In contrast, only one YSO candidate is associated with region C1, while no YSO candidate is found in region C3. This may indicate that region C2 is the most evolved region, followed by region C1 and C3. Figure 3 presents the intensity maps of 12 CO (1-0), 13 CO (1-0) (in contours), 12 CO (2-1), 12 CO (3-2), and 13 CO (3-2) (in contours) with different integrated velocity ranges. In this work, we refer to the emission integrated from −14 to −11 km s −1 as L1188a (see Figs. 3a, d, and g), the emission integrated from −11 to −8 km s −1 as the intermediate-velocity region (see Figs. 3b, e, and h), and the emission integrated from −8 to −5 km s −1 as L1188b (see Figs. 3c, f, and i). Thanks to the achieved high angular resolutions of 12-15 , 12 CO (2-1) and 12 CO (3-2) intensity maps show more fluffy structures than our previous 12 CO (1-0) maps. Short wispy and filamentary structures are ubiquitously present at the edges of L1188a, the intermediate-velocity region, and L1188b (see Figs. 3d-i); also, they appear to be randomly oriented. In the intermediate-velocity region, cavity-like structures are more evident than in previously obtained images of molecular line emission. YSO candidates seem to be spatially associated with L1188a or the intermediate-velocity region rather than L1188b (see . Given that the large-scale cloud structures have been investigated in 12 CO (1-0) and 13 CO (1-0) (Gong et al. 2017), we therefore concentrate on smaller structures in this work. In L1188b, there are at least two long filamentary structures, denoted as F1 and F2 in Fig. 3f, which are nearly parallel to each other and roughly separated by 0.8-1 pc on a linear scale. F1 appears to be filamentary, while the morphology of F2 displays an oblate "X"-like shape. Both of them are elongated and measure ∼2.7 pc in the northeast-southwest direction. Furthermore, they are roughly perpendicular to the long axis of L1188b, but nearly parallel to the long axis of L1188a. The crest of F1 is obtained from the CO (2-1) integrated intensity map of L1188b ( Fig. 3f) with the discrete persistent structure extractor (DisPerSe 6 ) algorithm (Sousbie 2011;Sousbie et al. 2011). Both the persistence and robustness thresholds were set to 3 K km s −1 to isolate filamentary structures. They were assembled to form a single filament when they form an angle smaller than 70 • . The resulting crest is shown as the red line in Fig. 3f.
Figures 4a-c present the position-velocity (PV) diagram of the CO lines along the crest of F1. One can clearly distinguish two velocity components with systemic velocities of about −11 and −8 km s −1 , which well match the systemic velocities of L1188a and L1188b (Gong et al. 2017). Furthermore, there are two bridging features at offsets of 3.5 and 9.5 that connect these two velocity components (indicated by the two green boxes in Figs. 4b and c), supporting that the two velocity components are associated with each other. Figure 4d shows the variation of the normalized integrated intensities of CO (2-1) and CO (3-2) along the crest of F1. This variation appears to be periodic with a sine wave-like shape, and one can identify at least four peaks indicated by the arrows in Fig. 4d. The four peaks correspond to offsets of 0.9, 1.7, 2.5, and 3.2 pc, respectively, suggesting that the peak-to-peak separations are nearly constant, that is, ∼0.8 pc. The maximum-to-minimum ratios are about 2.5-3.0 in this intensity profile.
Figures 5a-c show the PV diagrams of F2 in the CO lines. One can identify at least three velocity components with LSR velocities of −11.5, −10, and −7 km s −1 . The systemic velocities of L1188a and L1188b just lie within the gaps of the three velocity components. Furthermore, these three components are A115, page 4 of 16 Y. Gong et al.: Searching for further evidence for cloud-cloud collisions in L1188 In panels a-c, the contours represent the integrated intensity maps of 13 CO (1-0), and their integrated velocity ranges are the same as those of the corresponding CO maps. The contours start from 0.9 K km s −1 (5σ), and increase by 0.9 K km s −1 . Panels d-f and g-i: similar to panels a-c but for 12 CO (2-1) and 12 CO (3-2), respectively. In panels g-i, the contours represent the smoothed (55 ) integrated intensity maps of 13 CO (3-2), and start from 0.3 K km s −1 (5σ), and increase by 0.3 K km s −1 . Panels g-i have been scaled by a factor of two to match the color bar. In panels a-c, the green filled circles represent the YSO candidates from Szegedi-Elek et al. (2019) and the two yellow pentagrams represent the two 22 GHz water maser positions discussed in Sect. 3.3. The filamentary structures F1 and F2 are outlined in panel f. The beam size for CO (1-0), CO (2-1), and CO (3-2) are shown in the lower left corner of panels a, d, and g, respectively. In all panels, the (0, 0) offset corresponds to α J2000 = 22 h 18 m 00. s 44, δ J2000 = 61 • 49 03. 7. connected by bridging features over a linear scale of ∼1.1 pc (or 5 ), indicating that they are physically interacting. Consequently, the oblate X-like shape may be due to a superposition of the three components. Szegedi-Elek et al. (2019) identify a filamentary YSO distribution indicated in Fig. 1 which is nearly parallel to the long axis of L1188b and perpendicular to the long axis of L1188a, F1, and F2. Figure 6 presents a spatially complementary distribution between L1188a and L1188b. L1188a is sandwiched in L1188b, that is, the northeast end of L1188a matches the centric concave shape of L1188b very well (see the yellow dashed arc in the right panel of Fig. 6), while F1 matches the southern edge of L1188a. Despite the fact that complementary distributions have also been found in other cloud-cloud collision candidates associated with HII regions (e.g., Torii et al. 2017Torii et al. , 2018Fukui et al. 2018a;Fujita et al. 2019), this is the first time that such a feature is clearly present in a region without massive stars inside. Furthermore, enhanced 13 CO emission is curved in region C3 (see the yellow box in the right panel of Fig. 6), which is near the centric concave shape of L1188b. Toward region C3, 12 CO spectra look quite different from 13 CO spectra (see the top panel of Fig. 7). The 12 CO lines are asymmetric with enhanced redshifted emission with a A115, page 5 of 16 The contours start at 0.9 K (3σ) and increase by 0.9 K. The offsets are given with respect to the southwestern end of the PV cut. The systemic velocities of L1188a and L1188b are indicated by two blue vertical dashed lines. (b) Similar to a, but for 12 CO (2-1). (c) Similar to a, but for 12 CO (3-2). In panels b-c, the bridging features are indicated by the green dashed boxes. (d) Normalized integrated intensities of 12 CO (2-1) (blue points) and 12 CO (3-2) (red points) as a function of the offsets in the crest of F1. The offsets are calculated with respect to the southwestern end of F1.
slope that increases in intensity from lower to higher velocities, while the 13 CO lines appear to be Gaussian. The bottom panel of Fig. 7 shows a grid of the 12 CO (3-2) and 13 CO (3-2) line profiles. Most of the 12 CO (3-2) line profiles display two peaks, while 13 CO (3-2) lines peak toward the dips in the corresponding 12 CO (3-2) lines. Similar line profiles are also observed in 12 CO (1-0), 12 CO (2-1), and 13 CO (1-0). This is indicative of widespread 12 CO self-absorption in this region.
3.2. Discovery of a 1 pc-long arc Figure 8 shows the intensity maps of the three lowest J 12 CO transitions integrated from −20.5 to −12.5 km s −1 . This demonstrates that the most blueshifted component exhibits an arc-like filamentary structure which roughly measures ∼1 pc in length and ∼0.16 pc in width. This results in a high aspect ratio of ∼6. This arc is located at the northeastern end of F2, which is nearly perpendicular to the long axis of F2, and extends nearly parallel to the centric concave structure of L1188b. Furthermore, the radius of curvature is estimated to be about 0.6 pc for this arc, which is comparable to the centric concave shape of L1188b. Figure 9 shows the PV diagram along the arc, in which we can identify two velocity components centered at about −11 and −8 km s −1 , similar to the systemic velocities of L1188a and L1188b, at offsets of 2 -6 . The two velocity components are connected to each other by bridging features that extend across the whole arc structure. This indicates that the origin of the 1 pc-long arc could be related to the interaction between L1188a and L1188b. Moreover, weak low-velocity wing emission at velocities of <−12.5 km s −1 is also present over the offsets of 2 -6 and most of the emission can reach −16 km s −1 (indicated by the red dashed box in Fig. 9b). This wing emission is beyond the 13 CO velocity range, and extends ∼1 pc in length, indicating that this emission can be caused by large-scale shocks. The A115, page 6 of 16  In all panels, turquoise and red represent the intensity maps integrated from −14 to −11 km s −1 and from −8 to −5 km s −1 , respectively. In the right panel, the overlaid contours are the same as those in Fig. 3h. The curved region that shows enhanced 13 CO emission is indicated by the yellow box, while the centric concave shape of L1188b is indicated by the yellow dashed arc. In each panel, the two pentagrams represent the two 22 GHz water maser positions discussed in Sect. 3.3.
high-velocity wing emission is also uncovered, and is indicated by the three red arrows in Fig. 9c. We find high-velocity wing emission covering a linear scale of ∼0.1 pc, which is similar to the typical scale of molecular outflows in low-mass starforming regions (e.g., Arce et al. 2007). We also note that these wing emissions are only seen in the blueshifted velocity range, whereas the potential red-shifted emission may mix with the bridging features and the −8 km s −1 velocity component, making it undistinguishable in Fig. 9. At offsets of 6 -7 , one can find a velocity gradient of 3 km s −1 pc −1 .
Assuming that the CO wing emission integrated from −20.5 to −12.5 km s −1 is optically thin and that it follows local thermodynamic equilibrium (LTE), we are able to derive 12 CO column densities and rotational temperatures in the arc structure with the rotational diagram. Prior to the rotational diagram analysis, we convolved 12 CO (2-1) and 12 CO (3-2) data to the angular resolution of 55 , and all data were then linearly interpolated to the same grid. We used the following formula (e.g., Herbst & van Dishoeck 2009;Gong et al. 2015a): where k is the Boltzmann constant, W is the 12 CO integrated intensity, ν is the rest frequency, µ is the permanent dipole A115, page 7 of 16 A&A 632, A115 (2019) Fig. 7. Top: 12 CO (1-0), 13 CO (1-0), 12 CO (2-1), 12 CO (3-2), and 13 CO (3-2) spectra averaged over the region indicated by the yellow box in the right panel of Fig. 6. Bottom: grid of the 12 CO (3-2) (blue) and 13 CO (3-2) (black) line profiles in the yellow box indicated in Fig. 6. These spectra share the same velocity range between −15 and −3 km s −1 , and intensity range between −2 and 8 K. The 13 CO (3-2) spectra have been scaled by a factor of two. These spectra have been convolved to an angular resolution of 55 . moment, S is the transition's intrinsic strength, N tot is the total column density, T rot is the rotational temperature, Q is the partition function, and E u is the upper level energy of the transition. The values of Q and µ are taken from the CDMS 7 (Müller et al. 2005). A linear fit was performed toward every pixel to obtain their 12 CO rotational temperatures and column densities, and an example toward the peak in the arc structure is shown in Fig. 10a. In the following analysis, we only account for pixels with at least 5σ in all three 12 CO transitions. The resulting 12 CO rotational temperature and column density maps are shown in Figs. 10b and c. The derived rotational temperatures range from 8 to 12 K, while 12 CO column densities are within the range of (1.1-7.6) × 10 15 cm −2 . Assuming a typical 12 CO abundance of 8 × 10 −5 (Frerking et al. 1982) and the mean molecular weight of 2.8 per hydrogen molecule (e.g., Kauffmann et al. 2008), the total gas mass of the arc structure is about 0.34 M above the 5σ threshold. We also note that this arc is not detected in 13 CO lines. The 3σ upper limit of the 13 CO (1-0) integrated intensity is 1.2 K km s −1 , corresponding to a 13 CO column density of 6 × 10 14 cm −2 at an assumed excitation temperature of 10 K under LTE conditions. 7 https://zeus.ph1.uni-koeln.de/cdms/catalog

Water maser emission in L1188
Theoretical models have predicted that water masers can be created by cloud-cloud collisions (Tarter & Welch 1986), but direct observational evidence is still rare. We searched for 22 GHz water maser emission in a 360 × 75 area of region C1 close to the 1 pc-long shocked arc (see Figs. 8-10), and detected water maser emission in two positions. This is the first maser detection in L1188. Their spectra and spatial distributions are shown in Figs. 11a-d, and the two peak positions, M1 and M2, were mannually determined by the naked eyes. The observed maser properties are tabulated in Table 2.
The 22 GHz water maser spectral profiles are different from (quasi-)thermal line profiles of other molecules. The former ones display at least four peaks, whereas the latter ones only show two peaks at −11 and −9 km s −1 (see Figs. 11a and b). The peak velocities of the maser spectra are 5 km s −1 off the 13 CO (1-0) peak velocities. These peaks in the maser spectra can be divided into two groups separated by a velocity difference of >11 km s −1 , that is, a blueshifted group at a velocity range between −18 and −13 km s −1 , and a redshifted group at a velocity range between −2 and 0 km s −1 . The strongest maser component has an LSR velocity of about −15.2 km s −1 , which is blueshifted with respect to the systemic velocities of L1188a and L1188b that are traced by 13 CO. Its corresponding main beam brightness temperature is higher than the peak intensities of the three 12 CO lines. A comparison of the two linearly polarized components of these maser features shows no significant differences in line shape and flux density. Their isotropic H 2 O luminosities reach ∼2 × 10 −7 L , which is nearly an order of magnitude higher than those in nearby low-mass star-forming regions (Furuya et al. 2003), but lower than those in most high-mass star-forming regions (Urquhart et al. 2011).
Figures 11c and d show that the 22 GHz water maser emission mainly arises from two positions, M1 and M2, which are indicated by the two red crosses. They are separated by 48 , corresponding to ∼0.2 pc. Near M1, we find the Class I YSO WISEJ221729.87+615039.8 (details are described in Appendix A), which is seen in emission at 22 µm (see Figs. 11c and d). The separations between this YSO and the two maser positions are about 20 (or 0.08 pc) and 60 (or 0.23 pc), respectively. The offsets and the spectral profiles are similar to those found in other regions where masers are excited by YSO outflows (Furuya et al. 2001(Furuya et al. , 2003. This indicates that the detected masers could be excited by YSO outflows rather than cloud-cloud collisions. However, we cannot resolve outflows from this YSO in its northeastern direction where masers are located. This is likely due to insufficient sensitivities. Since previous distance measurements of L1188 (∼800 pc, see Sect. 1) still have large uncertainties (Gong et al. 2017;Szegedi-Elek et al. 2019), future trigonometric parallax measurements of the newly detected 22 GHz water masers with Very Long Baseline Interferometry (VLBI) measurements may provide distance estimates of a much higher quality.

Non-detection of SiO (2-1) and CH 3 OH (8 0 -7 1 A + ) emission
Since SiO (2-1) is a commonly used shock tracer (e.g., Martin-Pintado et al. 1992), its emission may trace the shock created by the cloud-cloud collision. The CH 3 OH (8 0 -7 1 A + ) line is a well-known Class I methanol maser transition that is widespread in the Milky Way (e.g., Yang et al. 2017) and frequently found to be associated with molecular outflows from YSOs. However,  it has also been found to occur in supernova remnant molecular cloud interactions (Frail 2011), and it is thought to be able to form via cloud-cloud collisions (e.g., Menten 1996;Leurini et al. 2016). We therefore searched for emission in the SiO and CH 3 OH lines in the region in which emission from other molecules is detected, but we were only able to establish an upper limit of 0.54 K km s −1 (3σ) for the intensity of SiO (2-1) integrated from −14 to −5 km s −1 . The 3σ upper limit for the SiO column density is ∼1.2 × 10 12 cm −2 by assuming LTE and an excitation temperature of 10 K. The upper limit of the main beam brightness temperature for CH 3 OH (8 0 -7 1 A + ) is 0.24 K (3σ), corresponding to a peak flux density of 9.8 Jy.
Here, we only took the pixels with at least 5σ detections into account. The line ratio maps of the observed region are shown in Fig. 12. It is important to note that R 12 2−1/1−0 varies from 0.7 to 1.2 with a median value of 0.9, suggesting a small variation in R 12 2−1/1−0 . This is likely because both 12 CO (1-0) and 12 CO (2-1) have very high opacities and nearly the same excitation temperature. Additionally, R 12 3−2/1−0 lies within a range of 0.2-0.7 with a median value of 0.5, generally lower than R 12 2−1/1−0 . This can be due to 12 CO (3-2)'s relatively lower opacity and/or its subthermal state. Furthermore, R 13 3−2/1−0 is found to range from 0.1 to 0.7. Since both 13 CO (3-2) and 13 CO (1-0) are likely optically thin, R 13 3−2/1−0 allows us to determine the excitation temperature in case of LTE. The observed ratios suggest an excitation temperature range of 6-10 K. However, the excitation temperature should be lower than the kinetic temperature since 13 CO (3-2) becomes subthermal (see discussions below). When the variation of 12 C/ 13 C isotopic ratios can be neglected, R 13/12 3−2 is sensitive to CO column densities and number densities. We find that the R 13/12 3−2 distribution is quite similar to that of the 13 CO emission; R 13/12 3−2 varies from 0.03 to 0.23. The highest R 13/12 3−2 corresponds to region C3, which shows enhanced 13 CO (3-2) emission and 12 CO (3-2) self-absorption (Figs. 6 and 7).
In order to derive physical properties, we employed a non-LTE analysis using the RADEX 8 code (van der Tak et al. 2007). The Einstein A coefficients are from the Leiden Atomic and Molecular Database (LAMDA 9 , Schöier et al. 2005), while collision rates are taken from Yang et al. (2010). An expanding spherical geometry, also known as large velocity gradient approximation (Sobolev 1960), has been assumed for the calculations, so the escape probability (β) is β = 1−exp(−τ) τ , where τ is the optical depth. The ortho-to-para ratio is fixed to three for H 2 . A total of 41 12 CO energy levels and 41 13 CO energy levels are included in the calculations. Since 12 CO (1-0) and 12 CO (2-1) tend to have very high opacities indicated by R 12 2−1/1−0 ∼1, we only used R 13 3−2/1−0 and R 13/12 3−2 to model physical properties by assuming that the three transitions can trace almost the same portion of molecular gas, following the prescription of previous studies (e.g., Nishimura et al. 2015). In order to fit the two ratios simultaneously, we assumed the isotopic ratio 12 C/ 13 C to be 70 (e.g., Gong et al. 2015b;Li et al. 2016), the 12 CO abundance to be 8 × 10 −5 (Frerking et al. 1982), and the velocity gradient d /dr to be 3 km s −1 pc −1 , which is based on the typical 13 CO line width and size of the cloud. During modeling, we considered a regular grid with kinetic temperature ranging from 5 to 100 K with a step of 1 K and the H 2 number density log[n(H 2 ) cm −3 ] varying from 2 to 6 with a step of 0.1. In order to obtain physical conditions for the observed region, we make use of the χ 2 analysis to evaluate different models toward every pixel. We note that χ 2 is defined as where R obs(i) and R model(i) represent the observed and modeled line ratios. Also, σ obs(i) represents the 1σ error in R obs(i) , which includes the flux calibration uncertainties and is calculated with the error propagation formula. The best fit is derived by minimizing χ 2 . The 1σ confidence range can be derived by χ 2 < 2.3 with two free parameters. Figure 13a presents an example of the modeling results toward the peak position of region C2 indicated by the cross in Figs. 12c and d. We can see that the kinetic temperature and H 2 number density are well constrained to be 18-25 K and 10 3.5−3.7 cm −3 , respectively, within the 1σ confidence range, while the best fit is 21 K for the kinetic temperature and 10 3.6 cm −3 for the H 2 number density. Our results also suggest that R 13/12 3−2 can be regarded as a reliable tracer to derive H 2 number densities because of its weak dependence on kinetic temperature. However, the assumed d /dr may have a large uncertainty. Since 13 CO (3-2) and 13 CO (1-0) are likely optically thin in the observed region, the d /dr values do not affect the modeled R 13 3−2/1−0 ratio. In order to quantify the impact of the variation of d /dr on R 13/12 3−2 , we modeled R 13/12 3−2 with a d /dr range of 1-10 km s −1 pc −1 and a fixed kinetic temperature of 21 K. The result is shown in Fig. 13b. This suggests that the derived H 2 number density becomes higher if one increases d /dr. Nevertheless, the derived H 2 number density only varies by a factor of two within a d /dr range of 1-10 km s −1 pc −1 .
Following the method described above, we determined the kinetic temperature and H 2 number density of each pixel with the line ratios R 13 3−2/1−0 and R 13/12 3−2 at a linear resolution of 0.2 pc, that is, an angular resolution of 55 . Their distributions are shown in Fig. 14. Figure 14a shows that most of the derived kinetic A115, page 11 of 16 A&A 632, A115 (2019) Fig. 12. Line ratio R 12 2−1/1−0 (a), R 12 3−2/1−0 (b), R 13 3−2/1−0 (c), and R 13/12 3−2 (d) distribution of the observed region. The overlaid contours are 13 CO (3-2) integrated intensities, which start from 0.9 K km s −1 (6σ) and increase by 0.9 K km s −1 . The beam sizes are shown in the lower left corner of panel a. In all panels, the (0, 0) offset corresponds to α J2000 = 22 h 18 m 00. s 44, δ J2000 = 61 • 49 03. 7. In panels c and d, the cross represents the peak position of region C2 that is analyzed in Fig. 13. temperatures lie in the range of 13-23 K with a median value of 18 K. These values are slightly higher than those (∼10 K) obtained in other nearby dark clouds (e.g., Benson & Myers 1989;Friesen et al. 2017;Keown et al. 2017), but they are much lower than in OMC-1, parts of which are likely heated by radiation from the Trapezium Cluster (Friesen et al. 2017;Tang et al. 2018). Although the pixels at the edges have higher kinetic temperature of >25 K, they have large uncertainties (see Fig. 14c) and are thus not included in our analysis. Higher kinetic temperatures are found within region C2 where star formation is more active than in the other regions. The high kinetic temperature could be due to feedback from YSOs. In contrast, H 2 number densities are better constrained to a range of 10 3 -10 3.6 cm −3 . The derived H 2 number densities are lower than the critical density (2 × 10 4 cm −3 ) of 13 CO (3-2) (Yang et al. 2010), demonstrating that this transition is subthermal throughout the entire observed region. Regions C1, C2, and C3 have similar H 2 number densities of ∼10 3.6 cm −3 . Given that star formation has already begun in regions C1 and C2, we suggest that region C3 is the next potential star-forming region because of similar physical conditions in the three regions. We also note that self-absorption in 12 CO (3-2) could affect the modeling results. In case of 12 CO (3-2) self-absorption, the observed R 13/12 3−2 is supposed to be higher than the modeled values without considering 12 CO (3-2) self-absorption. This will lead to an overestimate of the H 2 number density and an underestimate of kinetic temperature, indicating that region C3 tends to have H 2 number densities of <10 3.6 cm −3 and kinetic temperatures of 17 K.

Revisiting the cloud-cloud collision scenario for L1188
In this section, we assess the link between the newly observed features and the cloud-cloud collision scenario initially proposed by Gong et al. (2017). In Sect. 3.1, we not only confirm the presence of the two orthogonal molecular clouds, but also discover complementary distributions between L1188a and L1188b, which represent emission in distinct velocity ranges (see Fig. 6). Such complementary distributions are an expected outcome according to previous simulations of molecular cloud collisions (e.g., Habe & Ohta 1992;Anathpindika 2010;Takahira et al. 2014;Torii et al. 2017). Furthermore, L1188 has no massive stars inside, excluding the possibility that such a complementary The modeled line ratios, R 13 3−2/1−0 (red) and R 13/12 3−2 (blue), as a function of kinetic temperature and H 2 number density. The solid lines and the dotted lines represent the observed line ratios and their 1σ uncertainties for R 13 3−2/1−0 and R 13/12 3−2 . The shadowed region in the intersection represents χ 2 <2.3, corresponding to the 1σ range. (b) The modeled R 13/12 3−2 as a function of H 2 number density at a given kinetic temperature of 21 K. distribution is caused by the feedback of massive stars. Rather, it suggests a cloud-cloud collision in L1188 as its origin. Enhanced 13 CO emission found in the contact regions of L1188a and L1188b (see the right panel of Fig. 6) can result from a density enhancement caused by shocks expected during collisions. Loren (1976) points out that colliding molecular clouds could provide heating in the interface, which can lead to selfabsorption in optically thick lines due to temperature gradients. The observed widespread CO self-absorption (see Fig. 7) and high kinetic temperatures of >17 K (see discussions in Sect. 4.1) further support our proposed scenario.
The two parallel filamentary structures F1 and F2 are nearly perpendicular to the elongated axis of L1188b, and each of them has at least two velocity components. We therefore propose that the formation of F1 and F2 can also be attributed to the collision during which molecular gas is stripped off or dragged out. In such a case, there should be shear between different velocity components, which would lead to the Kelvin-Helmholtz (KH) instability. The KH instability could result in the periodic variation of intensities found in F1, which is similar to the scenario proposed to explain the Taurus striations (Heyer et al. 2016). For two incompressible fluids with densities ρ 1 and ρ 2 , and velocities, υ 1 and υ 2 , the maximum wavelength predicted by the KH instability (e.g., Chandrasekhar 1961) is where α 1 = ρ 1 /(ρ 1 + ρ 2 ), α 2 = ρ 2 /(ρ 1 + ρ 2 ), and g is the acceleration. The velocity difference is taken to be 2.5 km s −1 according to Figs. 4a-c. The density contrast (ρ 1 /ρ 2 ) of two fluids in the interface still has to be discussed. The integrated intensity contrast is found about three between the two velocity components in Fig. 4b. We assume the density contrast to be 3-10 at the initial contact stage. The acceleration is taken to be g = πGµm H N(H 2 ), where G is the gravitational constant, µ is the mean molecular weight per hydrogen molecule which is assumed to be 2.8, m H is the mass of the atomic hydrogen, and N(H 2 ) is the column density of the molecular hydrogen which is taken to be 2 × 10 21 cm −2 . This leads to a λ KH,max range of 7-24 pc. These maximum wavelengths are much higher than the observed periodic scale of 0.8 pc in F1 (see Fig. 4c). However, the observed regions differ from the ideal case (Chandrasekhar 1961). Hunter et al. (1997) incorporated thermodynamics into the KH instability by taking heating and CO cooling into account, and found that the perturbation wavelength is reduced to ∼0.9 pc for typical molecular regions (see their Sect. 3.4). The derived wavelength is similar to the observed periodic scale of 0.8 pc, indicating that such instability could occur in F1. The origin of the 1 pc-long shocked arc can also be attributed to the cloud-cloud collision. As shown in Sect. 3.2, we identified two velocity components connecting each other with bridging features across the whole arc. Such a velocity structure has been predicted by previous simulations of cloud-cloud collisions (Haworth et al. 2015a,b). Following the analyses for outflows (e.g., Arce et al. 2010), we estimate the momentum and kinetic energy of the blueshifted high-velocity wing emission to be 1.4 M km s −1 and 5.4 × 10 43 erg, respectively. According to previous theoretical predictions (e.g., Wu et al. 2018), such momentum and kinetic energy can be easily created by cloud-cloud collisions. Alternatively, the derived momentum and kinetic energy also lie in the momentum and kinetic energy range of outflows from low-mass YSOs (Arce et al. 2010). This indicates that even low-mass YSOs have the ability to inject the observed momentum and kinetic energy. However, the whole arc structure is nearly parallel to both the northeastern edge of L1188a and the centric concave shape of L1188b, and their linear scales are also similar. Such a morphology is more easily created by cloudcloud collisions. Nevertheless, we cannot exclude the possibility that the 1-pc shocked arc is driven by the outflow from YSOs.
The star-forming activities in C1-C3 support the sequential star formation from west to east reported by Szegedi-Elek et al. (2019). This can be alternatively explained by the cloudcloud collision scenario, in which L1188b moved from west to east. This scenario would readily result in the observed filamentary YSO distribution which is nearly parallel to the long axis of L1188b (see Fig. 1). Based on the proposed scenario, star formation is expected in the northwest, especially in the compressed region C3 that has similar physical conditions to the known star-forming regions C1 and C2 (see Sect. 4.1). Szegedi-Elek et al. (2019) suggest that star formation in L1188 started about 5 million years ago. Based on the assumption of the relative colliding speed of 5 km s −1 , Gong et al. (2017)   In panels c and d, the contours start from 0.9 K km s −1 (6σ) with a step of 0.9 K km s −1 . In all panels, the (0, 0) offset corresponds to α J2000 = 22 h 18 m 00. s 44, δ J2000 = 61 • 49 03. 7.
A115, page 14 of 16 Y. Gong et al.: Searching for further evidence for cloud-cloud collisions in L1188 growth time for KH instability (τ KH = λ KH,max υ 1 −υ 2 ρ 1 +ρ 2 (ρ 1 ρ 2 ) 1/2 , Drazin & Reid 1981) could be an independent estimate for the cloud-cloud collision timescale. By using the initial H 2 number densities of two clouds before the collision in Hunter et al. (1997), we reach a lower limit of 0.7 Myr for τ KH because the currently observed periodic length of 0.8 pc is likely shorter than λ KH,max . Hydrodynamical simulations have shown that a timescale of 2-5 Myr is needed to generate an arc structure (Takahira et al. 2014). The estimates point out that the collision might have happened about 0.7-5 Myr ago. This indicates that parts of YSOs in L1188 may have already formed before the collision, that is, parts of the star formation in L1188 are not caused by the collision. On the other hand, there are also stellar populations younger than 5 Myr. Therefore, the possibility of star formation triggered by the cloud-cloud collision should not be neglected.

Summary
In order to search for further observational evidence of cloudcloud collisions in L1188, we performed multiple molecular line observations toward the 20 × 20 intersection region of the two nearly orthogonal filamentary molecular clouds. The main results are summarized as follows.
1. Our 12 CO (2-1) and 12 CO/ 13 CO (3-2) observations provide a factor of approximately four better linear resolutions than previous CO observations. This leads to the discovery of a spatially complementary distribution between L1188a and L1188b. Enhanced 13 CO emission and 12 CO self-absorption is found toward their contacting regions. Furthermore, we found two filamentary structures F1 and F2, which are parallel to each other. Both of them have at least two velocity components that are connected with broad bridging features. At the most blueshifted side, we identified a 1 pc-long shocked arc ubiquitously showing 12 CO line wings. 2. In a 360 × 75 area of region C1 near the 1-pc shocked arc, we discover two positions showing 22 GHz water maser emission, which is the first maser detection in L1188. These masers are likely excited by the outflow from the Class I YSO WISEJ221729.87+615039.8. We also searched for thermal SiO emission and 95 GHz methanol maser emission, but none has been detected in the observed region. 3. We performed a non-LTE analysis of line ratios toward the observed region. This suggests that L1188 is characterized by kinetic temperatures of 13-23 K and H 2 number densities of 10 3 -10 3.6 cm −3 at a linear resolution of 0.2 pc. High kinetic temperatures of 20 K is found toward the most active starforming region, which is likely due to feedback from YSOs. 4. Compared with previous theoretical predictions and simulations, these new observational features can be readily explained by the scenario of cloud-cloud collisions. However, the feedback from low-mass young stellar objects may also make contributions to these observational features.  (Robitaille et al. 2006(Robitaille et al. , 2007. The solid black line indicates the best-fitting model. In order to study the nature of WISEJ221729.87+615039.8, we performed a spectral energy distribution (SED) fitting with the python SED fitter 15 (Robitaille et al. 2007) that uses radiative transfer models of YSOs in Robitaille et al. (2006). The modeled source distances are set to be within 700-1000 pc, while the modeled visual extinctions (A V ) range from 0 to 40. Figure A.1 presents the SED fitting result, which employed infrared data from 2MASS (Skrutskie et al. 2006), WISE (Wright et al. 2010), and AKARI (Kawada et al. 2007). This SED indicates that this YSO is a Class I object that is characterized by a 9.7 µm deep absorption of amorphous silicate in its circumstellar disk. From the best-fitting model, the distance is determined to be ∼788 pc and A V is 23.8, confirming that this YSO is deeply embedded in L1188. The luminosity and stellar mass are estimated to be 16 L and 1.2 M , suggesting a solar-type star. The inclination of the circumstellar disk is found to be about 57 • , but this can be underestimated because we only have an upper limit at 12 µm. Because of the potential outflow origin of the H 2 O masers (see Sect. 3.3), WISEJ221729.87+615039.8 may have a disk-jet system.