# Strong coupling between organic dye molecules and lattice modes of a
dielectric nanoparticle array

Rebecca Heilmann  Department of Applied Physics, Aalto University School of
Science, P.O. Box 15100, Aalto, FI-00076, Finland  Aaro I. Väkeväinen
Department of Applied Physics, Aalto University School of Science, P.O. Box
15100, Aalto, FI-00076, Finland  Jani-Petri Martikainen  Department of Applied
Physics, Aalto University School of Science, P.O. Box 15100, Aalto, FI-00076,
Finland  Päivi Törmä  Department of Applied Physics, Aalto University School
of Science, P.O. Box 15100, Aalto, FI-00076, Finland

###### Abstract

Plasmonic structures interacting with light provide electromagnetic resonances
which result in a high degree of local field confinement enabling the
enhancement of light-matter interaction. Plasmonic structures typically
consist of metals which, however, suffer from very high ohmic losses and
heating. High-index dielectrics, on the other hand, can serve as an
alternative material due to their low-dissipative nature and strong scattering
abilities. We study the optical properties of a system composed of all-
dielectric nanoparticle arrays covered with a film of organic dye molecules
(IR-792), and compare these dielectric arrays with metallic nanoparticle
arrays. We tune the light-matter interaction by changing the concentration in
the dye film and report the system to be in the strong coupling regime. We
observe a Rabi splitting between the surface lattice resonances (SLRs) of the
nanoparticle arrays and the absorption line of the dye molecules of up to 253
meV and 293 meV, for the dielectric and metallic nanoparticles, respectively.
The Rabi splitting depends linearly on the square root of the dye molecule
concentration , and we further assess how the Rabi splitting depends on the
film thickness for a low dye molecule concentration. For thinner films of
thicknesses up to 260 nm, we observe no visible Rabi splitting. However, a
Rabi splitting evolves at thicknesses from 540 nm to 990 nm. We perform
finite-difference time-domain simulations to analyze the near-field
enhancements for the dielectric and metallic nanoparticle arrays. The electric
fields are enhanced by a factor of 1200 and 400, close to the particles for
gold and amorphous silicon, respectively, and the modes extend over half a
micron around the particles for both materials.

Keywords: Strong coupling, dielectric nanoparticle arrays, Rabi splitting,
organic dye molecules, surface lattice resonances

##  1 Introduction

If an emitter and an optical mode of its local environment are coupled to each
other and exchange energy, hybrid states of the original modes can arise for a
sufficiently strong coupling. A system is in the strong coupling regime when
the coupling, or Rabi splitting, exceeds the linewidth of the hybrid modes.

Strong coupling has been observed between single emitters and high quality
cavity modes, and in various other systems, for a brief review see the
introduction of Ref. [1]. In this work, we focus on strong coupling in the
context of plasmonics [1] where the strong coupling typically involves
ensembles of emitters, although the single emitter regime has been approached
as well [2]. For ensembles of emitters, the Rabi splitting is in general
proportional to the square root of their concentration. In the strong coupling
regime of plasmonic systems, coherence phenomena can be observed [3, 4].
Plasmonic systems in the strong coupling regime have been reported to form
polariton lasers and condensates [5], exhibit enhanced conductivity [6] and
efficient second harmonic generation [7]. Light strongly coupled to molecular
transitions has also been reported to alter chemical reactions [8, 9].

An interesting plasmonic system where strong coupling has been observed and
utilized are metal nanoparticle arrays. In these arrays, the localized
plasmonic resonances of the nanoparticles hybridize with diffracted orders of
the periodic structure, forming lattice modes called surface lattice
resonances (SLRs) [10, 11, 12]. The SLR modes couple to the absorption line of
emitters, for example organic dye molecules, and reach the strong coupling
regime as evidenced by the square root dependence of the Rabi splitting on the
emitter concentration as reported by Väkeväinen et al. [13]. Strong coupling
in such arrays has already led to new coherence, lasing and condensation
phenomena [4, 5, 14, 15].

For the further exploitation of strong coupling in nanoparticle arrays, it
would be desirable to increase the lifetime of the system, and eventually make
them complementary metal-oxide-semiconductor (CMOS) compatible for integration
with present microelectronic components. The metal nanoparticles have a major
drawback: high absorption losses which lead to broad resonances and short
lifetimes. High-index dielectric materials on the other hand offer low-
dissipative properties and provide not only electric but also magnetic
resonances due to the excitation of circular displacement currents inside the
material [16].

Many plasmonic effects observed in metallic materials can be realized with
high-index dielectric materials. Optically induced resonances have been
realized with different shapes of dielectric nanoparticles, such as spheres,
spheroids, disks, rings and cylinders [17]. Dielectric nanoparticle arrays
have resonances analogous to the localized surface plasmon resonances and
diffracted orders, which form SLRs similar to those observed in metallic
nanoparticle arrays [18].

Silicon, in particular, has advantages such as its well established
fabrication technologies, CMOS-compability and the possibility to integrate it
into existing silicon optics, as well as being low-cost compared to
traditional plasmonic materials such as gold or silver [16]. In addition,
silicon nanoparticles are promising candidates for spectroscopic measurements
close to the particles due to their low heat generation [19]. Strong coupling
in silicon nanostructures has been predicted theoretically between
J-aggregates and silicon nanospheres [20] and demonstrated experimentally in
similar structures [21], as well as in systems of J-aggregates and silicon
nanodisks [22, 23].

In this work, we explore the strong coupling phenomena in dielectric
nanoparticle arrays. The cylindrical nanoparticles consist of amorphous
silicon and are arranged in square lattices. Organic dye molecules (IR-792)
act as emitters. Transmission measurements are conducted with the amorphous
silicon nanoarrays, and, for comparison, with similar gold nanoarrays. The
concentration of the dye molecules is varied and the Rabi splitting extracted.
To further investigate differences between amorphous silicon and gold
nanostructures we perform FDTD simulations using commercially available
software [24]. We find strong coupling between the SLRs of dielectric
nanoarrays and the absorption line of the emitters.

We observe an increase in the Rabi splitting with the square root of the dye
concentration. The amorphous silicon and gold nanoparticle arrays show similar
Rabi splittings of 253 meV and 293 meV, respectively. This makes amorphous
silicon nanoparticle arrays useful for applications as metal nanoparticle
arrays.

##  2 Materials and methods

We fabricate arrays of cylindrical gold and amorphous silicon nanoparticles.
The gold nanocylinders have a diameter of 100 nm and a height of 50 nm
arranged in a square lattice, while the amorphous silicon nanocylinders have a
diameter of 100 nm and a height of 100 nm with the same lattice dimensions. In
both cases the lattice periodicities are varied in 10 nm steps from 580 to 610
nm in x- and y-direction. Figure 1 A shows a scanning electron microscope
(SEM) picture of an amorphous silicon nanoarray.

The samples are fabricated on borosilicate substrates with electron-beam
lithography. To fabricate the gold samples, a resist layer of PMMA A4 is spin-
coated onto the substrate and 10 nm of aluminium is evaporated on top of the
resist to obtain a conductive layer for the electron beam exposure. The layout
for the nanoparticle arrays is patterned into the resist. The aluminium layer
is etched in AZ 351 B and the resist is developed in a IPA:MIBK (1:3)
solution. After this, 2 nm of titanium and 50 nm of gold are evaporated onto
the sample and a lift off is carried out in acetone. The titanium layer is
necessary to enhance the adhesion between the substrate and gold.

For the fabrication of the amorphous silicon samples, 100 nm of amorphous
silicon is deposited onto the substrate with chemical vapor deposition.
Following this, chromium nanoarrays are fabricated on top of the amorphous
silicon similarly to the gold nanoparticles. These chromium nanoarrays serve
as masks for the etching of amorphous silicon via reactive-ion-etching. The
reactive-ion-etching is conducted at a power of 200 W, a pressure of 10 mTorr,
with a SF6-flow of 15 sccm and a CHF3-flow of 50 sccm for 5 minutes. The
chromium mask is etched off afterwards in chromium wet etchant, so that only
the amorphous silicon nanoarrays are left on the sample.

The IR-792 molecules are dissolved into PMMA A3 (3% PMMA in anisol) with
different concentrations and spin-coated onto the samples at 3000 rpm for one
minute, which gives rise to a film thickness of approximately 150 nm. The
arrays are embedded in a symmetric refractive index environment by immersion
in index matching oil and addition of another borosilicate slide on top (n =
1.52).

Different film thicknesses of the dye films are obtained by changing the
viscosity of the PMMA mixture and varying the spin-coating speed. The higher
the amount of PMMA, the more viscous the mixture becomes. The dye molecules
are dissolved into PMMA-anisole mixtures with 3% and 6% PMMA in anisole and
spin-coated with speeds of 1000, 2000 and 3000 rpm. The film thicknesses are
measured with a profilometer.

The transmission measurement setup is depicted in Figure 1 B. The sample is
illuminated by a white light source, collected by a 0.3 NA objective and
guided to a spectrometer. The light excites in-plane propagating modes in
$k_{x}$\- and $k_{y}$-direction, each of which have two polarization
components. The modes with the E-field polarized perpendicular and parallel to
the propagation direction are transverse electric (TE) and transverse magnetic
(TM) modes, respectively. With the spectrometer, we only record the modes
propagating in $k_{y}$-direction and a linear polarizer is used to filter out
the TM modes. Hence, only the TE modes are visible in the measured extinction
spectra.

To characterize the Rabi splitting between the SLR and the absorption line of
the dye, we fit a coupled modes model to the measured extinction spectra.
First, a three-coupled-modes model is used to obtain the SLRs of the arrays
without the dye. The modes involved are the localized surface plasmon
resonances (LSPRs) and the two diffracted orders. Since only the upper
branches, i.e. the modes above the diffracted-order crossing, are involved in
the coupling, it is sufficient to only take these into consideration. After
obtaining the SLRs of the bare arrays, a two-coupled-modes model containing
the SLRs and the main absorption peak of the IR-792 is fitted to the measured
extinction spectra. The coupled modes model between the SLRs and the main
absorption peak of the IR-792 is given by the matrix

$M=\begin{pmatrix}E_{\mathrm{SLR}}-i\gamma_{\mathrm{SLR}}+s&\Omega/2\\\
\Omega/2&E_{\mathrm{dye}}-i\gamma_{\mathrm{dye}}\end{pmatrix},$ (1)

where $E$ and $\gamma$ are the energies and line widths of the uncoupled
modes, respectively, and $\Omega$ is the Rabi splitting. In addition, we
introduced a constant shift, $s$, to account for a shift in the refractive
index introduced by the dye molecule film. Here, the modes of the coupled
system are given by the eigenvalues of the matrix. We extract the peaks of the
SLRs in the measured extinction spectra for each k and apply a least-square-
sum model to minimize the difference between the eigenvalues of the matrix and
the extracted peaks from the measurements. Here, $\Omega$ and $s$ are used as
fitting parameters. Figure 1 C shows the absorption spectrum of IR-792 at a
concentration of 100 mM. The absorption spectra at the other concentrations
used are shown in the Supplementary Material, Figure 2.

The optical properties of the dielectric and gold nanoarrays are studied with
the finite-difference time-domain (FDTD) method by Lumerical’s commercially
available software. The arrays consist of cylinders in diameter and height as
in the experiments. The refractive indices for gold (“Au(Gold) - CRC”) and
titanium (”Ti(Titanium) - CRC”) are used from the material database provided
by Lumerical which are based on the CRC Handbook of Chemistry and Physics
[25]. The refractive index for amorphous silicon is added as a custom made
material according to the measured refractive index of amorphous silicon by Li
et al. [26]. The background refractive index is set to 1.52. To simulate the
array, periodic boundary conditions are used around the unit cell with
perfectly matched layers normal to the array to absorb the excitation and
scattered waves. The mesh size is set to 5 nm in all directions over the whole
simulation region. The excitation source is a plane wave source with a
wavelength range from 550 to 1000 nm with 225 frequency points. The light is
polarized orthogonally to the nanocylinder axis.

##  3 Results

