Spatial distribution of NH2D in massive star-forming regions

To understand the relation between NH$_2$D and its physical environment, we mapped ortho-NH$_2$D $1_{11}^s-1_{01}^a$ at 85.9 GHz toward 24 Galactic late-stage massive star-forming regions with Institut de Radioastronomie Millim$ \'e$trique (IRAM) 30-m telescope. Ortho-NH$_2$D $1_{11}^s-1_{01}^a$ was detected in 18 of 24 sources. Comparing with the distribution of H$^{13}$CN 1-0 as a dense gas tracer and radio recombination line H42$\alpha$, ortho-NH$_2$D $1_{11}^s-1_{01}^a$ present complex and diverse spatial distribution in these targets. 11 of the 18 targets, present a different distribution between ortho-NH$_2$D $1_{11}^s-1_{01}^a$ and H$^{13}$CN 1-0, while no significant difference between these two lines can be found in the other 7 sources, mainly due to limited spatial resolution and sensitivity. Moreover, with H42$\alpha$ tracing massive young stellar objects, ortho-NH$_2$D $1_{11}^s-1_{01}^a$ seems to show a relatively weak emission near the massive young stellar objects.


INTRODUCTION
Deuterium chemistry is an active process in molecular clouds (Caselli et al. 2008).The deuterated molecules, such as DCN, DCO + , N 2 D + and NH 2 D, can be formed at different evolution stages of star-forming (Caselli et al. 2008).Due to the relatively low zero-point energy, these deuterated molecules are generated more easily in molecular clouds (Herbst 1982).Therefore, the deuterium fractionation ( frac ), as the ratio between deuterated and normal molecule, is often higher than the D/H ratio of ≈ 1.5 × 10 −5 (Linsky et al. 1995;Oliveira et al. 2003).
As a single deuterated form molecule, NH 2 D was detected extensively in the molecular cloud. frac (NH 3 ) was observed to be ∼ 1 − 0.1 in the cold environment (e.g.Hatchell 2003;Crapsi et al. 2007;Harju et al. 2017) and ∼ 0.1 − 0.001 in the warm environment (e.g.Busquet et al. 2010;Pillai et al. 2011;Gong et al. 2015).Besides, with the model calculations, deuterium fractionation is sensitive to its environment changes, such as temperature and density (Roberts & Millar 2000a,b).
According to the high-resolution radio observation in late-stage massive star-forming regions, a series of the researches have reported the relation between NH 2 D and its physical environment.With the IRAM Plateau de Bure Interferometer (PdBI), Busquet et al. (2010) observed the ortho-NH 2 D 1  11 − 1  01 at 85.9263 GHz in the ultracompact HII region (UCHII region) IRAS 20293+3952, through which they reported that the  frac (NH 3 ) in young stellar objects (YSOs) was lower than in pre-protostellar cores, and the  frac (NH 3 ) could be related to the evolutionary stage in the pre-protostellar phase.With AtacamaLarge Millimeter/Submillimeter Array (ALMA) observation of Orion KL Hot Core and Compact Ridge for NH 2 D at 239.848 GHz, Neill et al. (2013) found that the deuterium fractionation in Orion KL Hot Core was lower than that in Compact Ridge, suggesting that the hot core originated in a slightly warmer environment result in low deuterium fractionation.Mills et al. (2018) reported that the  frac (NH 3 ) of Sgr B2(N2) is ≈ 10 times higher than that of Sgr B2(N1) with VLA observations of the NH 2 D 4  14 − 4  04 transitions at 25.0 GHz and ATCA observations of the NH 2 D 2  12 − 2  02 transitions at 49.9 GHz and the NH 2 D 3  13 − 3  03 transitions at 43.0 GHz, so that they attributed the abundance variations to the shock and fast warmup in Sgr B2(N2).However, the researches mentioned above are only based on one or two sources (Busquet et al. 2010;Neill et al. 2013;Mills et al. 2018).Consequently, a relatively large sample of massive star-forming regions is needed to understand the relation between NH 2 D and its physical environment.NH 2 D 1  11 − 1  01 at 110.1535 GHz with single point observations using Institut de Radioastronomie Millimétrique (IRAM) 30-m telescope toward 50 late-stage massive star-forming regions was reported with 72% detection rate (Li et al. 2022).With stronger emission (e.g.Shah & Wootten 2001;Wienen et al. 2021), ortho-NH 2 D 1  11 −1  01 at 85.9263 GHz is more suitable than para-NH 2 D 1  11 −1  01 at 110.1535 GHz to derive NH 2 D spatial distribution.Therefore, we mapped the ortho-NH 2 D 1  11 − 1  01 at 85.9263 GHz in late-stage high mass star-forming regions with IRAM 30-m to obtain spatial distribution information of NH 2 D.
With IRAM 30-m telescope observation, this work presents the results of ortho-NH 2 D 1  11 − 1  01 at 85.9263 GHz mapping toward a relatively large sample of 18 Galactic late-stage massive star-forming regions with known 6.7 GHz CH 3 OH masers.The observations are described in Section 2, the main results are reported in Section 3, a discussion is presented in Section 4, and a brief summary is given in Section 5.

