# Principal component analysis tomography in near-infrared integral field
spectroscopy of young stellar objects. I. Revisiting the high-mass protostar
W33A

F. Navarete1, A. Damineli1, J. E. Steiner1, R. D. Blum2
1 Universidade de São Paulo, Instituto de Astronomia, Geofísica e Ciências
Atmosféricas, Rua do Matão 1226,
Cidade Universitária São Paulo-SP, 05508-090, Brasil
2 NSF’s Optical–Infrared Astronomy Research Laboratory P.O. Box 26732, Tucson,
AZ 85719, USA
E-mail<EMAIL_ADDRESS>(FN)In Memoriam

(Accepted 2021 February 4. Received 2021 February 4; in original form 2020
November 29)

###### Abstract

W33A is a well-known example of a high-mass young stellar object showing
evidence of a circumstellar disc. We revisited the $K$-band NIFS/Gemini North
observations of the W33A protostar using principal components analysis
tomography and additional post-processing routines. Our results indicate the
presence of a compact rotating disc based on the kinematics of the CO
absorption features. The position-velocity diagram shows that the disc
exhibits a rotation curve with velocities that rapidly decrease for radii
larger than 0$\aas@@fstack{\prime\prime}$1 ($\sim$250 AU) from the central
source, suggesting a structure about four times more compact than previously
reported. We derived a dynamical mass of 10.0${}^{+4.1}_{-2.2}$
$\rm{M}_{\odot}$ for the “disc+protostar” system, about $\sim$33% smaller than
previously reported, but still compatible with high-mass protostar status. A
relatively compact H2 wind was identified at the base of the large-scale
outflow of W33A, with a mean visual extinction of $\sim$63 mag. By taking
advantage of supplementary near-infrared maps, we identified at least two
other point-like objects driving extended structures in the vicinity of W33A,
suggesting that multiple active protostars are located within the cloud. The
closest object (Source B) was also identified in the NIFS field of view as a
faint point-like object at a projected distance of $\sim$7,000 AU from W33A,
powering extended K-band continuum emission detected in the same field.
Another source (Source C) is driving a bipolar $\rm{H}_{2}$ jet aligned
perpendicular to the rotation axis of W33A.

###### keywords:

methods: statistical – techniques: imaging spectroscopy – stars: protostars –
stars: pre-main sequence – ISM: jets and outflows – ISM: individual: W33A

††pubyear: 2021††pagerange: Principal component analysis tomography in near-
infrared integral field spectroscopy of young stellar objects. I. Revisiting
the high-mass protostar W33A–4

##  1 Introduction

In the last few decades, integral field unit (IFU) observations have become
crucial to understand the kinematics of very complex extended objects over a
variety of scales, from protostellar objects and protoplanetary discs (e.g.
Murakawa et al. 2013; Alves et al. 2019) to galaxies (Kewley et al., 2019).
The powerful three-dimensional data cubes produced by IFU spectrographs are of
increasingly large size, containing an enormous amount of information and
noise that must be pulled from the cube by sophisticated means. For example, a
typical near-infrared data cube may have millions of pixels, and given the
complexity of the spatial and spectral information, often times much of the
information is hard to extract and may be simply ignored in the analysis in
preference of a restricted subset of the data (line maps and their products,
spectral analysis of bright regions, etc.). Frequently, a large fraction of
the data in these data cubes contains just instrumental or sky noise. Thus,
robust procedures to extract information from the large data sets are
increasingly required in astrophysics.

In this context, Steiner et al. (2009) developed a method for the use of the
Principal Component Analysis (PCA) Tomography in optical and near-infrared IFU
data cubes. PCA is a non-parametric technique, that is, it does not require
any assumption or expectation of the physical and geometrical parameters of
the data. PCA compresses the data expressed as a large set of correlated
variables (spatial and spectral pixels) in a small but optimal set of
uncorrelated and orthogonal variables (i.e. the eigenvectors and their spatial
projections), ordered by their relative importance to explain the initial data
set.

PCA tomography has been recently used in extra-galactic studies (e.g. Ricci et
al., 2011; Menezes et al., 2013; Ricci et al., 2014, 2015; Dahmer-Hahn et al.,
2019), but no application to stellar astrophysics cases has been made yet.
Thus, we applied this technique to the K-band IFU observations of the well-
known high-mass protostar W33A investigated by Davies et al. (2010) to check
if the PCA Tomography is able to recover their findings. W33A is a well-known
young stellar object (YSO) with indications of an active accretion process
ongoing. This protostar is associated with the Red MSX Source object
G012.9090$-$00.2607, with a bolometric luminosity of $\log(L_{\rm
bol}/L_{\odot})=4.51$ (Lumsden et al., 2013) and located at the distance of
2.4 kpc (Immer et al., 2013).

Davies et al. (2010) presented the study of this source through three-
dimensional spectroscopy in the $K$-band (2.2 µm) using NIFS/Gemini North IFU
observations. Those authors reported the identification of the key-structures
predicted by the disc-mediated accretion scenario (Shu et al., 1987): a
circumstellar disc traced by the CO features in absorption; a jet-like
structure traced by the $\rm{H}_{2}$ emission at 2.12 µm; and the cavity of
the jets close to their launching point traced by the spectro-astrometry of
the $\rm{Br\gamma}$ emission at milli-arcsecond scales. Those authors derived
the Keplerian mass of $\sim$15 $\rm{M}_{\odot}$ for the “protostar + disc”
system, assuming a distance of 3.8 kpc.

Given the importance of W33A as one of the rare cases of high-mass YSOs with
clear signs of active accretion, the K-band IFU observations of the high-mass
protostar W33A offer an excellent opportunity to verify the efficiency of PCA
Tomography on reproducing the results reported by Davies et al. (2010). We are
also interested to check how the results are sensitive to the additional post-
processing routines suggested by Steiner et al. (2009) (e.g. differential
atmospheric refraction correction, high-frequency noise filtering, and
instrumental fingerprint removal; see also Menezes et al. 2014), and if PCA
Tomography can provide additional information based on the correlation of the
spectral features present in the NIFS K-band observations.

This work is organised as follows. Section 2 presents the observations, the
post-processing routines of the data and the description of the PCA tomography
technique. Section 3 presents the results and interpretation of the data
through PCA tomography, the derivation of the physical parameters of the
source and the analysis of auxiliary near-infrared maps. Section 4 presents
the discussion of the data and comparison with recent findings in the
literature. Our conclusions are given in Section 5.

##  2 Observations and data processing

###  2.1 $K$-band IFU observations

The K-band ($\lambda$ $\approx$ 2.2 $\mu$m) IFU observations of W33A were
obtained with the Near-Infrared Integral Field Spectrometer (NIFS, McGregor et
al., 2003). NIFS has a full wavelength coverage in the K-band of about 4200 Å,
with a linear dispersion of $\Delta\lambda_{\rm disp}$ = 2.13 Å/pixel-1 and a
nominal spectral resolving power of $\lambda/\Delta\lambda_{\rm disp}$ = 5160
(Blum & McGregor, 2008), or $\approx$ 58 km s-1 in the velocity scale. NIFS
observes a 3″$\times$3″ field-of-view in a 60$\times$60 pixel2 array,
delivering a spatial sampling of 0$\aas@@fstack{\prime\prime}$050 pixel-1.

The observations were carried at the Gemini North 8-m telescope on the night
of May 25th, 2008. The raw data of W33A were downloaded from the Gemini
Science Archive (project GN-2008A-Q-44). The data include a set of calibration
frames required for processing the science observations: a set of flat-field
images taken with Quartz lamps on (flat-on) and off (flat-off); flat-field
images taken with a Ronchi mask (Ronchi flats), used to correct the spatial
distortion along the spectral dimension; a set of dark images, taken with the
same exposure time of the science and telluric sources; and an Argon-Neon lamp
spectrum (Ar-Ne) for wavelength calibration of the data. The A2 V double star
HR6798A+B was observed as the telluric standard star, used for correcting the
telluric absorption features, and for performing the flux-calibration of the
final data cubes.

The data were processed using the usual reduction procedure based on
IRAF/Gemini package. The telluric correction was applied to the data using a
one-dimensional template spectrum extracted from the telluric standard star.
The telluric $\rm{Br\gamma}$ feature was removed from the template by fitting
a Voigt profile. Then, the data cubes were constructed using the Gemini
NIFCUBE procedure, by interpolating the original data into a final spatial
pixel (spaxel) size of 0$\aas@@fstack{\prime\prime}$050, resulting in a set of
six science data cubes. The final spatial resolution achieved in the
observations is about 0$\aas@@fstack{\prime\prime}$14 $\pm$
0$\aas@@fstack{\prime\prime}$01, obtained as the mean FWHM of the telluric
standard double star.

Residual telluric emission features were identified as relatively strong
emission lines in all data cubes. Since these features often exhibit
significant variability in the FOV that might affect the results of the PCA
analysis, they were manually suppressed from the data by interpolating the
adjacent spectral pixels in the affected regions for each spaxel of the data
cube.

###  2.2 Supplementary near-IR maps

Additional $J$ ($\lambda$ $\approx$ 1.25 µm), H (1.64 µm), $\rm{H}_{2}$ (2.14
µm) and $K$ (2.22 µm) images of a 90″$\times$90″ FOV around the W33A protostar
were obtained from the UKIRT Widefield Infrared Survey for H2
(UWISH2)111http://astro.kent.ac.uk/uwish2/ (Froebrich et al., 2011). The large
scale maps have a pixel scale of 0$\aas@@fstack{\prime\prime}$2 pixel-1 and a
mean seeing of about 0$\aas@@fstack{\prime\prime}$96 at the K-band, offering a
complementary view of the stellar and nebular content in the vicinity of W33A.

###  2.3 Post-processing of the IFU data

The data cubes generated by the standard IRAF reduction procedure were further
processed using a set of routines written in IDL (Interactive Data Language).
The purpose of these additional procedures include $i)$ the atmospheric
refraction correction and combination of the data cubes, $ii)$ the flux-
calibration, and $iii)$ the high-frequency noise suppression from the data.

####  2.3.1 Differential Atmospheric Refraction Correction

The individual data cubes were corrected for the effects of the differential
atmospheric refraction (DAR) using the procedure described in Menezes et al.
(2014). Figure 1 compares the data before and after the DAR-correction based
on the centroid position of the point-like source of W33A as a function of the
wavelength for the entire spectral region covered by the K-band data cube. The
plot shows the rms of the source position reduces from $\sim$1.5 mas to about
0.5 mas, indicating the DAR-correction is a significant correction to the
displacement introduced by atmospheric effects at NIR wavelengths. The
procedure adopted for performing the DAR-correction of the data cubes also
takes into account any dithering offsets used for the observations, delivering
a set of aligned data cubes. The displacements introduced by the DAR effects
are significant for the analysis involving the entire continuum (e.g. see
Figs. 4 and 5 from Menezes et al. 2014), or different spectral features along
the wavelength range covered by the data (e.g. the $\rm{H}_{2}$ lines).
Indeed, Menezes et al. (2014) reported that the total displacement introduced
by the DAR effects can be larger than the typical NIFS spatial resolution from
$J$\- to K-band wavelengths. However, we do not expect a critical enhancement
of the results introduced by the DAR correction on the analysis of single
spectral lines with widths of a few tens to hundreds of Angstroms. However,
the single line analysis benefits from the post-processing routines aimed at
the suppression of the noise from the datacube, as further presented in Sect.
2.3.3 and 2.4.2.

A final data cube was obtained by median combining the six individual DAR-
corrected data cubes. We compared the residuals between the median-combined
and an average-combined data cube to check if the fluxes were preserved. We
noted that the largest deviations were about $\sim$ 1% of the flux in the
brightest regions of the FOV, closer to the W33A point source, and $\sim$ 5%
in low signal-to-noise regions, indicating no significant discrepancy between
the usage of the median or the average of the fluxes.

Figure 1: Centroid position of the W33A source as a function of the wavelength
(from 2.02 to 2.42 µm, from blue to red) of the same data cube before (top
panel) and after the atmospheric refraction correction (bottom). The (0,0)
position is indicated as a $\times$ symbol and the error bars corresponds to
the standard deviation of the position across the spectral dispersion.

####  2.3.2 Flux calibration of the data cubes

The median-combined data cube was flux-calibrated using a custom IDL routine
based on the STANDARD, SENSFUNC and CALIBRATE routines of the IRAF-NOAO
package. The three corrections performed by each routine are applied
simultaneously by computing the ‘calibration factor’ function ($C_{\lambda}$),
defined as:

$C_{\lambda}=2.5\log\left(\frac{O_{\lambda}}{t_{\rm{exp}}~{}B_{\lambda}(T_{\rm{eff}},m_{Ks})}\right)+A\,E_{\lambda}$
(1)

where, $O_{\lambda}$ corresponds to the observed telluric counts of the
observation (in ADU), $t_{\rm{exp}}$ is the exposure time of the observation
(in seconds), $B_{\lambda}(T_{\rm{eff}},m_{Ks})$ is the black-body flux (in
erg s-1 cm-2 $\mu$m-1) at a given effective temperature $T_{\rm{eff}}$, scaled
by the Ks-band magnitude of the telluric standard, $A$ is the average airmass
of the observations and $E_{\lambda}$ is the atmospheric extinction curve (in
mag airmass-1 units).

