# Near-Infrared Polarimetry and H2 emission toward Massive Young Stars:
Discovery of a Bipolar Outflow associated to S235 e2s3

R. Devaraj Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin
D02XF86, Ireland  A. Caratti o Garatti Istituto Nazionale di Astrofisica, 16
Salita Moiariello, Naples 80131, Italy  Dublin Institute for Advanced Studies,
31 Fitzwilliam Place, Dublin D02XF86, Ireland  L. K. Dewangan Physical
Research Laboratory, Navrangpura, Ahmedabad, Gujarat 380009, India  R.
Fedriani Department of Space, Earth and Environment, Chalmers University of
Technology, 412 96 Gothenburg, Sweden  T. P. Ray Dublin Institute for Advanced
Studies, 31 Fitzwilliam Place, Dublin D02XF86, Ireland  A. Luna Instituto
Nacional de Astrofísica, Óptica y Electrónica, Tonantzintla, Puebla 72840,
México

(Received 2022 September 23; Revised 2023 January 24; Accepted 2023 January
25)

###### Abstract

We present a near-infrared $H$ band polarimetric study toward the S235 e2s3
protostar, obtained using the POLICAN instrument on the 2.1m OAGH telescope.
The images reveal a bipolar outflow with a total length of about 0.5 pc. The
outflow nebulosity presents a high degree of linear polarization ($\sim$80%)
and reveals a centrosymmetric pattern with the polarization position angles.
The polarization characteristics suggest their origin to be single scattering
associated with dust in the outflow. Using multiwavelength archival data, we
performed spectral energy distribution (SED) fitting based on radiative
transfer models of turbulent core accretion theory. The best-fit SED model
indicated that the protostar has a mass of $6.8\pm 1.2$ $M_{\odot}$, with a
disk accretion rate of $3.6\pm 1.2\times 10^{-4}\,\mathrm{M_{\odot}\,yr^{-1}}$
and a total bolometric luminosity of $9.63\pm 2.1\times 10^{3}$ $L_{\odot}$.
Narrowband H2 ($2.12\,\mu$m) observations show shocked emission along the
bipolar lobes tracing the jet’s interaction with the surrounding medium. The
estimated H2 luminosity of the outflow is
$2.3_{-1.3}^{+3.5}\,\hbox{$L_{\odot}$}$, which matched the known power-law
correlation with the source bolometric luminosity, similar to other high-mass
outflows. The orientation of the bipolar outflow was found to be parallel to
the local magnetic field direction. The overall results assert the fact that
the S235 e2s3 source is a massive young star driving a highly collimated
bipolar outflow through disk accretion.

Uniﬁed Astronomy Thesaurus concepts: Massive stars (732); Polarimetry (1278);
Protostars (1302); Spectral energy distribution (2129); Star formation (1569);
Stellar jets (1607)

††journal: ApJ

##  1 Introduction

Figure 1: MIR three-color image of the S235 complex, constructed using Spitzer
MIPS and IRAC bands: $24\,\mu$m (red), $4.5\,\mu$m (green) and $3.6\,\mu$m
(blue) (Fazio et al., 2004). The left panel shows the larger S235 Main region
spanning $15\arcmin\times 14\arcmin$, corresponding to a physical scale of
$7.2\,\mathrm{pc}\times 6.7\,\mathrm{pc}$ at a distance of $1.65\,\rm{kpc}$.
The central ionizing O9.5V star for S235 Main is marked by the star symbol.
The right panel shows a zoomed-in view of the region centered on the S235 e2s3
source, spanning a $4\arcmin\times 4\arcmin$ FOV or a physical scale of about
$2\,\mathrm{pc}\times 2\,\mathrm{pc}$. The dashed white circle encompasses the
S235 e2s3 source.

Molecular outflows are one of the most prominent signposts of the star
formation process. They are intimately connected to accretion and ejection
phenomena of the embedded protostar (Pudritz & Norman , 1983; Ray et al.,
2007; Bally , 2016). The exact mechanism driving these outflows is unclear,
but according to models they can be due to a disk wind (Blandford & Payne ,
1982; Shu et al., 1994) or interactions between the disk field and
protostellar magnetosphere (Shang et al., 2007; Shu et al., 2007) or by
collimated flows from toroidal magnetic fields (Banerjee & Pudritz , 2006,
2007). A unifying factor in all these mechanisms is their dependence on
magnetic field interaction with the flows (Frank et al., 2014). Recent
observations suggest that the formation of intermediate- to high-mass stars
also proceeds via disk-mediated accretion as with their low-mass counterparts,
powering highly collimated outflows (Davis et al., 2004; Gredel et al., 2006;
Fedriani et al., 2018). The outflows hold insights about the associated
protostars and can reveal the nature of their origins through their
morphologies, kinematics, and physical properties. One way to the study the
interaction of these flows is through molecular H2 emission which is a
particularly good shock tracer as it is the primary coolant in the near-
infrared (NIR) and its emission can extend over large spatial scales (Caratti
o Garatti et al., 2006, 2008). Hence, observations of outflows with different
techniques are crucial in constraining the evolution of intermediate and high-
mass stars.

NIR imaging polarimetry offers a unique opportunity to study the morphology
and characteristics of outflows. The nebulosity associated to the outflows
presents high degrees of linear polarization that can only be attributed to
light scattering from dust particles (Werner et al., 1983; Bastien & Menard ,
1988). These observations can reveal information on the spatial distribution
of the scattering material and the inclination of the outflow (Tamura et al.,
1991; Jones et al., 2004; Meakin et al., 2005). The composition and size
distribution of the dust grains responsible for scattering can also be
investigated using polarimetry (e.g. Pendleton et al., 1990; Kim et al.,
1994). Generally, the young stars driving the most powerful outflows are
invisible at optical and NIR wavelengths as they are embedded in dust.
Polarimetric observations allow us to identify and locate the position of
these illuminating stars through their outflow polarization signatures (Tamura
et al., 1997; Weintraub et al., 2000).

The object of study in this work is ‘S235 e2s3’, which is located in the
eastern subregion of the S235 complex (see Figure 1). The S235 is a well
studied star-forming complex harboring an evolved H ii region: ‘S235 Main’
(Kirsanova et al., 2008; Dewangan & Anandarao, 2011; Dewangan et al., 2016).
The H ii region is ionized by a central massive star of O9.5V spectral type
(Georgelin et al., 1973). The region is estimated to be at a distance of
$1.65_{-0.10}^{+0.12}\,\rm kpc$ based on the Gaia DR2 (Brown et al., 2018)
parallax of the ionizing star. The molecular gas in the surrounding
environment has been traced with velocities ranging from
$-22\,\mathrm{kms^{-1}}$ to $-20\,\mathrm{kms^{-1}}$ and has been shown to
present a number of star forming clusters (Kirsanova et al., 2008, 2014;
Chavarría et al., 2014). Various studies have concluded that the clusters are
a result of positive stellar feedback driven by the expansion of the S235 Main
H ii region (Dewangan et al., 2016; Dewangan & Ojha , 2017; Devaraj et al.,
2021).

The source S235 e2s3 was identified by Dewangan & Anandarao (2011) through a
Spitzer IRAC (Fazio et al., 2004) mid-infrared (MIR) photometric study of the
S235 complex. They inferred it to be a class 0/I young stellar object (YSO)
lacking any NIR counterparts. Chavarría et al. (2014) used MIR color-color
diagrams and identified the source as a class I YSO. The object spatially
coincides less than $30\arcsec$ with the IRAS 05379+3550 source and is listed
in the Lumsden et al. (2013) catalog of massive young stars derived from the
MSX survey (Egan et al., 2003). Preliminary spectral energy distribution (SED)
modeling of the source by Dewangan & Anandarao (2011) revealed it to be a
high-mass protostellar object. Further studies of IRAC ratio maps
($4.5\,\mu$m/$3.6\,\mu$m) and narrowband H2 images (2.12$\,\mu$m) by Dewangan
et al. (2016) showed shocked and molecular emission surrounding the central
source. The maps also revealed an extended region with elongated structures,
indicating the presence of an outflow. These results suggest that the S235
e2s3 protostar presents active accretion and ejection processes. To understand
better the physical properties and the nature of the extended structures, we
have carried out a new NIR imaging polarimetric study toward the region of the
embedded protostar.