OBSERVATIONS
The 24 targets, selected from Reid et al. (2014), are late-stage massive star-forming regions with 6.7 GHz CH 3 OH masers with HC 3 N 12-11 higher than 1 K ( mb ) which is based on IRAM 30-m observations.The observations were performed with the Institut de Radioastronomie Millimétrique (IRAM) 30-m telescope at Pico Veleta, Spain during 2019 July, 2019 October, 2019 November, 2020 December and 2021 January.The on-the-fly (OTF) mode was used in our observations with the frequency coverage about 8 GHz bandwidth.The system temperatures were 80-150 K during our observation.We used the E090 receiver to cover 84 GHz to 92 GHz and the Fourier Transform Spectrometers (FTS) backend with a 195 kHz channel spacing, which corresponds to 0.68 km s −1 at 86 GHz.The map centers and map area are listed in table 1.
The antenna temperature ( * A ) is converted to the main beam brightness temperature ( mb ), using  mb =  * A •  eff / eff , the forward efficiency  eff is 0.95 and the beam efficiency  eff is 0.81.The used lines consist of ortho-NH 2 D 1  11 − 1  01 at 85.9263 GHz, H 13 CN 1-0 at 86.3387 GHz and H42 at 85.6884 GHz.Data reduction was conducted with GILDAS software 1 .The angular resolution of the IRAM 30-m telescope is about 29.3 arcsec at 84 GHz and about 26.7 arcsec at 92 GHz, while all sources were re-gridded to steps of 9 arcsecs, about 1/3 beam size.First-order baseline was used for all spectra.The typical rms are around 50 mK at 0.7 km s −1 velocity resolution.In the following analysis, we assume that the beam-filling factor is unity.
The detailed results of the three lines are given in Section 3.1-3.4,the ortho-NH 2 D 1  11 − 1  01 , H 13 CN 1-0 and H42 distribution details of each source are described in Section 3.5.

NH 2 D spatial distribution
The velocity integrated intensity of ortho-NH 2 D 1  11 −1  01 at 85.9263 GHz, observed with the IRAM-30m, are presented in Figure 1 -18 as red contour.The spectra of ortho-NH 2 D 1  11 − 1  01 at 85.9263 GHz are also presented in Figure 1 -18.For overlaying with H 13 CN 1-0 spectra, there are some of the sources in which several offsets are shown in Figure 1 -18.These several offsets are used to confirm the signal of ortho-NH 2 D 1  11 − 1  01 emission if this target has spatial structures, such as G081.75+00.59(see Figure 13).
Ortho-NH 2 D 1  11 − 1  01 has six hyperfine structures, with relative intensities listed in Table 2 (Tiné et al. 2000).However, the transitions of F=2-1, F=2-2, F=1-1 and F=1-2 for ortho-NH 2 D 1  11 − 1  01 are blended due to line broadening.For example, as shown in Figure 19, only three components can be found in G081.75+00.59for NH 2 D spectra.Thus, the fluxes of ortho-NH 2 D 1  11 − 1  01 were derived with the velocity integrated intensity of the emitting channels of its six hyperfine structures instead of Gaussian fitting.In the following analysis, we assume all NH 2 D spectra are optically thin.
The strongest velocity integrated intensity and distribution shape of ortho-NH 2 D for each source are listed in Table 4.The strongest ortho-NH 2 D emission is 6.95±0.14K km/s in G081.75+00.59and the weakest ortho-NH 2 D emission is 0.60±0.16and 0.60±0.11K km/s in G005.88-00.39 and G037.43+01.51,respectively.The median of ortho-NH 2 D emission is 2.27 K km/s.For the sources whose ortho-NH 2 D emission is lower than 1 K km/s, are considered as sources with weak NH 2 D emission.Conversely, we define strong NH 2 D emissions to be cases surpassing the 1 K km/s threshold.Accordingly, there are 4 sources characterized by weak NH 2 D emissions and 14 sources are classified into strong NH 2 D emissions.
With strong NH 2 D emission detected, 14 sources show complex and diverse distribution, which is different from other molecule tran-sitions, such as the dense gas tracer H 13 CN 1-0 (details described in the following).
For two hydrogen atoms in NH 2 D, different nuclear spin angular momentum can produce different states: ortho and para state.The ortho state has the largest nuclear spin statistical weight, and the para state has the largest nuclear spin statistical weight (Harju et al. 2017).Since normal radiative and collisional transitions can not change the spin of H nuclei, the transitions between ortho-and para-NH 2 D are forbidden (Ho & Townes 1983).Therefore, obtaining ortho-and para-NH 2 D could reflect an accurate abundance of NH 2 D.
The corresponding nuclear statistical weights are 27 and 9 for ortho-NH 2 D 1 The highest ratio is 4.73±1.99 in G035.20-01.73 and the lowest ratio is 0.89±0.44 in G005.88-00.39,both of them are with strong H42 emissions.In other sources, this ratio ranges from 1.26±0.07 to 2.95±0.12.Besides, the median of this ratio is 1.84.
2 https://cdms.astro.uni-koeln.de/classic/predictions/catalog/ 3.3 H 13 CN 1-0 maps result HCN lines can be used to trace dense gas (≳ 10 4 cm −3 ) because of high dipole moment and large Einstein A-coefficient (Evans et al. 2006).With the similar dipole moment and Einstein A-coefficient, lines of H 13 CN, as an isotope molecule of HCN, can be used as a dense gas tracer.So, H 13 CN lines can be more accurate to trace dense molecular gas, since they are mainly optically thin.Therefore, H 13 CN 1-0 is used as a dense gas tracer (Gao & Solomon 2004;Wu et al. 2010;Vasyunina et al. 2011).
In this work, the H 13 CN 1-0 emission was detected in all targets.H 13 CN 1-0 has three hyperfine structures.Similar to ortho-NH 2 D 1  11 − 1  01 , the velocity integrated intensity of H 13 CN 1-0 was derived from its three hyperfine structures with CLASS/GILDAS package.The mapping results and spectra of H 13 CN 1-0 are shown in Figure 1 -18 (black contour and gray scale).