Figure 2 shows the measured extinction spectra in $k$-space for an amorphous
silicon nanoarray without and with dye molecules at concentrations of 25 mM,
50 mM, 75 mM, 100 mM, 150 mM, and 200 mM. Figure 2 A clearly shows the
uncoupled SLRs of the array and Figure 2 C, although dye molecules at a
concentration of 25 mM are present, does not show any difference to the
uncoupled modes. As the dye concentration is increased further, see Figures 2
D-H, the system develops an anticrossing at the point where the energy of the
dye absorption maxima corresponds to the SLR energy. With increasing dye
concentration, the anticrossing or Rabi splitting increases and is expected to
increase linearly with the square root of the dye concentration [1, 11].
Figure 2 B shows a crosscut through Figure 2 A at $k_{y}=0$. The peak in the
spectrum is the $\Gamma$-point at 1.36 eV with a full width at half maximum of
9 meV. The corresponding measurements taken for the gold nanoparticle arrays
are shown in Figure 1 of the Supplementary Material.

To characterize the Rabi splitting, a coupled modes model is applied for gold
and amorphous silicon nanoarrays as described in the Materials and Methods
section. The Rabi splitting is calculated for gold and amorphous silicon
arrays with periodicities ranging from 580 nm to 610 nm in 10 nm steps.
Sufficiently away from the single particle resonance, the SLR modes for
different periodicities have similar density of states and thus couple in the
same way with the molecules. Therefore, we can use the standard deviation of
the values obtained for different periodicities, at the same concentrations,
in estimating the error limits for the Rabi splitting.

The strong coupling condition in the case of two modes with finite linewidths
is given by
$\Omega>\frac{\gamma_{\mathrm{dye}}}{2}+\frac{\gamma_{\mathrm{SLR}}}{2}$,
where $\Omega$ is the Rabi splitting and $\gamma_{\mathrm{dye}}$ and
$\gamma_{\mathrm{SLR}}$ are the linewidths of the uncoupled dye molecules and
SLRs, respectively [1]. It is evident from Figure 2 A that the uncoupled SLR
width stays constant for all $k$. The value for $\gamma_{\mathrm{SLR}}$ was
calculated to be 10 meV for both amorphous silicon and gold as obtained by
taking a crosscut at $k=0.5$ $\mu\mathrm{m}^{-1}$. To obtain
$\gamma_{\mathrm{dye}}$, the sum of two Gaussians is fitted to the measured
absorption spectra of the dye for each concentration. From this fit, the full-
width at half maximum of the main absorption peak is extracted. The obtained
values for $\Omega$ and $\gamma_{\mathrm{dye}}$ are listed in Table 1. The dye
linewidth can be alternatively determined by fitting one Gaussian to the
absorption spectra, see corresponding results in the Supplementary Material
(Table 1 and Figure 3); while the numerical values of the Rabi splitting and
the linewidths are slightly changed, the qualitative behaviour is the same.
The strong coupling condition is met for all concentrations except 25 mM.
Therefore, and since there is no visible splitting for a concentration of 25
mM (see Figure 2 C), it is concluded that the system with such a low
concentration is not in the strong coupling regime (splitting below the
linewidth). Hence, we exclude this particular data point in Figure 3 A.

Figure 3 A shows the Rabi splitting as a function of the square root of the
molecular concentration, $\sqrt{c}$, and a linear dependence is clearly
evident. Both systems, gold and emitters as well as amorphous silicon and
emitters, are in the strong coupling regime. In general, the Rabi splitting is
larger for the gold than for the amorphous silicon nanoparticle arrays.

Now, one has to keep in mind that quite thin films were used for the above
results. If the film covers the whole mode volume of the plasmonic system,
then the concentration determines the Rabi splitting. However, if the mode
volume extends beyond the film thickness, then by increasing the thickness,
one introduces more emitters into the mode volume and the Rabi splitting may
increase with the film thickness also for the same concentration. Therefore,
by varying the film thickness, one can obtain indirect information about the
extent of the mode volume which is in general a difficult question to address
in plasmonics [27, 28].

To investigate this issue, we measure different dye layer thicknesses with 25
mM dye concentration, as shown in Figure 3 B. With thicker dye layers, the
system evolves into the strong coupling regime. It is clear that the splitting
increases with the dye layer thickness, suggesting that the amount of emitters
participating in the interaction with the SLR increases. The Rabi splitting
arises for a 530 nm thick dye layer but its size does not increase much when
increasing the layer thickness further. The saturation of the Rabi splitting
as a function of dye layer thickness implies that the SLR mode volume, despite
having high field enhancement near the nanoparticles, reaches to some extent
up to 700 nm from the array plane. It has been predicted theoretically [29]
and observed in silver hole arrays coupled to CdSe nanoplatelets [30] that the
Rabi splitting first increases linearly and then saturates with increasing
film thickness. The surface plasmon polariton modes decay exponentially from
the metal, and additional emitters cannot couple to the modes beyond the
confinement length. Since the amorphous silicon particles show in general
similar results as the gold particles, we conclude that the modes in the
dielectric nanoparticle arrays behave similarly to the modes of the metal
nanoparticle structures.

To investigate the near field distribution of the arrays and their possible
connection to strong coupling, the structures are simulated using the FDTD-
approach. The results are shown in Figure 4 and Figure 5 for amorphous silicon
and gold, respectively. For both materials, the electric field distribution is
shown at the first and second SLR energies, which correspond to the first and
second diffracted order, as crosscuts through the center of the nanoparticles
in xy- and yz-planes. The first SLR is at 916 nm for gold and at 913 nm for
amorphous silicon, while the second SLR is at 646 nm for gold and at 647 nm
for amorphous silicon, see Figures 4 A and 5 A. In the case of gold, the
second SLR overlaps with the single particle resonance and therefore appears
as a dip in the spectrum. The single particle resonances are located at 690 nm
and 570 nm for gold and amorphous silicon, respectively. Figures 4 B and 5 B
show the measured extinction spectra for comparison to the simulations. In the
case of amorphous silicon, the measured SLRs are located at 914 nm and 654 nm.
For gold, the first SLR is located at 913 nm. The second SLR overlaps with the
single particle resonance, which is at 691 nm. In general, the simulations and
measurements agree very well.

The local electric field intensity for the amorphous silicon is enhanced by a
factor of 400 for the first SLR and by a factor of 80 for the second SLR. In
the case of gold, the local electric field is enhanced by a factor of 1200 for
the first SLR and by a factor of 70 close to the second SLR at 690 nm. The
field enhancements are therefore larger for gold than for amorphous silicon
for this set of parameters. Note that it is difficult to compare the field
enhancement systematically since particles which are of the same size, but
made of different materials will have their single particle resonances at
different energies. Nonetheless, the difference in the near field enhancement
is probably the reason for the larger values for the Rabi splitting in case of
the gold nanoparticle arrays.

The crosscuts through the electric field magnitude in the xy-plane (see
Figures 4 C, E and 5 C, E) and yz-plane (see Figures 4 D, F and 5 D, F) are
shown in a 1.2 $\mu$m range around the particle. In the case of gold
nanoparticles, the field enhancement concentrates on the edges of the
particle, while the case of amorphous silicon, the electric field has a
substantial component inside the particle. Figures 4 D, F and 5 D, F show that
the fields also extend well over half a micron, despite field enhancement
around the particles. This finding, supports the observation that strong
coupling of the SLRs to the dye molecules is possible in thicker films even at
25 mM dye concentration. The simulations were run without the titanium disk
beneath the gold nanoparticles, as we found that the titanium does not
significantly affect the results (not shown).

##  4 Discussion

In summary, we have studied light-matter interaction between all-dielectric
nanoparticle arrays consisting of amorphous silicon and organic dye molecules,
and compared them to a similar system of gold nanoparticle arrays. We found
that the SLRs of the dielectric nanoarrays form strong coupling with the
emitters, which has, to our knowledge, not been observed before. We show a
linear dependence of the Rabi splitting on the dye molecule concentration. We
found Rabi splittings up to 253 eV and 293 eV for amorphous silicon and gold,
respectively. In our work we used square arrays of cylinders. Recently Todisco
et al.[31] demonstrated that further increases in Rabi splittings can be
achieved by using disk dimers with narrow gap in between. This is due to
stronger field localization between disks when incident light polarization is
along the dimer. It would be interesting to explore similar effects in our
systems. We also found that the Rabi splitting which was not observed at a
concentration of 25 mM and a film thickness of 150 nm can be observed at film
thicknesses of around 540 nm. The simulated near-field enhancements were up to
a factor of 1200 close to the particles in the case of gold nanoparticles and
400 for amorphous silicon nanoparticles. The simulations show that the modes
extend considerably beyond the particles in the substrate normal direction.

The simulations and measurements generally agree very well. The SLRs in the
measurements are red-shifted by approximately 3 nm with respect to the SLRs in
the FDTD simulations. This difference is most likely due to small
irregularities in the fabricated arrays and a change in the refractive index
introduced by the presence of the dye molecules and the PMMA film.

In our system, the resonances of amorphous silicon nanoarrays show similar
behaviour in the strong coupling regime and provide similar linewidths and
quality factors for the resonant modes as the gold nanoarrays. However, it has
been proposed that dielectric nanostructures could provide much narrower
linewidths due to their non-dissipative nature [19]. The linewidths are
strongly dependent on size and shape of the particles [32] as well as on the
location of the SLR mode relative to the single particle resonance [33, 26].
The linewidth is strongly related to the lifetime of the mode and plays
therefore an important role in coherence phenomena such as lasing. Strongly
coupled systems of dielectric nanoarrays can therefore lead to new
developments and improvements in lasing applications, for example toward
polariton lasers and condensates with lower thresholds.

##  5 Acknowledgements

Funding: This work was supported by the Academy of Finland under project
numbers 303351, 307419, 327293, 318987 (QuantERA project RouTe), and 318937
(PROFI), by Centre for Quantum Engineering (CQE) at Aalto University, and by
the European Research Council (ERC-2013-AdG-340748-CODE). Part of the research
was performed at the Micronova Nanofabrication Centre, supported by Aalto
University. Competing interests: The authors declare no competing financial
interests. Data and materials availability: All data needed to evaluate the
conclusions in the paper are present in the paper.