Figure 2: NIR gray-scale image of the total intensity emission obtained in $H$
band from the polarimetric observations. The image was constructed by average-
combining a set of 18 images resulting in a total integration time of 2 hrs.
The left panel image spans a $4\arcmin\times 4\arcmin$ FOV of the POLICAN
instrument. The right panel shows a zoomed-in view of the S235 e2s3 outflow,
spanning a $1.2\arcmin\times 1.2\arcmin$ field equivalent to a physical scale
of about 0.57 pc. The two bipolar lobes are marked NW and SE. The total length
of the bipolar nebulosity ($\sim 100,000$ au) is shown as a red dashed line.
The line also corresponds to the position angle of the outflow axis
($137^{\circ}$). In both the panels the location of the S235 e2s3 protostar is
marked by a purple cross.

The paper is organized as follows: in Section 2, we describe details of the
observations and data reduction techniques. In Section 3, we present
polarization maps and show the spatial distribution of the outflow. In Section
4, we discuss the results about the morphology of the outflow, the
polarization properties, magnetic fields, and the H2 luminosity. Finally in
Section 5, we summarize the work with concluding remarks.

##  2 Observations and Data Sets

###  2.1 NIR polarization observations

NIR linear polarimetric observations were obtained in $H$ band ($1.6\,\mu$m)
using the POLICAN instrument (Devaraj et al., 2018b) attached to the CANICA
camera (Carrasco et al., 2017) on the $2.1\,\mathrm{m}$ OAGH telescope in
Cananea, Sonora, Mexico. The camera has a HgCdTe detector with
1024$\times$1024 pixels and a pixel-scale of $0.32\,\mathrm{arcsec}$. The
field of view (FOV) achieved with the instrument is $4\arcmin\times 4\arcmin$.
The observations were acquired on the 26th of September 2017, centered on the
S235 e2s3 source with coordinates $\alpha_{2000}$ = 05h41m22.705s,
$\delta_{2000}$ = +35$\arcdeg$52$\arcmin$08.93$\arcsec$ (Dewangan & Anandarao,
2011). The average atmospheric seeing during the observations was about
$0\arcsec.9$. The angular size of the $4\arcmin$ FOV corresponds to a physical
scale of about $2\,\mathrm{pc}$ at a distance of $1.65\,\mathrm{kpc}$.

Polarimetric imaging was performed through combination of a rotating half-wave
plate (HWP) as the modulator and a fixed wire-grid polarizer as the analyser.
In order to measure the linear polarization, the HWP is rotated to four
position angles $(0^{\circ},22$.$5^{\circ},45^{\circ},67$.$5^{\circ})$ to
complete one set of images. Each set of images were obtained at 18 sky
dithered positions to estimate the sky contribution and boost the signal-to-
noise ratio when the images were combined. The exposure times were set to
$100\,\mathrm{s}$ per image, totaling an integration time of $2\,\mathrm{hrs}$
for the 72 images (18$\times$4 HWP).

Image reduction consisted of dark subtraction, flat-fielding, and sky removal.
After basic reduction, the dithered images were aligned and combined for each
HWP angle. They were then astrometry-corrected using data from 2MASS
(Skrutskie et al., 2006). The pixel intensities resulting from the four
combined HWP images were used to compute the Stokes parameters $I$, $Q$, and
$U$ for the field (see Devaraj et al., 2018b).

The polarization degree ($P$) and position angle ($PA$) were computed from the
Stokes parameters, respectively, as follows: $P=\sqrt{Q^{2}+U^{2}}/I$,
$PA=\frac{1}{2}\tan^{-1}\left(\frac{U}{Q}\right)$. The polarization
uncertainty ($\sigma_{P}$) was also computed from the corresponding Stokes
errors. Since $P$ is positively biased from the quadrature combination of $Q$
and $U$, the de-biased degree of polarization ($P^{{}^{\prime}}$; also
expressed as percentage polarization) was estimated using a Ricean correction
as $P^{{}^{\prime}}=\sqrt{P^{2}-\sigma_{P}^{2}}$ (Wardle & Kronberg, 1974).
The foreground polarization contribution from interstellar dust toward the
S235 region is minimal at $\sim 0.3\%$ (Devaraj et al., 2021). Hence, the
removal of foreground polarization is neglected in the analysis. Full
descriptions of the calibration procedures and polarimetric analyses are given
in Devaraj et al. (2018a, b). The final polarization values at each pixel were
stored as an image array and were binned for every 3$\times$3 pixels to create
the polarization maps for the observed region.

###  2.2 Archival data sets

Additional archival data sets, ranging from NIR to radio wavelengths, were
used in our study. MIR images were obtained from Spitzer IRAC and MIPS data
sets through guaranteed time observations (GTOs) (Fazio et al., 2004). IRAC
images (at 3.6, 4.5, 5.8, and $8\,\mu$m) have an angular resolution of $\sim
2\arcsec$ and the MIPS $24\,\mu$m images have $6\arcsec$ resolution.
Narrowband H2 ($\nu=1-0$ S(1); $\lambda/\Delta\lambda=2.122\,\mu$m/0.032
$\mu$m) imaging data were obtained from the extended H2 emission survey
(Navarete et al., 2015), conducted using WIRcam at the CFHT. The survey also
provided $K$ band continuum images, which were used to get the final
continuum-subtracted H2 map.

For the purpose of SED modeling, NIR $K$ band photometric data were taken from
2MASS (Skrutskie et al., 2006). MIR photometric fluxes were taken from the
WISE all sky release (3.4, 4.6, 12, $22\,\mu$m) (Wright et al., 2010), MSX
Galactic plane survey (8.3, 12.1, 14.7, 21.3$\,\mu$m) (Egan et al., 2003) and
AKARI point source catalog (Abrahamyan et al., 2011). Far-IR and millimeter
fluxes (at 350$\mu$m and 1.1 mm) were taken from the extended Bolocam Galactic
Plane Survey, version 2.1 (Aguirre et al., 2011; Ginsburg et al., 2013;
Merello et al., 2015). James Clark MAxwell Telescope (JCMT) Scuba observations
were used for the flux at 850$\,\mu$m (Klein et al., 2005).

##  3 Analyses and Results

###  3.1 NIR total intensity and Polarization Maps

The NIR $H$ band total intensity image from the observations is shown in
Figure 2. The image is centered on the S235 e2s3 source (indicated by a purple
cross) and extends to a size of $4\arcmin\times 4\arcmin$, equivalent to a
physical scale of $2\,\mathrm{pc}\times 2\,\mathrm{pc}$. The most prominent
result seen in the image is the two large extended structures of nebulosity
appearing to trace an outflow. The nebulosity clearly shows its distinction as
a collimated outflow and differs from the typical IR reflection nebulosity due
to a dusty envelope around a compact H ii region. A zoomed-in view of the
region toward the nebulosity is shown in the right panel of Figure 2. The two
structures are diametrically opposite to each other and form the morphology of
bipolar outflow lobes. One lobe is oriented north-west, whereas the other is
pointed to the south-east. These are marked as NW and SE, respectively. The
nebulosity in the NW lobe is much closer to the central source when compared
to the SE lobe. We determined a rough value for the position angle of the
outflow axis by bisecting the outflow lobes around their spatial distribution.
The position angle was estimated to be around $137^{\circ}$ (measured from the
north-up to east-left direction). The red dashed line in Figure 2 is shown
parallel to the outflow axis. There is no significant NIR emission in the
central region between the lobes. This could be due to the fact that the
central source is deeply embedded in circumstellar material and not visible in
our $H$ band observations.

The polarimetric imaging results are presented in Figure 3 as maps of the
degree of polarization and position angles. The maps have been constructed
using polarization values estimated at each pixel as described in Section 2.1.
The left panel in Figure 3 shows a gray-scale image of the degree of
polarization with the color bar at the top indicating the range of values.
Darker regions correspond to higher polarization and lighter regions for lower
polarization. The polarization levels in the region are remarkably high and
reach a maximum of about 80$\%$. The bipolar nebulosity is the most prominent
structure seen in the polarized flux. The regions associated to the nebulosity
show high degree of polarization across both lobes, averaging about 70$\%$. A
few foreground/background stars toward the outflow are revealed in the
polarization map that are not visible in the total intensity image in Figure
2. The stars have significantly lower polarizations of less than 10%. There is
no significant polarization emission at the location of the central source. A
streamer-like elongated structure is visible from the source center toward the
SE lobe with polarization values of up to 30%. This could be from dust in or
around the jet from the young star. In general, the high polarization levels
in the nebulosity show excellent correlation, spatially matching the
morphology of the outflow seen in Figure 2.