H42𝛼 maps result
The young massive stars, with spectral types of B or O, can produce ultra-violet photons, which can ionize the surrounding interstellar medium and create HII regions.Consequently, the abundantly ionized hydrogens exist in HII regions and can recombine with electrons to produce the radio recombination lines (RRLs).Therefore, the RRLs can be used to trace young massive star-forming and HII regions.RRL H42 was detected in nine sources (see also in figure 1 -18, the black contour and gray scale).Without turbulence, the full width at half maxima (FWHM) of RRLs is 19.1 km s −1 at 8000 K (Peters, Longmore, & Dullemond 2012).Therefore, in the case of a faint RRLs signal observed in G031.28+00.06, the spectral resolution was smoothed to 1.4 kms −1 to confirm whether it is H42 emission or not in the present work.Among these nine targets, a completely different spatial distribution between ortho-NH 2 D and H42 is presented in five targets, while the spatial distribution of ortho-NH 2 D and H42 have overlapping regions in the rest four targets.
A strong H42 emission presents in G005.88-00.39(see Figure 1), while the ortho-NH 2 D were marginally detected and almost unresolved in the same position.There are two NH 2 D emission features and a relatively weak H42 structure with associated NH 2 D emission in G011.91-00.61(see Figure 2).G049.48-00.36 and G049.48-00.38 are two adjacent sources (see Figure 10 and Figure 11).Both of them, the ortho-NH 2 D and H42 present different spatial distributions.However, the ortho-NH 2 D emission peaks at the H42 emission peak in G049., the ortho-NH 2 D is weaker than that in G049.48-00.38.

Individual target detail
In our observation, most sources show different spatial distributions between ortho-NH 2 D and H 13 CN 1-0 emission, or between ortho-NH 2 D and H42 emission, while a few sources show similar spatial distributions.For example, the ortho-NH 2 D and H 13 CN 1-0 emission present similar spatial distribution in G011.91-00.61.The details each source are described in following.

G005.88-00.39
G005.88-00.39 is a bright expanding shell-like UCHII region (Acord, Churchwell, & Wood 1998).According to near-infrared NACO-VLT observations, Feldt et al. (2003) suggested G005.88-00.39contains a massive young O5 V star.Furthermore, this target is one of the most powerful outflows in Galaxy (Harvey & Forveille 1988;Zapata et al. 2020).In the present work, strong ( ∫  mb d ≥40 K km s −1 ) H 13 CN 1-0 and H42 emission were detected in G005.88-00.39,and their velocity integrated intensity maps(see figure 1).However, ortho-NH 2 D at 85.9 GHz was only observed at the source center with a weak emission, which is the weakest among our samples.Therefore, the ortho-NH 2 D may be marginally detected in this source and it was plotted from 2.Moreover, the ortho-NH 2 D velocity integrated intensity mapping presents less than a beam, so the spatial structure might not be distinguished in G005.88-00.39.Due to the weak ortho-NH 2 D emission and indistinguishable spatial structure, we can not confirm whether the distribution of ortho-NH 2 D is similar to other molecules or not.
Although the intensity of H 13 CN 1-0 is similar in these two structures, the ortho-NH 2 D appears weaker in the northeast structure where H42 is detected.

G023.44-00.18
The strong 6.7 GHz CH 3 OH maser has been detected in the massive star-forming region G023.44-00.18(Walsh et al. 1998).With Submillimeter Array (SMA) observation, Ren et al. (2011) reported that G023.44-00.18presented two dust cores that were north and south arranged, the south one has active bipolar outflow.In addition, they suggested that these two dust cores were in the pre-UCHII evolutionary stage.By observing the ortho-NH 2 D 1  11 − 1  01 at 85.9263 GHz toward PdBI, Zhang et al. (2020) found several NH 2 D cores in G023.44-00.18with D frac ∼0.1.But these cores are offset by two dust cores which are reported in Ren et al. (2011).However, because of the spatial resolution limitation, both ortho-NH 2 D and H 13 CN 1-0 emission can not be distinguished by two cores in our observation, and H42 was not detected (see Figure 4).Nevertheless, in the northwest of G023.44-00.18, the obvious NH 2 D signal was detected with weak H 13 CN 1-0 emission.