We adopted the Ks-band flux of the 2MASS catalogue (Skrutskie et al., 2006) to
derive the flux-calibration of the data cubes. Since the HR6798A+B system is
unresolved at the 2MASS catalogue but resolved in the NIFS observations, the
construction of $O_{\lambda}$ was performed by computing the integrated
spectrum of the resolved components of the telluric standard star over the
spatial pixels (spaxels) with average flux above a 3-$\sigma$ threshold above
the value of the sky background. Then, the continuum emission was modelled by
fitting a polynomial function of order $n$ $\leq$ 4, using the outlier-
resistant IDL algorithm robust_poly_fit. Next,
$B_{\lambda}(T_{\rm{eff}},m_{Ks})$ was constructed by following a two-step
procedure. First, the shape of the black-body spectrum was obtained by using
the effective temperature corresponding to the spectral type of the telluric
standard (A2 V), taken from Tokunaga (2000). Then, the black-body emission was
scaled using the Ks-band flux of the telluric star (Ks = 5.87 mag). The 2MASS
Ks-band magnitudes ($m_{\rm Ks}$) were converted to fluxes ($F_{\rm Ks}$)
using Eq 2:

$m_{Ks}=-2.5\log\left(\frac{F_{Ks}}{F_{Ks,0}}\right)$ (2)

where $F_{Ks,0}$ = 4.28 $\cdot$ 10-7 erg s-1 cm-2 µm-1 is the zero-point flux
at the Ks-band (Cohen et al., 2003).

Finally, the additive term of Eq. (1) requires a model for the atmospheric
extinction as a function of the wavelength. Thus, we adopted the mean
extinction through the broad K-band filter for Mauna Kea, $E_{\lambda}$ =
0.033 mag/airmass, reported by Tokunaga et al. (2002).

After constructing $C_{\lambda}$, each spaxel of the data cube (here,
represented as $D_{xy\lambda}$) was divided by the exposure time of the
observations, and was extinction-corrected and flux-calibrated as indicated by
Eq. (3):

$\log(D_{fcal,xy\lambda})={0.4\cdot
C_{\lambda}}\cdot\log\left(\frac{D_{xy\lambda}}{t_{\rm{exp}}}\right)$ (3)

The results of the flux calibration were checked using the integrated telluric
standard spectrum as follows. First, the total flux of the telluric standard
star was integrated over the bandwidth of the 2MASS Ks-band filter
($\lambda_{c}$ = 2.159 µm, $\Delta\lambda$ = 0.262 µm, Cohen et al., 2003).
Then, the corresponding Ks-band magnitude was evaluated using Eq. (2), and
compared with the 2MASS Ks-band magnitude of the telluric star (5.87 mag).
Their difference, $\Delta K_{\rm s}$ = $-$0.03 mag, corresponds to a Ks-band
flux ratio of $\approx$1.028 between the NIFS and 2MASS fluxes, or an error of
about 3 per cent.

Figure 2 shows the integrated emission of the flux calibrated W33A data,
showing the point like object at the upper central region of the FOV
(0$\aas@@fstack{\prime\prime}$0, 1$\aas@@fstack{\prime\prime}$2) and the
associated K-band continuum nebulosity towards the South.

Figure 2: Integrated flux of the K-band NIFS observations of W33A. The image
is shown in the logarithm scale indicated by the colour bar in the right. The
contours are placed at the 1, 10 an 50% of the peak intensity of the map. The
W33A point source is located around the (0$\aas@@fstack{\prime\prime}$0,
1$\aas@@fstack{\prime\prime}$2) position in the FOV.

####  2.3.3 High-frequency spatial noise filtering

NIFS data are provided in a grid of $x$, $y$ and $\lambda$ coordinates, where
the astrophysical information is often mixed with a complex noise structure,
including the photon noise, high-frequency noise and systematic instrumental
noise that is often associated with astronomical data cubes (see, e.g.,
Menezes et al. 2014). Although the photon noise is a well-known and inherent
source of noise in astronomical observations, the nature of the last two
components are in general unknown and may be introduced during the data
acquisition or data processing (e.g. the resampling of the original
0$\aas@@fstack{\prime\prime}$103 $\times$ 0$\aas@@fstack{\prime\prime}$43 to
the 0$\aas@@fstack{\prime\prime}$050 $\times$ 0$\aas@@fstack{\prime\prime}$050
spatial pixels). Despite the unknown origin of the noise in IFU observations,
a considerable fraction of it can be successfully suppressed from the data
based on the application of a low-pass Butterworth filter (explained below)
and the PCA Tomography (Sect. 2.4) on the data.

In the Fourier domain, the spatial information can be expressed in terms of
frequencies ($\nu_{\rm x}$, $\nu_{\rm y}$), allowing us to identify and
suppress high-frequency noise using an appropriate low-pass band filter. Thus,
we used a two-dimensional Butterworth filter in the spatial frequency
dimension to suppress the high-frequency noise present in the data cube,
preserving the frequencies containing the astrophysical information. The
details on the construction of the Butterworth filter can be found in Sect. 5
from Menezes et al. (2014). We used a Butterworth filter with order $n$ = 2
and cut-off frequencies defined in terms of the Nyquist frequency of the data
($R_{\rm Ny}$), corresponding to the highest-order sampled by the Fast Fourier
Transform (FFT). Due to the different size of the original spaxels in the
$x$\- and $y$-direction of the NIFS data (0$\aas@@fstack{\prime\prime}$103 and
0$\aas@@fstack{\prime\prime}$043, respectively), the cut-off frequencies must
be slightly different in each direction. The determination of the cut-off
frequency is a crucial step and must be performed carefully: if overestimated,
it partially corrects the high-frequency noise that should be suppressed; if
underestimated, it filters frequencies containing the astrophysical signal,
causing an enlargement of the Point Spread Function (PSF).

For W33A observations, frequencies greater than $f_{\rm cutoff,x}$ = 0.5
$R_{\rm Ny}$ and $f_{\rm cutoff,y}$ = 0.7 $R_{\rm Ny}$ were filtered. Figure 3
presents the mean spatial FFT of the W33A data cube before and after the
Butterworth filtering. The cut-off frequencies define an ellipsoid in the
Fourier space, indicated by the white contour on top of the mean FFT maps. We
note that almost all the high-frequency features shown in the top panel of
Fig. 3 were suppressed by the Butterworth filtering as indicated in the bottom
panel of the figure.

Figure 3: Mean Fourier Transform of the spatial dimensions of the data cube
before (top panel) and after the Butterworth filtering (bottom). The white
ellipsoid indicates the cut-off frequency used for constructing the
Butterworth filter to suppress the high-frequency noise of the data.

A simple subtraction of the original cube and the filtered cube can be used to
check if the Butterworth filtering was successfully applied. Figure 4 exhibits
the residuals between the average of the filtered and the original data cubes.
Since the noise is proportional to the signal-to-noise ratio of the data, the
residuals are larger towards the point-like source in the FOV. Despite this,
the residuals exhibit a noisy pattern over the entire FOV. Even though most of
the residuals correspond to only a fraction of the noise of the datacube, the
noise suppression makes a significant difference for the analysis of weak
spectral features (e.g. the CO absorption features at 2.3 µm).

Figure 4: Difference between the NIFS data cube before and after the spatial
Butterworth filtering. The residuals are shown in units of $\sigma$, in a red-
to-blue colour-scale to highlight the details of the noise structure. The
contours corresponds to the 1, 10 and 50% of the peak flux in the original
data cube.

###  2.4 The PCA Tomography

This Section presents a brief description of the Principal Component Analysis
(PCA) Tomography222the routines can be downloaded in
http://www.astro.iag.usp.br/~pcatomography/ technique and its interpretation.
Further details can be found in (Steiner et al., 2009) and references therein.

The PCA is a statistical procedure that transforms a set of multi-dimensional
variables which are mutually dependent (i.e. correlated) into a set of
orthogonal and linearly uncorrelated variables. The PCA Tomography is applied
to the data cubes in order to quantify and extract information, ordered by the
fraction of the variance of these new variables: the first of them explains
the highest amount of the data variance, and the higher-order variables
explain successively smaller amounts of the data variance. In the transformed
coordinate system, the spectral pixels of the data cube correspond to the
variables, and the spaxels are the observables. The new set of uncorrelated
variables – the eigenspectra – are obtained as a function of the wavelength.
The projection of the observables (that is, the spaxels) on the eigenspectra
are images, referred as tomograms.

The new coordinates of the data are called Principal Components (PCs), and
each PC is orthogonal (i.e. uncorrelated) to the others. Each principal
component (PC) is composed of: $i)$ an eigenvalue, $A_{i}$, which indicates
the relevance of the PC in explaining the information contained in the
original data cube; $ii)$ a tomogram, $T_{i}(x,y)$, corresponding to the
projection of the PC on the spatial dimensions of the image; and $iii)$ an
eigenspectrum, $E_{i}(\lambda)$, which is the projection of the PC on the
spectral direction. These three components must be analysed together to
properly evaluate the information contained in each PC. The maximum number of
PCs corresponds to the number of the spectral elements of a given data cube.

The PCA also works as a linear dimensional reduction technique by compressing
the data, expressed initially as a large set of correlated variables, into a
small but optimal set of uncorrelated variables, ordered by their eigenvalues.
Hence, most of the information of the original dataset can be explained with a
relatively small number of PCs. The analysis of the relative value of $A_{i}$
as a function of the PC number – the so-called Scree test – allows us to
determine the best number of PCs that explains the initial data (see Sect.
2.4.1).

Due to the orthogonality of the PCs, the original data cube can be
reconstructed by considering a smaller set of PCs with insignificant losses of
its astrophysical information. This procedure allows us to suppress a
significant amount of noise which is often associated with high-order PCs (see
Sect. 2.4.2).

####  2.4.1 The Eigenvalues and the Scree Test

Each Principal Component is ordered by its eigenvalue (largest to smallest),
$A_{i}$, which indicates the relevance of the PC in explain the information
contained in the data cube. In other words, the eigenvalue is a direct measure
of the variance contained in the PC. The fractional variance of each principal
component $i$, $V_{i}$, is expressed as the ratio between its eigenvalue
$A_{i}$ and the sum of the eigenvalues of all the PCs, given by:

$V_{i}=\frac{A_{i}}{\displaystyle\sum\limits_{i=1}^{N}A_{i}}$ (4)

Table 1 lists the fractional and cumulative variance for selected Principal
Components of the PCA Tomography analysis of W33A. We note that the first PC
contains about 99.8% of the total variance of the data cube, while higher-
order PCs exhibit successively smaller fractional variance values. For W33A,
the maximum number of components is equal to 1782, which is the number of
spectral pixels of the data cube.

Table 1: Scree test of the PCA tomography for the K-band data cube of W33A. Ni | Variance | Cumulative
---|---|---
| (%) | Variance (%)
1 | 99.80360 | 99.8036
2 | 0.18817 | 99.9918
3 | 0.00253 | 99.9943
4 | 0.00109 | 99.9954
5 | 0.00042 | 99.9958
6 | 0.00025 | 99.9960
7 | 0.00021 | 99.9962
8 | 0.00013 | 99.9964
9 | 0.00008 | 99.9965
… |  |
260 | 5.06$\cdot$10-6 | 99.999
… |  |
1782 | 4.61$\cdot$10-10 | 100.00

The analysis of the fractional variance alone does not provide a useful way to
determine how many PCs contains astrophysical information in the datacube.
This diagnosis can be assessed by constructing the Scree test of the PCA,
shown in Fig. 5, which presents the distribution of the fractional variance as
a function of the order of the PC. The Scree test shows a rapid decrease of
the variance as a function of the number of the PC, exhibiting two distinct
regions: the signal-dominated PCs (from PC1 to PC5) from which the variances
are significantly larger than the noise level, and the noise-dominated PCs
(PC6 and onward). The noise level was estimated by performing a linear fit of
$V_{i}$ by $i$ for the components dominated by noise (PC $\geq$ 10).

The analysis and interpretation of the tomograms and eigenspectra of the
Principal components are presented in Sect. 3.3.

Figure 5: Scree test showing the variance of each Principal Component (PC)
versus the Principal Component order. The variance of each PC is presented as
a percentile of the total variance of the data cube (see Eq. (4)). A dashed
line indicates the noise level, based on the linear fit of the data associated
with PCs of order 10 or higher.

####  2.4.2 Noise suppression and reconstruction of the data cube

The first 5 PCs contain about $\approx$99.996% of the input variance of the
data cube (Table 1) and the individual variance of these PCs are above the
noise level (see Fig. 5). Thus, we reconstructed the datacube using the first
five PCs, which led to the suppression of a small but important amount of
noise from the initial data cube.

Figure 6 exhibits the spatial and spectral residuals between the original and
the reconstructed data cube. Cuts of the residuals in the spatial and spectral
dimension exhibit a large-scale noise pattern distributed over the whole FOV,
with relatively stronger intensity closer to the bright point-like source.
These noise patterns are often referred as “instrumental fingerprints" (see
details in Sect. 6 from Menezes et al., 2014) and they correspond to a
systematic noise structure within IFU data cubes. As expected, despite the
presence of the instrumental fingerprints in the residuals, no coherent
spatial or spectral structures related to the astrophysical object were
identified.

Figure 6: Residuals between the original W33A data cube and the reconstructed
data cube using the first 5 PCs. The top panels present the spatial
distribution of the residuals of individual spectral channels indicated by the
vertical filled red lines of the mean spectrum shown on the bottom panel.