Figure 3: Polarization maps of the S235 e2s3 bipolar outflow. The left panel
shows the de-biased degree of polarization. A gray-scale color bar at the top
shows the percentage range of the degree of polarization. The right panel
shows the polarization position angles represented by red vectors. The lengths
of the vectors correspond to the degree of polarization. The polarization
angles display a centrosymmetric pattern around the S235 e2s3 source, which is
marked by a purple cross in both panels.

The distribution of polarization position angles is shown as red vectors in
the right panel of Figure 3. Only the regions associated to the nebulosity
having a polarization level threshold above 3$\sigma$ are plotted. The most
notable feature in the map is the centrosymmetric pattern of the polarization
vectors. The pattern is equally distributed in both bipolar lobes. This
feature together with the very high degree of polarization is indicative that
the nebulosity is being illuminated by a single source at its center. The
light from the illuminating source is singly scattered at near right angles
into our line of sight. Such scattering mechanism has been associated with
regions with low optical depth (Gledhill et al., 2001; Gledhill , 2005). The
polarization at the central source location (near the purple cross) is faint
and the polarization vectors seem to appear perpendicular to the axis of the
outflow. This polarization pattern could originate from the dusty disk around
the young star (see Section 4.1 for more details). However, there is no direct
indication of the presence of a disk, as our observations lack the resolution
to resolve one. Overall, the results of the polarization maps with the
centrosymmetric pattern suggest that the NW and SE lobes belong to the same
region and emanate from a single source.

###  3.2 Narrowband H2 image

Observations in NIR narrowband H2 emission trace shock-heated gas in molecular
outflows and are useful to investigate the ejection processes in young stars
(Navarete et al., 2015; Kim et al., 2018). In Figure 4 we present a three-
color composite image of the outflow region using NIR continuum-subtracted H2
emission at 2.12$\,\mu$m (green), NIR $H$ band at 1.6$\,\mu$m (red), and
Spitzer 4.5$\,\mu$m (blue). The displayed H2 emission is $3\sigma$ above the
background. We identified several H2 emission features toward the nebulosity
which are composed of discrete ﬂows/shocks. These are named as 1a, 2a, 2b, and
2c in Figure 4. The most remarkable result is that the H2 emission is bipolar
and very closely coincides with the axis of the outflow in the plane of the
sky. The shocked emission is narrow and well-collimated within the opening
angles of the outflow. The extent of H2 emission is similar to the dusty
nebulosity and reaches to $\sim 19\arcsec$ or about 31,000 au from the center.
The maximum cross sectional width is $\sim 6\arcsec$ or about 9,900 au.

The NW outflow apparently consists of single shocked H2 emission area at its
tip (marked as 1a). The SE region has two distinct elongated emissions. One
very close to the central source (marked as 2a) and the other located before
the SE dusty lobe (marked as 2b). The 2b H2 emission coincides exceptionally
well with the boundary of the dusty nebulosity. This region particularly
reveals a scenario consisting of a bow shock and the different layers
associated with it. There is a small emission region located away from the
southern tip of the nebulosity (marked as 2c). This emission does not coincide
directly with the outflow axis, but is within $10\arcsec$ of the outflow
region. The positional offset in the 2c emission could be due to several
reasons, such as projection effects or precession of the jet. We do not find
any further H2 emission related to the outflow beyond our FOV when comparing
the large-scale H2 map of the region (Dewangan et al., 2016). We note that in
Figure 4 toward the west of the nebulosity, there is a long vertical feature
with diffuse H2 emission that is not spatially associated to the outflow. This
feature is due to UV-excited emission (fluorescence) from the S235 main
ionizing star (see Figure 1 for details). The overall morphology of the H2
emission, along with the NIR dusty nebulosity and the position of the central
source (blue in Figure 4) exhibit a well-defined ejection process occurring
with the S235 e2s3 source.

Figure 4: Three-color composite map of the outflow region, constructed using
$H$ band 1.6$\,\mu$m (red), H2 emission 2.12$\,\mu$m (green), Spitzer
4.5$\,\mu$m (blue). The image spans $1.2\arcmin\times 1.2\arcmin$,
corresponding to a projected physical scale of $0.57\,\mathrm{pc}\times
0.57\,\mathrm{pc}$. The H2 emission in green with its multiple components are
marked as 1a, 2a, 2b, and 2c. The S235 e2s3 source is seen in blue from the
Spitzer 4.5$\,\mu$m image and is outlined by a dashed circle.

###  3.3 SED modeling

=13mm  Facility/ Wavelength Flux Densitya Comments Instrument
$\lambda\,(\mu\rm{m})$ (Jy) 2MASS K 2.2 0.46 (0.08) $\times 10^{-3}$ upper
limit WISE: W1 3.4 7.87 (0.37) $\times 10^{-3}$ upper limit WISE: W2 4.6 34.2
(0.8) $\times 10^{-3}$ upper limit Spitzer/Ch3 5.8 51.3 (0.1) $\times 10^{-3}$
upper limit Spitzer/Ch4 8.0 105 (0.04) $\times 10^{-3}$ MSX/A 8.3 0.20 (0.01)
WISE: W3 11.6 0.62 (0.01) MSX/C 12.13 0.96 (0.11) MSX/D 14.6 1.16 (0.1)
AKARI/IRC-L 18 4.50 (0.17) MSX/E 21.3 5.29 (0.34) WISE: W4 22.1 11.18 (0.2)
AKARI/FIS-N60 65 178 (8.13) AKARI/FIS-L 140 230 (105) CSO/Bolocam 350 19.95
(0.38) JCMT/SCUBA 850 1.38 (0.52) CSO/Bolocam 1100 1.04 (0.12)

Table 1: Fluxes from NIR to radio wavelengths used for the SED modeling.

The SED modeling of the S235 e2s3 source was carried out to derive its
physical parameters and properties of the mass accreting disk/envelope.
Dewangan & Anandarao (2011) previously carried out preliminary SED analysis
for the S235 e2s3 source with limited data points using Robitalle et al.
(2007) SED fitting. However, the Robitalle et al. (2007) SED models were
developed mainly for low-mass protostars without considering massive star
formation properties such as high accretion rates and high-mass surface
densities. De Buizer et al. (2017) studied a sample of massive protostars from
the SOFIA survey and showed that Robitalle et al. (2007) models typically
result in high-mass infall rates but about 100 times smaller disk accretion
rates, leading to smaller bolometric luminosities and a larger protostellar
masses.

Hence, here we use SED fitting tool ‘SedFitter’ that is part of the python
package `sedcreator`111https://github.com/fedriani/sedcreator implemented by
Fedriani et al. (2023). The SED fitting code is developed using radiative
transfer grid models of intermediate and high-mass stars from Zhang & Tan
(2018), based on the turbulent core accretion theory of massive star formation
(McKee & Tan, 2003). The initial conditions in the massive star formation
model are pressurized, dense, massive cores embedded in high-mass surface
density clump environments. These massive stars are assumed to form from
preassembled massive prestellar cores supported by internal pressure from a
combination of magnetic fields and turbulence. In the model grid there are
three main parameters that set the physical properties: the initial mass of
the core ($M_{c}$), the mean mass surface density of the clump environment
($\Sigma_{\mathrm{cl}}$), and the protostellar mass ($M_{*}$), which defines
the location along an evolutionary track from a given initial condition. The
properties of protostellar cores, including the protostar, disk, infall
envelope, outflow, and their evolution, are calculated self-consistently from
the given initial conditions. In addition to the three main physical
parameters, there are also two secondary parameters: the inclination angle of
the outflow to the line of sight ($\theta_{\mathrm{view}}$) and the level of
foreground extinction ($A_{V}$). Thus, there are five parameters
($M_{c},\Sigma_{\mathrm{cl}},M_{*},\theta_{\mathrm{view}},$ and $A_{V}$) that
give rise to a protostellar SED. The observations are then fitted to the model
grid by varying these 5 parameters and calculating the $\chi^{2}$ (see Eq. 1
from De Buizer et al., 2017). Other properties, such as accretion rate, infall
envelope mass, outflow cavity opening angle, and disk size, are prescribed for
a given set of values of $M_{c},\Sigma_{\mathrm{cl}}$, and $M_{*}$.