G031.28+00.06
G031.28+00.06 is located in W43, which is a molecular cloud complex (Nguyen Luong et al. 2011).As is shown in figure 5, ortho-NH 2 D, H 13 CN 1-0 and H42 were detected in G031.28+00.06.The ortho-NH 2 D emission presents two structures in this source, while H 13 CN 1-0 and H42 (see figure 5) present one structure.Compared with other targets, the ortho-NH 2 D and H42 emission in G031.28+00.06 are relatively weak.Furthermore, two NH 2 D structures seem to be isolated by the H42 structure.Besides, the ortho-NH 2 D and H 13 CN present a different distribution in this source.

G035.20-01.73
G035.20-01.73(W 48) contains an UCHII region (Wood & Churchwell 1989).The east to west evolutionary gradient was found in the HII region (Rygl et al. 2014).We obtained the H 13 CN 1-0, distributed from east to west with two strong peaks and a weak peak (see figure 8).Pillai et al. (2011) present the distribution of ortho-NH 2 D at 85.9 GHz oriented west to east in G035.20-01.73toward PdBI observations.Combined with the 3.5 mm dust emission map, they found the dust emission peak is not the brightest NH 2 D peak in this target.H42 is observed in G035.20-01.73,which is near from the center peak of H 13 CN 1-0 distribution (see figure 8).The ortho-NH 2 D was also observed with strong emission (see figure 8).The ortho-NH 2 D, away from the H42, is located in the west of G035.20-01.73 and shows a similar distribution to the result in Pillai et al. (2011).Besides, the NH 2 D peak offsets H 13 CN 1-0 peak, while the NH 2 D distribution is similar to the 3.5 mm dust emission in Pillai et al. (2011).
3.5.10G049.48-00.36G049.48-00.36(W 51 IRS2) is one of the massive protocluster candidates in giant molecular cloud W51 (Ginsburg et al. 2012).W 51 is a complex structure molecular cloud with several velocity components (Scoville & Solomon 1973).In the present work, the ortho-NH 2 D and H 13 CN 1-0, H42 (see figure 10) maps were obtained in G049.48-00.36.At the map center, the structures of the ortho-NH 2 D, H 13 CN 1-0 and H42 distribution are close.Although part of ortho-NH 2 D emission is detected near one of the H42 structures, most ortho-NH 2 D emission distributions deviate from the H42 emission.With strong emission, the ortho-NH 2 D, H 13 CN 1-0 and H42 were also detected in the southeast where G049.48-00.38 is located..5.11 G049.48-00.38G049.48-00.38 (W 51 M) is adjacent to G049.48-00.36 and it is also one of the massive protocluster candidates in W 51 (Ginsburg et al. 2012).As figure 11 shows, strong ortho-NH 2 D and H 13 CN 1-0, H42 were observed in this work (see figure 10).The ortho-NH 2 D distribution is similar to H 13 CN 1-0, while the ortho-NH 2 D and H 13 CN 1-0 emission offset H42 peak.With the PdBI observations, Vastel et al. (2017) reported that ortho-NH 2 D was detected in two marginally resolved sources, while we have detected only one here due to the spatial resolution limitation.
3.5.13G081.75+00.59G081.75+00.59(DR 21) is a giant star-forming complex in Cygnus X molecular cloud complex (Dickel, Dickel, & Wilson 1978;Kirby 2009).It had been mapped to a 4 pc length ridge by Herschel observations (Motte et al. 2007;Hennemann et al. 2012) and contains three OH 1665 MHz emissions regions DR 21M, DR 21OH and W 75S (Koley et al. 2021).With IRAM observation, Christensen et al. (2023) reported integrated column density map of NH 2 D, DCN, DNC and DCO + toward Cygnus-X.In our previous work, this target has been detected the strongest velocity integrated intensity of para-NH 2 D 1  11 − 1  01 at 110.1535 GHz with IRAM-30m and the highest deuterated fraction of NH 3 (D frac =0.044) (Li et al. 2022).In the present work, we mapped DR 21 (OH) and W 75S in DR 21 and we describe them separately in the following two parts.
DR 21OH contains several high mass protostars (Mangum, Wootten, & Mundy 1991).It is younger than DR 21M because it was undetected in the radio continuum emission (Davis et al. 2007).In the present work, we only obtained the north part in DR 21OH.The ortho-NH 2 D and H 13 CN 1-0 emissions are strong in this region (see figure 13).The H 13 CN 1-0 presents a strong peak in the southernmost toward our map and ortho-NH 2 D offset the H 13 CN peak.