##  3 Results

###  3.1 The large-scale near-infrared view of W33A

The analysis of the large-scale K-band emission of W33A was presented by
Davies et al. (2010, hereafter D10), exhibiting the embedded protostellar
object W33A and its extended K-band nebulosity (see their Fig. 1). We took
advantage of the available $JHK$ and $\rm{H}_{2}$ maps of the UWISH2 survey
(Froebrich et al., 2011) to present a more complete view of the infrared
content associated with W33A, including any large-scale $\rm{H}_{2}$ emission
at 2.12 µm, a clear signature of star-forming activity for a broad range of
bolometric luminosity values (e.g. Varricatt et al., 2010; Navarete et al.,
2015).

Figure 7 presents a false-colour $JHK$ map of a 1.5 arcmin2 region around W33A
from the UWISH2 data, overlaid by contours of the continuum-subtracted
$\rm{H}_{2}$ emission. The top panel, shows the K-band reflection nebula
illuminated by W33A, oriented in the SE-NW direction. The weakness of the NW
lobe compared to the SE is then due to the orientation of the system combined
with a strong extinction gradient from the (SE) blue- to the (NW) red-shifted
lobe of the K-band outflow (in Sect. 3.4.2, we show that the $\rm{H}_{2}$
emission identified is blue-shifted in the NIFS FOV, and the presumptive red-
shifted lobe is off the NIFS FOV to the NW). The NIFS FOV is indicated by the
white box, the position of the K-band compact source is indicated by the black
$\times$ symbol, and the position angle (PA = 150∘) of the disc reported by
D10 is shown by the blue arrows, which is aligned with the K-band nebulosity.

The white contours on top of the NIR map indicate the continuum-subtracted
$\rm{H}_{2}$ emission. The $\rm{H}_{2}$ contours indicate no clear extended
emission aligned with the main K-band nebulosity of W33A. Instead, they are
tracing a set of $\rm{H}_{2}$ knots roughly aligned in the E-W direction, as
previously identified by Lee et al. (2013) (see their Fig. 2). At least six
bright $\rm{H}_{2}$ knots (labelled from “a” to “f”) were identified as
arising from a bipolar jet. The jet has a projected length of $\sim$50″,
corresponding to a linear size of 0.58 parsecs at the distance of W33A (2.4
kpc). The large-scale $\rm{H}_{2}$ jet is offset by a few arcseconds to the S
direction of W33A, suggesting that at least another active protostar is
located in the vicinity of W33A.

Figure 7: Top panel: False-colour $JHK$ map (blue: $J$, green: $H$, red: K)
around W33A from the UWISH2 data. The position of the NIFS map is indicated by
the white box. W33A is indicated by the black $\times$ symbol, corresponding
to the (0,0) position in the map (RA = 18:14:39.53, Decl. = $-$17:52:00.0).
The rotation axis of the disc reported by Davies et al. (2010) is indicated by
the blue arrow. The white contours indicate the emission of the continuum-
subtracted $\rm{H}_{2}$ image at the levels $\sigma^{n}$ (where $n$ =
$1,2,3...$). Knots of $\rm{H}_{2}$ emission are labelled in white, and the
white dashed line connects the knots labelled as “d” and “e”. Bottom panel: A
closer view of the $JHK$ map in a 24″$\times$24″box around the W33A. The RGB
map is overlaid by contours of each near-infrared filter (blue: $J$, green:
$H$, red: K) to highlight the details closer to the brightest region of the
main nebula. The $\rm{H}_{2}$ knots “d” and “e” are indicated by the white
contours, together with a dashed white line connecting them. Additional point-
like sources are indicated by black $\times$ symbols and labelled as B and C,
respectively.

To better constrain the origin of the $\rm{H}_{2}$ jet, a zoomed in view of a
24″ $\times$ 24″ region around W33A is shown in the bottom panel of Fig. 7.
Contours of the $J$-, $H$-, and $K$-band filters (blue, green, and red,
respectively) are overlaid to highlight the structures around the W33A source,
including two additional point-like sources located a few arcseconds to the S
of W33A (labelled as B and C), and an extended $H$-band emission oriented to
the NW direction. Source B (RA = 18:14:39.584, Dec. = $-$17:52:02.71) is
likely driving an extended $H$-band emission oriented to the NW direction,
overlapping the $K$-band nebular emission arising from W33A. Source C (RA =
18:14:39.606, Dec. = $-$17:52:06.06) is observed as a relatively weak compact
emission in the $H$-band, located at the central position of the line
connecting the $\rm{H}_{2}$ knots “d” and “e”, suggesting this object is
actually launching the large-scale $\rm{H}_{2}$ knots indicated by the white
contours. Despite the fact that Sources B and C do not exhibit unequivocal
K-band counterparts in the large-scale maps, any K-band emission associated
with Source B and its nebulosity might still be detected in the NIFS FOV.

In Fig. 8 we present the NIFS FOV overlaid by contours of the $JHK$ filters.
Thanks to the relatively small plate scale of the NIR maps
(0$\aas@@fstack{\prime\prime}$2 pixel-1), we were able to extract maps with
the same size as the NIFS FOV, allowing us to compare the high-resolution NIFS
data (0$\aas@@fstack{\prime\prime}$15) with coarser resolution data taken at
shorter wavelengths. The map exhibits the point-source of W33A at the top
central region and the main K-band nebulosity to the South. The comparison
between the black NIFS K-band contours and the contours of the $H$\- and
$K$-band (green and red, respectively) suggests that the main nebulosity may
be a superposition of distinct structures within the NIFS FOV: it is partially
arising from W33A at the upper region of the FOV, while the K-band nebula is
better correlated with the extended $H$-band emission (green contours) towards
the South.

The relatively blue colour of the main nebulosity when compared to W33A also
favours that the extended emission is brighter at shorter wavelengths, as
expected for the extended $H$-band emission detected in the NIFS FOV. At the
bottom left (SE direction from W33A), a faint K-band point-like source peaks
at the brightest and central region of the $H$-band contours, matching the
expected K-band counterpart of Source B identified in the $JHK$ map (see
bottom panel of Fig. 7).

A detailed inspection in the vicinity of W33A confirms that multiple
(proto)stars are located at the line-of-sight of the main K-band nebulosity of
W33A. The complementary information from the $JHK$ maps indicates that both
the compact emission arising from Source B and its extended $H$-band emission
might contribute to the extended structures identified in the NIFS FOV, as we
cannot be sure whether such emission comes from the scattered light in the
outflow cavity or from source B. In addition, the $J-$band counterpart seems
to be a chance alignment of an optically visible star as indicated by DSS
images (not shown in the present work).

Figure 8: False colour composite map of the NIFS field-of-view (blue: 2.05 µm,
green: 2.20 µm, red: 2.35 µm) overlaid by the $JHK$ contours (in blue, green
and red, respectively). The $JHK$ contours are placed at 25, 40, 75, and 90%
of the peak intensity of each map. For comparison, the black contours
delineate the K-band emission of the NIFS observations. The main structures
are labelled on the image.

###  3.2 The NIFS K-band spectra of the bright compact source and the nebular
emission in W33A

In this section, we briefly analyse the overall properties of the K-band NIFS
observations of W33A to guide the reader through the PCA Tomography analysis
presented in Sect. 3.3.

The left panel of Figure 9 presents a composite RGB map of the NIFS FOV based
on the integrated emission of a 0.4 µm band-pass centred at 2.05 (blue), 2.12
(green), and 2.35 µm (red). The red and blue contours correspond to the
regions used for integrating the spectrum of the W33A protostar and the
extended nebulosity, respectively, shown in the right panel. The spectra shown
in the right panel were integrated over the regions indicated by the contours:
the W33A protostar (red), the main nebulosity (blue) and the spaxels
corresponding to Source B (in green). We excluded the contribution of Source B
(the corresponding spaxels are flagged by the green contour in the left panel)
from the nebular emission. The K-band spectrum of Source B (shown as green) is
about $\sim$100 times fainter than the extended nebula (in blue), and does not
significantly affect the integrated spectrum of the brighter region. The
spectrum of Source B will be further discussed in Sect. 3.4.1.

The integrated spectra indicate that the compact source shows a steeper
increase in flux towards larger wavelengths when compared to the nebulosity,
as expected for enshrouded protostellar objects. In addition, the nebulosity
is brighter than the protostar at short wavelengths, suggesting a relatively
strong local reddening effect by the dust towards W33A, that decreases at
larger distances from the source. This may be due to the viewing geometry.
Since the compact source is likely seen through the obscuring disc/torus, it
appears more reddened than the cavity of the outflow (i.e. the nebulosity) we
see in scattered light. The low-$J$ CO absorption lines support this picture
since they are strong toward the point source and missing from the nebular
spectrum as expected if dense cool gas surrounds the point source.

Both the protostar and the nebulosity shows the presence of the
$\rm{Br\gamma}$ feature in emission at 2.166 µm, the first four CO bandhead
features in emission at $\lambda$ 2.29 µm, and the R- and P-branches of the CO
absorption features between 2.32 and 2.37 µm. A weak emission of the
$\rm{Na\,\textsc{i}}$ doubled at 2.21 µm is also observed in both spectra. At
least five ro-vibrational transitions of the $\rm{H}_{2}$ molecule are only
detected towards the extended nebula, which is also roughly oriented in the
direction of an extended $\rm{H}_{2}$emission driven by W33A as reported by
D10.

Figure 9: Left panel: False-colour RGB image of the NIFS FOV (left, blue:
2.05µm, green: 2.12µm, red:2.35µm) using the flux-calibrated data cube. The
contours indicate the regions where the integrated K-band spectra (right
panel) were extracted: the W33A protostar (in red), the main nebulosity (in
blue) and the compact emission corresponding to Source B (in green).. Right:
Integrated K-band spectra of each region indicated in the FOV. The spectrum of
Source B (in green) was multiplied by a factor of 50. No nebular contribution
was subtracted. The main spectral transitions present in the data are
indicated by the vertical dashed lines: the $\rm{H}_{2}$ transitions, the
$\rm{Br\gamma}$ feature at 2.16 µm, the $\rm{Na\,\textsc{i}}$ doublet at 2.21
µm, the CO bandhead features in emission at $\lambda$ 2.29 µm, and the CO
absorption features at $\lambda$ $\sim$ 2.32 µm.

###  3.3 The PCA Tomography of the $K$-band NIFS data cube

In this section, we present the analysis and interpretation of the first five
Principal Components of the PCA tomography applied to the K-band data cube of
W33A. The tomograms and eigenspectra associated to these PCs are shown in
Figs. 10 to 15.

#### PC1