In the current Zhang & Tan (2018) model grid, the main parameters are sampled
as follows: $M_{c}$ is sampled at 10, 20, 30, 40, 50, 60, 80, 100, 120, 160,
200, 240, 320, 400 and 480 M⊙ and $\Sigma_{\mathrm{cl}}$ is sampled at 0.10,
0.316, 1.0 and 3.16 $\mathrm{g\,cm^{-2}}$, which makes a total of 60
evolutionary tracks. For each track, $M_{*}$ is sampled at 0.5, 1, 2, 4, 8,
12, 16, 24, 32, 48, 64, 96, 128 and 160 M⊙. In the end, this yields a total of
432 physical models that have difference combinations of
$M_{c},\Sigma_{\mathrm{cl}},M_{*}$. Then, for each physical model there are 20
viewing angles sampled uniformly at
$\cos\theta_{\mathrm{view}}=0.975,0.925,\cdots,0.025$, i.e., equally
distributed from 1 (face-on) to 0 (edge-on). Therefore, a total of $432\times
20=8640$ SEDs are considered. The SED models are also convolved with the
transmission profiles of the instrument filters to simulate the fluxes
detected in the observational bands of the various instruments. Finally, the
visual extinction is constrained from 0 to a user-defined value $A_{V,{\rm
max}}$. For our SED fit, we provided $A_{V,{\rm max}}=60$ mag based on $A_{V}$
estimates of the S235 e2s3 clump (Devaraj et al., 2021). A distance of 1.65
kpc was adopted for the analysis. Table 1 summarizes the reliable fluxes used
from the NIR to radio wavelengths for the SED modeling. Note that the data at
$\lesssim 8\,{\rm\mu m}$ are treated as upper limits for the SED fitting since
short-wavelength fluxes at $\lesssim 8\,{\rm\mu m}$ are affected by polycyclic
aromatic hydrocarbons (PAHs) and thermal emission from very small grains, and
these effects have not been included in the radiative transfer models (see
Zhang & Tan, 2018).

Figure 5: SED of the S235 e2s3 source using the Fedriani et al. (2023) fitting
tool. The filled red circles are observed fluxes taken from archival data sets
(see Table 1), with upper limits represented with downward-facing arrows. The
gray lines represent the ‘good’ models that satisfy the condition
$\chi^{2}<2\chi_{\mathrm{min}}^{2}$ (gray line shades decrease with increasing
$\chi^{2}$). The black solid line corresponds to the best-fit model. The
derived best-fit model physical parameters of the source are listed in Table
2.

Figure 5 shows results of the SED modeling with observed fluxes over
wavelength. Only models that satisfy the criterion of
$\chi^{2}<2\chi_{\mathrm{min}}^{2}$, which are considered ‘good’ models (see
Fedriani et al., 2023, for more details), were selected in the analysis. A
total of 10 models satisfied this criterion and are shown as gray lines in
Figure 5 with their shades decreasing with increasing $\chi^{2}$. In order to
obtain the optimum values from the SED fitting, we take the average of the 10
‘good’ SED model parameters (Fedriani et al., 2023). These values then
represent the best-fit model for the source, which is shown as a solid black
line in Figure 5.

The computed values from the best-fit model suggest that S235 e2s3 is a
protostar with a stellar mass of $6.8\pm 1.2\,\hbox{$M_{\odot}$}$ and with a
total bolometric luminosity ($L_{\mathrm{bol}}$) of $9.63\pm 2.1\times
10^{3}\,\hbox{$L_{\odot}$}$. The isotropic bolometric luminosity
($L_{\mathrm{bol,iso}}$) i.e. without correction for the inclination and
extinction is $1.68\pm 0.74\times 10^{3}\,\hbox{$L_{\odot}$}$. Other
parameters estimated indicate the star has a disk accretion rate of
3.6$\times$10-4 M${}_{\odot}\,\mathrm{yr^{-1}}$ and a foreground extinction
$A_{V}$ of $\sim 22$ mag. Table 2 summarizes the values of the best-fit model
into two sections highlighting the SED dependent and derived parameters.

Table 2: SED fitting results for S235 e2s3

Parameter | Best-fit valuesb
---|---
Core mass, $M_{\mathrm{c}}$ ($M_{\odot}$) | 18.6 (17.4–19.8)
Mass surface density, $\Sigma_{\mathrm{cl}}$ ($\mathrm{g\,cm^{-2}}$) | 3.16 (2.16–4.16)
Stellar mass, $M_{*}$ ($M_{\odot})$ | 6.8 (5.6–8.0)
Inclination angle, $\theta_{\mathrm{view}}$ (∘) | 78 (66–87)
Foreground extinction, $A_{\mathrm{V}}$ (mag) | 22.5 (12.8–32.2)
Stellar age, $T_{\mathrm{now}}$ (yr) | 3.1 (2.0–4.2)$\times 10^{4}$
Stellar temperature, $T_{*}$ (K) | 7095 (6078–8112)
Core radius, $R_{\mathrm{c}}$ (au) | 3670 (2500–4840)
Envelope mass, $M_{\mathrm{env}}$ ($M_{\odot}$) | 5.4 (4.1–6.8)
Disk mass, $M_{\mathrm{disk}}$ ($M_{\odot}$) | 2.3 (1.5–3.2)
Disk accretion rate, $\dot{M}_{\mathrm{disk}}$ ($M_{\odot}\mathrm{\,yr^{-1}}$) | 3.6 (2.4–4.8)$\times 10^{-4}$
Bolometric luminosity, $L_{\mathrm{bol}}$ ($L_{\odot}$) | 9.63 (7.5–11.7)$\times 10^{3}$
Isotropic Bol. luminosity, $L_{\mathrm{bol,iso}}$ ($L_{\odot}$) | 1.68 (0.94–2.4)$\times 10^{3}$
Outflow opening angle, $\theta_{\mathrm{w,esc}}$ (∘) | 34 (28–40)

  * b

Mean values for $\chi^{2}<2\chi_{\mathrm{min}}^{2}$. The values given in the
parenthesis are the ranges in the parameter values for the ‘good’ models (see
text).

##  4 Discussion

NIR polarimetric observations have revealed a bipolar nebulosity around the
S235 e2s3 source. The polarization maps show a centrosymmetric pattern due to
dust scattering. Molecular H2 emission is found to display extended features
along the nebulosity, tracing shocked regions of the jet. All this evidence
illustrates that the S235 e2s3 source is dynamic and presents a disk/outflow
system. In this section, we discuss the outflow morphology, polarization
characteristics, orientation of the outflow with respect to the local magnetic
field, and the properties of the molecular H2 emission.

###  4.1 Outflow morphology and polarization characteristics

The most important results seen in our observations is the discovery of the
bipolar nebulosity of an outflow (see Figure 2 and 3). The presence of bipolar
outflow structures in the NIR continuum around young stars is very rare since
normally they are very extincted or most of the dust in the outflow is
destroyed due to the energetics of the protostar jets. Conversely, bipolar
structures have been observed in several cases around post-AGB stars due to
their slow mass loss leading to dusty nebulosities (Weintraub et al., 2000;
Serrano Bernal et al., 2020). The visibility of structures in our NIR total
intensity image could be due to longer integration time which is limited in
other surveys such as 2MASS (Skrutskie et al., 2006) and UKIDSS (Lawrence et
al., 2007).

The total projected length of the bipolar nebulosity spans a size of about
100,000 au ($\sim$0.5 pc), indicating a large extent in its dust distribution.
The nebulosity in individual lobes has an extension of about 33,000 au.
Several studies of outflows from embedded Class 0/I high-mass YSOs using
molecular observations have shown that they have lengths greater than 1 pc
(Navarete et al., 2015; Bally , 2016). However, in our case the NIR continuum
observations reveal the dust in the outflow and it is not the best tracer to
constrain the full length of the outflow.