## References

  * [1] P Törmä and W L Barnes.  Strong coupling between surface plasmon polaritons and emitters: a review.  Reports on Progress in Physics, 78(1):013901, January 2015.
  * [2] Vasyl G Kravets, Andrei V Kabashin, William L Barnes, and Alexander N Grigorenko.  Plasmonic surface lattice resonances: a review of properties and applications.  Chemical reviews, 118(12):5912–5951, 2018.
  * [3] S Aberra Guebrou, C Symonds, E. Homeyer, et al.  Coherent emission from a disordered organic semiconductor induced by strong coupling with surface plasmons.  Physical Review Letters, 108(6):066401, 2012.
  * [4] L. Shi, T. K. Hakala, H. T. Rekola, J.-P. Martikainen, R. J. Moerland, and P. Törmä.  Spatial coherence properties of organic molecules coupled to plasmonic surface lattice resonances in the weak and strong coupling regimes.  Physical Review Letters, 112(15), April 2014.  arXiv: 1404.4160.
  * [5] Aaro I Väkeväinen, Antti J Moilanen, Marek Nečada, Tommi K Hakala, Konstantinos S Daskalakis, and Päivi Törmä.  Sub-picosecond thermalization dynamics in condensation of strongly coupled lattice plasmons.  arXiv preprint arXiv:1905.07609, 2019.
  * [6] E Orgiu, J George, JA Hutchison, et al.  Conductivity in organic semiconductors hybridized with the vacuum field.  Nature Materials, 14(11):1123, 2015.
  * [7] Thibault Chervy, Jialiang Xu, Yulong Duan, et al.  High-efficiency second-harmonic generation from hybrid light-matter states.  Nano letters, 16(12):7352–7356, 2016.
  * [8] Raphael F Ribeiro, Luis A Martínez-Martínez, Matthew Du, Jorge Campos-Gonzalez-Angulo, and Joel Yuen-Zhou.  Polariton chemistry: controlling molecular dynamics with optical cavities.  Chemical science, 9(30):6325–6339, 2018.
  * [9] Manuel Hertzog, Mao Wang, Jürgen Mony, and Karl Börjesson.  Strong light–matter interactions: a new direction within chemistry.  Chemical Society Reviews, 48(3):937–961, 2019.
  * [10] FJ Garcia De Abajo.  Colloquium: Light scattering by particle and hole arrays.  Reviews of Modern Physics, 79(4):1267, 2007.
  * [11] Robert J. Moerland, Tommi K. Hakala, Jani-Petri Martikainen, Heikki T. Rekola, Aaro I. Väkeväinen, and Päivi Törmä.  Strong coupling between organic molecules and plasmonic nanostructures.  In Sergey I. Bozhevolnyi, Luis Martin-Moreno, and Francisco Garcia-Vidal, editors, Quantum Plasmonics, pages 121–150, Cham, 2017. Springer International Publishing.
  * [12] Weijia Wang, Mohammad Ramezani, Aaro I Väkeväinen, Päivi Törmä, Jaime Gómez Rivas, and Teri W Odom.  The rich photonic world of plasmonic nanoparticle arrays.  Materials today, 21(3):303–314, 2018.
  * [13] A. I. Väkeväinen, R. J. Moerland, H. T. Rekola, et al.  Plasmonic surface lattice resonances at the strong coupling regime.  Nano Letters, 14(4):1721–1727, April 2014.
  * [14] Mohammad Ramezani, Alexei Halpin, Antonio I Fernández-Domínguez, et al.  Plasmon-exciton-polariton lasing.  Optica, 4(1):31–37, 2017.
  * [15] Milena De Giorgi, Mohammad Ramezani, Francesco Todisco, et al.  Interaction and coherence of a plasmon–exciton polariton condensate.  Acs Photonics, 5(9):3666–3672, 2018.
  * [16] Isabelle Staude and Jörg Schilling.  Metamaterial-inspired silicon nanophotonics.  Nature Photonics, 11(5):274, 2017.
  * [17] Arseniy I. Kuznetsov, Andrey E. Miroshnichenko, Mark L. Brongersma, Yuri S. Kivshar, and Boris Luk’yanchuk.  Optically resonant dielectric nanostructures.  Science, 354(6314):aag2472, November 2016.
  * [18] Andrey B Evlyukhin, Carsten Reinhardt, Andreas Seidel, Boris S Luk’yanchuk, and Boris N Chichkov.  Optical response features of si-nanoparticle arrays.  Physical Review B, 82(4):045404, 2010.
  * [19] Alex Krasnok, Martín Caldarola, Nicolas Bonod, and Andrea Alú.  Spectroscopy and biosensing with optically resonant dielectric nanostructures.  Advanced optical materials, 6(5), 2018.
  * [20] Christos Tserkezis, PAD Gonçalves, Christian Wolff, Francesco Todisco, Kurt Busch, and N Asger Mortensen.  Mie excitons: Understanding strong coupling in dielectric nanoparticles.  Physical Review B, 98(15):155439, 2018.
  * [21] Hao Wang, Yanlin Ke, Ningsheng Xu, et al.  Resonance coupling in silicon nanosphere–j-aggregate heterostructures.  Nano letters, 16(11):6886–6895, 2016.
  * [22] Shao-Ding Liu, Jin-Li Fan, Wen-Jie Wang, Jing-Dong Chen, and Zhi-Hui Chen.  Resonance coupling between molecular excitons and nonradiating anapole modes in silicon nanodisk-j-aggregate heterostructures.  ACS Photonics, 5(4):1628–1639, 2018.
  * [23] Francesco Todisco, Radu Malureanu, Christian Wolff, et al.  Magnetic and electric mie-exciton polaritons in silicon nanodisks.  arXiv preprint arXiv:1906.09898, 2019.
  * [24] Lumerical Inc.  https://www.lumerical.com/products/.
  * [25] David R Lide.  CRC handbook of chemistry and physics: a ready-reference book of chemical and physical data.  CRC press, 1995.
  * [26] Jiaqi Li, Niels Verellen, and Pol Van Dorpe.  Engineering electric and magnetic dipole coupling in arrays of dielectric nanoparticles.  Journal of Applied Physics, 123(8):083101, February 2018.
  * [27] R Ruppin.  Electromagnetic energy density in a dispersive and absorptive material.  Physics letters A, 299(2-3):309–312, 2002.
  * [28] A Femius Koenderink.  On the use of purcell factors for plasmon antennas.  Optics letters, 35(24):4208–4210, 2010.
  * [29] A González-Tudela, PA Huidobro, Luis Martín-Moreno, C Tejedor, and FJ García-Vidal.  Theory of strong coupling between quantum emitters and propagating surface plasmons.  Physical review letters, 110(12):126801, 2013.
  * [30] Jan M Winkler, Freddy T Rabouw, Aurelio A Rossinelli, et al.  Room-temperature strong coupling of cdse nanoplatelets and plasmonic hole arrays.  Nano letters, 19(1):108–115, 2018.
  * [31] Francesco Todisco, Marco Esposito, Simone Panaro, Milena De Giorgi, Lorenzo Dominici, Dario Ballarini, Antonio I. Fernández-Domínguez, Vittorianna Tasco, Massimo Cuscunà, Adriana Passaseo, Cristian Ciracì, Giuseppe Gigli, and Daniele Sanvitto.  Toward Cavity Quantum Electrodynamics with Hybrid Photon Gap-Plasmon States.  ACS Nano, 10(12):11360–11368, 2016.
  * [32] Jon A Schuller and Mark L Brongersma.  General properties of dielectric optical antennas.  Optics express, 17(26):24084–24095, 2009.
  * [33] Andrey G Nikitin, Andrei V Kabashin, and Hervé Dallaporta.  Plasmonic resonances in diffractive arrays of gold nanoantennas: near and far field effects.  Optics express, 20(25):27941–27952, 2012.

Table 1: Measured linewidth of the dye molecules, $\gamma_{\mathrm{dye}}$, the strong coupling condition, $(\gamma_{\mathrm{dye}}+\gamma_{\mathrm{SLR}})/2$, and Rabi splittings of amorphous silicon and gold, $\Omega_{\mathrm{a-Si}}$ and $\Omega_{\mathrm{Au}}$, for different concentrations, $c$, for a film thickness of 150 nm. The molecule linewidth is obtained by fitting two Gaussians to the measured absorption spectrum (see Figure 1), and choosing the underlying Gaussian for the main absorption peak. c [mM] |  $\gamma_{\mathrm{dye}}$ [eV] |  $(\gamma_{\mathrm{dye}}+\gamma_{\mathrm{SLR}})/2$ [eV] |  $\Omega_{\mathrm{a-Si}}$ [eV] |  $\Omega_{\mathrm{Au}}$ [eV]
---|---|---|---|---
25 | 0.10 | 0.06 | - | -
50 | 0.11 | 0.06 | 0.13 | 0.13
75 | 0.12 | 0.07 | 0.16 | 0.19
100 | 0.11 | 0.06 | 0.17 | 0.20
150 | 0.12 | 0.07 | 0.21 | 0.23
200 | 0.13 | 0.07 | 0.25 | 0.29

Figure 1: (A) SEM picture of an amorphous silicon sample before chromium
etching. (B) Transmission measurement setup and schematic of the sample with
the dye. (C) Absorption spectrum of IR-792 with a concentration of 100 mM
solved in PMMA A3. We fit two Gaussians to the spectrum, shown by the pink and
blue curves. The pink curve corresponds to the main absorption peak.

Figure 2: Dispersion measurements of a nanoarray of amorphous silicon
nanoparticles with a periodicity of 600 nm in x- and y-direction, diameter of
100 nm and dye concentrations of (A) 0 mM (C) 25 mM, (D) 50 mM, (E) 75 mM, (F)
100 mM, (G) 150 mM and (H) 200 mM. The film thickness is 150 nm. Maximum
extinction is given in yellow. (B) shows a crosscut through (A) at $k_{y}=0$.
In (A) the red dashed lines depict the uncoupled diffracted orders of the
lattice and the black solid lines the SLRs obtained by the coupled mode
fitting of the extinction maxima. In (C), (D), (E), (F), (G) and (H) the
horizontal red dashed line depicts the absorption line of the dye molecules
and the diagonal red lines correspond to the fitted SLRs (upper branch) from
Figure 2 (A). The solid black lines depict the coupled mode fitting of the
extinction maxima. It is clearly visible that the Rabi splitting increases
with increasing dye molecule concentration.

Figure 3: Dependence of the Rabi splitting on the square root of the
concentration $\sqrt{c}$ for a film thickness of 150 nm (A) and Rabi splitting
for $c$=25 mM as a function of different film thicknesses (B). In (A), the
data point at $c=25$ mM is excluded since the splitting is smaller than the
linewidth and thus not measurable.

Figure 4: FDTD simulation results for amorphous silicon. (A) Simulated
transmission spectrum and (B) measured spectra at $k_{y}$=0 for comparison.
Electric field intensity distribution at resonance at 646 nm crosscuts in the
xy-plane (C) and xz-plane (D). Electric field intensity distribution at
resonance at 913 nm crosscuts in the xy-plane (E) and xz-plane (F).

Figure 5: FDTD simulation results for gold. (A) Simulated transmission
spectrum and (B) measured spectra at $k_{y}$=0 for comparison. Electric field
intensity distribution at resonance at 690 nm crosscuts in the xy-plane (C)
and xz-plane (D). Electric field intensity distribution at resonance at 916 nm
crosscuts in the xy-plane (E) and xz-plane (F). As it is difficult to see from
the color scale, we note that the field does not decay to negligible values
within the shown $z$-range.

## Supplementary Material

Table 1: Measured linewidth of the dye molecules, $\gamma_{\mathrm{dye}}$, the strong coupling condition, $(\gamma_{\mathrm{dye}}+\gamma_{\mathrm{SLR}})/2$, and Rabi splittings of amorphous silicon and gold, $\Omega_{\mathrm{a-Si}}$ and $\Omega_{\mathrm{Au}}$, for different concentrations, $c$, for a film thickness of 150 nm. The molecule linewidth is obtained by fitting one Gaussian to the measured absorption spectrum. c [mM] |  $\gamma_{\mathrm{dye}}$ [eV] |  $(\gamma_{\mathrm{dye}}+\gamma_{\mathrm{SLR}})/2$ [eV] |  $\Omega_{\mathrm{a-Si}}$ [eV] |  $\Omega_{\mathrm{Au}}$ [eV]
---|---|---|---|---
25 | 0.14 | 0.08 | - | -
50 | 0.17 | 0.09 | 0.18 | 0.19
75 | 0.20 | 0.11 | 0.24 | 0.25
100 | 0.20 | 0.11 | 0.24 | 0.27
150 | 0.23 | 0.12 | 0.30 | 0.31
200 | 0.28 | 0.15 | 0.36 | 0.39
Figure 1: Dispersion measurements of a nanoarray of gold nanoparticles with a
periodicity of 600 nm in x- and y-direction, diameter of 100 nm and dye
concentrations of (A) 0 mM (C) 25 mM, (D) 50 mM, (E) 75 mM, (F) 100 mM, (G)
150 mM and (H) 200 mM. The film thickness is 150 nm. Maximum extinction is
given in yellow. (B) shows a crosscut through (A) at $k_{y}=0$. In (A) the red
dashed lines depict the uncoupled diffracted orders of the lattice and the
black solid lines the SLRs obtained by the coupled mode fitting of the
extinction maxima. In (C), (D), (E), (F), (G) and (H) the horizontal red
dashed line depicts the absorption line of the dye molecules and the diagonal
red lines correspond to the fitted SLRs (upper branch) from Figure 2 (A). The
solid black lines depict the coupled mode fitting of the extinction maxima. As
can be seen the Rabi splitting increases with increasing dye molecule
concentration.

Figure 2: Absorption spectra of IR-792 solved in PMMA A3 with dye
concentrations of (A) 25 mM, (B) 50 mM, (C) 75 mM, (D) 100 mM, (E) 150 mM and
(F) 200 mM. We fit two Gaussians to the spectrum, shown by the pink and blue
curves. The pink curve corresponds to the main absorption peak. Figure 3:
Dependence of the Rabi splitting on the square root of the concentration
$\sqrt{c}$ for a film thickness of 150 nm. The data point $c=25$ mM is
excluded since the splitting is smaller than the linewidth and thus not
measurable. For this fit, the molecule linewidth is obtained by fitting one
Gaussian to the measured absorption spectrum. For corresponding values, see
Table 1.