G109.87+02.11
G109.87+02.11(Cep A) is an active massive star-forming region, and it has been detected in several YSOs (Hughes & Wouterloot 1984).The brightest radio source in G109.87+02.11 is HW2 (where locates at (0 ′′ ,0 ′′ ) in our maps), associated with an early-B spectral type zeroage main-sequence (ZAMS) star (Hughes & Wouterloot 1984).With IRAM 30-m observations, Li et al. (2017) observed strong ortho-NH 2 D emission in the northeast of G109.87+02.11.As is shown in figure 15, the H 13 CN 1-0 and strong ortho-NH 2 D emission were detected in G109.87+02.11.The distribution of ortho-NH 2 D is presented in the northeast of G109.87+02.11,and it is similar to the result of Li et al. (2017).However, the ortho-NH 2 D is dislocation with the strongest H 13 CN 1-0 emission.The strongest H 13 CN 1-0 emission is situated in proximity to HW2, and the ortho-NH 2 D emissions appearing to be isolated by HW2.Besides, in this target, the strong ortho-NH 2 D was detected in the northeast of G109.87+02.11and H 13 CN was detected relatively weakly.

G111.54+00.77
G111.54+00.77(NGC 7538) is a massive star-forming region which locates in a giant molecular cloud complex (Ungerechts, Umbanhowar, & Thaddeus 2000).With Herschel observation, Fallscheer et al. (2013) reported that G111.54+00.77has several high-mass dense clumps candidates which may produce several intermediate-to highmass stars in the future.A young O star with 25 M ⊙ embedded in NGC 7538 IRS 1 has been observed with a widely distributed molecular outflow (Qiu, Zhang, & Menten 2011).In G111.54+00.77, the widely distributed ortho-NH 2 D and H 13 CN 1-0 distribution were detected (see figure 16), and H42 emission was detected in the north where is NGC 7538 IRS 1.A small part of ortho-NH 2 D emission is near to the H42 emission while a large part appears in the south of G111.54+00.77.In the west of G111.54+00.77, the ortho-NH 2 D emission is strong while H 13 CN 1-0 emission is weak.

G121.29+00.65
Toward SMA observations, Juárez et al. (2019) reported that G121.29+00.65 (L 1287) have a high fragmentation level and a bipolar outflow is oriented northeast to southwest direction.G121.29+00.65 have been detected 6.7 GHz CH 3 OH, implying that there is forming a massive star (Rygl et al. 2010).In G121.29+00.65, the distribution of H 13 CN 1-0 was obtained (see figure 17), similar to the distribution of HCN 1-0 which is reported by Yang et al. (1991).The distribution of ortho-NH 2 D is perpendicular to H 13 CN 1-0 and the widely distributed ortho-NH 2 D is oriented northwest to the southeast.The distribution of ortho-NH 2 D is also perpendicular to the outflow reported in Juárez et al. (2019).

G188.94+00.88
It was reported that G188.94+00.88(S 252) includes two massive protostars and is near to the Sh 2-247 which is a complex of HII region and OB stars (Minier et al. 2005).Several very weak ortho-NH 2 D fragments were detected in G188.94+00.88(see figure 18).These ortho-NH 2 D emissions are marginally detected in this target.Similar to G005.88-00.39, the spatial structure might also not be distinguished in G188.94+00.88due to the beam size.The main ortho-NH 2 D emission region is presented in the south of H 13 CN 1-0 distribution.In G188.94+00.88,H42 was not detected.
Without asymmetrically distributed spatial structures detected and unresolved sources, only G081.75+00.59presents a relatively similar distribution between NH 2 D and H 13 CN 1-0.There may exist contingencies in this distribution.
As described in Section 3.3, H 13 CN 1-0 is used as a dense gas tracer (Gao & Solomon 2004;Wu et al. 2010;Vasyunina et al. 2011) in the present work to analyze the distribution of NH 2 D and dense gas in massive star-forming late stage.Crapsi et al. (2007) performed observations of para-NH 2 D 1  11 − 1  01 at 110.1535 GHz toward low-mass starless core L 1544 with PdBI and they found that NH 2 D intensity peaks were at 1.2 mm continuum dust peak.Therefore, they suggested that NH 2 D is a good tracer of the high density gas in the low-mass starless core.However, toward the gas-phase model, Roberts & Millar (2000b) simulated that gas density seems to have no effect on D frac (NH 3 ).
In present work, NH 2 D seems to be enhanced in G031.28+00.06,G035.19-00.74,G081.87+00.78,G109.87+02.11and G121.29+00.65,at relatively low gas density part (see Figure 5, Figure 7, Figure 14, Figure 15 and Figure 17, respectively).Therefore, such results indicate that the gas density may not be an important physical parameter to influence the abundance of NH 2 D, which is opposite to that of Crapsi et al. (2007) and agrees with the gas-phase model of Roberts & Millar (2000b).Compared with other physical parameters (such as temperature), the gas density might not be the important physical parameter to affect the abundance of NH 2 D in the late-stage massive star-forming regions.But without any information of NH 3 , we can not confirm that the D frac (NH 3 ) is affected by dense gas.More detail should be obtained with the information of NH 3 .