In general, the overall spatial and spectral information from PC1 roughly
correspond to projections of the integrated properties of the data cube (i.e.
its mean in the spatial and spectral dimensions). However, these projections
are not given in positive or negative fluxes, but rather in correlations and
anti-correlations. Such a strong resemblance to the data is expected since PC1
contains the largest relative variance of the data cube (explaining roughly
99.8% of the data; see Table 1.

The tomogram of PC1 exhibits a peak at the position of the bright point-like
source corresponding to the W33A protostar. The main nebulosity is also
identified in the tomogram, exhibiting larger correlation at the vicinity of
W33A and decreasing towards the South. The correlation between the compact and
the extended emission (i.e. the fact that both appear as positive variances in
the tomogram) indicates that both structures shares the same spectral
information shown in the eigenspectrum, shown in the top right panel of Fig.
10.

The eigenspectrum shows that the variance of the continuum emission increases
towards larger wavelengths. Such behaviour (in the flux spectrum) is often
associated with thermal emission of a relatively hot source (dust), as
expected for a protostar. Discrete spectral features such as the
$\rm{Br\gamma}$, the CO bandhead features, and a relatively weak Na i doublet
are also observed with positive variance and, thus, correlated with the
continuum emission. In addition, the low-$J$ CO absorption lines are
identified as absorption features.

Figure 10: The tomogram (left panel) and the corresponding eigenspectrum
(right) of the Principal Component 1 of the PCA Tomography of W33A. The
variance of the tomogram is indicated in a rainbow colour scale, where black
indicates weak correlation (min), and red indicates strong correlations (max).
The contour-levels are chosen to highlight the extended emission and the
point-like source. The eigenspectrum exhibits the variance ($E_{1}$) as a
function of wavelength. The position of spectral features guiding the analysis
of the data are labelled and indicated by vertical dashed lines.

#### PC2

The higher-order PCs present the correlation and anti-correlation between the
spectral features in the eigenspectra, and their association with the
structures probed by the tomograms. Thus, for better highlighting the variance
gradient between distinct structures, the tomograms are shown in a red-to-blue
scale, where correlated (positive variance) and anti-correlated regions
(negative variance) are shown in blue and red, respectively.

In the top left panel of Fig. 11, the tomogram of PC2 exhibits the extended
nebulosity with positive variance (shown in blue), anti-correlated with a
ring-like structure around W33A exhibiting negative variance values (in red).
The eigenspectrum shows continuum and discrete spectral features that are
related to the structures observed in the tomogram. The continuum emission
shows positive variance values for $\lambda$ $\lesssim$ 2.32 µm, and is
negative for longer wavelengths. Combining the spectral and spatial
information, we interpret the continuum emission as dominated by the main
nebulosity at short wavelengths, while the structure around the compact source
is brighter at longer wavelengths. This behaviour can also be interpreted in
terms of the local extinction of the dust: the reddening is stronger towards
the compact source than that observed in the line-of-sight of the nebulosity.
This interpretation is also consistent with the RGB image shown in Fig. 8.

The top right and bottom panels of Fig. 11 also exhibits the H2 transitions,
$\rm{Br\gamma}$ and $\rm{Na\,\textsc{i}}$ lines correlated with the blue
continuum, indicating that these lines are produced in the main nebulosity and
absent in the compact source. Davies et al. (2010) presented spectra of the
individual bright regions shown in the nebulosity of Figure 11 (see their
Figs. 3 and 10). The three bright regions in the cavity from PC2 correspond to
the regions defined as B, C and D by D10: B is closest to the point source and
shows weaker low-$J$ CO absorption strength than the point source, and C and D
show yet weaker absorption of the narrow CO transitions. This is consistent
with the spots showing a reflected spectrum of the point source.

The narrow CO absorption lines detected at $\lambda$ $\gtrsim$ 2.32 µm in the
eigenspectrum appear as anti-correlated features, just as the continuum
emission at longer wavelengths. The CO absorption features are produced in a
relatively dense and cold environment, where the CO molecules are shielded
from the radiation of the protostar. These physical parameters are often
observed in the outer regions of circumstellar discs, consistent with their
association with the anti-correlated structure around the W33A protostar. The
absence of the CO bandhead features in the eigenspectrum of PC2 is also
remarkable. These features are produced at the inner region of the
circumstellar disc, at temperatures of around 2,000 K (Blum et al., 2004),
unresolved at the spatial resolution of the NIFS observations
($\sim$0$\aas@@fstack{\prime\prime}$15, or 360 AU at 2.4 kpc). Indeed, the PC2
tomogram exhibits a ring-like structure around the W33A protostar, with
variance values dropping to zero towards its central region, which is
consistent with the unresolved origin of the CO bandhead emission.

Figure 11: The tomogram (top left) and the corresponding eigenspectrum (top
right) of the Principal Component 2 of the PCA Tomography of W33A. In the left
panel, the contours of the PC1 tomogram are overlaid as black curves. For
orders $n$ $\geq$ 2, the tomogram is shown in a red-to-blue colour scale
indicated in the top right, where blue indicates the strongest correlations
and red indicates the strongest anti-correlations. In the right panel, the
eigenspectrum associated with PC1 is overlaid as the blue curve. The
horizontal dashed line indicates variance equal to zero. A detailed view of
the variances of selected spectral lines of the eigenspectrum are shown in the
bottom panels (from left to right: the H2 line at 2.12 µm, the $\rm{Br\gamma}$
feature at 2.16 µm, the CO (2-0) bandhead emission at 2.29 µm, the CO
absorption features at 2.32 µm, and the $\rm{Na\,\textsc{i}}$ doubled at 2.20
µm).

#### PC3

The tomogram of PC3 exhibits a spatially resolved structure around the compact
source, with coherent changes between positive (blue) and negative variance
values (red) at the central position of the point-like object. In extra-
galactic objects, this pattern is often identified in rotating structures,
such stellar discs or torus around super-massive black holes (e.g. Ricci et
al., 2014). In the context of protostellar objects, it is likely arising from
a rotating circumstellar disc, as identified by Davies et al. (2010) when
analysing the CO absorption features (see their Fig. 14). The rotation axis of
the structure has a position angle of about $-$45∘ (from N to E).

The eigenspectrum of PC3 exhibits no variance of the continuum, but shows a
series of narrow features for $\lambda$ $\gtrsim$ 2.3 µm, likely arising from
the low-J CO lines (see details in the bottom panel of Fig. 12). The blue-
shifted components are correlated and so, associated with the positive
variance values in the tomogram (shown as blue) while the red-shifted
components are associated with negative values (shown as red), indicating the
disc is rotating from NE to SW direction.

A similar kinematic pattern is also observed for the $\rm{Br\gamma}$ feature
(also indicated in the bottom panel of the figure), suggesting that the
rotating disc exhibits an ionised component (see discussion in Section 4.4).

Figure 13 shows a closer view of the rotating structure identified in the
tomogram. The rotation axis of the disc is observed with a position angle of
PA = 140 $\pm$ 10∘. A radial profile of the variance along the disc plane is
presented in the bottom panel of the figure, showing that the disc kinematic
signature is only observed within the inner
$\pm$0$\aas@@fstack{\prime\prime}$4 region around W33A (about $\sim$3 times
larger than the resolution of the observations). Even though the CO absorption
features are observed in the extended nebulosity, the principal component 3 is
only showing the regions exhibiting ordered kinematics on these lines,
restricted to the disc. This analysis suggests that the ordered kinematics are
observed in a much more compact region than previously reported by D10.

Figure 12: The tomogram (top left) and the corresponding eigenspectrum (top
right) of the Principal Component 3 of the PCA Tomography of W33A. For
details, see Figs. 10 and 11.
Figure 13: A detailed view of the disc-like structure probed by the PC3
tomogram. Top panel: The normalised variance of the tomogram is shown in a
divergent red-to-blue colour-scale. The position of the compact source is
indicated by the black contours. The red arrow and blue dashed lines indicate
the position angle of the rotation axis of the disc and its error,
respectively, estimated as PA = 140 $\pm$ 10∘. The dashed black line shows the
direction from which the radial variance profile was extracted. Bottom: the
radial variance profile extracted along the plane of the disc. The points and
error bars correspond to the mean variance and the 1-$\sigma$ error of a
3-pixel width region perpendicular to the sampled spatial direction. The
dashed vertical lines indicate the region where the absolute value of the
variance is larger than zero ($|r|$ $\leq$ 0$\aas@@fstack{\prime\prime}$2).

#### PC4

The tomogram of PC4 exhibits a correlated structure roughly perpendicular to
the rotating disc identified in PC3 (see Fig. 12).

The eigenspectrum of PC4 shows both continuum and discrete features, as
observed in PC2 (Fig. 11). The continuum emission between 2.17 and 2.35 µm is
anti-correlated with the emission at shorter and at longer wavelengths. We
interpret this result as if the same region is associated with at least three
distinct processes: $i)$ the dust content is reflecting the radiation emitted
by the protostar at short wavelengths; $ii)$ the thermal emission of the dust
is observed at longer wavelengths; $iii)$ the same region exhibits the CO
bandhead features in emission, and the profiles are quite narrow when compared
to the emission observed in PC1. This may indicate that we observe this
emission moving perpendicularly to the line-of-sight.

In addition, the eigenspectrum of PC4 shows a strong $\rm{Br\gamma}$ feature
in emission, suggesting the presence of partially-ionised gas that could also
contribute as a source of scattering electrons producing the reflected light
at shorter wavelengths.

Figure 14: The tomogram (top left) and the corresponding eigenspectrum (top
right) of the Principal Component 4 of the PCA Tomography of W33A. For
details, Figs. 10 and 11.

#### PC5

The tomogram of PC5 exhibits an extended emission with a different morphology
to the main nebulosity of W33A (shown by the black contours), suggesting that
this structure is either produced by a distinct mechanism or exhibits distinct
physical conditions.

Indeed, the eigenvector exhibits at least six $\rm{H}_{2}$ transitions
correlated with the extended emission shown in the tomogram, indicating that
PC5 is either tracing the base of a disc wind, or the cavity of the large-
scale outflow shown in Fig. 8, observed with a position angle (PA) of 145∘
(from N to E), measured along the brightest and more elongated region of the
molecular emission. In addition, the spatial distribution of the $\rm{H}_{2}$
emission probed by the tomogram of PC5 is similar to the continuum-subtracted
map of the $\rm{H}_{2}$ (1-0) S(1) transition at 2.12 µm presented by D10 (see
their Fig. 15). The most significant difference between the tomogram and the
line map is that the first considers the emission of all the $\rm{H}_{2}$
transitions available in the spectral range (producing a higher signal-to-
noise map) while the second method is based only on the intensity of a single
line (associated with a relatively larger noise).

In addition, the $\rm{Br\gamma}$ feature and CO bandhead emission are anti-
correlated with the H2 lines, indicating they arise from the compact region at
the launching point of the jet (indicated as red in the tomogram). In general,
the CO features are produced in a dense region (e.g. the circumstellar disc)
while the $\rm{H}_{2}$ is excited through shocks in a more diffuse medium.
Therefore, PC5 also shows the spatial separation of regions with distinct
density regimes.

Figure 15: The tomogram (top left) and the corresponding eigenspectrum (top
right) of the Principal Component 5 of the PCA Tomography of W33A. For
details, Figs. 10 and 11.

###  3.4 Derivation of the physical parameters of the structures associated
with the W33A protostar

The following results were obtained by analysing the data cube reconstructed
using the first five Principal Components of the PCA Tomography (Sect. 2.4.2).
We also compare our findings with those reported by D10 to check whether the
post-processing of the NIFS data cubes led to a significant improvement on the
determination of the physical parameters.

Figure 16 presents a composite RGB image of the NIFS FOV similar to the one
presented in Fig. 9. The spectra shown in the right panel were integrated over
the regions indicated by the contours: the compact K-band object (red), the
disc wind/cavity of the outflow probed by the molecular $\rm{H}_{2}$ emission
(green) and the main nebulosity (blue).

Figure 16: Left panel: False-colour RGB image of the NIFS FOV (left, blue:
2.05µm, green: 2.2µm, red:2.35µm) recreated using PCs 1–5. The contours
indicate the regions where the integrated spectra were extracted. Right:
integrated K-band spectra of the W33A point-like source (in red), the main
nebulosity (in blue), and the $\rm{H}_{2}$ wind (in green).

The integrated spectra shows that the main nebulosity is brighter than the
protostar at shorter wavelengths, while the embedded object gets brighter for
$\lambda$ $\gtrsim$ 2.35 µm. This behaviour was also observed in the analysis
of PC2 (see Fig. 11), showing that the continuum emission is correlated with
the nebulosity at short wavelengths due to the reflection of the radiation
from the embedded protostar, while it is correlated with the compact source
for longer wavelengths likely due to local extinction effects.

In addition, the spectra of the main nebulosity (in blue) and the $\rm{H}_{2}$
wind (in green) exhibit at least five ro-vibrational transitions of the
$\rm{H}_{2}$ molecule. The $\rm{H}_{2}$ emission observed in the spectrum of
the main nebulosity is likely due to the overlap between this region and the
$\rm{H}_{2}$ wind. The most energetic transition corresponds to the H2 2–1
S(1) line at 2.2477 µm, with upper-level energy of 12,550 K. The absence of
transitions with larger energy indicates that the physical conditions of the
$\rm{H}_{2}$ gas (i.e. temperature and column density) are not sufficient to
produce more energetic transitions. We further investigate the physical
parameters of the molecular $\rm{H}_{2}$ emission in Sect. 3.4.2.

####  3.4.1 The nature of the YSO candidates in the vicinity of W33A

We performed the photometry of the UKIDSS $JHK$ images using the IRAF/DAOPhot
package to better characterise the protostellar object candidates in the
vicinity of W33A (Sources B and C). We used the 1.4′ $\times$1.4′ FOV shown in
the top panel of Fig. 7, containing enough point-like sources to obtain a good
photometric calibration by cross-matching with the UKIDSS catalogue.

First, we obtained an average map of the individual $J$, $H$, and $K$ images
to identify the positions of the sources that appears in, at least, one of the
three filters. We obtained a list of 298 point-like objects with peak
intensities above a 3-$\sigma$ threshold. We performed the photometry of the
individual images, using an aperture with radii of $r$ = 1″, and a ring with
inner and outer radius of 1$\aas@@fstack{\prime\prime}$6 and
2$\aas@@fstack{\prime\prime}$0 to estimate the sky background. The fluxes were
converted to instrumental magnitudes, and the photometric calibration was done
by cross-matching the sources detected in all the three bands ($J,H$ and $K$)
with the UKIDSS catalogue, leading to a list of 132 objects. A final list
containing the photometry of 287 objects, was obtained after removing the
sources detected only in a single filter, excluding, for example, any bright
K-band knot associated with the large-scale outflow in the FOV from the
initial list.