# Methanol and water maser observations separate disc and outflow sources in
IRAS 19410+2336

Darwish, M. S.1,2, Edris, K. A.3, Richards, A. M. S.4, Etoka, S.4, Saad, M.
S.1,2, Beheary, M. M.3, Fuller, G. A.4
1Astronomy Department, National Research Institute of Astronomy and Geophysics
(NRIAG), 11421 Helwan, Cairo, Egypt.
2Kottamia Center of Scientific Excellence in Astronomy and Space Science
(KCScE, STDF No. 5217, ASRT), Cairo, Egypt.
3Astronomy and Meteorology Department, Faculty of Science, Al-Azhar
University, Cairo, Egypt.
4Jodrell Bank Centre for Astrophysics, Department of Physics $\&$ Astronomy,
The University of Manchester, M13 9PL, UK.
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We investigate the kinematics of high mass protostellar objects within the
high mass star forming region IRAS 19410+2336. We performed high angular
resolution observations of 6.7-GHz methanol and 22 GHz water masers using the
MERLIN (Multi-Element Radio Linked Interferometer Network) and e-MERLIN
interferometers. The 6.7-GHz methanol maser emission line was detected within
the $\sim$ 16–27 km s-1 velocity range with a peak flux density $\sim$50 Jy.
The maser spots are spread over $\sim$1.3 arcsec on the sky, corresponding to
$\sim$2800 au at a distance of 2.16 kpc. These are the first astrometric
measurements at 6.7 GHz in IRAS 19410+2336. The 22-GHz water maser line was
imaged in 2005 and 2019 (the latter with good astrometry). Its velocities
range from 13 to $\sim$29 km s-1. The peak flux density was found to be 18.7
Jy and 13.487 Jy in 2005, and 2019, respectively. The distribution of the
water maser components is up to 165 mas, $\sim$350 au at 2.16 kpc. We find
that the Eastern methanol masers most probably trace outflows from the region
of millimetre source mm1. The water masers to the West lie in a disc (flared
or interacting with outflow/infall) around another more evolved millimetre
source (13-s). The maser distribution suggests that the disc lies at an angle
of 60∘ or more to the plane of the sky and the observed line of sight
velocities then suggest an enclosed mass between 44 M⊙ and as little as 11 M⊙
if the disc is edge-on. The Western methanol masers may be infalling.

###### keywords:

Stars: formation – stars: massive – stars: individual: IRAS 19410+2336 –
masers – ISM

††pubyear: 2019††pagerange: Methanol and water maser observations separate
disc and outflow sources in IRAS 19410+2336–Methanol and water maser
observations separate disc and outflow sources in IRAS 19410+2336

##  1 Introduction

The ambiguity of how massive stars were formed makes studying their formation
one of the most critical current issues in astrophysics. Both observation and
theory face many challenges (Beuther et al., 2007; Zinnecker & Yorke, 2007).
The challenge from the observational side is due to several reasons: Firstly,
they form at large distances ($>1$ kpc, with the exceptions of the Orion
(Hirota et al., 2007) and Cepheus A (Torrelles et al., 2012) regions).
Secondly their formation happens on much shorter time scales than their low
mass counterparts. Furthermore, high-mass star formation takes place inside
high density clouds, obscuring our view of the process. Theoretically, the
fundamental problem in high-mass star formation is known as the "radiation
pressure problem". That is, the radiation pressure from the forming
protostellar object is predicted to stop the accretion process when the star
reaches a limit $\sim$10M$\odot$ (Wolfire & Cassinelli, 1987). Since more
massive stars do exist, the theoretical models must deal with this challenge
to explain how the mass build up proceeds.

So far, the most widely accepted model for high-mass star formation is based
on a ’scaled-up’ version of the low-mass model, but with a few alterations
(e.g. Keto & Zhang, 2010, and references their in). This model suggested that
in order for stars to gain enough mass to be considered high-mass stars, they
must undergo much higher accretion rates through a disc which allows the
radiation pressure to be overcome (Tan & McKee, 2004; Krumholz et al., 2009).

Interstellar masers are believed to be one of the best tracers of massive star
forming regions (Ellingsen et al., 2007) and the study of many maser
transitions has shed a lot of light on the high-mass star formation (HMSF)
process. Masers are used to characterize circumstellar discs, bipolar outflows
(Stahler et al., 2000), magnetic fields (e.g. Cohen, 1984; McCall et al.,
1990; Edris et al., 2005; Etoka et al., 2012) and internal motions (e.g. Reid
et al., 2009; Sanna et al., 2010; Sugiyama et al., 2011). Ellingsen et al.
(2007) provide a rudimentary time line for classifying the evolutionary stage
of young high-mass stars.

H2O water and CH3OH class II methanol masers are believed to be closely
associated with the early stages of massive star formation (Caswell, 1997;
Szymczak & Gérard, 2004). It is known for quite some time now that class-II
6.7-GHz methanol masers are commonly found associated with ground state OH
masers within $\simeq$1 arcsec radius (Caswell, 1996). High angular resolution
imaging has revealed that this association can reach $\sim$100 mas radius but
close association (i.e., down to a $\simeq$10 mas radius) is rare (Etoka et
al., 2005) due to the combination of two factors: the important role played by
grain-mantle evaporation in the production of both OH and CH3OH during HMSF
(Hartquist et al., 1995) and the expected evolutionary offset in the
appearance and duration of these two types of maser emission, where the
results of large surveys (Breen et al., 2010a) and (Breen et al., 2010b)
suggest that although OH and methanol can be found in the same regions in
agreement with the models of Cragg et al. (2002), methanol masers appear first
and disappear while OH and water are still detectable.

Another well known association is that of water and Class-II 6.7-GHz methanol,
with their exclusive association (i.e., no sign of OH maser emission) thought
to be tracing earlier stages of the star formation process (Ellingsen et al.,
2007; Breen et al., 2010a). The spatial correlation between H2O and CH3OH
masers towards a number of sources has been studied by Beuther et al. (2002a),
who concluded that although the typical separation between H2O and CH3OH
masers is 0.03 pc there is no evidence for closer associations. They felt that
this was unlikely, given that the maser species are collisionally and
radiatively pumped, respectively. Breen et al. (2018) compared the few-hundred
sources in the HOPS (H2O Southern Galactic Plane) and MMB (Methanol Multibeam)
surveys and found that a majority of associated methanol and water masers
occurred within 1$\aas@@fstack{\prime\prime}$0, or 0.024 pc (at a typical
distance of 5 kpc), comparable to the combined astrometric accuracy.

Additionally, these masers can be detected originating from deep inside the
dusty molecular envelope of high-mass stars in early stages of formation,
making them a sensitive probe for discovering stars in their embryonic state.

Considering these issues, we present high angular resolution observations of
H2O and CH3OH masers toward the high mass star forming region IRAS 19410
+2336. The previous observations carried out toward the IRAS 19410+2336 region
suggest that it is a young massive star forming region in an early
evolutionary state with a complex morphology (Qiu et al., 2008, and references
therein). Kurtz et al. (1994) reported that there is no Ultracompact H II
region detected at 8.4 GHz, while Sridharan et al. (2002) detected weak 3.6 cm
continuum emission of 1 mJy.

Beuther et al. (2002b), using the IRAM 30-m MAMBO bolometer (pixel spacing
$\sim 5\aas@@fstack{\prime\prime}5$), found a $\lambda$ 1.2-mm source mm1 at
$\alpha$(J2000) = 19h 43m 11s.24, $\delta$(J2000) = 23∘ 44′
$3\aas@@fstack{\prime\prime}4$. The total astrometric accuracy of the peak is
$0\aas@@fstack{\prime\prime}5$ or better but the source is extended over
$>5^{\prime\prime}$. One of the $\lambda$ 3-mm sources mapped by Rodón et al.
(2012) using the IRAM interferometer has a similar extent and position,
$\alpha$(J2000) = 19h 43m 11s.206, $\delta$(J2000) = 23∘ 44′
$3\aas@@fstack{\prime\prime}29$. A number of their detections at $\lambda$ 1.2
mm fall in this region, the closest being 13-s at $\alpha$(J2000) = 19h 43m
11s.206, $\delta$(J2000) = 23∘ 44′ $3\aas@@fstack{\prime\prime}29$, with an
astrometric error estimated at $0\aas@@fstack{\prime\prime}18$.

The 6.7-GHz methanol maser toward IRAS 19410+2336 was first discovered by
Menten (1991). He found 6.7-GHz emission in the velocity range from 15 to 28
km s-1 with a peak flux density of 103 Jy. MacLeod et al. (1993) reported a
similar detection and also monitored the 12-GHz transition, which had a peak
flux density of 10.6 Jy at 27 km s-1. The first recorded imaging of this
region at 22 GHz, in 1979, was by Lada et al. (1981). They reported that the
water maser had a velocity range from 12 to 20.6 km s-1 with a peak flux
density of 21 Jy.

The distance to this source is estimated to be 2.16 kpc (Xu et al., 2009),
while the total infrared luminosity is around 104 ${L_{\odot}}$ (Martín-
Hernández et al., 2008). The average temperature for the whole region
estimated to be 40$\pm$15 K, while the temperature of the dense cores is
80$\pm$40 K (Rodón et al., 2012).

To investigate what the CH3OH and H2O masers trace and how they are
distributed and correlated with other tracers, IRAS 19410+2336 has been
observed using MERLIN and e-MERLIN at high angular resolution. The details of
the observations and data reduction are presented in Sec. 2 while the results
are introduced in Sec.3. The discussion is given in Sec. 4, and conclusions
drawn in Sec. 5.

Table 1: MERLIN and e-MERLIN observational parameters for methanol and water masers. Observational parameters | CH3OH | H2O
---|---|---
Date of observation | 31-Dec-2004 | 01-Jan-2005 | 20-March-2019
No. antenna | Five antennas | Five antennas | Five antennas
Field centre ($\alpha$, $\delta$ J2000) | 19${}^{h}\\!$ 43${}^{m}\\!$ 11s.247, 23${}^{\circ}\\!$ 44${}^{\prime}\\!$ $03\aas@@fstack{\prime\prime}32$ | 19${}^{h}\\!$ 43${}^{m}\\!$ 11s.184, 23${}^{\circ}\\!$ 44${}^{\prime}\\!$ $02\aas@@fstack{\prime\prime}961$
Rest frequency (MHz) | 6668.518 | 22235.079 | 22235.079
No. of frequency channels | 511 | 255 | 512
Total band width (MHz) | 0.5 | 4 | 4
Bandpass calibrator | 0552+398 | 3C273 | 3C84
Phase calibrator | 1932+204A | 1932+204A | J1946+2300
Restoring beam (mas$\times$mas) | 70$\times$45 | 20$\times$20 | 29$\times$28
rms in quiet channel (Jy beam-1) | 0.0346 | 0.0707 | 0.050

##  2 Observations and data reduction

Our initial observations of IRAS 19410+2336 were carried out using the MERLIN
interferometer. Observations in dual-polarization (LL,RR) of the Class-II
6.7-GHz methanol and 22-GHz water maser emission were performed in December
2004 and January 2005, respectively. In both cases, the target data were
adjusted to fixed velocity with respect to the LSR (Local Standard of Rest) in
the correlator. The data were extracted from the MERLIN archive and converted
to FITS at Jodrell Bank Center for Astrophysics (JBCA) using local software
(dprogs) and the AIPS (Astronomical Image processing system) software package.
CASA (the Common Astronomy Software Application) software package version
5.1.1 was used in order to calibrate and reduce the data. The observational
parameters for IRAS 19410+2336 are listed in Table 1. The observations and
data reduction followed normal MERLIN procedures, see the MERLIN User Guide
(Diamond et al., 2003).

Subsequent observations of the 22-GHz water maser emission were performed with
eMERLIN in March 2019; the data were reduced manually in CASA following a
strategy similar to that described in
_https://github.com/e-merlin/eMERLIN_CASA_pipeline_.

###  2.1 CH3OH masers