The general morphology of the outflow lobes is symmetric extending in the NW-
SE direction. The symmetry of the outflow and its bipolar nature has few
implications. Firstly, it is reasonable to assume that the environment around
the protostar is less turbulent, thus allowing for the outflow to propagate
unperturbed and form a high degree of symmetry. This also indicates that the
outflow is likely driven by an isolated protostar and is not located around
clustered star formation. Secondly, a quiet low turbulent environment might
imply particularly smooth and undisturbed accretion from the core onto the
protostellar disk and the star, thus favoring a very regular and uniform
occurrence of accretion. Although these assumptions are simplified, it
presents an overview of the conditions around the outflow.

The outflow morphology is also useful to infer the inclination of the star
system with respect to the line of sight of the observer. Since the outflow is
bipolar, with the lobes being similar in size and seen easily separated
spatially in the sky, we can presume that the outflow is nearly perpendicular
to the observer. The derived inclination angle from the SED fitting in Section
3.3 corroborates this view. The polarization levels in our observations reach
a very high degree of about $70-80\%$ in the outflow nebulosity. Only few
observations (for e.g. from Hubble Space Telescope) have reported similarly
high values (Sahai et al., 1998; Meakin et al., 2005). Such high levels of
polarization can only be measured if the scattered light from the dust
particles are nearly at right angles to the line of sight of the observer.
Hence, the scattering region should be perpendicular to the line of sight,
therefore indicating the bipolar outflow lies approximately on the plane of
sky.

Figure 6: Magnetic field map toward the S235 e2s3 outflow region obtained from
NIR background starlight polarimetry (Devaraj et al., 2021). The polarization
position angles are shown as red vectors and they trace the magnetic field
orientation in the plane of sky. The lengths of the vectors indicate the
polarization percentage and a reference vector length is shown at the bottom
left. The right panel shows a histogram of the polarization position angles
representing the distribution of magnetic field orientations. The dotted red
dashed line represents the position angle of the outflow axis.

The polarization angles demonstrate centrosymmetric pattern indicative of
single scattering mechanism in an optically thin medium. Gledhill  (2005) used
polarimetric models of bipolar nebulosities to derive dust scattering
properties. They found that objects with low optical depth ($\tau$) in the
range $0.1<\tau<1$ tend to produce highly symmetric polarization similar to
our observations. The high polarization levels are also indicative that the
dust grains in the nebulosity are smaller ($<0.3\,\mu$m) than the scattering
wavelength and could be made up of graphites or silicates (Shure et al., 1995;
Lucas & Roche , 1998). At the region toward the central source, we see a small
area with polarization angles perpendicular to the outflow axis. The optical
depth toward the center is generally higher as these regions are extincted in
the observations. The circumstellar envelope of the accreting protostar could
be dominant in the central region, since a dense dust core has been identified
at this position at millimeter wavelengths (Dewangan et al., 2016; Devaraj et
al., 2021). According to the dust scattering models of Whitney & Hartmann
(1993); Fischer et al. (1994), a polarization disk can be expected for
scattered light emerging from an optically thick disk around the illuminating
source. This particular case of disk polarization occurs since light passing
through the disk will be scattered at least twice (above and below the plane)
and loses the simple centrosymmetric polarization pattern produced by single
scattering. Hence, it is likely that the disk of S235 e2s3 is edge on to the
observer, producing the observed polarization angles at the center. Again,
this is consistent with the outflow being in the plane of the sky.

It is expected that the driving source of a bipolar outflow is likely at the
center of the two lobes. However, this identification may not be
straightforward when several infrared sources are visible in the field which
can contribute to the observed outflow structures. Hence, the identification
and position of the driving source needs to be confirmed beyond what may be
inferred based on the simple morphological properties and spatial
distribution. We carried out polarization centroiding (Weintraub et al., 2000)
to identify the location of the illuminating source. This technique uses the
scattered light to trace back the source responsible for scattering. The
centrosymmetric polarization vectors in our images were rotated by
$90^{\circ}$ and traced to identify their point of intersection, where the
expected illuminating star is located. All the polarization vectors in our
observations intersect within $3\arcsec$ of the outflow central region. This
coincides with the exact position of the S235 e2s3 source as revealed from MIR
observations. Hence, it is clear that our protostar is illuminating and
driving the bipolar outflow.

###  4.2 Outflow orientation with respect to local magnetic field

Theoretical developments in star formation studies suggest that the collapse
of a cloud may proceed mainly along the magnetic field orientation, since the
partly ionized gas can flow easily along the field lines (Mouschovias , 1976;
Shu, Adams & Lizano , 1987). At smaller scales when the collapse of the core
is predominantly along the magnetic field lines, it is expected that the
magnetic braking of the rotating core is more efficient when the spin axis of
the core is parallel to the magnetic field (Mouschovias & Paleologou , 1980).
This scenario would result in the formation of a disk perpendicular to the
field lines. Hence, one expects the field lines to be drawn into an hourglass
shape and the outflows oriented parallel to the magnetic field (Galli & Shu ,
1993).

A number of observations toward young stars have found outflows parallel to
the ambient magnetic field. These include the outflow from the Mon R2 core
(Bally & Lada , 1983), HH 80-81 jet (Carrasco-Gonzalez et al., 2010), IRAS
18089-1732 high-mass system (Beuther et al., 2010), and few low-mass Class 0
protostars (Chapman et al., 2013). However, more recent statistical studies of
outflows using millimeter observations do not provide a consistent picture.
Hull et al. (2013) studied the alignment of 30 cores ($\sim 1000$ au) using
dust polarization maps from the TADPOL survey (Hull et al., 2014) and found
that the magnetic field lines are not tightly aligned with the outflow axis.
Rather, they find the data to be consistent with scenarios where outflows and
magnetic fields are if anything preferentially perpendicular.

In Figure 6 we present a local magnetic field map toward the S235 e2s3 region
to compare their relative orientations with respect to the bipolar outflow.
The magnetic field directions are obtained from NIR background starlight
polarimetry based on dichroic extinction from dust grains aligned to the
magnetic field. The complete polarimetric analysis and selection of bona fide
stars for tracing magnetic field directions is presented in Devaraj et al.
(2021). We used only the sample of stars that were associated within the 2 pc
FOV of the outflow region. The results in the left panel of Figure 6 show that
the global magnetic field orientations in the region are uniform. The
morphology of the outflow also appears to be along the magnetic field lines.
To compare the observations quantitatively we have shown a histogram plot of
the polarization position angles in the right panel of Figure 6. The histogram
distribution peaks at values of about $145-150^{\circ}$. The position angle of
the outflow axis determined earlier ($137^{\circ}$) is shown as a red dashed
line. Both the outflow axis and the magnetic field lines have orientations
within $10^{\circ}$ and indicate that the outflow is parallel to the ambient
magnetic field. This result in isolation supports the general theory of core
collapse along the field lines.

Jones & Amini  (2003) found that there is a rough correlation between the
outflow axis and the magnetic field orientation at the farther edge of the
outflow nebulosity. They determined that when the outflow lies in the plane of
sky, the magnetic field tends to be parallel to the outflow axis. This
particularly matches our observations as our outflow is in the plane of sky
and the field lines are traced in the surrounding medium. Matsumoto et al.
(2006) studied the alignment of outflows with magnetic fields in cloud cores
through numerical simulations. They found that the outflow tends to be aligned
with the large-scale ($>5000$ au ) magnetic field if the magnetic field in the
core is larger than $80\,\mu G$. Devaraj et al. (2021) studied the magnetic
field strength in the core of S235 e2s3 and found it to be subcritical with a
value of $74\pm 14\,\mu G$. This value including the uncertainties is close to
the threshold value of $80\,\mu G$ and may well indicate the observations
match the predictions.