Spatial distribution of ortho-NH 2 D and massive star formation
To analyze whether NH 2 D is affected by the formed massive young stellar objects or not, we compare the distribution between ortho-NH 2 D and H42.In the present work, ortho-NH 2 D and HII regions present a different distribution in most targets.Both ortho-NH 2 D and H42 emissions were detected in nine targets.In five sources, G015.03-00.67,G031.28+00.06,G035.20-01.73,G075.76+00.33 and G111.54+00.77, the NH 2 D emission are away from HII regions (see Figure 3, Figure 5, Figure 8, Figure 12 and Figure 16, respectively).In these five HII regions, the NH 2 D may be depleted during the massive stars forming.For G049.48-00.36 and G049.48-00.38, the relatively strong ortho-NH 2 D emission was detected in HII regions (see Figure 10 and Figure 11).The structure in W 51 is complex (Scoville & Solomon 1973), we can not confirm if that the ortho-NH 2 D and H42 emission come from the same structure.But it is obvious that ortho-NH 2 D emission is away from H42 emission, especially the strong H42 emission and the center of the HII region.
The ortho-NH 2 D presents two structures in G011.91-00.61(see Figure 2).One is the northeast structure which is similar to the distribution of H42 emission; another is the southwest structure without H42 emission.However, the ortho-NH 2 D emissions in the southwest structure are stronger and more extensive than emissions in the northeast structure.The reason may be that the massive star-forming stage in the northeast part is later than the southwest structure, and as a result, NH 2 D in the northeast part is much more depleted than NH 2 D in the southwest part.
In cold molecular gas, there is an active process of deuterated molecules forming (Caselli & Ceccarelli 2012).And during this process, several neutral species are frozen-out onto dust grains and molecules are deuterated on the surface of dust grains (Caselli & Ceccarelli 2012).Moreover, according to the simulation, Caselli et al. (2008) suggested that the deuterium fractionation of H + 3 (which means the abundance ratio of H + 3 to H 2 D + ) dropped rapidly when the gas kinetic temperature is above ∼ 15 K.The formation process of NH 2 D is related to H 2 D + (Emprechtinger et al. 2009).It may indicate that the  frac (NH 3 ) decreases with temperature increasing when the temperature is relatively high.Besides, with PdBI observation toward the massive star-forming region IRAS 20293+3952, Busquet et al. (2010) reported that they found strong NH 2 D emission in starless cores and did not detect NH 2 D near YSOs.Therefore, according to previous research of observation and theory (e.g.Caselli et al. 2008;Emprechtinger et al. 2009;Busquet et al. 2010;Caselli & Ceccarelli 2012), the abundance of NH 2 D should decrease with the evolution of massive star formation.
Without H42 detected, these nine sources present both strong and weak NH 2 D emission.In two sources, G037.43+01.51 and G188.94+00.88, the velocity integrated intensity of ortho-NH 2 D is lower than 1 K km/s (see Figure 9 and Figure 18).In G037.43+01.51, the NH 3 column density has been derived with VLA observation (Lu et al. 2014).Compared with other sources, this target presents low column density, such as G034.39+00.22whose column density is one order of magnitude higher than G037.43+01.51(Lu et al. 2014).Similarly, the NH 3 column density of G188.94+00.88 is also lower than other sources, such as G081.75+00.59whose column density is one order of magnitude higher than G188.94+00.88 with Effeslsberg observation (Jĳina, Myers, & Adams 1999).With low NH 3 column density, the NH 2 D might generate relatively less in these two sources.
The NH 2 D seems to be depleted near the HII regions which indicates a later evolution stage.Furthermore, in the region without H42 emission, the strong ortho-NH 2 D emission is detected.Therefore, combining with previous researches (e.g.Caselli et al. 2008;Emprechtinger et al. 2009;Busquet et al. 2010;Caselli & Ceccarelli 2012), we suggest that the abundance of NH 2 D may be affected by the evolution of massive star formation.

NH 2 D enhancement
With the Green Bank Telescope and the K-Band Focal Plane Array observation, Urquhart et al. (2015) have mapped the NH 3 (1,1) and NH 3 (2,2) inversion emission toward 66 massive star-forming regions, identified by the Red Midcourse Space Experiment Source survey (Lumsden et al. 2013).These 66 massive star-forming regions include several of our targets, such as G035.19-00.74(see figure 20).In G035.19-00.74,strong NH 2 D emissions present a northwest to southeast distribution, while the northwest and the southeast regions present weak H 13 CN 1-0 emission (see figure 7), these regions might be devoid of dense molecular gases.However, in the weak H 13 CN 1-0 region, NH 3 (1,1) and NH 3 (2,2) were detected.Therefore, the molecular gases that exist in these regions and their density is relatively low.
The obvious NH 2 D enrichment is presented in the southeast and northwest of G035.19-00.74.In particular, in the northwest of G035.19-00.74, the very high deuterium fractionation of NH 3 is detected.Compared with the H 13 CN emission in G035.19-00.74, the deuterium fractionation of NH 3 is either not affected by gas density.Moreover, combining with the ammonia data, we could obtain the distribution of the deuterium fractionation in the following work.