The 6.7-GHz methanol line was observed in a spectral bandwidth of 0.5 MHz
corresponding to 30 km s-1 velocity range with a channel separation of 0.04 km
s-1. The source 0552+398 was observed as a bandpass calibrator. At this
frequency the flux density for 0552+398 was set to be 4.6 Jy based on previous
scaling using 3C286, (Baars et al., 1977), allowing for the resolution of
MERLIN.

As a phase reference calibrator, the compact Quasar 1932+204A was observed
through the observations at $\alpha$(J2000) = 19h 35m 10s.473, $\delta$(J2000)
= 20∘ 31′ $54\aas@@fstack{\prime\prime}154$ with a cycle of phase
reference:target 1.7:5.3 min, in a useful bandwidth of 13 MHz overlapping the
line band.

0552+398 was observed in both bandwidths and used to measure the phase offset
between configurations. This correction, along with the bandpass table and the
phase reference solutions for phase and amplitude, were applied to the target.
The brightest target channel, at $V_{\mathrm{LSR}}$ 26.9 km s-1 was used for
self-calibration and the solutions were applied to all channels. The absolute
position errors for 6.7-GHz methanol masers are 12 mas in each of RA and Dec,
derived from the error contributions described by Diamond et al. (2003).

The position used for astrometry was measured before self-calibration but any
changes were negligible. The CLEAN algorithm as implemented in the CASA
package "tclean" was used for cleaning and de-convolution of the image cube.
This was exported as FITS to allow component fitting in AIPS. We edited the
velocity labeling to correct inconsistencies between CASA and AIPS in order to
ensure that the data for analysis were correctly labeled in velocities with
respect to the LSR. Some parameters of the final cube are given in Table 1.
Table 2 lists the absolute positions as well as the corresponding velocities
of the brightest feature for both methanol and water masers.

In order to fit the components through the different channels in a cube image
we used the AIPS task SAD (Search And Destroy). The positions of the maser
components were determined via 2D Gaussian fitting above 3$\sigma$ threshold
in each channel map. The components were required to appear in at least three
consecutive channels with positions within 20 mas of their occurrence in each
consecutive channel, such groups forming spectral features. The peak
intensities of these features with their relative positions and velocities
were consider to represent the group. The spectra in Figs. 1 and 2 show the
component flux densities for each separate group. In most cases the spectral
profiles are approximately Gaussian although some are blended.

###  2.2 H2O masers

The 22 GHz water maser line was observed in 2005 with a spectral bandwidth of
4 MHz corresponding to 50 km s-1 velocity range with a channel separation of
0.21 km s-1.

In order to calibrate the variation of instrumental gain and phase across the
spectral bandpass, the source 3C273 was selected to be a bandpass and flux
calibrator. The flux density for 3C273 was set to be 24 Jy (flux monitoring
kindly supplied by Tersranta, Metsahovi, private communication). Observations
of 1932+204A were again carried out for phase referencing but unfortunately it
was too faint and far away to be useful at 22 GHz. The brightest maser channel
was therefore used for self-calibration and the solutions applied to all
channels, but this means that there is no improvement of the astrometry over
the input observing position.

A short e-MERLIN observation was made in 2019 using a spectral bandwidth of 4
MHz in 512 channels, giving a channel separation of 0.105 km s-1, at constant
frequency. The upgraded array allows greater sensitivity due to new receiver
systems, a broader continuum bandwidth thanks to optical fibers and a new,
more flexible correlator. Thus, we could also observe calibration sources at
125 MHz bandwidth. 3C84 was used as a bandpass and flux scale calibrator
(assuming a flux density of 20 Jy) and also to calculate the phase offset
between the 4- and 125-MHz configurations. J1925+2106, 4.8∘ from the target,
flux density about 2 Jy beam-1, was observed every 20 min as time-dependent
amplitude calibrator. J1946+2300 was used as the phase reference source, at a
position of $\alpha$(J2000) = 19h 46m 06s.2510, $\delta$(J2000) = 23∘ 00′
$04\aas@@fstack{\prime\prime}414$ (a separation of 0.99∘ from the target) with
a flux density of 0.07 Jy beam-1.

We applied the corrections derived from these calibrators to IRAS 19410+2336.
The IRAS 19410+2336 dataset was then transformed to constant velocity in the
LSR frame. The best few hours of data were used to image a channel at 23.7 km
s-1 which contained a bright peak at $\alpha$(2000)= 19h 43m 11s.184,
$\delta$(2000)= 23∘ 44′ $03\aas@@fstack{\prime\prime}008$, at a signal-to-
noise ratio of 23, giving a stochastic position error of 2 mas.

The largest quantifiable contribution to position uncertainty is the phase
change in the offset between phase-reference and target, contributing 4 mas.
Allowing for other errors, due to antenna position uncertainties and
calibration errors, the astrometric accuracy of this position is better than
10 mas. This peak was then used for self-calibration, applying the corrections
to all line data, and a data cube made with a restoring beam of 29
mas$\times$28 mas, rms $\sim 50$ mJy except in the brightest channels which
were dynamic range limited.

Other aspects of data reduction were similar to that for methanol and “tclean"
was used to make a final data cube; some parameters are given in Table 1. We
fitted components to the emission and grouped the spots into features as
described in Section 2.1, with a positional match requirement of $\sim$20 mas.

Table 2: The absolute positions and the radial velocities of the brightest maser feature at epochs 2005 (methanol) and 2019 (water). The positions leading terms are $\alpha$(J2000) = 19h 43m, $\delta$(J2000) = 23∘ 44′. Type of maser |  $\alpha$(J2000) |  $\delta$(J2000) |  $V_{\mathrm{LSR}}$ (km s-1)
---|---|---|---
CH3OH | 11.238$\pm$0s.013 | 3.297$\pm$0′′.015 | 26.548
H2O | 11.1840$\pm$0s.1 | 2.961$\pm$0′′.1 | 26.420

##  3 Results

The main aim of this work was to understand better the kinematics of the high
mass star forming region IRAS 19410+2336 as well as to investigate the
possibility of determining the position of the embedded protostellar objects.
The 6.7-GHz CH3OH and 22 GHz water masers were detected using MERLIN and
e-MERLIN. In addition to the distribution of each maser type, the positional
and velocity relationships between the two masers species were examined.
Further details about each maser type are discussed below.

Figure 1: Panel (a) shows the spectra of H2O maser associated with IRAS
19410+2336, the numbers represent the features listed in table 3, the blue
spectra is for the 2005 and the cyan one represents the 2019 epoch, while (b)
represents the distribution of the H2O maser components toward the IRAS
19410+2336 region for 2005 epoch (plus signs) and 2019 epoch (circles), the
size of the symbol is proportional to the log of their intensity.

###  3.1 Water masers

The H2O masers in the IRAS 19410+2336 region were observed at two epochs, 2005
and 2019. Five spectral features were identified in 2005 while four features
were found in 2019. Table 3 lists the peak intensities, positions and
velocities for both water maser epochs. Figure 1 (a) represents the spectral
components observed in both epochs. In 2005, the velocity range covered by the
emission was 12–29 km s-1, grouped into 5 features with the highest flux
density determined to be 18.7 Jy and the faintest feature peak 3.304 Jy.

In 2019 we identified four features with peaks in the range 18.44 – 28.76 km
s-1, with peak flux densities ranging from 1.6 to 13.5 Jy. Features 2005 a and
d have counterparts close in velocity in 2019 f and i, at around 19 and 28 km
s-1. 2005 features b and c lie between a and f but are slightly less red-
shifted whilst 2019 features g and h lie between d and i and are slightly less
blue-shifted. The most blue-shifted feature, e, was only detected by us in
2005.

Figure 1 (b) shows that the water maser components are distributed within an
area about 165 mas across, $\sim$350 au at 2.16 kpc. The masers detected in
2019 are about 90 mas (200 au) in angular extent, also with the more blue-
shifted emission to the North. The mid point of the extreme velocity peaks is
$\alpha$(J2000) = 19h 43m 11s.1840, $\delta$(J2000) = 23∘ 44′
$2\aas@@fstack{\prime\prime}961$.

At all epochs, the 22 GHz H2O maser peaks fall within a velocity range between
12–29 km s-1 with a maximum angular extent of 200 mas. The more blue-shifted
features are always the most northerly. Only the 2019 data have
astrometrically accurate positions. We assumed that the masers are centred on
a common centre of excitation, at an intermediate position and velocity. We
aligned the 2005 data with the 2019 data by assuming that the red- and blue-
shifted clumps were offset symmetrically about a common centre at both epochs.

Table 3: The properties of the detected H2O maser features toward IRAS 19410+2336 with relative position errors, the positions leading terms are $\alpha$(J2000) = 19h 43m, $\delta$(J2000) = +23∘ 44′. The top and bottom halves are 2005 and 2019 data, respectively. Epoch | No. | Vel.(km s-1) | Flux(Jy) |  $\alpha$ (s) | Error(ms) |  $\delta$ (′′) | Error(mas)
---|---|---|---|---|---|---|---
2005 | a | 28.11 | 3.304 $\pm$ 0.070 | 11.1807 | 0.016 | 2.879 | 0.236
| b | 27.06 | 4.588 $\pm$ 0.079 | 11.1807 | 0.012 | 2.881 | 0.188
| c | 26.42 | 18.708 $\pm$ 0.078 | 11.1812 | 0.002 | 2.890 | 0.045
| d | 19.89 | 3.743 $\pm$ 0.079 | 11.1873 | 0.014 | 3.043 | 0.227
| e | 13.36 | 10.013 $\pm$ 0.119 | 11.1868 | 0.008 | 3.035 | 0.129
2019 | f | 28.76 | 1.784 $\pm$ 0.102 | 11.1851 | 0.036 | 2.915 | 0.5
| g | 24.76 | 2.287 $\pm$ 0.162 | 11.1836 | 0.094 | 3.006 | 0.7
| h | 23.39 | 13.487 $\pm$ 0.205 | 11.1844 | 0.007 | 3.008 | 0.1
| i | 18.44 | 1.699 $\pm$ 0.135 | 11.1829 | 0.072 | 3.007 | 1.0

The results, especially the disappearance of the 13 km s-1 feature in 2019,
confirm the suggestion by Caswell et al. (1995); Szymczak et al. (2000) that
the region is strongly variable. We discuss in Section 4 whether the other
features at similar velocities can be related between epochs.

###  3.2 Methanol masers

Our 6.7-GHz masers observations represent the first astrometric measurements,
since the positions quoted by Beuther et al. (2002a) are derived from Minier
et al. (2000) and are actually for the 12 GHz masers. We detected 20 6.7-GHz
methanol maser features over a velocity range of 16.3 – 27.6 km s-1, with a
complex spatial distribution. Figure 3 (upper) illustrates the spectra of the
6.7-GHz masers, where the brightest feature (1) $\sim$50 Jy has the maximum
velocity $\sim$27 km s-1, while the faintest peak is $\sim$0.388 Jy at
velocity 20.46 km s-1. Fig. 3 (lower) shows that our 6.7-GHz methanol masers
are extended over $\sim 1\aas@@fstack{\prime\prime}3$, with the greatest
elongation E-W. The brightest emission, features 1 and 18, are separated by
$\sim$$0.^{\prime\prime}8$ at a position angle of $\sim$60∘ (from N through
E), with the most red-shifted emission to the East. The weaker, mostly blue-
shifted methanol emission covers a range of position angles extending to
$\sim$110∘. These orientations are close to the position angles of the
outflows identified as South C and A by Beuther et al. (2003).

Weaker emission at intermediate velocities is more scattered, mostly but not
exclusively associated with the locations of the peaks. The components within
a 50 mas region around $\alpha$(J2000) = 19h 43m 11s.25, $\delta$(J2000) = 23∘
44′ $03\aas@@fstack{\prime\prime}3$ cover the greatest localized methanol
maser velocity span, of 19–27.5 km s-1.

The parameters of the 6.7 GHz CH3OH maser spectral features are listed in
Table 4, while Figure 4 (a and b) represents the distribution of the methanol
and water masers components throughout the region, color-coded by velocity.