One important aspect that should be considered in the context of the S235 e2s3
outflow is the impact of the massive star in the S235 complex. Since the
outflow is at the boundary of the expanding shell-like H ii region, it is
possible that the turbulent shock fronts and ionizing radiation from the
massive star influence core collapse at its edges. In this scenario the
magnetic fields are also dragged along the shell. This has been investigated
by many studies (Kirsanova et al., 2014; Dewangan et al., 2016; Devaraj et
al., 2021) suggesting triggered star formation in this region. Hence, it is
also possible that the outflow axis being parallel to the magnetic field in
our case is due to the impact of the massive star. To draw a definite
conclusion on such sources, it will be useful to carry out a comprehensive
study of the alignment of outflows with magnetic fields in intermediate to
high-mass stars and compare the results with their magnetic field strengths
and outflow inclination angles.

###  4.3 H2 Luminosity and shocked emission properties

Many studies have revealed H2 outflows on scales of thousands of astronomical
units to parsecs around intermediate and high-mass YSOs (Davis et al., 2004;
Varricatt et al., 2010; Caratti o Garatti et al., 2015). These observations
show that the outflows are well-collimated and present similarities with
outﬂows from low-mass YSOs. A key parameter in these studies is the H2
luminosity and its relation to the source luminosity, which directly links the
outflow properties with the driving source properties.

We studied the S235 e2s3 outflow H2 luminosity and compared it with the source
properties. The H2 luminosity was derived by estimating the H2 1-0 S(1)
emission line flux ($F_{2.12}$) from the continuum-subtracted image. First, we
calculated the total sum of pixel counts for each of the H2 emission features
(1a, 2a, 2b, and 2c) indicated in Figure 4. A factor of 1.10 was used to
multiply to the sum in order to compensate for the H2 line flux that is
subtracted from the K band image (e.g. Lee et al., 2019). We get the total sum
of H2 emission to be 50,310 counts for an area of 272 arcsec2. After zero-
point correction and calibration, this resulted in a total 2.12$\,\mu$m line
flux of $F_{2.12}=2.28\times 10^{-16}$ ${\rm W\,m}^{-2}$

Next, the 2.12$\,\mu$m H2 luminosity ($L_{2.12}$) was estimated using the
relation

$L_{2.12}=4\pi{D^{2}}F_{2.12}\times 10^{0.4A_{2.12}}$ (1)

where $A_{2.12}$ is the foreground extinction at 2.12$\,\mu$m adopted as
$2.3\pm 1$ mag (see Section 3.3; $A_{2.12}\approx 0.1A_{V}$) and $D$ is the
distance ($1.65\,\rm{kpc}$) to the source.

The total rovibrational H2 luminosity ($L_{\rm H_{2}}$) can then be
approximated from $L_{2.12}$ by assuming that the shocked emitting gas is in
local thermal equilibrium (LTE). Typically, H2 jets from low-mass YSOs have an
average temperature of $\sim 2000\,\rm{K}$, and the $L_{2.12}$ is about $10\%$
of the total $L_{\rm H_{2}}$ (Caratti o Garatti et al., 2006). However, more
massive jets driven by massive YSOs usually have average temperature of
$2500-3000\,\rm{K}$ (Caratti o Garatti et al., 2015), higher than low-mass
jets. Therefore, $L_{2.12}$ is about $7\%$ of the total $L_{\rm H_{2}}$. This
yields a total H2 luminosity of the S235 e2s3 outflow to be $L_{\rm
H_{2}}=2.3_{-1.3}^{+3.5}\,\hbox{$L_{\odot}$}$.

This estimation has many uncertainties as we are using a general value of the
extinction toward individual H2 emission features, and assuming an approximate
ratio of the total H2 to H2.12 luminosity. Nevertheless the estimated H2
luminosity places reasonable constraints on the properties of the jet.

Figure 7: Log ($L_{\rm H_{2}}$) vs. Log ($L_{\rm bol}$) from Caratti o Garatti
et al. (2015), with results combined from previous works (Caratti o Garatti et
al., 2006, 2008), as labeled in the upper left corner. The red dashed line
indicates the fit for data points from low-mass YSOs (Caratti o Garatti et
al., 2006). The blue dashed line shows the best linear fit resulting from the
sample of high-mass YSOs. The black dashed-dotted line shows the best linear
fit of high-mass YSOs, excluding the five high-mass outliers (see Caratti o
Garatti et al., 2015, for more details). The S235 e2s3 source is represented
by a filled purple circle.

The H2 luminosities of outflows driven by low-mass YSOs are typically lower
than the luminosity of the S235 e2s3 outflow. For example, the $L_{\rm H_{2}}$
in Caratti o Garatti et al. (2006) study of 23 protostellar jets driven by
low-mass YSOs ranges from 0.007 to 0.76 $L_{\odot}$. Similarly, the outﬂows
detected in Serpens/Aquila from the UWISH2 survey (Ioannidis & Froebrich ,
2012) show observed $L_{\rm H_{2}}$ ranging from 0.01 to 1.0 $L_{\odot}$.
These studies are indicative that the luminosity of the S235 e2s3 outflow
corresponds to that of a higher-mass star.

Several studies show that a tight correlation exists between outflow
energetics and YSO properties. Caratti o Garatti et al. (2006) derived an
empirical relationship between $L_{\rm H_{2}}$ and $L_{\mathrm{bol}}$ for a
large sample of Class 0 and Class I low-mass YSOs. They found a power-law
correlation of $L_{\rm H_{2}}\propto L^{0.58}_{\mathrm{bol}}$. Caratti o
Garatti et al. (2008) extended this study to include a high-mass YSO: IRAS
20126+4104 and found that the same correlation exists for this source. This
relationship was then investigated for a sample of 18 massive jets from
intermediate and high-mass YSOs by Caratti o Garatti et al. (2015). They
inferred that the $L_{\rm H_{2}}$ and $L_{\mathrm{bol}}$ relationship behaves
similarly as the low-mass sample and the correlation was found to be $L_{\rm
H_{2}}\propto L^{0.62}_{\mathrm{bol}}$.

Figure 8: Plot of the length of the bipolar jet as a function of the
bolometric luminosity of its associated source. Values from Navarete et al.
(2015) are shown as filled green circles while those from Varricatt et al.
(2010) are indicated by opened red circles. Dashed lines mark the boundary of
Log ($L/\hbox{$L_{\odot}$}$=2.5), at which there is a change in the upper
limit of $l_{\rm proj}$. The values of the S235 e2s3 source are represented by
the filled purple circle.

In Figure 7 we show the plot of $L_{\rm H_{2}}$ versus $L_{\mathrm{bol}}$ from
all the aforementioned studies. The obtained best linear fits are indicated by
different colored dashed lines. The $L_{\rm H_{2}}$ and $L_{\mathrm{bol}}$ of
S235 e2s3 are also shown by filled purple circle. We find that the values of
S235 e2s3 is distributed close to other high-mass YSOs, matching the $L_{\rm
H_{2}}$ versus $L_{\mathrm{bol}}$ relationship. An overarching trend seen from
these observations is that the H2 luminosity increases with the bolometric
luminosity of the source, which at the protostellar stage, accounts for most
or all of the source’s bolometric luminosity. This suggests that the mass loss
is mainly driven by the accretion power, which grows as the mass of the
central source increases. In addition, it is also evident that the $L_{\rm
H_{2}}$ versus $L_{\mathrm{bol}}$ is uniform across both high and low-mass
protostars. Therefore, this strongly suggests that the outflows from massive
YSOs are similar in nature to those from low-mass YSOs (i.e., produced by
disk-mediated accretion).

The comparison between the morphology of the outflow and the intrinsic
characteristics of the driving source is another important aspect to consider
in the evolution of protostars. The projected length ($l_{\rm proj}$) and
collimation factor ($R$ = length/width) are useful parameters in this context.
In Section 3.2 we find that the S235 e2s3 jet has a $l_{\rm proj}$ of 31,000
au (i.e. the distance from the source center to the farthest emission of the
jet) and a width of 9900 au. The $l_{\rm proj}$ value is generally a lower
limit for the real length of the structure as the inclination angle
$\theta_{\mathrm{view}}$ of the system affects the jet’s appearance on the
plane of sky. The real length of the jet can be estimated as $l_{\rm
real}=l_{\rm proj}/\sin(\theta_{\mathrm{view}})$ which leads to $l_{\rm
real}=31700$ au for our outflow. With these values known, the collimation
factor of the S235 e2s3 jet was estimated to be $R=3.2$.