SUMMARY AND CONCLUSION REMARKS
Using of the IRAM-30m, we carried the ortho-NH 2 D 1  11 − 1  01 at 85.9263 GHz mapping toward a sample of 24 late stage massive starforming.The ortho-NH 2 D 1  11 − 1  01 at 85.9263 GHz was detected in 18 targets.H 13 CN 1-0 at 86.3387 GHz was also detected in these 18 targets, while H42 at 85.6884 GHz was detected in nine targets.
11 of 18 sources, present a different distribution between ortho-NH 2 D 1  11 − 1  01 and H 13 CN 1-0 with asymmetrically and resolvable distributed spatial structure.For the other 7 sources, no significant difference between these two lines can be found.This is mainly due to the limited spatial resolution and sensitivity.
Our main results include: 1.The ortho-NH 2 D 1  11 − 1  01 emissions present a complex distribution, which is different from other molecule transitions, such as H 13 CN 1-0. 2. Most of the targets show the different distribution between ortho-NH 2 D 1  11 − 1  01 and H 13 CN 1-0.Compared with other physical parameters (such as temperature), therefore, the dense gas (or gas density) may not be an important physical parameter to affect the NH 2 D enhancement in these targets.

Figure 1 .Figure 2 .
Figure1.(a): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (gray scale and black contour) in G005.88-00.39.The contour levels start at 2 in steps of 1 for NH 2 D, while the contour levels start at 18 in steps of 36 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour) in G005.88-00.39.The contour levels start at 2 in steps of 1 for NH 2 D, while the contour levels start at 7 in steps of 24 for H42.The gray scale starts at 3. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G005.88-00.39.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (0 ′′ ,0 ′′ ).(d): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G005.88-00.39.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (0 ′′ ,0 ′′ ).

Figure 3 .
Figure 3. (a): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (Gray scale and black contour) in G015.03-00.67.The contour levels start at 3 in steps of 2 for NH 2 D, while the contour levels start at 18 in steps of 12 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour) in G015.03-00.67.The contour levels start at 3 in steps of 2 for NH 2 D, while the contour levels start at 6 in steps of 15 for H42.The gray scale starts at 3. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G015.03-00.67.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (45 ′′ ,36 ′′ ).The black box in the picture indicates the range of flux integration of ortho-NH 2 D at 85.9263 GHz.(d): Spectra of ortho-NH 2 D at 85.9263GHz and H 13 CN 1-0 in G015.03-00.67.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (-9 ′′ ,-18 ′′ ).(e): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G015.03-00.67.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (-54 ′′ ,-108 ′′ ).

Figure 4 .Figure 5 .
Figure 4. (a): NH 2 D at 85.9263  GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (Gray scale and black contour) in G023.44-00.18.The contour levels start at 5 in steps of 3 for NH 2 D, while the contour levels start at 12 in steps of 12 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour, while the H42 is not detected) in G023.44-00.18.The contour levels start at 5 in steps of 3 for NH 2 D. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G023.44-00.18.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (0 ′′ ,0 ′′ ).The black box in the picture indicates the range of flux integration of ortho-NH 2 D at 85.9263 GHz.(d): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G023.44-00.18.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (0 ′′ ,-9 ′′ ).

Figure 6 .
Figure6.(a): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (Gray scale and black contour) in G034.39+00.22.The contour levels start at 5 in steps of 3 for NH 2 D, while the contour levels start at 9 in steps of 9 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour, while the H42 is not detected) in G034.39+00.22.The contour levels start at 5 in steps of 3 for NH 2 D. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G034.39+00.22.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (0 ′′ ,-45 ′′ ).The black box in the picture indicates the range of flux integration of ortho-NH 2 D at 85.9263 GHz.(d): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G034.39+00.22.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (-18 ′′ ,-27 ′′ ).(e): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G034.39+00.22.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (-9 ′′ ,-9 ′′ ).

Figure 7 .
Figure 7. (a): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (Gray scale and black contour) in G035.19-00.74.The contour levels start at 5 in steps of 5 for NH 2 D, while the contour levels start at 21 in steps of 12 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour, while the H42 is not detected) in G035.19-00.74.The contour levels start at 5 in steps of 5 for NH 2 D. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G035.19-00.74.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (0 ′′ ,0 ′′ ).The black box in the picture indicates the range of flux integration of ortho-NH 2 D at 85.9263 GHz.(d): Spectra of ortho-NH 2 D at 85.9263GHz and H 13 CN 1-0 in G035.19-00.74.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (0 ′′ ,0 ′′ ).(e): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G035.19-00.74.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (-45 ′′ ,45 ′′ ).(f): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G035.19-00.74.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (45 ′′ ,-45 ′′ ).

Figure 8 .Figure 9 .Figure 10 .Figure 11 .
Figure 8. (a): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (Gray scale and black contour) in G035.20-01.73.The contour levels start at 5 in steps of 5 for NH 2 D, while the contour levels start at 9 in steps of 12 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour) in G035.20-01.73.The contour levels start at 5 in steps of 5 for NH 2 D, while the contour levels start at 12 in steps of 12 for H42.The gray scale starts at 3. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G035.20-01.73.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (0 ′′ ,0 ′′ ).The black box in the picture indicates the range of flux integration of ortho-NH 2 D at 85.9263 GHz.(d): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G035.20-01.73.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (-54 ′′ ,0 ′′ ).