Table 4: The properties of the detected CH3OH maser features toward IRAS 19410+2336 with position errors. The position leading terms are $\alpha$(J2000) = 19h 43m, $\delta$(J2000) = +23∘ 44′. No. |  Vel.(km s-1) | Flux(Jy) |  $\alpha$ (s) | Error(ms) |  $\delta$ (′′) | Error(mas)
---|---|---|---|---|---|---
1 | 26.92 |  49.183 $\pm$ 0.079 | 11.2464 | 0.007 | 3.298 | 0.1
2 | 24.46 |  10.165 $\pm$ 0.070 | 11.2500 | 0.007 | 3.288 | 0.1
3 | 23.67 |  5.763 $\pm$ 0.066 | 11.2500 | 0.014 | 3.296 | 0.2
4 | 21.65 |  0.514 $\pm$ 0.077 | 11.2524 | 0.211 | 3.113 | 3.0
5 | 21.43 |  3.257 $\pm$ 0.066 | 11.1958 | 0.021 | 2.916 | 0.3
6 | 20.60 |  1.221 $\pm$ 0.078 | 11.1958 | 0.087 | 2.934 | 1.3
7 | 20.07 |  9.683 $\pm$ 0.075 | 11.1956 | 0.007 | 2.942 | 0.1
8 | 19.67 |  19.902 $\pm$ 0.084 | 11.1950 | 0.007 | 2.951 | 0.1
9 | 21.65 |  0.496 $\pm$ 0.090 | 11.2533 | 0.335 | 3.312 | 4.1
10 | 20.99 |  0.494 $\pm$ 0.103 | 11.2716 | 0.364 | 3.263 | 6.7
11 | 20.46 |  0.388 $\pm$ 0.060 | 11.2519 | 0.160 | 3.328 | 1.5
12 | 20.42 |  0.555 $\pm$ 0.077 | 11.2514 | 0.211 | 3.283 | 2.9
13 | 20.07 |  2.880 $\pm$ 0.078 | 11.2744 | 0.029 | 3.398 | 0.4
14 | 19.72 |  2.949 $\pm$ 0.087 | 11.1874 | 0.036 | 3.138 | 0.5
15 | 19.23 |  1.033 $\pm$ 0.122 | 11.2477 | 0.233 | 3.294 | 4.0
16 | 18.88 |  4.842 $\pm$ 0.069 | 11.2714 | 0.014 | 3.243 | 0.2
17 | 18.97 |  9.960 $\pm$ 0.076 | 11.2688 | 0.007 | 3.218 | 0.1
18 | 16.95 |  30.107 $\pm$ 0.106 | 11.1960 | 0.007 | 2.943 | 0.1
19 | 16.82 |  9.814 $\pm$ 0.148 | 11.2672 | 0.021 | 3.250 | 0.5
20 | 17.08 |  19.015 $\pm$ 0.111 | 11.2673 | 0.007 | 3.208 | 0.1

Figure 2: The integrated spectra of the 6.7 GHz CH3OH and 22 GHz water masers
for both epochs 2005 and 2019

.

Figure 3: The upper panel represents the spectra of 6.7-GHz methanol maser
emission, the numbers refer to the features listed in Table 4. Numbers in
black, dark and mid grey correspond to full, dotted and dashed lines. The
lower panel shows the distribution of the CH3OH maser components toward the
IRAS 19410+2336 region, the size of the symbol is proportional to the log of
their intensities.

It also shows that the methanol maser spots are spread over $\sim
1\aas@@fstack{\prime\prime}3$ corresponding to $\sim$2808 au at a distance
2.16 kpc. The radial velocities were found to be in the range of 16.38 km s-1
to 27.53 km s-1.

The spatial correlation between water and methanol is closest for the West
part of the blue-shifted methanol masers, see Figure 4(b), but the closest
angular separation is $>0\aas@@fstack{\prime\prime}1$. The red-shifted
methanol maser components, Figure 4(a), are offset by
$>0\aas@@fstack{\prime\prime}5$ from the water masers, with more blue-shifted
methanol still further East. The flux densities for the methanol peaks range
from 0.514 Jy to 49.183 Jy, where the brightest maser is the most red-shifted,
located in between the more blue-shifted clumps.

Figure 4: The distribution of the 6.7 GHz CH3OH(crosses) and H2O maser
components denoted by plus signs and circles for the 2005 and 2019 epochs
respectively toward the IRAS 19410+2336 region. The size of the symbol is
proportional to the log of the maser component intensity, colour-coded by the
LSR velocities. Panels (a) and (b) represent the Eastern and Western regions
where masers were imaged, respectively.

##  4 Discussion

The IRAS 19410+2336 was previously observed by Caswell et al. (1995) where 6.7
GHz and 12 GHz masers were detected in velocity ranges 14 to 27 km s-1 and 16
to 28 km s-1 respectively, while the peak flux densities were found to be 42
Jy and 15.8 Jy, respectively. Szymczak et al. (2000), using single dish
observations, found that the strongest feature of 34 Jy occurred at 17.3 km
s-1 and the spectrum had changed considerably compared with the previous
observations (see Figure 1 in Szymczak et al. 2000).

Figure 2 shows the integrated spectra of 6.7-GHz methanol and 22-GHz water
masers. The water masers appears to cover a wider range in LSR velocity than
the methanol but the methanol peaks are brighter (during our observations).
Figure 6, right, plots our maser observations, denoting methanol by crosses
and water masers by plus signs (2005) and circles (2019). The size of methanol
and water feature symbols is proportional to the log of their intensities,
color-coded by velocity. Methanol masers at 12 GHz are shown by triangles
(Minier et al., 2000). The c3mm-sII and 13-s mm continuum peaks are shown by
solid stars (position accuracies
$0\aas@@fstack{\prime\prime}5,0\aas@@fstack{\prime\prime}2$) (Rodón et al.,
2012). The mm1 1.2-mm continuum peak is shown by an unfilled star, (position
accuracy 0$\aas@@fstack{\prime\prime}$5) Beuther et al. (2003). The peak of
the thermal SiO clump is shown by a square (measured in a beam size of 4”)
(Widmann et al., 2016). This shows that the 12-GHz masers observed by Minier
et al. (2000), have an angular separation of $\sim$
$0\aas@@fstack{\prime\prime}7$ and appear to be spatially correlated (within
astrometric uncertainties) with the brightest 6.7 GHz peaks at similar
velocities. It is likely that they are tracing the same phenomena around the
embedded protostellar object/s (see later).

Turning to the 22 GHz water masers, Figure 1 (b) shows that the 2019
components have a smaller separation and are rotated with respect to the 2005
data. The displacement of each feature is about 70 mas or 150 au, which would
imply a proper motion speed of 36 km s-1 if the clumps are physically the same
at both epochs. This seems unlikely for two reasons. Firstly, although it
could signify infall, or rotation of a NE-SW disc, these mechanisms would have
to be very non-radial or warped. Secondly, the speed is far greater than the
$V_{\mathrm{lsr}}$ span, and would imply an unreasonably large enclosed mass,
see discussion in Section 4.2 para 4.

Beuther et al. (2002a) used 1999 VLA observations to measure the positions of
the H2O maser peaks. The peaks at $V_{\mathrm{LSR}}$ 27 and 17 km s-1 were
located at $\alpha$(J2000) = 19h 43m 11.2s, $\delta$(J2000) = +23∘ 44′
$03\aas@@fstack{\prime\prime}0$ and $\alpha$(J2000) = 19h 43m 11.2s,
$\delta$(J2000) = +23∘ 44′ $03\aas@@fstack{\prime\prime}1$, respectively, with
an absolute astrometric error estimated at $1^{\prime\prime}$ and a relative
accuracy of $0\aas@@fstack{\prime\prime}1$. This gives a separation of (0,
100) $\pm$ 100 mas, with the more blue-shifted peak to the north, as in our
data. The mid point of these peaks is at (11.2s,
$03\aas@@fstack{\prime\prime}05$) in the coordinates of Figure 5. This is
offset by (220, 44) mas from our 2019 mid point at (11s.1840
$2\aas@@fstack{\prime\prime}961$), which is well within the astrometric
uncertainty of the VLA data. It seems reasonable to align the mid-point of
Beuther et al. (2002a)’s data with our observations, as in Figures 5 and 6.

The 22 GHz maser toward IRAS 19410+2336 was first resolved by Lada et al.
(1981), as previously mentioned. Their position (in J2000) is over 10 arcsec
from recent measurements, probably due to the difficulties in determing
astrometric positions from fringe-rate analysis in 1981\. The maser velocities
cover a range overlapping our 2005 observations, at 12.0, 18.7 and 20.6 km
s-1. The 12 km s-1 feature was brightest, at $\sim$21 Jy, the others, of a few
Jy, being offset from this by (–94, –40) and (–64, –20) mas, uncertainties
$\leq$10 mas. We aligned the average position of these features with the
average position of the features we detected at similar velocities in 2019.

Comparing the observations of Lada et al. (1981) and Beuther et al. (2002a)
with our data shows that emission around 18–20 km s-1 was detected at all
epochs. The most blue-shifted feature, around 13 km s-1, was detected in 1979
and 2005 but not in 1999 nor 2019, and the red-shifted features were detected
in 1999, 2005 and 2019. We note that 22 GHz single-dish observations by
Sridharan et al. (2002) found a peak of 110 Jy, and single-dish observations
by Xi et al. (2015) detected emission between 3–33 km s-1 with a peak $\sim$28
Jy but in both cases it is possible that emission from other parts of the
region is in the beam. If the features at similar velocities are from the same
gas at different epochs (which seems unlikely, see Section 4.2), their flux
densities vary roughly 2-fold between epochs; if they are from different parts
of IRAS 19410+2336 then its water maser variability is even more pronounced.
The flux variability could be caused by the variability of infrared sources
associated with this region.

Our 2005 and 2019 astrometric observations, with accuracies of 12 and 10 mas
respectively, locate most of the 22-GHz H2O masers firmly on the Western edge
of the 6.7 and 12-GHz methanol maser distribution.

Figure 5:  Comparison between previous and current water maser imaging towards
IRAS 19410+2336. The X are for Lada et al. (1981) VLBI measurements, the
squares are for Beuther et al. (2002a) measurements, while the crosses (+) and
the circles are for our MERLIN and eMERLIN measurements respectively. The
epochs are aligned as outlined in Sections 3.1 and 4.

In conclusion we can claim that for Lada et al. (1981) observations the
astrometry was uncertain, being based on only a single baseline, so it is
possible that the emission detected by Lada et al. (1981) and by Beuther et
al. (2002a) as well as ours came from a similar region. The flux variability
could be caused by variable infrared source that associated with this region.

Class I methanol masers at 95 and 44 GHz have been detected towards this
region by Val’tts et al. (1995) and Litovchenko et al. (2011). The peak flux
density at 95 GHz was estimated to be 3.9 Jy at a velocity of 22.59 km s-1,
while Litovchenko et al. (2011) reported three maser features with the
brightest feature at 17 Jy at a peak velocity of 23.6 km s-1. The presence of
Class I methanol masers associated with shocked material in protostellar
outflows (see Beuther et al. 2003) in addition to Class II methanol masers is
a further evidence for tracing very early stages of massive star formation
(Cyganowski et al., 2009) and Leurini et al. (2016, and references therein).

The OH maser was first reported by Braz & Sivagnanam (1987) with flux density
1.1 Jy and corresponding radial velocity 20.9 km s-1. Edris et al. (2007),
using the Greenbank single dish telescope, detected peaks at flux densities of
1.7 Jy and 0.7 Jy in velocity ranges of 19 to 22 km s-1 and 20 to 23 km s-1
respectively, showing modest variability. The position was $\alpha$(J2000) =
19h 43m 12s.2 $\pm 1^{s}.1$, $\delta$(J2000) = 23∘ 44′ $3^{\prime\prime}$ $\pm
12^{\prime\prime}$ As noted in Section 1, the presence of OH masers implies a
later evolutionary stage than does methanol class II and H2O without OH. This
suggests that different prototars in this region are at different evolutionary
stages, but high-resolution OH imaging is needed to confirm which masers are
associated with which objects.