Studies of high-velocity outflows from massive YSOs have shown collimation
factors of $R=2.05$ as compared to $R=2.8$ for low-mass stars (Wu et al.,
2004), indicating a weak tendency that outflows associated with massive stars
are less collimated than those from low-mass stars (Richer et al., 2000).
However, S235 e2s3 has a large collimation factor suggesting that massive YSOs
can also drive highly collimated outflows.

Navarete et al. (2015) carried out an H2 survey on a large sample of massive
YSOs and compared the distribution of jet projected lengths as a function of
source bolometric luminosity. They also included sources from the studies of
Varricatt et al. (2010). Figure 8 shows a plot of $l_{\rm proj}$ versus
$L_{\mathrm{bol}}$ from both these studies. We have also included the values
of the S235 e2s3 outflow represented by the filled purple circle. The length
of the S235 e2s3 jet is seen to be lower than the average length of sources
with similar $L_{\mathrm{bol}}$. This indicates that jet lengths can vary
enormously for a given stellar luminosity. The overall distribution of data
points in the plot shows a large scatter, although we also find a weak trend
that the largest outflows are associated with high-luminosity sources. This
may suggest that generally luminous protostars drive longer jets, but it is
not the main ingredient in determining the outflow length. Instead, the
outflow lengths can be dependent on other factors such as the evolutionary
stage of the driving source. A comparison of $l_{\rm proj}$ versus
$L_{\mathrm{bol}}$ for different classes of YSOs will be valuable in this
context.

##  5 Conclusions

We have presented a highly collimated bipolar outflow discovered from the S235
e2s3 protostar using NIR polarimetric observations. The outflow reveals a rare
bipolar dusty nebulosity in the NIR continuum images. Narrowband H2
observations show shocked emission along the outflow tracing the bipolar jet’s
interaction with the ambient material. Based on the observations and using
archival data, we have studied the properties of the outflow and derived the
physical parameters of the protostar. The main results of the study are
summarized as follows.

1\. The bipolar outflow presents a high degree of polarization ($\sim$80%) and
reveals a centrosymmetric pattern in the polarization position angles. This
suggests that the outflow nebulosity is illuminated by a single source at its
center whose light is seen singly scattered at near right angles to our line
of sight. Based on this, the outflow was inferred to be perpendicular to the
observer. These characteristics also indicate that the S235 e2s3 protostar is
the source responsible for driving the outflow.

2\. SED fitting of the source using models of turbulent core accretion theory
suggest that the S235 e2s3 protostar has a mass of $6.8\pm
1.2\,\hbox{$M_{\odot}$}$ with a total bolometric luminosity
($L_{\mathrm{bol}}$) of $9.63\pm 2.1\times 10^{3}\,\hbox{$L_{\odot}$}$. The
protostar has a disk accretion rate of 3.6$\times$10-4
M${}_{\odot}\mathrm{\,yr^{-1}}$ and an age of approximately $3.1\times 10^{4}$
yr. All these parameters are consistent with the values obtained for other
massive protostars.

3\. The orientation of the outflow axis was found to be parallel to the
ambient magnetic field direction, conforming to the expectations from
magnetically induced core collapse theory. However, since the outflow is at
the boundary of an H ii region, feedback from the massive star might also
influence the orientation of the outflow.

4\. Narrowband H2 ($2.12\,\mu m$) observations show multiple molecular
emission features along the bipolar jet. These were identified as regions of
H2 flows, shocks, or knots. We estimated the total H2 luminosity of the jet to
be $L_{\rm H_{2}}=2.3_{-1.3}^{+3.5}\,\hbox{$L_{\odot}$}$. This value, along
with the $L_{\rm bol}$, closely matched the relationship between H2 luminosity
and bolometric luminosity obtained for low and high-mass YSOs (Caratti o
Garatti et al., 2015). These results support the theory that jets from high-
mass YSOs are similarly launched to those from low-mass YSOs via disk
accretion.

5\. The total length of the outflow spans about 0.5 pc, with its collimation
factor being $R=3.2$. This indicates massive YSOs can also drive highly
collimated flows. Comparison of the jet length with its bolometric luminosity
does not indicate any relationship between those parameters. Hence, the size
of the jet generally does not scale with the mass of the protostar.

We suggest that follow-up spectroscopic observations of the protostar along
with longer wavelength observations of the outflow will be very useful to
understand better the dynamics and properties of this massive star
disk/outflow system.

This research has been supported by the European Research Council advanced
grant H2020-ER-2016-ADG-743029 under the European Union’s Horizon 2020
Research and Innovation program. D.R. acknowledges INAOE and CONACyT-Mexico
for SNI grant (CVU 555629). A.L. acknowledges financial support by CONACyT-
Mexico through project CB-A1S54450. A.C.G. has been supported by PRIN-INAF-
MAIN-STREAM 2017 “Protoplanetary disks seen through the eyes of new-generation
instruments” and by PRIN-INAF 2019 “Spectroscopically tracing the disk
dispersal evolution (STRADE)”. R.F. has received funding from the European
Union’s Horizon 2020 research and innovation program under the Marie
Skłodowska-Curie grant agreement No 101032092. We thank all the OAGH staff for
their help in conducting observations with POLICAN. The authors acknowledge
the use of SAOImage DS9 software which is developed with the funding from the
Chandra X-ray Science Center, the High Energy Astrophysics Science Archive
Center and JWST Mission office at Space Telescope Science Institute.