The $JHK$ magnitudes of the W33A, Source B and Source C are listed in Table 2,
and the full list containing the 287 sources is available as on-line material
(The first ten rows are listed in Table 4. W33A was not detected in the $J$\-
and $H$-band and, thus, the reported values correspond to upper limits, as
well as the $J$-band magnitude of Source C. In addition, the reported $J$-band
magnitude of Source B corresponds to a foreground star in the line-of-sight of
the object detected in the $H$\- and $K$-bands.

Table 2: $JHK$ photometry of the protostellar object candidates in the
vicinity of W33A.

Source | J (mag) | H (mag) | K (mag)
---|---|---|---
W33A | 20.77∗ | 15.05∗ | 9.51 $\pm$ 0.19
Source B | –+ | 13.76 $\pm$ 0.09 | 10.30 $\pm$ 0.20
Source C | 20.77∗ | 16.31 $\pm$ 0.11 | 12.69 $\pm$ 0.23

Notes: $\ast$ indicates upper limit due to the non-detection of a point-like
source at the given position in the corresponding filter. $+$ a foreground
source lying very close to the object exhibits $J$ = 15.45 $\pm$ 0.08 mag.

Fig. 17 presents the $K$ versus $H-K$ colour-magnitude diagram (top panel) and
the $J-H$ versus $H-K$ diagram (bottom) of the sources within the 1.4′
$\times$1.4′ FOV around W33A. The colour-magnitude diagram indicates that the
three sources (W33A, Source B and Source C) are relatively bright K-band
objects ($K$ $<$ 13 mag), exhibiting $H-K$ colour indices consistent with
embedded young objects ($H-K$ $>$ 3 mag). W33A is the most reddened source,
with $H-K$ = 5.54 $\pm$ 0.21 mag, followed by Source C ($H-K$ = 3.62 $\pm$
0.26 mag) and Source B ($H-K$ = 3.46 $\pm$ 0.22 mag). These objects are
located among the brightest K-band sources of a sequence of reddened objects,
confirming that they are deeply embedded within their circumstellar
environments as expected for protostars.

The bottom panel of Fig. 17 presents a colour-colour diagram based on the
$J-H$ and $H-K$ colour indices. Given that the $J$-band magnitude of the three
sources corresponds to upper limits of their real values, their position in
the $y$-axis of the diagram correspond to lower limits. Despite of that, the
three sources are consistent with reddened objects, corroborating the young
nature of these objects.

Figure 17: $K$ versus $H-K$ colour-magnitude diagram (top panel) and the $J-H$
versus $H-K$ colour-colour diagram (bottom) of the 1.4′ $\times$ 1.4′ region
around W33A. W33A is indicated as the red diamond symbol, Source B as the blue
triangle and Source C as the green squared symbol. The photometric errors are
shown as error bars and the arrows indicates lower limits of the actual values
(see text).

In Fig. 18, we present the spectrum of W33A (shown as the red curve) and
Source B (in blue), integrated over a 0$\aas@@fstack{\prime\prime}$15 radius,
each with the extended nebular emission subtracted using an annular region
around each compact sources, with inner and outer radii of
0$\aas@@fstack{\prime\prime}$20 and 0$\aas@@fstack{\prime\prime}$30,
respectively (indicated by the white contours in Fig. 18). We note that W33A
has a steeper continuum due to the emission of the dust, while Source B has a
blue continuum. For comparison, we also show the integrated spectrum of Source
B extracted using the original data cube (shown as the green curve), before
the post-processing of the NIFS data cube and its reconstruction (see details
in Sect. 2.3). The comparison between the spectrum of Source B before and
after the reconstruction of the data cube shows the effectiveness of the noise
suppression due to the post-processing procedures applied to the data cube. In
addition, the K-band spectrum of Source B exhibits relatively weak emission of
the $\rm{Br\gamma}$ and the $\rm{H}_{2}$ transition at 2.12 µm. These features
are typically observed in YSOs, corroborating the photometric results from
Fig. 17.

Figure 18: Left panel: False-colour image of the NIFS FOV integrated from 2.03
to 2.04 µm. The point-like source is saturated to enhance the contrast between
W33A and the fainter, K-band counterpart of Source B at the SE edge of the
datacube. The black contours indicate the regions where the integrated spectra
were extracted. The nebular emission, evaluated over the regions indicated by
the white contours, was subtracted from the spectra. Right: integrated K-band
spectra of the W33A point-like source (in red), and the Source B before
(green) and after the post-processing routines applied to the NIFS datacube
(in blue). The spectra of Source B were multiplied by a factor of 400.

####  3.4.2 The physical parameters of the molecular hydrogen gas

In Figure 19 we present the velocity and the dispersion map of the
$\rm{H}_{2}$ (1–0) S(1) transition at 2.12 µm, the brightest $\rm{H}_{2}$
transition identified in the K-band NIFS data cube. The radial velocity was
corrected by the rest velocity of W33A ($V_{\rm LSR}$ = 36.7 km s-1, Lumsden
et al. 2013). The top panel of Fig. 19 shows that the bulk of the $\rm{H}_{2}$
emission is located at a radial velocity of $-$56 km s-1, showing no
significant velocity gradient in the FOV. In addition, no substantial change
in the full-width at zero intensity (FWZI) of the line profile is observed
towards the emitting region, exhibiting values around 110 km s-1 (middle
panel). The bottom panel of Fig. 19 shows the high-velocity of the
$\rm{H}_{2}$ emission profile, evaluated as the terminal velocity at the blue
side of the $\rm{H}_{2}$ profile (estimated from the FWZI parameter). The
velocity map exhibits relatively higher velocity values close to the W33A
protostar ($\sim$ $-$130 km s-1) and is roughly constant at large distances
($\sim$ $-$113 km s-1).

Figure 19: Velocity map of the $\rm{H}_{2}$ (1–0) S(1) transition at 2.12 µm.
The top panel presents the centroid velocity of the emission,the middle panel
shows the full-width at zero intensity of the line profile, and the bottom
panel indicates the terminal velocity at the blue side of the $\rm{H}_{2}$
emission. The black contours corresponds to 10%, 25%, 50% and 75% of the peak
intensity of the $\rm{H}_{2}$ emission.

We further measured the fluxes of the individual $\rm{H}_{2}$ transitions in
the integrated spectrum of the $\rm{H}_{2}$ cavity (shown as the green region
in Fig. 16) to construct a ro-vibrational Boltzmann diagram of the
$\rm{H}_{2}$ molecule and, thus, to evaluate the physical parameters of the
jet. Table 3 summarises the fluxes and the physical properties of each
$\rm{H}_{2}$ transition. Figure 20 presents the Boltzmann diagram of the
$\rm{H}_{2}$ emission, constructed using the methodology adopted by Caratti o
Garatti et al. (2006, 2015) and briefly explained below.

Table 3: H2 transitions identified in the spectrum of the cavity of the
outflow.

Transition | Flux | $\sigma_{F}$ | $\lambda$ | $g$ | $E_{\rm upper}$ | $A$
---|---|---|---|---|---|---
| (10-19 erg cm-2 s-1) | (µm) |  | (K) | (10-7 s-1)
1-0 S(2) | 8.51 | 0.36 | 2.0338 | 9 | 7584 | 3.98
1-0 S(1) | 30.19 | 0.40 | 2.1218 | 21 | 6956 | 3.47
1-0 S(0) | 9.09 | 0.36 | 2.2235 | 5 | 6471 | 2.53
2-1 S(1) | 5.16 | 0.40 | 2.2477 | 21 | 12550 | 4.98
1-0 Q(1) | 58.49 | 0.44 | 2.4066 | 9 | 6149 | 4.29

Notes: the energy levels are from Dabrowski (1984), and the Einstein
coefficients are from Turner et al. (1977).

Figure 20: Boltzmann diagram of the $\rm{H}_{2}$ gas associated with W33A. The
dashed line indicates the best fit to the data. The corresponding temperature,
total column density and reddening values are indicated in the upper-right
corner.

The upper-level column density of each $\rm{H}_{2}$ transition
($N_{\upsilon,J}$, where $\upsilon$ and $J$ denotes the upper vibrational and
rotational levels, respectively) is evaluated as:

$N_{\upsilon,J}=\frac{4\pi\Omega}{h}\frac{F_{\upsilon,J}}{A_{\upsilon,J}\nu_{\upsilon,J}g_{\upsilon,J}}$
(5)

where $h$ is the Planck constant (in erg s), $\Omega$ is the solid angle of
the emitting region (in sr), $F_{\upsilon,J}$ is the integrated intensity of
the transition (in erg cm-2 s-1 sr-1), $A_{\upsilon,J}$ corresponds to the
Einstein coefficient of the transition (in s-1), $\nu$ is the frequency (in
Hz), and $g_{\upsilon,J}$ is the statistical weight of the transition. The
upper-level column density of each $\rm{H}_{2}$ transition is plotted against
their upper-level energy ($E_{\upsilon,J}$, in K). Assuming the gas is in
Local Thermodynamic Equilibrium (i.e. follows the Boltzman distribution), a
linear fit to the data in a log-normal scale (i.e. $\log(y)$ =
$\alpha\,x+\beta$) provides an estimate of the excitation temperature
($T_{\rm{exc}}$ = $-1/\alpha$), and the total column density [$N_{\rm{tot}}$ =
$Q(T_{\rm{ex}})\exp(\beta)$] of the $\rm{H}_{2}$ gas. The K-band extinction
($A_{\rm{K}}$) can also be inferred by choosing the best $A_{\rm{K}}$ value
that maximises the correlation between the quantities shown in the Boltzmann
diagram (for details, see Caratti o Garatti et al. 2015).

The best-fit to the Boltzmann diagram shows that the molecular gas is
associated with a $T_{\rm{exc}}$ value of (2.3$\pm$0.1)$\cdot$103 K, and a
total $\rm{H}_{2}$ column density of $\log$($N_{\rm{tot}}$/cm-2) =
18.8$\pm$0.1 dex. The K-band reddening is about 4.20$\pm$0.1 mag, which is
equivalent to an $A_{V}$ of 63$\pm$2 mag assuming the $A_{\rm V}/A_{\rm K}$
ratio of 14.95 from Damineli et al. (2016).

The origin of the $\rm{H}_{2}$ emission can be inferred using the ratio of
$\rm{H}_{2}$ transitions such as the (1–0) S(1)/(2–1) S(1). Based on the
fluxes reported in Table 3 and correcting them for reddening, we obtained a
(1–0) S(1)/(2–1) S(1) ratio of 9.3 $\pm$ 1.4. The observed ratio is close to
the theoretical value of $\sim$10 for thermal excitation of the $\rm{H}_{2}$
emission through shocks, and significantly larger than the predicted value of
1.8 for the case of radiative excitation through UV photons from the protostar
(Black & Dalgarno, 1976).

By taking advantage of the spatial information available in the NIFS data, we
implemented a pixel-to-pixel evaluation of the physical parameters of the
$\rm{H}_{2}$ gas, as presented in Fig. 21. Spatial maps of $A_{\rm{K}}$,
$T_{\rm{exc}}$ and $N_{\rm{tot}}$ were obtained by constructing and fitting
the Boltzmann diagram of the $\rm{H}_{2}$ transitions for the spaxels
exhibiting fluxes above a 3-$\sigma$ threshold for all the five $\rm{H}_{2}$
transitions listed in Table 3. In addition, the mass of the $\rm{H}_{2}$ gas
($M_{H2}$) was computed based on the total column density of the $\rm{H}_{2}$
molecules as:

$M_{\rm H2}(x,y)=2\mu m_{H}A_{\rm pix}N_{\rm tot}(x,y)$ (6)

where the (x,y) indices correspond to the position of a given spaxel of the
data cube, $\mu$ = 1.24 is the mean atomic weight of the $\rm{H}_{2}$, $m_{\rm
H}$ is the proton mass (in g) and $A_{\rm pix}$ corresponds to the projected
area of the spaxel (in cm2).

Figure 21: Pixel-to-pixel derivation of the $\rm{H}_{2}$ parameters for those
pixels exhibiting fluxes above a 3-$\sigma$ threshold. From top left to bottom
right: the $A_{\rm{K}}$, $T_{\rm{exc}}$, column density and mass of the
$\rm{H}_{2}$ gas. The black contours indicates the integrated H2 (1–0) S(1)
flux. The red contours indicate the position of the point-like source and its
extended continuum emission. The range of each parameter was chosen in order
to improve the visualisation of the maps. The mean value of the $A_{\rm{K}}$
and $T_{\rm{exc}}$, and the total column density and mass of the $\rm{H}_{2}$
gas are shown at the bottom of their respective panels.

The top left panel of Fig. 21 presents the distribution of the K-band
extinction over the $\rm{H}_{2}$ emission associated with the cavity of the
outflow. The median K-band extinction of the $\rm{H}_{2}$ emission is
$A_{\rm{K}}$ = 4.1$\pm$0.2 mag, corresponding to an $A_{\rm{V}}$ of $\sim$61
$\pm$ 3 mag (Damineli et al., 2016).

The map of the excitation temperature of the gas is presented in the top right
panel of Fig. 21. The median $T_{\rm{exc}}$ value is (2.2 $\pm$ 0.1) $\cdot$
103 K. The total column density of the $\rm{H}_{2}$ gas is presented in the
bottom left panel. The distribution of $N_{\rm{tot}}$ is similar to the flux
of the H2 (1–0) S(1) line at 2.12 µm, indicated by the contours, suggesting a
linear relationship between the brightness of this transition and the
logarithm of the column density. The total column density of the gas is about
$\log$($N_{\rm{tot}}$)$\sim$ 18.7$\pm$0.3 dex. The bottom right panel of Fig
21 shows the mass of the $\rm{H}_{2}$ gas, computed using Eq. (6). The total
mass of the emitting $\rm{H}_{2}$region is about (5.4$\pm$0.1) $\cdot$ 10-5
$\rm{M}_{\odot}$.

A rough estimate of the dynamical timescale of the cavity of the outflow,
probed by the $\rm{H}_{2}$ emission, can be evaluated as

$t_{\rm dyn}=\frac{\ell_{\rm H2}}{|v_{\rm
H2}|}=\frac{\cos{(i)}}{\sin(i)}\frac{\ell_{\rm proj}}{|v_{\rm
rad}|}\approx\frac{4.74}{\tan(i)}\frac{\ell_{\rm proj}\,(\rm{AU})}{|v_{\rm
rad}|\,(\rm{km\,s^{-1}})}\,\rm{yr}$ (7)

where $\ell_{\rm proj}$ is the projected length of the emitting region in the
plane of the sky ($\ell_{\rm proj}$ = $\ell_{\rm H2}\cdot\sin(i)$), $v_{\rm
rad}$ is its radial velocity of the emission ($v_{\rm rad}=v_{\rm
H2}\cdot\cos(i)$), and $i$ is the inclination angle of the structure ($i$ =
30∘). The $\rm{H}_{2}$ emission extends up to 1$\aas@@fstack{\prime\prime}$65
from the driving source, corresponding to a projected linear scale of 3960 AU
at the distance of 2.4 kpc. The emission peaks at 2.1216 µm, corresponding to
a velocity of about $-$28.3 km s-1. Considering the rest velocity of W33A,
$v_{lsr}$ = 36.7 km s-1 (Lumsden et al., 2013), the projected radial velocity
of the $\rm{H}_{2}$ emission is about $-$65 km s-1, in agreement with the
analysis presented in Fig. 19. Using Eq. (7) and propagating the error on each
input parameters333We considered a typical $\pm$20∘ on the inclination angle,
$\pm$29 km s-1 on the velocity (half of the velocity resolution of the NIFS
data), and an error of 10% on the projected length, we obtained a dynamic
timescale of about (5.0 $\pm$ 2.4)$\cdot$102 years, suggesting the
$\rm{H}_{2}$emission arises from a relatively young structure. We further
obtained a rough estimate of the mass-loss rate of the ejected gas as
$\dot{M}_{\rm H2}$ = $M_{\rm H2}$/$t_{\rm dyn}$ $\sim$ (1.1$\pm$0.5) $\cdot$
10-7 $\rm{M}_{\odot}$ yr-1.

####  3.4.3 Spectro-astrometry of the Brackett-Gamma feature

According to the PCA results, the $\rm{Br\gamma}$ emission is observed towards
the W33A protostar (PC1) and the extended nebulosity to the south (PC1 and
PC2), the rotating disc associated with W33A (PC3), and compact emitting
regions closer to the W33A protostar (PC4 and PC5).

Davies et al. (2010) reported that the $\rm{Br\gamma}$ emission was probing
the cavity of the large-scale K-band outflow at very small scales using the
spectro-astrometry of the $\rm{Br\gamma}$ line, a technique that combines the
spectral and spatial information of the data cube to resolve kinematic
structures at sub-pixel scales with impressive accuracy (i.e. one tenth of a
milliarcsecond). Thus, we adopted the same methodology presented by D10 to
check if the post-processing of the NIFS data led to a significant improvement
of the data at sub-pixel scales.

The spectro-astrometry of the $\rm{Br\gamma}$ feature of W33A is presented in
Fig. 22. The left panel presents the offset of the centroid in the right
ascension (black curve) and declination axes (red curve) as a function of the
wavelength (top panel), together with the $\rm{Br\gamma}$ flux at the same
wavelength scale (bottom). The right panel presents the spatial variation of
the centroid position as a function of the velocity, matching the colour scale
of the $\rm{Br\gamma}$ profile shown in the left panel. The positional error
of 0.12 mas, indicated in the top left corner, was based on the centroid
position offset of the nearby continuum around the $\rm{Br\gamma}$ line
(corresponding to the spectral regions shown in red).

Figure 22: Spectro-astrometry of the $\rm{Br\gamma}$ feature. Top left panel:
the centroid position in the right ascension (RA, in black) and declination
axis (Decl., in red) as a function of the wavelength. The baseline is
indicated by a horizontal grey line. The peak of the $\rm{Br\gamma}$ feature
(at 0 km s-1) is indicated by a vertical line. Bottom left: the normalised
$\rm{Br\gamma}$ flux. The coloured region indicates the spectral interval
shown on the right panel. The spectral regions used to estimate the positional
error of the centroid (about $\sim$0.12 mas in each direction) are indicated
in red. Right panel: Spatial variation of the $\rm{Br\gamma}$ centroid
position. The velocity with respect to the peak of the emission is indicated
by the horizontal colour bar shown in the bottom region of the plot. The
positional error is shown in the left top corner of the plot. The dashed line
indicates the position-angle (PA) of the red and blue high-velocity
components.

The blue and red high-velocity components (corresponding to velocities of
$-$300 to +300 km s-1) are offset by roughly $\sim$ 0.55 mas, from the south-
east to the north-west direction, corresponding to $\approx$ 1.3 AU at the
distance of 2.4 kpc. The displacement between the blue and red high-velocity
components roughly follows a linear distribution with a position angle (PA) of
154$\pm$26∘ (from N to E), which is compatible with the results reported by
D10(156$\pm$25∘444equivalent to the value reported by the authors, PA =
113$\pm$25∘, oriented from E to N).

####  3.4.4 Kinematics of the CO absorption features

We used the Penalized Pixel-Fitting algorithm (pPXF Cappellari & Emsellem,
2004) to derive the velocity map of the CO absorption features observed in the
spectrum of W33A. The pPXF uses a template spectrum to fit and evaluate the
radial velocity at each spaxel in the data cube. The procedure starts with an
initial guess of the velocity and width parameters and solves a Gauss-Hermite
function based on the template spectrum to fit the observed one. The residuals
between the model and the observed spectra are perturbed using a penalised
function and the best model is chosen based on a $\chi^{2}$ minimisation
procedure, delivering a set of parameters (V, $\sigma$, $h_{3}$, … $h_{m}$)
and their uncertainties for each spaxel of the datacube. For this work, we
focused the further analysis on the interpretation of the radial velocity
measurements only. The brightest spaxel of the continuum emission at $\sim$2.3
µm (at the central position of the point-source object) was adopted as the
template spectrum (corresponding to the rest-velocity at 0 km s-1) to derive
the kinematics of the CO absorption lines using the spectral region ranging
from 2.32 to 2.37 µm, containing both the low-J CO features of the R and the P
branches. To spatially sample the rotating structure, the calculations were
performed only for the spaxels within a 0$\aas@@fstack{\prime\prime}$5 radii
around the central position of W33A.

Figure 23 presents the results of the pPXF analysis. The left panel indicates
the map of radial velocities overlaid with contours of the continuum emission.
The map indicates a clear gradient between negative and positive radial
velocities very similar to the pattern observed in the tomogram of PC3 (see
Figs. 12 and 13), suggesting that the kinematics arise from a rotating disc as
reported by D10. The rotating axis of the structure is aligned with a position
angle of about 150 $\pm$ 10∘ (from N to E), compatible with the orientation of
the rotating axis observed in the tomogram of PC3 (see Fig. 13). The velocity
gradient arises from a relatively compact region of
$\sim$0$\aas@@fstack{\prime\prime}$5, corresponding to a projected length of
about 1,200 AU at the distance of 2.4 kpc.

Figure 23: Top panel: Radial velocity map of the CO absorption features
towards W33A, derived with the pPXF fitting. The black contours are placed at
the 3${}^{n}-\sigma$ level ($n$ = $0,1,2,...$) of the integrated continuum
emission at $\sim$2.3 µm. The velocity scale is presented on the right side of
the map. The dashed black line indicates the position angle of 150∘ (from N to
E) of the rotation axis of the disc, following the velocity gradient between
the blue- and the red-shifted disc components. The filled black line is
perpendicular to the disc axis, indicating the direction where the position-
velocity diagram was extracted. Middle: The position-velocity diagram for the
CO absorption features. The black circles and their errors corresponds the
measurements using the pPXF analysis. The best Keplerian rotation curve for
$M\sin(i)=8.7$ $\rm{M}_{\odot}$ is indicated by the solid blue line. The
dotted and dashed blue lines indicate the Keplerian curves considering the
upper and lower 1-$\sigma$ errors on the mass, respectively. Bottom panel: The
cumulative distribution function (CDF) of the $M\sin(i)$ value evaluated for
10,000 Monte-Carlo simulations is shown as the solid blue curve. The relative
histogram of $M\sin(i)$ values is shown as the solid black line. The solid
vertical line corresponds to the median value and the dashed lines are placed
at the 16 and 84% confidence values, corresponding to the 1-$\sigma$ error
bars.

We used the velocity map delivered by the pPXF analysis to extract a position-
velocity (PV) diagram perpendicularly to the rotation axis of the disc. The
results are presented in the middle panel of Fig. 23, exhibiting the mean
value and the standard deviation of the velocity at each position, evaluated
within 3-pixel wide bins across the extraction axis (indicated by the filled
black line in the top panel of Fig. 23). The PV diagram indicates a clear
separation between the blue- and red-shifted components (with peak velocities
around $+$ 5 and $-$ 4 km s-1, respectively) separated by roughly
0$\aas@@fstack{\prime\prime}$2\. The best Keplerian model, shown as the solid
blue curve, indicates that the total mass of the “protostar+disc” system is
about 8.7${}^{+1.8}_{-1.6}$ $\rm{M}_{\odot}$. We note that the velocity falls
more rapidly than any Keplerian model at larger radii. The mass and its
standard deviation were evaluated using a Monte-Carlo method (see results in
the bottom panel of Fig. 23), varying the input data randomly within their
errors. We used a total of 10,000 simulations to achieve the reported values.
Assuming the inclination angle of $i$ = 60∘ adopted by D10 (estimated from
interferometric mid-IR observations from de Wit et al., 2010), the
inclination-corrected mass of the “protostar+disc” system in W33A is
10.0${}^{+2.1}_{-1.8}$ $\rm{M}_{\odot}$.

##  4 Discussion

This work presents a re-analysis of the $K$-band IFU observations of the well-
known high-mass protostar W33A with an eye to further extend the analysis
presented originally by Davies et al. (D10) through application of a Principal
Components Analysis.

First, we performed the flux-calibration of the data cube followed by the
differential atmospheric refraction (DAR) correction. The DAR effects
correspond to a slight positional offset of the source as a function of the
wavelength (see Fig. 1). The proper correction of such effects is necessary to
properly compare the spatial information at distinct wavelengths in the data
cube, and to perform measurements at a very small fraction of a pixel as
presented in this work; see also D10.

Then, we performed the PCA tomography of the data cube, and found that most of
the astrophysical information of the data cube is well-explained by the first
five Principal Components (PCs): PC1 shows that the embedded compact source
associated with W33A dominates the emission in the K-band; PC2 presents the
differences between the compact and the extended emission regions, mostly
showing an anti-correlation introduced by reddening effects; PC3 exhibits the
kinematics of a rotating circumstellar disc around W33A, probed by the low-$J$
CO absorption lines and exhibiting a contribution of ionised gas; PC4 probes
the cavity of the jets, showing a roughly perpendicular structure aligned to
the rotating axis of the disc, which is likely reflecting the emission of the
inner region of the disc and the protostar itself; and PC5 exhibits the anti-
correlation between the extended $\rm{H}_{2}$ emission and the $\rm{Br\gamma}$
and CO emission arising from the compact object.

Along with the PCA results, the analysis of the near-infrared $JHK$ maps of
W33A allowed us to obtain an overall interpretation of the environment
associated with the protostar and its surroundings, as illustrated in Fig. 24.

Figure 24: Interpretation of the near-infrared $JHK$ maps and NIFS K-band
observations of W33A. The illustration roughly follows the same orientation of
the observations (North is up and East is to the left). The structures are
labelled and shown in different colours, as indicated in the top left corner.

The data cube was reconstructed using the first five PCs, suppressing a
considerable fraction of noise associated with PC6 and higher-orders PCs (see
Scree test in Fig. 5). Then, we followed a similar analysis as that performed
by D10 to validate our results using the reconstructed data cube. In general,
we found that the post-processing routines and noise suppression of the data
cube led a significant improvement of the results, as discussed below.

###  4.1 Multiple IR sources within W33A

The comparison of the large-scale NIR maps with the NIFS FOV in Fig.s 7 and 8
shows that different sources are likely contributing to the observed extended
K-band emission detected in the field.

In Fig. 11, PC2 shows the anti-correlation between the circumstellar region of
W33A (probed mostly by the resolved CO absorption features) and the extended
emission mostly probed by the continuum and $\rm{Br\gamma}$ emission. The
tomogram of PC2 exhibits a clear separation between the extended emission
closer to W33A and the larger structure towards the S.

A point-like source is observed at the SE edge of the NIFS FOV, consistent
with the position of the Source B (see Fig. 7). At the distance of W33A, the
linear separation of W33A and Source B is about $\sim$7000 AU. In addition,
the detection of a large-scale bipolar jet in the E-W direction associated
with the $H$-band compact Source C (see Fig. 7) indicates that multiple
sources are located within this region. The $JHK$ photometry and the colour-
magnitude and colour-colour diagrams presented in Fig. 17 allowed us to
confirm that Source B and C are compatible with young, protostellar objects.

###  4.2 The mass of W33A and the compact circumstellar disc

We obtained a robust estimate of the total mass of the “disc+protostar” system
and its error based on the pPXF analysis of the CO absorption features at 2.35
µm (see Fig. 23). The best fit of the position-velocity diagram is consistent
with a Keplerian rotation structure around a total mass of M$\sin(i)$ =
8.7${}^{+1.8}_{-1.6}$ $\rm{M}_{\odot}$, assuming a distance of 2.4 kpc.
Assuming the inclination angle of 60∘ as adopted by D10, the mass is
10.0${}^{+2.1}_{-1.8}$ $\rm{M}_{\odot}$. If considering a typical $\pm$20∘
uncertainty on the inclination angle, the upper-limit of the mass increases by
a factor of two, 10.0${}^{+4.1}_{-2.2}$ $\rm{M}_{\odot}$, indicating that the
inclination uncertainty is a significant source of error on the total mass.
The inclination angle value of 60∘ adopted in D10 and in this work was derived
from the analysis of interferometric mid-IR observations from de Wit et al.
(2010). Those authors argued that the derived inclination angle of the system
is consistent with the fact that W33A exhibits no strong $H-$band emission
arising from the central star (as one would expect for a less inclined
system), and it is also in agreement with the relatively weak large-scale
K-band redshifted outflow lobe (Fig. 7). Ilee et al. (2013) reported a more
accurate determination of the inclination angle of the W33A system as
$i(^{\circ})$ = 37${}^{+16}_{-11}$ based on the modelling of the CO bandhead
emission at 2.3 µm using high spectral resolving power observations
(R$\sim$30,000). More recently, Fedriani et al. (2020) analysed the K-band
emission of the high-mass YSO IRAS 11101$-$5829\. Those authors reported that
the CO bandhead emission exhibits low velocities, likely arising from a face-
on rotating disc (i.e. high inclination angle). However, the geometry of
$\rm{H}_{2}$ the jet driven by the IRAS source requires that the disc must be
roughly edge-on, with low inclination angle. They conclude that most of the CO
emission is observed as reflected and scattered light in the outflow cavity
walls, and therefore the inclination of the system can be incorrectly
interpreted based on the modelling of the CO bandhead emission. In addition,
the identification of the low-$J$ CO features in absorption in the spectrum of
W33A (e.g. see Fig. 9) is also indicative of an inclined disc and, therefore,
not compatible with such a small inclination value as that reported by Ilee et
al. (2013).

Davies et al. (2010) derived the total mass of W33A as
15${}^{+5}_{-3}$$\rm{M}_{\odot}$, considering a distance of 3.8 kpc. Scaling
their mass to the revised distance adopted in our work, 2.4 kpc, the total
mass of W33A is 9.5${}^{+3.2}_{-1.9}$ $\rm{M}_{\odot}$, which is in good
agreement with the mass we obtained with the pPXF analysis. The pPXF method
allows the determination of the kinematics of very faint spectral features,
such as the CO absorption features, and also provides a robust estimate of the
errors, an important extension of the D10 result.

The rotation curve based on the CO absorption lines reported by D10 was
derived by smoothing the original CO absorption velocity centroids using a
5-pixel boxcar filtering. This procedure suggests this component is a
relatively extended structure with an angular size of about 1″ around the
source (see their Fig. 14), corresponding to a radius of $\sim$2,000 AU from
the protostar. Our results show a slightly more ordered view of the absorbing
gas and indicate that the velocities rapidly decrease for radii larger than
0$\aas@@fstack{\prime\prime}$1 ($\sim$500 AU) from the central position. The
velocities deviate from Keplerian outside the inner
0$\aas@@fstack{\prime\prime}$2 region of the PV diagram shown in Fig. 23
(corresponding to a projected length of $\sim$500 AU at the distance of the
source). These results suggest the circumstellar disc associated with W33A is
a compact structure, similar to the one observed towards the high-mass YSO
IRAS 13481-6124 by Kraus et al. (2010). Based on interferometric infrared
observations, those authors reported a $\sim$130 AU flared disc around a 18
$\rm{M}_{\odot}$ protostellar object lying within the IRAS source.

Recently, Maud et al. (2017) presented the analysis of interferometric ALMA
observations of the W33A protostar at Band 6 and 7 (230 and 345 GHz,
respectively) with angular resolutions of
$\sim$0$\aas@@fstack{\prime\prime}$2, similar to that of the NIFS data
(0$\aas@@fstack{\prime\prime}$15). Those authors have found that the kinematic
signatures associated with W33A are inconsistent with a Keplerian disc at the
angular resolution of their data. However, they did not exclude the
possibility of an underlying compact disc with Keplerian rotation at smaller
scales, such as confirmed by the NIFS K-band IFU observations. ALMA
observations obtained at larger baseline lengths (e.g. 16 km) are, therefore,
necessary to disentangle the structure and kinematics of the innermost region
of W33A.

###  4.3 A compact wind traced by the molecular hydrogen gas

The near-IR maps from Fig. 7 shows no evidence for a large-scale $\rm{H}_{2}$
jet associated with W33A, as one would expect for an active protostar (e.g.
Varricatt et al. 2010; Navarete et al. 2015). At longer wavelengths, recent
ALMA observations of W33A (Maud et al., 2017) have shown that high-velocity
emission of CO and SiO lines are probing a relatively compact jet-like
structure with angular sizes of $\sim$0$\aas@@fstack{\prime\prime}$5\. Despite
of that, no evidence for a collimated jet-like structure was observed in the
NIFS FOV, but a wide-angle $\rm{H}_{2}$ polar emission driven by the
protostar. The nature of the $\rm{H}_{2}$ emission can be traced by the other
structures associated with W33A. The orientation of the sub-pixel structure
probed by the spectro-astrometry of the $\rm{Br\gamma}$ feature (PA =
154$\pm$26∘, see Fig. 22) is compatible with the disc rotation axis
(150$\pm$10∘, see Fig. 23), suggesting that the ionised Hydrogen is tracing
the base of the jet/wind in high-mass YSOs as also seen in IRAS 13481-6124
through near-IR interferometry (Caratti o Garatti et al., 2016). In addition,
the orientation of the $\rm{H}_{2}$ emission (PA $\sim$ 145∘, see Fig. 15) is
also compatible with the ionised gas and the disc axis, favouring the
interpretation of a wide wind traced by the $\rm{H}_{2}$ emission rather than
the emission of scattered light in the cavity of the large-scale outflow. This
interpretation is also in agreement with the physical parameters of the
$\rm{H}_{2}$ emission derived through the $\rm{H}_{2}$ (1–0) S(1)/(2–1) S(1)
line ratio and the ro-vibrational analysis presented in Sect. 3.4.2,
confirming the thermal origin of the emission.

The physical parameters of the $\rm{H}_{2}$ emission were derived by
constructing the Boltzmann diagram of the different $\rm{H}_{2}$ transitions
of the entire $\rm{H}_{2}$ emitting region (see Fig. 20), and also by
evaluating the parameters for each spaxel in the data cube (see Fig. 21).
Based on the pixel-to-pixel analysis, the $\rm{H}_{2}$ emission exhibits a
mean $A_{\rm{K}}$ reddening of 4.3 $\pm$ 0.2 mag, an excitation temperature of
about 2.2$\times$103 K, and a total $\rm{H}_{2}$ column density of 1018.7±0.2
cm-2. The total mass of the $\rm{H}_{2}$ emission was estimated to be
M${}_{\rm H_{2}}$ = (4.6$\pm$1.0) $\cdot$ 10-5 $\rm{M}_{\odot}$. In general,
there is a good agreement between the physical parameters evaluated from both
methods (i.e. the integrated spectrum and the pixel-to-pixel analysis). The
spatially resolved analysis also shows the clumpy character of the
$\rm{H}_{2}$ emission. The $A_{\rm{K}}$ map shown in the top left panel of
Fig. 21 indicates the extinction towards the $\rm{H}_{2}$ jet is strong and
variable, with values ranging from $\sim$1.5 up to 5.0 magnitudes in the
K-band. The $A_{\rm{K}}$ values are larger towards the positions where the
$\rm{H}_{2}$ emission is stronger indicating dense clumps of gas and dust
perhaps along the cavity walls of the large-scale outflow (recall the emission
is consistent with shocked gas). But the conclusion is that the $\rm{H}_{2}$
is probing a wind.

We further compared the physical parameters of the $\rm{H}_{2}$ emission of
W33A using the typical values observed toward large-scale $\rm{H}_{2}$ jets of
high-mass protostars from Caratti o Garatti et al. (2015). The total column
density of the $\rm{H}_{2}$ wind driven by W33A is consistent with the range
of column density values (1017 to 1020 cm-2) reported by those authors. In
addition, they also derived the mean visual extinction ($A_{V}$) of the
$\rm{H}_{2}$ jets ranging from 1 to 50 mag, with a median of $A_{V}$ = 15 mag.
The upper limit of the visual extinction values reported by them is compatible
with the median $A_{V}$ value for the $\rm{H}_{2}$ wind associated with W33A,
$A_{V}$ = 63$\pm$2 mag.

Finally, we cannot infer the $A_{\rm{K}}$ towards the W33A protostar based on
the $\rm{H}_{2}$ analysis due to the absence of molecular gas emission closer
to the protostar. Alternatively, the $A_{\rm{K}}$ value of the point-like
source can be evaluated based on its K-band excess ($E_{H-K}$). Assuming the
extinction law from Damineli et al. (2016), the $A_{\rm{K}}$ is obtained as
$A_{K}=1.3\cdot E_{H-K}$. For W33A, the observed $H-Ks$ colour index is about
4.0 mag. Considering the intrinsic $H-Ks$ colour index for a B2 V star
($\sim$10 $\rm{M}_{\odot}$), $(H-Ks)_{0}$ = $-$0.05 mag (Wegner, 2003), we
obtain an extinction in the $K_{s}$-band of $A_{Ks}$ $\approx$ 5.3 mag for the
point-like source (or $A_{V}$ $\approx$ 80 mag). This value is in agreement
with the upper-limit of the $A_{\rm{K}}$ values associated with the
$\rm{H}_{2}$emission, shown in Fig. 21.

###  4.4 Different structures associated with ionised gas in W33A

The $\rm{Br\gamma}$ emission arises from gas at typical electron temperatures
of 8,000-10,000 K and it is observed in distinct structures of W33A,
suggesting that the hydrogen recombination plays a significant role in cooling
the gas. The Principal Component 1 (Fig. 10) indicates the bulk of the
$\rm{Br\gamma}$ emission arises at the central source of W33A. The source of
ionised gas could be an unresolved ultra-compact $\rm{H\,\textsc{ii}}$ region
around the protostar, in agreement with the absence of 5 GHz radio continuum
emission towards the protostar (Purcell et al., 2013).

In Principal Component 2 (Fig. 11), the $\rm{Br\gamma}$ emission mostly arises
from the the extended emission observed in the FOV, either originated by
Source B or W33A itself.

In PC3 (Fig. 12), the $\rm{Br\gamma}$ feature exhibits a similar kinematic
pattern as the CO absorption features. Our interpretation is that the CO
features are probing a compact and rotating circumstellar disc around W33A
(see Fig. 23), and the surface of the disc is partially ionised, producing a
disc wind, similar to that observed in Herbig Ae/Be stars (Tambovtseva et al.,
2016) but perhaps reaching larger radii owing to the more massive central
object. Although the analysis of the $\rm{H}_{2}$ emission does not indicate
the presence of any photo-ionised region in the extended wind, we cannot rule
out the presence of photo-evaporating wind in $\sim$100 AU scales closer to
the surface of the underlying cold disc, where the $\rm{H}_{2}$ emission is
weak probably because the molecules are being dissociated. As the wind moves
away from the central region, the $\rm{H}_{2}$ molecules are excited but not
dissociated, producing the observed extended structure as observed in e.g.
Fig. 21. This interpretation is complementary to the spectro-astrometric
result shown in Sect. 3.4.3, which shows the ionised gas tracing the cavity of
the outflow at milli-arcsecond scales, while PC3 shows evidences for the
presence of ionised gas at $\sim$100 AU scales, in the outer regions of the
circumstellar disc. The NIFS resolution is not sufficient to disentangle the
likely ionised surface from the underlying cold disc probed by the CO
features. Thus, as we pointed out above in Sect. 4.2, follow-up observations
at higher spatial resolutions (e.g. using the very long baseline of the ALMA
array) focusing on the investigation of the circumstellar environment may help
provide a better and more accurate view of the structure on small scales.

The PC4 in Fig. 14 shows that strong $\rm{Br\gamma}$ emission is observed very
close to the W33A protostar, while the spectro-astrometry of the
$\rm{Br\gamma}$ feature (Fig. 22) shows the innermost region of the
$\rm{H}_{2}$ jet cavity at milli-arcsecond scales ($\sim$1.3 AU), with a
Position Angle of 115${}^{\circ}\pm$26∘. Our results are in agreement with
those reported by D10, suggesting that the spectro-astrometry of the
$\rm{Br\gamma}$ feature is likely probing the base of the bipolar wind/jet
traced by the $\rm{H}_{2}$ emission. We also observe that the structure
identified in Fig. 22 is less affected by noise than that presented in D10,
indicating that the post-processing of the NIFS data cubes and noise
suppression plays a significant role on the analysis of sub-pixel scale
structures within the data.

##  5 Conclusions

We presented a reanalysis of the $K$-band NIFS/Gemini North observations of
the protostar W33A, first done by D10, based on Principal Component Analysis
Tomography. The PCS technique is well adapted to the analysis of complex data
sets like the one analysed here.

  1. 1.

The PCA tomography was able to recover the spectral and spatial information of
structures associated with the three key-ingredients of the disc-mediated
accretion scenario of the high-mass star formation process: the cavity of the
large-scale outflow probed by the ionised hydrogen, a rotating disc-like
structure probed by the CO absorption features, and a wide-angle wind traced
by the $\rm{H}_{2}$ emission. All these structures were present in the first
five Principal Components of the PCA Tomography. In addition, the technique
also reveals structures with different kinematics (e.g. the rotation of the
disc and the extended wind structure) in a much more efficient way than
investigating each spectral feature using traditional methods.

  2. 2.

The ionised gas ($\rm{Br\gamma}$) contributes to every structure associated
with the W33A protostar: although the bulk of the ionised gas is observed at
the central source as expected for a high-luminosity protostar, it also tracks
the inner region of the jet and its cavity, as well as the surface of the
disc.

The comparison of the structures identified in the NIFS FOV with large-scale
near-infrared maps around W33A led us to the following conclusions:

  1. 3.

The propagation of the $\rm{H}_{2}$ jet identified in the NIFS FOV is
consistent with the orientation of the $K$-band reflection nebula identified
in the $JHK$ maps. However, no evidence for a large-scale $\rm{H}_{2}$ jet
driven by W33A was found.

  2. 4.

A point-like object offset by $\sim$3″to the SE of W33A was identified in the
NIFS FOV. This object (Source B) has a relatively blue spectrum, which is
consistent with the bright $J$-Band source identified in the $JHK$ maps.
Source B is also associated with a $H$-band extended emission oriented to the
NW direction that overlaps the emission from W33A, contributing to the
observed $K$-band continuum extended emission observed in the NIFS FOV.

  3. 5.

A third point-like source (Source C), offset by 6″to the S direction of W33A,
was identified in the large-scale near-infrared maps. Source C is likely the
driving source of a set of $\rm{H}_{2}$ knots oriented in the E-W direction.
The identification of sources B and C revealed that W33A is located in a
region with multiple (proto)stellar objects.

We reconstructed the data cube using only the first five Principal Components
which represents 99.996% of the variance and from this data cube we derive the
following conclusions:

  1. 6.

First, we derived the physical parameters of the $\rm{H}_{2}$ jet and found
that it is associated with an excitation temperature of (2.2 $\pm$ 0.1)
$\cdot$ 103 K, and a mean $\rm{H}_{2}$ column density of
$\log$($N_{\rm{tot}}$/cm-2) = 18.7$\pm$0.2. The mean $K$-band extinction is
about 4.1$\pm$0.2 mag (corresponding to an $A_{\rm V}$ of $\sim$61 mag), and
the mass of the $\rm{H}_{2}$ jet is about $\log$(MH2/$\rm{M}_{\odot}$) =
$-$4.3$\pm$0.3.

  2. 7.

We reanalysed the spectro-astrometry of the $\rm{Br\gamma}$ emission,
confirming that the ionised gas is tracing the cavity of the jets at sub-pixel
scales. The positional offset between the blue- and red-shifted components is
$\sim$0$\aas@@fstack{\prime\prime}$55, about three times smaller than
previously reported ($\sim$1$\aas@@fstack{\prime\prime}$5), corresponding to a
linear size of $\sim$1.3 AU at 2.4 kpc.

  3. 8.

We measured the kinematics of the CO absorption lines to infer the dynamical
mass of the W33A protostar. The analysis of the position-velocity diagram of
the CO absorption features indicates a truncated rotating structure, with
Keplerian-like rotation in the inner 0$\aas@@fstack{\prime\prime}$2 radii. We
used a Keplerian model to evaluate the mass of the “protostar + disc” system
as 10.0${}^{+4.1}_{-2.2}$ $\rm{M}_{\odot}$ assuming a distance of 2.4 kpc. The
mass is about $\sim$33% smaller than previously reported by Davies et al.
(2010), placing W33A closer to the lower mass-limit for a high-mass protostar.

  4. 9.

The above-mentioned results were obtained thanks to the application of the PCA
Tomography towards a bona fide high-mass protostar, showing the technique is
encouraging to be applied to other high-mass young stellar object candidates.

## Acknowledgements

FN thanks to Fundação de Amparo à Pesquisa do Estado de São Paulo - FAPESP for
support through proc. 2017/18191-8. AD acknowledges FAPESP for support through
proc. 2019/02029-2. The authors would like to thank the valuable comments and
suggestions made by the anonymous referee that helped to improve the
manuscript.

## Data availability

The data sets were derived from observations available in the public domain:
http://www1.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/gsa/ (under the Gemini program
GN-2008A-Q-44).