Figure 6 (right) shows the positions of other emission detected in this
region.The thermal SiO (Widmann et al., 2016) was observed with a 4-arcsec
beam, larger than the total maser extent. Millimetre source mm1, position
accuracy $0.^{\prime\prime}5$ (Beuther et al., 2003) is located close to the
center of the methanol maser distribution whilst 13-s, position accuracy
$0.^{\prime\prime}13$ (Rodón et al., 2012) is close to the water masers. We
suggest that the Eastern methanol masers (Figure 3 (a)) are associated with
thermal SiO and a less evolved protostar than the Western methanol and water
masers. This is explained in the following Sections 4.1 and 4.2.

###  4.1 Outflow activity

Outflow activity in high-mass star forming regions is predicted by different
massive stars formation scenarios. The competitive accretion model suggests
that outflows in high mass star forming regions should be less collimated than
those produced by low mass star forming regions (Tan et al., 2014).

On the other hand, the core accretion model, which is considered as a scaled-
up version of low-mass star formation, suggests well-collimated outflows
similar to those associated with low mass protostars (Widmann et al., 2016,
and references therein). According to Arce et al. (2007) the opening angles
for well-collimated outflows in massive star forming regions are between 25∘
and 30∘.

The outflow activity in IRAS 19410+2336 has been studied extensively by
Beuther et al. (2003); Martín-Hernández et al. (2008) and Widmann et al.
(2016). Beuther et al. (2003) observed CO (J = 1-0) in addition to the 2.6 mm
continuum using the IRAM 30 m telescope and the Plateau de Bure
interferometer, and the H2 line at 2.12 $\mu$m, using the Calar Alto 3.5 m
telescope. These results support the presence of well-collimated outflows
similar to those found in low-mass star formation. They concluded that the
region contains a strong outflow activity with at least seven outflows. They
also report that the kinematics of the individual outflows resemble those in
low mass-star forming regions.

The line widths of the SiO detected by Widmann et al. (2016) in the region
(typically 18 km s-1) along with a Gaussian profile, but with an additional
signal in the line wings, strongly suggests that SiO traces dynamical gas in
jets and/or outflows (Widmann et al., 2016).

The methanol masers show the most red-shifted emission near the center of the
distribution with more blue-shifted emission extending to the SE and SW. The
peak position of SiO, shown in Fig. 6 was measured with a $4^{\prime\prime}$
beam, is thought to trace an outflow. Its location points to a close
association with the Eastern methanol maser features. It is therefore likely
that these maser features probe the same outflow as the SiO. The mid-point of
the red- and blue-shifted methanol masers is centred at about $\alpha$(J2000)
= 19h 43m 11.26s, $\delta$(J2000) = +23∘ 44′ $03\aas@@fstack{\prime\prime}3$).
This is within the $0\aas@@fstack{\prime\prime}5$ position uncertainty of the
mm1 peak (Beuther et al., 2002b). Beuther et al. (2003) suggest that mm1
produces multiple outflows. The maser distribution could trace ongoing outflow
activity, in a similar direction to H2 feature (2), associated with outflow C
by Beuther et al. (2003) (Figure 6). The distribution of the masers with
respect to their mid-point suggests outflow opening angles of about 40∘ and
20∘. This results agrees with the well-collimated outflow criteria in Arce et
al. (2007).

de Villiers et al. (2015) presented an evolutionary sequence based on the
association of 6.7-GHz methanol maser with outflows in the high-mass star
forming region. This sequence suggested that: Firstly, the outflow appears and
grows without any presence of methanol maser. Then, both outflow and methanol
masers are present and could be observed. Later, methanol maser emission
disappears but the outflow remains detectable and the UC H II region appears.
Eventually, the outflow switches off while the UC H II region remains.

Although de Villiers et al. (2015) didn’t test the relation between water
masers and class II methanol masers in the presence of outflows, Ellingsen et
al. (2007) concluded that methanol class II masers in high mass star forming
regions are tracing an earlier evolutionary phase than water masers.

Figure 6:  Right hand side: Species in IRAS 19410+2336 detected close to the
masers imaged by (e-)MERLIN. Methanol masers (crosses) and water masers (plus
signs and circles) (this work). The size of methanol and water feature symbols
is proportional to the log of their intensities, color-coded by velocity. 12
GHz methanol masers (triangles) (Minier et al., 2000). c3mm-sII and 13-s mm
continuum peaks (solid stars, position accuracies
$0\aas@@fstack{\prime\prime}5,0\aas@@fstack{\prime\prime}2$) (Rodón et al.,
2012). mm1 1.2-mm continuum peak (unfilled star, position accuracy
0$\aas@@fstack{\prime\prime}$5) (Beuther et al., 2003). Peak of the $\sim
4^{\prime\prime}$ thermal SiO clump (square) (Widmann et al., 2016). The left
hand side represents CO (1-0) data from Beuther et al. (2003) which appears in
the background of the map while the ellipsis and arrows represent the proposed
outflow directions from Beuther et al. (2003). The maser area is marked by the
green square.

Bearing in mind the scheme of de Villiers et al. (2015) and also the maser
based evolutionary sequence of Ellingsen et al. (2007), we suggest that the
embedded protostar associated with the Eastern methanol masers (Figure 6) is
in an early stage of evolution. These masers lie within the fitted extent of
the 3-mm source detected by Rodón et al. (2012). The closest 1.2-mm source to
the Eastern masers is the peak of mm1 detected by Beuther et al. (2002b). We
suggest that the Western methanol and water masers are associated with a more
evolved 1.2-mm source, 13-s, which is located close to the centre of the water
maser distribution (within the combined position errors of
$0\aas@@fstack{\prime\prime}2$) (Section 4.2). Rodón et al. (2012) estimated
the temperature of 13-s at 90 K, the hottest source in the region at about
twice the temperature assumed by Beuther et al. (2003) for mm1. 13-s is also
associated with a higher gas column density at 25 1023cm-2, at the location
where Martín-Hernández et al. (2008) and Qiu et al. (2008) detected a near-
infrared(NIR) and mid-infrared(MIR) source.

mm1 and 13-s are separated by $0\aas@@fstack{\prime\prime}8$, greater than
their combined, nominal astrometric errors of $0\aas@@fstack{\prime\prime}55$.
However, we note that mm1, and also the source c3mm-s11 detected at 3-mm
wavelength by Beuther et al. (2002b), were observed at lower resolution than
13-s. Each of mm1 and c3mm-s11 cover at least $\sim$5 arcsec, including
contributions from other sources outside the region containing the masers
reported here. If the source associated with the Eastern methanol masers is at
an early stage of evolution, this may explain why is too faint to be detected
as a distinct 1.2-mm source between the Eastern red- and blue-shifted methanol
features.

If the water masers we detect toward the IRAS 19410+2336 region were
associated with a bipolar outflow from the same origin, s-13, the jets must be
precessing or have an opening angle of $\sim$85∘, and the masers seen at
different epochs are in most cases either different ejecta, or interactions
with different parts of the surrounding medium. Thus, the outflow could be
characterized as poorly collimated. Arce et al. (2007) found that an opening
angle of $>50$∘ is usually associated with more evolved sources than those
producing well-collimated outflows. It is possible, therefore, that the water
masers are associated with an outflow from 13-s, poorly collimated due to age
or the competative accretion mechanism, but this is not consistent with other
observations of this region and we consider the alternative possibility that
they instead trace a disc.

###  4.2 Disc and infall

The presence of circumstellar gas discs around young stars was inferred to be
very common by Beckwith & Sargent (1993). The 6.7 GHz methanol maser commonly
traces the circumstellar gas disc or outflow of high-mass star-forming regions
(Fujisawa et al., 2014, and references therein). Observations by Norris et al.
(1998) added support to the hypothesis that the 6.7 GHz methanol masers lie
within rotating edge-on circumstellar discs.

The combined distribution of the methanol masers we mapped in this region of
IRAS 19410+2336 does not appear to be a disc, as the most red-shifted emission
has blue-shifted emission both to the East and to the West. As explained in
Sec. 4.1, the Eastern masers may trace an outflow from mm1. Bartkiewicz et al.
(2009) found that at high resolution, 6.7 GHz methanol masers trace a variety
of morphologies including infall, and it is possible that the Western masers
close to 13-s may be a sign of infall.

The closest methanol and water masers are clumps 14 and d, separated by 96
mas, equivalent to 205 au, 0.001 pc. Their velocities only are 0.17 km s-1
apart. Although this is much closer than the separations measured at lower
resolution by Breen et al. (2018) and Beuther et al. (2002a), we note that, as
discussed by Beuther et al. (2002a), the different pumping mechanisms mean
that they are likely to trace different protostellar phenomena (e.g
Bartkiewicz et al., 2012; Beuther et al., 2002a; Sanna et al., 2010), although
these can be associated with the same exciting source. The methanol masers are
separated from the water masers by more than half the total size of the water
maser area, the combined astrometric error is much smaller (16 mas) and the
physical separation may be greater than 205 au if the masing clumps are at
different distances along the line of sight. Nonetheless, the masers have
appeared at different locations at different epochs. Given the variability of
the masers and the excitation of emission at different locations at different
epochs, as well as proper motions, it is possible that methanol and water
masers could appear closer in future observations.

Section 4.1 explains that, whilst it is possible that the water masers trace
an outflow, this is not consistent with other studies which suggest that an
outflow from the central, high-mass YSO would be more tightly collimated.
Although the 22-GHz water maser is not widely known to trace discs or rotating
material around protostellar objects, this is not unknown, e.g. Imai et al.
(2006) found strong evidence for water masers tracing rotating material (i.e.
a disc) through monitoring the high-mass star forming region G192.16-3.84.
Figure 5 illustrates all published H2O masers imaged towards the IRAS
19410+2336 region, which lead us to suggest that they are tracing rotating
gas. At all epochs, these masers show a consistent offset between more blue-
shifted emission to the North and red-shifted to the South (Figure 5). The
maximum angular extent is equivalent to $\sim$350 au at 2.16 kpc, and the mid
point is close to the mm source 13-s (Rodón et al., 2012), within its
astrometric uncertainty. If the position change of the MERLIN – e-MERLIN maser
features is due to physical rotation, this corresponds to $V_{\mathrm{pm}}$ =
36 km s-1 in a tilted disc of approximately 175 au radius. The
$V_{\mathrm{lsr}}$ velocities span $\sim$15 km -1, or a line-of-sight velocity
of magnitude 7.5 km s-1 with respect to the mid-point, giving a total orbital
velocity $\sim$37 km -1. However, this would correspond to an enclosed mass
$\sim$270 M⊙, which is far greater than the mass estimate for 13-s of 8.1 M⊙
by Rodón et al. (2012). Moreover, although the Beuther et al. (2002a) maser
positions have large astrometric errors, the position angle between the red-
and blue-shifted emission has a smaller uncertainty, and is inconsistent with
rotation in a consistent direction. This makes it more likely that the water
maser position changes are due to different parts of the disc being excited at
the different epochs.

Thus, if we assume that the maser positions observed in 2005 represent the
maximum extent of the disc, and deduce from the distribution of masers at all
epochs that the angle of inclination with respect to the plane of the sky is
$\sim 60^{\circ}$, then the rotation velocity is $\sim$15 km s-1. This
corresponds to a more feasible enclosed mass of $\sim$44 M⊙. The main
uncertainties are due to the extent and nature of the disc, which could be
warped or flared, and the possibility of infall or outflow as well as
rotation. If the disc is close to edge-on and the apparent E–W extensions are
due to flaring or ablation, the enclosed mass could be as low as 11 M⊙, close
to the estimate by Rodón et al. (2012).

##  5 Conclusions

Methanol class II 6.7-GHz and water 22-GHz masers were observed toward IRAS
19410+2336 region at high angular resolution using MERLIN and e-MERLIN. The
astrometric uncertainty is 12–10 mas, the most accurate measurements yet made
in this region.

The CH3OH maser features span a velocity range of about 16 to 27 km s-1 with a
peak flux density $\sim$50 Jy. They are spread over a region $\sim
1\aas@@fstack{\prime\prime}3$ corresponding to $\sim$2808 au at a distance
2.16 kpc (Figure 4). They are separated into two groups
$0\aas@@fstack{\prime\prime}8$ apart, which might indicate that our methanol
masers are associated with more than one protostar in this region. The Eastern
group of CH3OH masers have an elongated distribution with a red-blue offset
similar to the direction of ‘C’ (see arrows in Figure 6) proposed previously
by Beuther et al. (2003). The apparent centre of the Eastern methanol maser
outflow is at about $\alpha$(J2000) = 19h 43m 11.26s $\delta$(J2000) = +23∘
44′ $03\aas@@fstack{\prime\prime}3$ which is within the mm1 peak position
uncertainty of $\sim 0\aas@@fstack{\prime\prime}5$.