## References

  * Abrahamyan et al. (2011) Abrahamyan, H. V., Mickaelian, A. M, & Knyazyan, A. V., 2015, A&C, 10, 99
  * Aguirre et al. (2011) Aguirre, J. E., Ginsburg, A. G., Dunham, M. K., et al. 2011, ApJS, 192, 4
  * Banerjee & Pudritz (2007) Banerjee, R., Pudritz, R. E., 2007, ApJ, 660, 479
  * Banerjee & Pudritz (2006) Banerjee, R., Pudritz, R. E., 2006, ApJ, 641, 949
  * Bally & Lada (1983) Bally, J., & Lada, C. J., 1983, ApJ, 265, 824
  * Bally (2016) Bally, J., 2016, ARA&A, 54, 491
  * Blandford & Payne (1982) Blandford, R. D., Payne, D. G., 1982, MNRAS, 199, 883
  * Bastien & Menard (1988) Bastien P., Menard F., 1988, ApJ, 326, 334
  * Beuther et al. (2010) Beuther, H., Vlemmings, W. H. T., Rao, R., van der Tak, F. F. S., 2010, ApJ, 724, L113
  * Brown et al. (2018) Brown, A. G. A., Vallenari, A., Prusti, T., et al. 2018, A&A, 616, A1
  * Caratti o Garatti et al. (2006) Caratti o Garatti, A, Giannini, T., Nisini, B., Lorenzetti, D., 2006, A&A, 449, 1077
  * Caratti o Garatti et al. (2008) Caratti o Garatti, A, Froebrich, D., Eisloffel, J., Giannini, T. & Nisini, B., 2008, A&A, 485, 137
  * Caratti o Garatti et al. (2015) Caratti o Garatti, A., Stecklum, B., Linz, H., Garcia Lopez, R. & Sanna, 2015, A&A, 573, A82
  * Carrasco-Gonzalez et al. (2010) Carrasco-Gonzalez C., Rodriguez L. F., Anglada G. et al., 2010, Science, 330, 1209
  * Chapman et al. (2013) Chapman, N., Davidson, J. , Goldsmith, P. F., Houde, M. et al., 2013, ApJ, 770, 151
  * Chavarría et al. (2014) Chavarría L., Allen L., Brunt C., et al. 2014, MNRAS, 439, 3719
  * Carrasco et al. (2017) Carrasco, L., Hernández Utrera, O., Vázquez, S., et al. 2017, Rev. Mexicana Astron. Astrofis., 53, 497
  * Davis et al. (2004) Davis, C. J., Varricatt, W. P., Todd, S. P., & Ramsay Howat, S. K. 2004, A&A, 425, 981
  * De Buizer et al. (2017) De Buizer, J. M., Liu, M., Tan, J. C., et al., 2017, ApJ, 843, 33
  * Devaraj et al. (2018a) Devaraj, R., Mayya, Y. D., Carrasco, L., & Luna, A., 2018a, PASP, 130, 055001
  * Devaraj et al. (2018b) Devaraj, R., Luna, A., Carrasco, L. et al., 2018b, PASP, 130, 055002
  * Devaraj et al. (2021) Devaraj, R., Clemens, D. P., Dewangan, L. K., Luna, A., et al., 2021, ApJ, 911, 81
  * Dewangan & Anandarao (2011) Dewangan, L. K., & Anandarao, B. G. 2011, MNRAS, 414, 1526
  * Dewangan et al. (2016) Dewangan, L. K., Ojha, D. K., Luna, A., et al. 2016, ApJ, 819, 66
  * Dewangan & Ojha (2017) Dewangan, L. K. & Ojha, D. K., 2017, ApJ, 849, 66
  * Egan et al. (2003) Egan, M. P., Price, S. D., & Kraemer, K. E. 2003, BAAS, 35, 1301
  * Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10
  * Fedriani et al. (2018) Fedriani, R., Caratti o Garatti, A., Coffey, D., Garcia Lopez, R., et al., 2018, , A&A, 616, A126
  * Fedriani et al. (2023) Fedriani, R., Tan, J. C., Telkamp, Z., Zhang, Y., et al., 2023, ApJ, 942, 7
  * Frank et al. (2014) Frank, A., Ray, T. P., Cabrit, S., Hartigan, P., Arce, H. G., et al., 2014, PPVI conf., Univ. of Arizona Pres, 451
  * Fischer et al. (1994) Fischer, O., Henning, T., Yorke, H. W., 1994, A&A, 284, 187
  * Galli & Shu (1993) Galli, D., Shu, F. H., 1993, ApJ, 417, 220
  * Georgelin et al. (1973) Georgelin, Y. M., Georgelin, Y. P., & Roux S. 1973, A&A, 25, 337
  * Ginsburg et al. (2013) Ginsburg, A., Glenn, J., Rosolowsky, E., et al. 2013, ApJS, 208, 14
  * Gledhill et al. (2001) Gledhill T. M., Chrysostomou A., Hough J., Yates J., 2001, MNRAS, 322, 321
  * Gledhill (2005) Gledhill T. M., 2005, MNRAS, 356, 883
  * Gredel et al. (2006) Gredel, R. 2006, A&A, 457, 157
  * Jones & Amini (2003) Jones, T. J., & Amini, H., 2003, AJ, 125, 1418
  * Jones et al. (2004) Jones, T. J., Woodward, C. E., & Kelley, M. S. 2004, ApJ, 128, 2448
  * Hull et al. (2013) Hull, C. L. H., Plambeck, R. L., Bolatto, A. D., et al. 2013, ApJ, 768, 159
  * Hull et al. (2014) Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, ApJS, 213, 13
  * Ioannidis & Froebrich (2012) Ioannidis, G., & Froebrich, D. 2012, MNRAS, 425, 1380
  * Kim et al. (1994) Kim, S., Martin, P. G., & Hendry, P. D., 1994, ApJ, 422, 164
  * Kim et al. (2018) Kim, H-J, Koo, B-C, Pyo, T-S, & Davis, C. J., 2018 ApJ, 863, 74
  * Klein et al. (2005) Klein, R., Posselt, B., Schreyer, K., Forbrich, J., & Henning, T. 2005, ApJS, 161, 361
  * Kirsanova et al. (2008) Kirsanova, M. S., Sobolev, A. M., Thomasson, M., et al. 2008, MNRAS, 388, 729
  * Kirsanova et al. (2014) Kirsanova, M. S., Wiebe, D. S., Sobolev, A. M., et al. 2014, MNRAS, 437, 1593
  * Lawrence et al. (2007) Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599
  * Lee et al. (2019) Lee Y.-H., Koo B.-C., Lee J.-J., Burton, M. & Ryder, S., 2019, AJ, 157, 123
  * Lucas & Roche (1998) Lucas P. W., Roche P. F., 1998, MNRAS, 299, 699
  * Lumsden et al. (2013) Lumsden, S. L., Hoare, M. G., Urquhart, J. S., et al. 2013, ApJS, 208, 11
  * Matsumoto et al. (2006) Matsumoto, T., Nakazato, T., & Tomisaka, K. 2006, ApJ, 637, 105
  * McKee & Tan (2003) McKee, C. F. & Tan, J. C. 2003, ApJ, 585, 850
  * Meakin et al. (2005) Meakin, C. A., Hines, D. C., & Thompson, R. I., 2005, ApJ, 634, 1146
  * Merello et al. (2015) Merello, M., Evans II, N. J., Shirley, Y. L., Rosolowsky, E., et al., 2015, ApJS, 218, 1
  * Mouschovias (1976) Mouschovias, T. C. 1976, ApJ, 207, 141
  * Mouschovias & Paleologou (1980) Mouschovias, T. C., & Paleologou, E., V. 1980, ApJ, 237, 877
  * Navarete et al. (2015) Navarete, F., Damineli, A., Barbosa, C. L., & Blum, R. D. 2015, MNRAS, 450, 4364
  * Pendleton et al. (1990) Pendleton, Y., Tielens, A. G. G. M., & Werner, M. W. 1990, ApJ, 349, 107
  * Pudritz & Norman (1983) Pudritz, R. E., & Norman, C. A. 1983, ApJ, 274, 677
  * Ray et al. (2007) Ray, T., Dougados, C., Bacciotti, F., et al. 2007, PPV conf., Univ. of Arizona Press, 231
  * Richer et al. (2000) Richer J. S., Shepherd D. S., Cabrit S., Bachiller R., Churchwell E., 2000, PPIV conf., Univ. of Arizona Press, 867
  * Robitalle et al. (2007) Robitaille T. P., Whitney B. A., Indebetouw R., Wood K., 2007, , ApJS, 169, 328
  * Sahai et al. (1998) Sahai, R., Hines, D. C., Kastner, J. H., Weintraub, D., et al. 1998, ApJ, 492, L163
  * Serrano Bernal et al. (2020) Serrano Bernal. E. O, Sabin, L., Luna, A., Devaraj, R., et al., 2020, MNRAS, 495, 2599
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Shang et al. (2007) Shang H., Li Z.-Y., Hirano N., 2007, PPV conf., Univ. of Arizona Press, 261
  * Shu, Adams & Lizano (1987) Shu, F. H., Adams F. C., Lizano S. 1987, ARA&A, 25, 23
  * Shu et al. (1994) Shu F., Najita J., Ostriker E., Wilkin F., Ruden S., Lizano S., 1994, ApJ, 429, 781
  * Shu et al. (2007) Shu F. H., Galli D., Lizano S., Glassgold A. E., Diamond P. H., 2007, ApJ, 665, 535
  * Shure et al. (1995) Shure, M., Sellgren, K., Jones, T. J., Klebe, D., 1995, AJ, 109, 721
  * Tamura et al. (1991) Tamura, M., Gatley, I., Joyce, R. R., et al. 1997, ApJ, 378, 611
  * Tamura et al. (1997) Tamura, M., Hough, J. H., Chrysostomou, A., et al. 1997, MNRAS, 287, 894
  * Varricatt et al. (2010) Varricatt W. P., Davis C. J., Ramsay S., Todd S. P., 2010, MNRAS, 404, 661
  * Wardle & Kronberg (1974) Wardle, J. F. C., & Kronberg, P. P. 1974, ApJ, 194, 249
  * Werner et al. (1983) Werner, M. W., Capps, R. W., & Dinerstein, H. L. 1983, ApJ, 265, L13
  * Weintraub et al. (2000) Weintraub, D. A., Kastner, J. H., Hines, D. C., Sahai, R., 2000, ApJ, 531, 401
  * Whitney & Hartmann (1993) Whitney B. A., Hartmann, L., 1993 ApJ, 402, 605
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, et al. 2010, AJ, 140, 1868
  * Wu et al. (2004) Wu Y., Wei Y., Zhao M., Shi Y., Yu W., Qin S., Huang M., 2004, A&A, 426, 503
  * Zhang & Tan (2018) Zhang, Y. & Tan, J. C. 2018, ApJ, 853, 18