## References

  * Alves et al. (2019) Alves F. O., Caselli P., Girart J. M., Segura-Cox D., Franco G. A. P., Schmiedeke A., Zhao B., 2019, Science, 366, 90
  * Black & Dalgarno (1976) Black J. H., Dalgarno A., 1976, ApJ, 203, 132
  * Blum & McGregor (2008) Blum R. D., McGregor P. J., 2008, AJ, 135, 1708
  * Blum et al. (2004) Blum R. D., Barbosa C. L., Damineli A., Conti P. S., Ridgway S., 2004, ApJ, 617, 1167
  * Cappellari & Emsellem (2004) Cappellari M., Emsellem E., 2004, PASP, 116, 138
  * 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. (2015) Caratti o Garatti A., Stecklum B., Linz H., Garcia Lopez R., Sanna A., 2015, A&A, 573, A82
  * Caratti o Garatti et al. (2016) Caratti o Garatti A., et al., 2016, A&A, 589, L4
  * Cohen et al. (2003) Cohen M., Wheaton W. A., Megeath S. T., 2003, AJ, 126, 1090
  * Dabrowski (1984) Dabrowski I., 1984, Canadian Journal of Physics, 62, 1639
  * Dahmer-Hahn et al. (2019) Dahmer-Hahn L. G., et al., 2019, MNRAS, 489, 5653
  * Damineli et al. (2016) Damineli A., Almeida L. A., Blum R. D., Damineli D. S. C., Navarete F., Rubinho M. S., Teodoro M., 2016, MNRAS, 463, 2653
  * Davies et al. (2010) Davies B., Lumsden S. L., Hoare M. G., Oudmaijer R. D., de Wit W.-J., 2010, MNRAS, 402, 1504
  * Fedriani et al. (2020) Fedriani R., et al., 2020, A&A, 633, A128
  * Froebrich et al. (2011) Froebrich D., et al., 2011, MNRAS, 413, 480
  * Ilee et al. (2013) Ilee J. D., et al., 2013, MNRAS, 429, 2960
  * Immer et al. (2013) Immer K., Reid M. J., Menten K. M., Brunthaler A., Dame T. M., 2013, A&A, 553, A117
  * Kewley et al. (2019) Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, ARA&A, 57, 511
  * Kraus et al. (2010) Kraus S., et al., 2010, Nature, 466, 339
  * Lee et al. (2013) Lee H.-T., et al., 2013, ApJS, 208, 23
  * Lumsden et al. (2013) Lumsden S. L., Hoare M. G., Urquhart J. S., Oudmaijer R. D., Davies B., Mottram J. C., Cooper H. D. B., Moore T. J. T., 2013, ApJS, 208, 11
  * Maud et al. (2017) Maud L. T., Hoare M. G., Galván-Madrid R., Zhang Q., de Wit W. J., Keto E., Johnston K. G., Pineda J. E., 2017, MNRAS, 467, L120
  * McGregor et al. (2003) McGregor P. J., et al., 2003, in Iye M., Moorwood A. F. M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4841, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. pp 1581–1591, doi:10.1117/12.459448
  * Menezes et al. (2013) Menezes R. B., Steiner J. E., Ricci T. V., 2013, ApJ, 765, L40
  * Menezes et al. (2014) Menezes R. B., Steiner J. E., Ricci T. V., 2014, MNRAS, 438, 2597
  * Murakawa et al. (2013) Murakawa K., Lumsden S. L., Oudmaijer R. D., Davies B., Wheelwright H. E., Hoare M. G., Ilee J. D., 2013, MNRAS, 436, 511
  * Navarete et al. (2015) Navarete F., Damineli A., Barbosa C. L., Blum R. D., 2015, MNRAS, 450, 4364
  * Purcell et al. (2013) Purcell C. R., et al., 2013, ApJS, 205, 1
  * Ricci et al. (2011) Ricci T. V., Steiner J. E., Menezes R. B., 2011, ApJ, 734, L10
  * Ricci et al. (2014) Ricci T. V., Steiner J. E., Menezes R. B., 2014, MNRAS, 440, 2419
  * Ricci et al. (2015) Ricci T. V., Steiner J. E., Giansante L., 2015, A&A, 576, A58
  * Shu et al. (1987) Shu F. H., Adams F. C., Lizano S., 1987, ARA&A, 25, 23
  * Skrutskie et al. (2006) Skrutskie M. F., Cutri R. M., Stiening R., Weinberg M. D., et al., 2006, AJ, 131, 1163
  * Steiner et al. (2009) Steiner J. E., Menezes R. B., Ricci T. V., Oliveira A. S., 2009, MNRAS, 395, 64
  * Tambovtseva et al. (2016) Tambovtseva L. V., Grinin V. P., Weigelt G., 2016, A&A, 590, A97
  * Tokunaga (2000) Tokunaga A. T., 2000, Infrared Astronomy. Springer, p. 143
  * Tokunaga et al. (2002) Tokunaga A. T., Simons D. A., Vacca W. D., 2002, PASP, 114, 180
  * Turner et al. (1977) Turner J., Kirby-Docken K., Dalgarno A., 1977, ApJS, 35, 281
  * Varricatt et al. (2010) Varricatt W. P., Davis C. J., Ramsay S., Todd S. P., 2010, MNRAS, 404, 661
  * Wegner (2003) Wegner W., 2003, Astronomische Nachrichten, 324, 219
  * de Wit et al. (2010) de Wit W. J., Hoare M. G., Oudmaijer R. D., Lumsden S. L., 2010, A&A, 515, A45+