The H2O masers, observed in 2005 and 2019, span a velocity range from 13 to
$\sim$29 km s-1. The peak flux density was found to be 18.7 Jy and 13.487 in
2005, and 2019, respectively. They are distributed over 165 mas ($\sim$350 au
at 2.16 kpc). The H2O masers may be associated with rotating material (a
tilted disc) around a protostellar object centered at $\alpha$(J2000) = 19h
43m 11.184s $\delta$(J2000) = +23∘ 44′ $2\aas@@fstack{\prime\prime}96$ which
coincides with the mm source 13-s within the combined position accuracy of
$\sim$$0.^{\prime\prime}2$ (Rodón et al., 2012). The estimated mass for this
protostar is 11–44 M⊙, the lower limit being close to the estimate by Rodón et
al. (2012). The Western methanol masers are all slightly blue-shifted, and if
they are also associated with 13-s it is possible that they are in infall.

This scenario is consistent with the multiple outflows identified by Beuther
et al. (2003) in the IRAS 19410+2336 region, and with the development of
multiple protostars. Our findings also imply that the Eastern masers, where no
H2O is detected, are associated with a less evolved protostellar object than
the Western masers. The angular separation between the two putative
protostellar objects is $\sim$$0.^{\prime\prime}36$, while the physical
separation at 2.61 kpc would be at least $\sim$780 au.

Eventually, high angular resolution observations of other maser species in the
region are needed to investigate which phenomena they trace, their association
with the embedded protostars and provide clues to their evolutionary stages in
the massive star forming region IRAS 19410+2336. Observations of dust and
thermal lines at higher resolution, ideally allowing a position and
astrometric accuracy of 100 mas or better, are needed to confirm the
associations between discs or outflows and the protostars (as in the ALMA
observations by Sanna et al. (2019) of G023.01-00.41). ALMA – or NOEMA, given
the Declination of IRAS 19410+2336 – could make this possible.

## acknowledgements

We warmly thank Prof. R. Battye, M. Gray and R. Beswick, and the rest of the
e-MERLIN team for guidance in reducing these data. We also remember the
important role of the late Dr Jim Cohen in the initiation of this project. We
thank the referee for very helpful comments which have improved this paper.
e-MERLIN is the UK radio interferometer array, operated by the University of
Manchester on behalf of STFC. We acknowledge the use of MERLIN archival data
as well as NASA’s Astrophysics Data System Service. M.Darwish would like to
acknowledge the Science and Technology Development Fund (STDF) N5217, Academy
of Scientific Research and Technology (ASRT), Cairo, Egypt and Kottamia Center
of Scientific Excellence for Astronomy and Space Sciences (KCSEASSc), National
Research Institute of Astronomy and Geophysics. Our sincere thanks to H.
Beuther and F. Widmann for their helpful discussion.

## References

  * Arce et al. (2007) Arce H. G., Shepherd D., Gueth F., Lee C.-F., Bachiller R., Rosen A., Beuther H., 2007, Protostars and Planets V, pp 245–260
  * Baars et al. (1977) Baars J. W. M., Genzel R., Pauliny-Toth I. I. K., Witzel A., 1977, A&A, 61, 99
  * Bartkiewicz et al. (2009) Bartkiewicz A., Szymczak M., van Langevelde H. J., Richards A. M. S., Pihlström Y. M., 2009, A&A, 502, 155
  * Bartkiewicz et al. (2012) Bartkiewicz A., Szymczak M., van Langevelde H. J., 2012, A&A, 541, A72
  * Beckwith & Sargent (1993) Beckwith S. V. W., Sargent A. I., 1993, in Levy E. H., Lunine J. I., eds, Protostars and Planets III. pp 521–541
  * Beuther et al. (2002a) Beuther H., Walsh A., Schilke P., Sridharan T. K., Menten K. M., Wyrowski F., 2002a, A&A, 390, 289
  * Beuther et al. (2002b) Beuther H., Schilke P., Menten K. M., Motte F., Sridharan T. K., Wyrowski F., 2002b, ApJ, 566, 945
  * Beuther et al. (2003) Beuther H., Schilke P., Stanke T., 2003, A&A, 408, 601
  * Beuther et al. (2007) Beuther H., Churchwell E. B., McKee C. F., Tan J. C., 2007, Protostars and Planets V, pp 165–180
  * Braz & Sivagnanam (1987) Braz M. A., Sivagnanam P., 1987, A&A, 181, 19
  * Breen et al. (2010a) Breen S. L., Ellingsen S. P., Caswell J. L., Lewis B. E., 2010a, MNRAS, 401, 2219
  * Breen et al. (2010b) Breen S. L., Caswell J. L., Ellingsen S. P., Phillips C. J., 2010b, MNRAS, 406, 1487
  * Breen et al. (2018) Breen S. L., et al., 2018, MNRAS, 474, 3898
  * Caswell (1996) Caswell J. L., 1996, MNRAS, 279, 79
  * Caswell (1997) Caswell J. L., 1997, MNRAS, 289, 203
  * Caswell et al. (1995) Caswell J. L., Vaile R. A., Ellingsen S. P., Whiteoak J. B., Norris R. P., 1995, MNRAS, 272, 96
  * Cohen (1984) Cohen R. J., 1984, The Observatory, 104, 125
  * Cragg et al. (2002) Cragg D. M., Sobolev A. M., Godfrey P. D., 2002, MNRAS, 331, 521
  * Cyganowski et al. (2009) Cyganowski C. J., Brogan C. L., Hunter T. R., Churchwell E., 2009, ApJ, 702, 1615
  * Diamond et al. (2003) Diamond P. J., Garrington S. T., Gunn A. G., Leahy J. P., McDonald A., Muxlow T. W. B., Richards A. M. S., Thomasson P., 2003, Technical report, The MERLIN User Guide. Jodrell Bank Observatory, UK., http://www.merlin.ac.uk/user_guide/MUG.ps.gz
  * Edris et al. (2005) Edris K. A., Fuller G. A., Cohen R. J., Etoka S., 2005, A&A, 434, 213
  * Edris et al. (2007) Edris K. A., Fuller G. A., Cohen R. J., 2007, A&A, 465, 865
  * Ellingsen et al. (2007) Ellingsen S. P., Voronkov M. A., Cragg D. M., Sobolev A. M., Breen S. L., Godfrey P. D., 2007, in Chapman J. M., Baan W. A., eds, IAU Symposium Vol. 242, Astrophysical Masers and their Environments. pp 213–217 (arXiv:0705.2906), doi:10.1017/S1743921307012999
  * Etoka et al. (2005) Etoka S., Cohen R. J., Gray M. D., 2005, MNRAS, 360, 1162
  * Etoka et al. (2012) Etoka S., Gray M. D., Fuller G. A., 2012, MNRAS, 423, 647
  * Fujisawa et al. (2014) Fujisawa K., et al., 2014, PASJ, 66, 31
  * Hartquist et al. (1995) Hartquist T. W., Menten K. M., Lepp S., Dalgarno A., 1995, MNRAS, 272, 184
  * Hirota et al. (2007) Hirota T., et al., 2007, PASJ, 59, 897
  * Imai et al. (2006) Imai H., Omodaka T., Hirota T., Umemoto T., Sorai K., Kondo T., 2006, PASJ, 58, 883
  * Keto & Zhang (2010) Keto E., Zhang Q., 2010, MNRAS, 406, 102
  * Krumholz et al. (2009) Krumholz M. R., Klein R. I., McKee C. F., Offner S. S. R., Cunningham A. J., 2009, Science, 323, 754
  * Kurtz et al. (1994) Kurtz S., Churchwell E., Wood D. O. S., Myers P., 1994, in American Astronomical Society Meeting Abstracts #184. p. 907
  * Lada et al. (1981) Lada C. J., Blitz L., Reid M. J., Moran J. M., 1981, ApJ, 243, 769
  * Leurini et al. (2016) Leurini S., Menten K. M., Walmsley C. M., 2016, A&A, 592, A31
  * Litovchenko et al. (2011) Litovchenko I. D., Alakoz A. V., Val’Tts I. E., Larionov G. M., 2011, Astronomy Reports, 55, 1086
  * MacLeod et al. (1993) MacLeod G. C., Gaylard M. J., Kemball A. J., 1993, MNRAS, 262, 343
  * Martín-Hernández et al. (2008) Martín-Hernández N. L., Bik A., Puga E., Nürnberger D. E. A., Bronfman L., 2008, A&A, 489, 229
  * McCall et al. (1990) McCall M. L., Richer M. G., Visvanathan N., 1990, ApJ, 357, 502
  * Menten (1991) Menten K. M., 1991, ApJ, 380, L75
  * Minier et al. (2000) Minier V., Booth R. S., Conway J. E., 2000, A&A, 362, 1093
  * Norris et al. (1998) Norris R. P., et al., 1998, ApJ, 508, 275
  * Qiu et al. (2008) Qiu K., et al., 2008, ApJ, 685, 1005
  * Reid et al. (2009) Reid M. J., Menten K. M., Zheng X. W., Brunthaler A., Xu Y., 2009, ApJ, 705, 1548
  * Rodón et al. (2012) Rodón J. A., Beuther H., Schilke P., 2012, A&A, 545, A51
  * Sanna et al. (2010) Sanna A., Moscadelli L., Cesaroni R., Tarchi A., Furuya R. S., Goddi C., 2010, A&A, 517, A71
  * Sanna et al. (2019) Sanna A., et al., 2019, A&A, 623, A77
  * Sridharan et al. (2002) Sridharan T. K., Beuther H., Schilke P., Menten K. M., Wyrowski F., 2002, ApJ, 566, 931
  * Stahler et al. (2000) Stahler S. W., Palla F., Ho P. T. P., 2000, Protostars and Planets IV, pp 327–352
  * Sugiyama et al. (2011) Sugiyama K., et al., 2011, PASJ, 63, 53
  * Szymczak & Gérard (2004) Szymczak M., Gérard E., 2004, A&A, 414, 235
  * Szymczak et al. (2000) Szymczak M., Hrynek G., Kus A. J., 2000, A&AS, 143, 269
  * Tan & McKee (2004) Tan J. C., McKee C. F., 2004, in Lamers H. J. G. L. M., Smith L. J., Nota A., eds, Astronomical Society of the Pacific Conference Series Vol. 322, The Formation and Evolution of Massive Young Star Clusters. p. 263 (arXiv:astro-ph/0403498)
  * Tan et al. (2014) Tan J. C., Beltrán M. T., Caselli P., Fontani F., Fuente A., Krumholz M. R., McKee C. F., Stolte A., 2014, Protostars and Planets VI, pp 149–172
  * Torrelles et al. (2012) Torrelles J. M., Patel N. A., Curiel S., Gómez J. F., Anglada G., Estalella R., 2012, arXiv e-prints, p. arXiv:1204.3808
  * Val’tts et al. (1995) Val’tts I. E., Dzyura A. M., Kalenskii S. V., Slysh V. I., Bus R., Vinnberg A., 1995, Astronomy Reports, 39, 18
  * Widmann et al. (2016) Widmann F., Beuther H., Schilke P., Stanke T., 2016, A&A, 589, A29
  * Wolfire & Cassinelli (1987) Wolfire M. G., Cassinelli J. P., 1987, ApJ, 319, 850
  * Xi et al. (2015) Xi H., Zhou J., Esimbek J., Wu G., He Y., Ji W., Tang X., 2015, MNRAS, 453, 4203
  * Xu et al. (2009) Xu Y., Voronkov M. A., Pandian J. D., Li J. J., Sobolev A. M., Brunthaler A., Ritter B., Menten K. M., 2009, A&A, 507, 1117
  * Zinnecker & Yorke (2007) Zinnecker H., Yorke H. W., 2007, ARA&A, 45, 481
  * de Villiers et al. (2015) de Villiers H. M., et al., 2015, MNRAS, 449, 119