Figure 12 .Figure 13 .
Figure12.(a): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (Gray scale and black contour) in G075.76+00.33.The contour levels start at 3 in steps of 2 for NH 2 D, while the contour levels start at 15 in steps of 12 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour) in G075.76+00.33.The contour levels start at 3 in steps of 2 for NH 2 D, while the contour levels start at 5 in steps of 5 for H42.The gray scale starts at 3. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G075.76+00.33.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (0 ′′ ,0 ′′ ).The black box in the picture indicates the range of flux integration of ortho-NH 2 D at 85.9263 GHz.(d): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G075.76+00.33.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (-18 ′′ ,-27 ′′ ).

Figure 14 .
Figure14.(a): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (Gray scale and black contour) in G081.87+00.78.The contour levels start at 5 in steps of 2 for NH 2 D, while the contour levels start at 30 in steps of 15 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour, while the H42 is not detected) in G081.87+00.78.The contour levels start at 5 in steps of 2 for NH 2 D. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G081.87+00.78.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (0 ′′ ,0 ′′ ).The black box in the picture indicates the range of flux integration of ortho-NH 2 D at 85.9263 GHz.(d): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G081.87+00.78.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (0 ′′ ,0 ′′ ).(e): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G081.87+00.78.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (+36 ′′ ,0 ′′ ).

Figure 15 .Figure 16 .
Figure15.(a): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (Gray scale and black contour) in G109.87+02.11.The contour levels start at 3 in steps of 3 for NH 2 D, while the contour levels start at 12 in steps of 7 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour, while the H42 is not detected) in G109.87+02.11.The contour levels start at 3 in steps of 3 for NH 2 D. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G109.87+02.11.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (0 ′′ ,0 ′′ ).The black box in the picture indicates the range of flux integration of ortho-NH 2 D at 85.9263 GHz.(d): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G109.87+02.11.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (+51 ′′ ,+18 ′′ ).

Figure 17 .
Figure 17.(a): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (Gray scale and black contour) in G121.29+00.65.The contour levels start at 3 in steps of 2 for NH 2 D, while the contour levels start at 5 in steps of 6 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour, while the H42 is not detected) in G121.29+00.65.The contour levels start at 3 in steps of 1 for NH 2 D. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G121.29+00.65.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (0 ′′ ,0 ′′ ).The black box in the picture indicates the range of flux integration of ortho-NH 2 D at 85.9263 GHz.(d): Spectra of ortho-NH 2 D at 85.9263GHz and H 13 CN 1-0 in G121.29+00.65.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (0 ′′ ,0 ′′ ).(e): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G121.29+00.65.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (-27 ′′ ,9 ′′ ).(f): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G121.29+00.65.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (18 ′′ ,-18 ′′ ).

Figure 18 .
Figure18.(a): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H 13 CN 1-0 velocity integrated intensity image (Gray scale and black contour) in G188.94+00.88.The contour levels start at 3 in steps of 1 for NH 2 D, while the contour levels start at 5 in steps of 2 for H 13 CN 1-0.The gray scale starts at 3. (b): NH 2 D at 85.9263 GHz velocity integrated intensity contour (red contour) overlaid on H42 velocity integrated intensity image (Gray scale and black contour, while the H42 is not detected) in G188.94+00.88.The contour levels start at 3 in steps of 1 for NH 2 D. (c): Spectra of ortho-NH 2 D at 85.9263 GHz and para-NH 2 D at 110.1535 GHz in G188.94+00.88.The para-NH 2 D was observed by IRAM-30m with position-switching mode.The red spectra is ortho-NH 2 D and the black is para-NH 2 D. The offset for two spectra is (0 ′′ ,0 ′′ ).The black box in the picture indicates the range of flux integration of ortho-NH 2 D at 85.9263 GHz.(d): Spectra of ortho-NH 2 D at 85.9263 GHz and H 13 CN 1-0 in G188.94+00.88.The red spectra is ortho-NH 2 D and the black is H 13 CN 1-0.The offset for two spectra is (0 ′′ ,-18 ′′ ).

Figure 20 .
Figure 20.(a): NH 2 D at 85.926 GHz velocity integrated intensity contour (red contour) overlaid on NH 3 (1,1) main group velocity integrated intensity image (Gray scale and black contour) in G035.19-00.74.The contour levels start at 5 in steps of 5 for NH 2 D, while the contour levels start at 90 in steps of 54 for NH 3 (1,1) main group.The gray scale starts at 3. (b): NH 2 D at 85.926 GHz velocity integrated intensity contour (red contour) overlaid on NH 3 (2,2) main group velocity integrated intensity image (Gray scale and black contour) in G035.19-00.74.The contour levels start at 5 in steps of 5 for NH 2 D, while the contour levels start at 54 in steps of 27 for NH 3 (2,2) main group.The gray scale starts at 3.

Table 1 :
Source List