##  Appendix A Aperture photometry of the UKIDSS JHK images

The $JHK$ magnitudes of the W33A, Source B and Source C are listed in Table 2.
The full list containing the 287 sources is available online as supplementary
material and in the CDS.

Table 4: $JHK$ Photometry of the sources (first 10 rows). The full table is
available online as supplementary material and in the CDS.

ID | RA | Dec | UKIDSS | J | $\sigma_{J}$ | H | $\sigma_{H}$ | K | $\sigma_{K}$ | $J-H$ | $\sigma_{J-H}$ | $J-K$ | $\sigma_{J-K}$ | $H-K$ | $\sigma_{H-K}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| (J2000) | (J2000) | Designation | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag)
1a | 273.66458 | -17.86669 | J181439.51-175200.9 | 20.77∗ | – | 15.05∗ | – | 9.51 | 0.19 | 5.72 | 0.64 | 11.26 | 0.66 | 5.54 | 0.21
2b | 273.66486 | -17.86736 | J181439.58-175202.6 | 15.45+ | 0.08 | 13.76 | 0.09 | 10.30 | 0.20 | 1.70 | 0.13 | 5.16 | 0.22 | 3.46 | 0.22
3c | 273.66495 | -17.86827 | – | 20.77∗ | – | 16.31 | 0.11 | 12.69 | 0.23 | 4.46 | 0.64 | 8.08 | 0.68 | 3.62 | 0.26
4 | 273.66562 | -17.86503 | J181439.75-175154.1 | 11.37 | 0.07 | 10.99 | 0.08 | 10.29 | 0.20 | 0.38 | 0.11 | 1.08 | 0.21 | 0.70 | 0.21
5 | 273.66211 | -17.85791 | J181438.90-175128.3 | 13.24 | 0.08 | 12.63 | 0.09 | 12.35 | 0.21 | 0.61 | 0.12 | 0.89 | 0.22 | 0.28 | 0.23
6 | 273.67337 | -17.87352 | J181441.61-175224.6 | 13.35 | 0.08 | 13.07 | 0.09 | 12.93 | 0.21 | 0.28 | 0.12 | 0.42 | 0.23 | 0.14 | 0.23
7 | 273.67474 | -17.86697 | J181441.94-175201.1 | 13.57 | 0.08 | 13.18 | 0.09 | 12.99 | 0.21 | 0.39 | 0.12 | 0.59 | 0.23 | 0.20 | 0.23
8 | 273.66714 | -17.88072 | J181440.11-175250.5 | 14.11 | 0.08 | 13.14 | 0.09 | 12.63 | 0.21 | 0.96 | 0.12 | 1.48 | 0.22 | 0.51 | 0.23
9 | 273.66238 | -17.87348 | J181438.97-175224.5 | 14.10 | 0.08 | 13.30 | 0.09 | 12.88 | 0.21 | 0.81 | 0.12 | 1.22 | 0.23 | 0.42 | 0.23
10 | 273.66678 | -17.86504 | J181440.03-175154.1 | 14.29 | 0.08 | 13.40 | 0.09 | 12.95 | 0.21 | 0.89 | 0.12 | 1.35 | 0.23 | 0.46 | 0.23

Notes: the columns are as follows: (1) ID of the source (a: W33A; b: Source B;
c: Source C); (2) Right ascension (J2000, in degrees); (3) Declination (J2000,
in degrees); (4) Name of the source or its designation in the UKIDSS
catalogue; (5) $J$-band magnitude; (6) $J$-band magnitude error; (7) $H$-band
magnitude; (8) $H$-band magnitude error; (9) $K$-band magnitude; (10) $K$-band
magnitude error; (11) $J-H$ colour index (in mag); (12) $J-H$ colour index
error (in mag); (13) $J-K$ colour index (in mag); (14) $J-K$ colour index
error (in mag); (15) $H-K$ colour index (in mag); (16) $H-K$ colour index
error (in mag). An $\ast$ symbol indicates upper limit due to the non-
detection of a point-like source at the given position in the corresponding
filter. $+$ indicates the $J$-band flux of a foreground source lying very
close to the targeted object.