\section{Availability:} All code is available from the first author upon request.
\end{abstract}

\section{Introduction}

Copy number alterations, including both germline variants (CNVs) and somatic copy number aberrations (CNAs) are associated with disease \cite{pollack99genomewide}. Somatic aberrations (CNAs) are particularly important for tumourigenesis, contributing to genomic instability and gene deregulation. For example, oncogene activation by gene amplification or tumour suppressor loss as a result of gene deletion (Fig.\ref{dna-mrna}) can cause transcriptional changes and contribute to the pathogenesis of cancer.
On the other hand, expression can be influenced by the proximal genes within a several Mb window (\emph{cis}-acting), but can also be affected by a remote CNA throughout the genome (\emph{trans}-acting), as depicted in Fig.\ref{dna-mrna}. One of the challenges in cancer genomics is to characterize such intermediate phenotypic changes caused by both {\it cis-} and {\it trans-} CNAs that ultimately lead to tumour progression.

\begin{figure}[tb]
\centering
\includegraphics[width=.5\textwidth]{Figure1}
\caption{ Expression abundance level can be affected by both {\it cis}- and {\it trans}- acting copy number alterations. integrating gene expression data and gene dosage information measured by microarray or sequencing technologies. }
\label{dna-mrna}
\end{figure}

\textbf{Detecting {\it cis}- and {\it trans}-regulatory effects}
An integrative analysis of CNA and expression has the potential to uncover regulatory relationships between CNA and gene expression (Fig.\ref{dna-mrna}). The availability of genome-wide copy number and gene expression data facilitates the discovery of such regulatory maps, but few studies have quantitatively assayed the effect of both {\it cis-} and {\it trans-} copy number on the expression abundance in the same primary tumour sample. Given the inherent noise in such data and the large-scale of genomic information, the detection of both {\it cis}- and {\it trans}-acting effects between DNA copy number and mRNA on a genome-wide basis is a difficult task.  Such effects are small and infrequent \cite{sladek06elucidating}, which further imposes computational challenges. Detection of  the driving alterations amidst numerous random passenger events is complicated by the frequency of co-occurring events such as co-amplifications and co-deletions. Thus despite considerable efforts, the location of many relevant CNAs that result in gene disregulation remain elusive.

\textbf{Previous research}
This area of research is of key importance since many somatic alterations could contribute to tumorigenesis and their identities are critical for the development of specific anti-cancer agents. Not surprisingly, the association between DNA copy number and mRNA expression data has gained considerable attention. One of the first papers to combine gene expression and DNA copy number data by  \cite{pollack99genomewide} concluded that 62\% of highly amplified genes demonstrated moderate to high expression levels. Others have similarly identified genes that exhibit gain/loss and simultaneous over/under-expression in cancer \cite{myllykangas08integrated}, but it has been noted that chromosomal amplifications do not necessarily result in the over-expression of the genes located in the same cytoband  \cite{platzer02silence}.

In summary, recent studies have focused on the local effects of copy number on expression, and few have considered {\it trans}-acting effects on a global scale. Correlation analysis remains in widespread use for these types of comparisons due to its efficiency \cite{bussey06integrating}, and has more recently been used to query {\it trans}-effects on a genome-wide level \cite{lee08integrative}.  Several more sophisticated tools have also been proposed.  For example,  \cite{menezes09integrated} made use of gene sets instead of individual genes in regression models in their efforts to integrate the two data types. However, simple regression has limited power in face of the immense amount of genomic data. Moreover, their focus is on {\it cis}- events and therefore cannot incorporate remote interactions.

\textbf{Our approach}
In order to find copy-number driven expression, we propose a new framework for reconstructing a genome-wide regulatory map on a genomic and transcriptomic level. The objective of the proposed framework is to set up an interaction network between CNA and expression changes (arrows in Fig.\ref{dna-mrna}). By interpreting global gene expression changes in the context of CNA, the network has the potential to uncover both {\it cis}- and {\it trans}-regulatory elements
Our approach combines several key ideas: First, we employ an $L_1$-regularized regression model which allows simultaneous feature selection and model fitting which generalizes well in high-dimensions. Second, we introduce a combination of global and local search strategies to optimally set the regularization parameter in our model. Third, we quantitatively measure each transcript's  copy-number dependence.

\section{Method}
We assume copy number data of $p$ DNA regions $X= (x_1, x_2, ..., x_p)$  and mRNA expression data of $q$ genes $Y= (y_1, y_2, ..., y_q)$, both available from $n$ tumour samples.

The framework for investigating CN trans- effect is built on a linear regression method with $L_1$-regularization (section~\ref{method-lasso}). The strength of regularization is chosen adaptively for each gene by a combined global and local search over the space of regularization parameters (Section~\ref{method-search}). The results are incorporated into a quantitative measure for copy-number dependence by comparing variances of the model residuals with and without a locus of interest (Section~\ref{method-granger}).

\subsection{$L_1$-constrained regression selects most influential genomic loci for each transcript}\label{method-lasso}
We use $L_1$ constrained regression (a.k.a the {\it Lasso} \cite{lasso}) to derive the causal effects of genome-wide copy number data, as the predictor variables, on gene expression profiles, as the response variables. $L-1$-regression minimizes the sum of squared errors between response and prediction, while keeping a bound on the sum of the absolute values of the regression coefficients. By penalizing the absolute sum of coefficients, $L_1$-constrained regression effectively shrinks non-significant coefficients to zero. In contrast to this, $L_2$-constrained regression (aka \emph{ridge regression}) keeps a bound on the sum of the squared values of the regression coefficients and generally results in small, but non-zero, regression coefficients. Hence,  $L_1$-regression is a prominent method in high-dimensional studies for its efficiency in performing sparse statistical inference on thousands of variables \cite{tibshirani97lasso}. Efficient computational methods exist to infer the entire path of variable inclusions over different choices of regularization strength \cite{lars}.

\textbf{In our scenario:} For gene $i$  with expression profile $y_i$, let $X \beta_i$ be the candidate prediction model with regression coefficients $\beta_i = (\beta_{i1}, \beta_{i2}, ..., \beta_{in})$. The Lasso model is then given by
\begin{equation}\label{EqLasso}
\hat{\beta}_i = \argmin_{\beta_i} \   \left\|\  y_i - X \beta_i \  \right \|^2_2 +\lambda_i  \left\| \beta_{i} \right\|_1 ,
\end{equation}
where $\| x \|_2 = \sum x_i^2$ indicates the $L_2$ norm measuring the distance between model prediction, and response and $\| x \|_1 = \sum |x_i|$ indicates the $L_1$ norm used to constrain the regression coefficients.

The $\lambda_i$ in Eq.(\ref{EqLasso})  is a regularization parameter controlling the sparsity and  strength of regularization.
By varying $\lambda_i$,  $L_1$-regression/Lasso adaptively includes predictors in the model (see Figure~\ref{lambda}). The subscript ``$i$" emphasizes that the optimization problem in Eq.(\ref{EqLasso}) has to be solved for every gene. In particular, for each gene a regularization parameter $\lambda_i$ has to be chosen. In the next section, we propose to utilize an efficient combination of a global and local search strategy.

\begin{figure}[t]
\centering
\includegraphics[width=.48\textwidth]{MYClassopath}
\caption{The  Lasso path along a range of $\lambda$ values for  {\it MYC}, a well known oncogene,  in the experimental dataset. Each line corresponds to a regression coefficient measuring the influence of copy number variation at a particular locus (names on the right) on {\it MYC} expression. Regularization strength decreases with $\lambda$ from the left to the right. This can be seen by more and more coefficients hitting the zero line and vanishing from the model.}
\label{lambda}
\end{figure}

\subsection{Adaptive regularization by combined global and local search}\label{method-search}

The regularization parameter $\lambda$ is crucial since it directly determines how many predictor vectors are to be included into the model.
We perform a two-step procedure for $\lambda$ selection that combines a global search over a discrete set of candidate values with a local search around the optimal candidate.
For the global search the likelihood function is evaluated by cross-validation, for the local search Brent's minimization without derivatives \cite{brent73} implemented in R package ``penalized'' \cite{geoman09l1} is used.
Brent's method combines steps of golden section (a fast method based on binary divisions of the search space) with parabolic interpolation to efficiently zoom in onto the maximum.
\begin{enumerate}
\item[Step 1] Global search: We perform 10-fold cross-validation across a wide range of candidate $\lambda$ values (grey lines in Fig.\ref{loglike}). Log-likelihood scores (black dots) are computed for each of these values to provide a point of focus for the next step.
\item[Step 2] Local search: At the point of the maximum log-likelihood of Step 1 (black line in Fig.\ref{loglike}), we perform a local search (grey box) with Brent's optimization to locate the optimal value of $\lambda$ (dashed line).
\end{enumerate}
\begin{figure}[!h]
\centering
\includegraphics[width=.45\textwidth]{Figure3}
\caption{$\lambda$ selection procedure base on the log-likelihood function. The grey lines mark the global search points, and the two black lines (one solid, one dashed) mark the optimal selection points of $\lambda$ by Step 1 and 2.}
\label{loglike}
\end{figure}

Having chosen an optimal $\lambda$ value and computing the Lasso solution of Eq.(\ref{EqLasso}) for each gene, we obtain a coefficient matrix $B = (\hat{\beta}_1, \hat{\beta}_2, \ldots, \hat{\beta}_n)$. Because of the $L_1$-constraint in the Lasso, this matrix will be sparsely populated and many coefficients will be zero. Matrix $B$ represents the regulatory relationships  between copy number changes and transcriptional changes that was indicated by arrows in Fig.\ref{dna-mrna}.

\textbf{A technical note} Lasso assumes the predictor variables to be linearly independent. However, DNA copy number data obtained from high resolution microarrays can exhibit high correlation between probes. Therefore, copy number data should be first tested for independence and, if necessary, the data should be merged with an appropriate method such as those implemented in R package ``CGHregion'' \cite{cghregions} to produce regions with independent signatures. We implemented our framework in R with the packages ``penalized"  and ``CGHregions".

\begin{figure}[!b]
 \sffamily
 \begin{minipage}{.48\textwidth}
\textbf{A.}  Precision curves

\centerline{\includegraphics[width=.95\textwidth]{precisionCurves}}
\end{minipage}
\hfill
 \begin{minipage}{.48\textwidth}
\textbf{B.} Recall curves

\centerline{ \includegraphics[width=.95\textwidth]{recallCurves}}
\end{minipage}

\caption{Performance of three methods in the experiments of network structure learning, data is simulated with varied number of samples $n$ and number of predictor variables $p$ and number of response variables $q, p=q$. }
    \label{structure}
\end{figure}

\subsection{A quantitative measure of copy-number dependence}\label{method-granger}

The proposed framework allows us to quantitatively  measure the dependency between gene expression and copy number, both in {\it cis} and in {\it trans}.  The underlying assumption is that the predictive power of a model can be measured by the proportion of variance explained, i.e., the ratio of the variance of modelling residuals with or without incorporating  one or more predictors.  More formally, given two models --one {\it with} and one {\it without} incorporating a set of predictors-- the negative natural logarithm of the variance ratio between the prediction residuals,
\begin{equation}\label{EqGranger}
\text{score}_i \ = \ - \ln \left(\frac{\sigma^2_{with}}{\sigma^2_{without}}\right),
\end{equation}
can be used as a score indicating the dependency of the response variable on the predictor(s).
Therefore, by checking if the variance of residuals is reduced after incorporating the copy-number predictors,  the dependency of expression responses on copy number can be quantitatively measured.

Furthermore, comparisons can be conducted with copy-number {\it cis}- and {\it trans}-, indicating whether a gene is significantly influenced by the {\it trans}- effect.
Formally, we do this by comparing a sequence of nested models that successively incorporate more and more copy number predictors: from no predictor at all, to using only copy number changes in  {\it cis}, and finally to using both  {\it cis}- and  {\it trans}- copy number changes.
\begin{eqnarray}
y_i &=&y_i \label{EqPrediction1}\\
y_i &= &x_i\beta_i  + \epsilon_i \label{EqPrediction2}\\
y_i &=&x_i\beta_i  + \sum_{j, j \neq i} x_j \beta_j + \epsilon_i' \label{EqPrediction3}
\end{eqnarray}
Eq.(\ref{EqGranger}) then allows us to estimate the effect of adding more predictors by measuring the difference of variance explained between any two equations in Eq.(\ref{EqPrediction1})-(\ref{EqPrediction3}).
For example, the dependency score of  expression of a gene $i$' on copy number can be calculated based on Eq.(\ref{EqPrediction1}) and Eq.(\ref{EqPrediction3})
\begin{equation}\label{EqCopyNumberDependent}
\text{score}_i^{both} \ = \ -\ln \left(\frac{\sigma^2_{\epsilon_i'}}{\sigma^2_{y_i}}\right).
\end{equation}
And the  dependency score of expression $y_i$ on the copy number {\it trans}-effect conditioned on the {\it cis}- effect based on Eq.(\ref{EqPrediction2}) and Eq.(\ref{EqPrediction3}) is
\begin{equation}\label{EqTransDependent}
\text{score}_i^{trans} \ = \ -\ln \left(\frac{\sigma^2_{\epsilon_i'}}{\sigma^2_{\epsilon_i}}\right).
\end{equation}
If the score is significant, there is likely to be a copy number {\it trans}-effect on gene $i$. Note that by incorporating the  {\it cis}- predictor $x_i$ in both equation, we remove the {\it cis}-effect before exploring the {\it trans}-effect $ \sum_{j, j \neq i} x_j$, which is equivalent to a conditional independence test.

\section{Experiments on simulated networks} \label{simulation}

Simulated data was generated to test the performance of the $L_1$ method on network structure inference. Let $X\sim N(0, 1)$ be a matrix with $p$ Gaussian variables each of $n$ samples, and $T$ be a sparse matrix of $p$ by $q$ generated from randomly sampling 50 non-zero coefficients from a distribution $N(0, 1)$. From $X$ and $T$, a response vector $Y = XT$  with $q$ variables is generated. Both $X$ and $Y$ are added with centered Gaussian noise $N(0, \sigma^2_\epsilon)$. The objective is to infer $T$ from $X$ and $Y$.

Two popular network reconstruction methods in the literature that are applicable to the small $n$, large $p$ problem are compared with the $L_1$ method. They are the Gaussian Graphical Models  \cite{schafer05empirical} and Partial Least Square \cite{vasyl08reconstruction}. However, to remove any ambiguity that could be caused by the different numbers of inferred relations, the number of relations, i.e. total number of non-zero coefficients in $T$ are known to all methods.

The performance of the methods in comparison is judged by whether the non-zero coefficients in $T$ can be correctly detected. Data is simulated with $p$, the number of variables in $X$,  varying from 100 to 2000 and sample size $n$ varying from 50 to 500. The number $q$ of variables in $Y$ is equal to $q$. For each setting of $p$ and $n$, datasets are generated with different noise levels $\sigma^2_\epsilon \in \{0.05, 0.1, 0.2\}$. For the result, precision and recall statistics for the inference results of all variables in $Y$ are computed and averaged across noise levels.

The precision and recall curves, grouped by values of $p$, are plotted in Fig.\ref{structure}. The advantage of the $L_1$ method can be clearly seen:   in all scenarios precision and recall are higher than for the competitors. The difference is most pronounced for large values of $p$, which is especially important for real-world applications. In summary the simulations show that the proposed method outperforms other state-of-the-art methods in its accuracy to reconstruct regulatory networks.

\section{A genome-wide map of $cis$- and $trans$-regulation in breast cancer} \label{results}

To demonstrate the effectiveness of the proposed method in a real-world application we used a breast cancer dataset with 89 samples analyzed using both mRNA expression and DNA copy number microarrays \cite{chin06genomic}.

\subsection{A quantitative assessment of prediction accuracy} \label{experiment}

While the simulations showed the accuracy of our framework in network reconstruction compared to other network models,  on real data we can not repeat this analysis since the true network is unknown. However, a complementary question that we can answer is how well our model predicts gene expression from copy-number alterations compared to other state-of-the-art regression models. We compare $L_1$-constrained regression with three other well known prediction models: Random Forests \cite{randomforests}, SVM \cite{svm}, and Recursive Partitioning \cite{recursivepartitioning}. Each method uses copy number data to predict expression data and accuracy is measured in terms of the mean squared error of the prediction.

Models were trained with a subset of samples in the \cite{chin06genomic} data and then the fitted models were used to predict the rest of samples.  In detail: 100 expression profiles are randomly selected from the Chin data. For each of the expression profile, models were fitted with 7/8 (78) of the samples and then the mean squared errors of predicting the test set 1/8 (11) of the samples.
The Random Forests method was trained with 500 trees so that the computational time is on the same scale as that of other methods. SVM was trained with a linear kernel.

\begin{figure}
\includegraphics[width=.45\textwidth]{predictiveBoxplots}
\caption{Mean squared errors of predictions by four methods on Chin's dataset, the dashed line shows the median of mean squared errors of Lasso. The average of the mean squared errors for four methods are 0.526, 0.528, 0.635, 0.707. }
\label{Fig.predictive}
\end{figure}

Results in Fig.\ref{Fig.predictive} show that the $L_1$-constrained regression is not out-performed by any competitor. Additionally, it provides a sparser solution than the other models: On average $L_1$-regression used less than 7 features, while Recursive Partitioning used more than 32 and SVM and Random Forests use all of them.  Thus, the $L_1$-model combines sparsity with high prediction accuracy. This is important since the selected features directly point to the relevant copy number alterations leading to expression changes. This experiment proves the proposed framework's potential with applications on copy number and expression data.

\subsection{The global impact of copy-number change on gene expression}

The proposed framework was used to study the impact of copy-number aberrations on every gene expression profile in the Chin dataset.
Using Eq.(\ref{EqCopyNumberDependent}) we can rank genes based on the copy-number dependent score, as shown in Table~\ref{TableDependentGenes}.
Several known oncogenes are present such as {\it ERBB2} and its neighbors, {\it STARD1} and {\it MGC9753}, as well as {\it ESR1}. For example, {\it ERBB2} is amplified in approximately 20\% of the breast cancer cases for which concomitant over-expression is also typically observed. Many other genes with less well characterized functions were also identified and merit additional investigation, e.g. TFF1 which is a stable secretory protein expressed in the gastric mucosa.

\begin{table}[!htp]
\caption{Top 20 genes ranked by the copy-number dependent score}
\centering
\begin{tabular}{|c|c|c||c|c|c|} \hline
Gene &Chr  &  Score & Gene &Chr &  Score \\ \hline
TFF1 & 21 &   3.31 &MGC9753 & 17 &   2.38 \\
NAT1 &  8 &  2.73  &PROML1 &  4 &  2.33\\
TFF3 & 21 &   2.63 &CEACAM6 & 19 &  2.31\\
AGR2 &  7 &   2.61 &SCGB2A2 & 11 &  2.31\\
GABRP &  5 &  2.53 &CPB1 &  3 & 2.27 \\
MGC9753 & 17 &   2.5 &MLPH &  2 & 2.26\\
ESR1 &  6 &  2.49 &GRB7 & 17 &   2.26 \\
S100A7 &  1 &  2.44 &SFRP1 &  8 &  2.26 \\
PIP &  7 &  2.39 &ERBB2 & 17 &   2.24\\
HNF3A & 14 &  2.38 &CEACAM5 & 19 & 2.21\\
\hline
\end{tabular}
\label{TableDependentGenes}\end{table}

Expressions driven by the same CNAs region are also enriched with cancer-specific functions.  We test whether the downstream targets of these DNA regions share similar functions or motifs, since genes with similar functions are often co-regulated. Gene Ontology (GO) \cite{GO} enrichment analysis was performed on the downstream targets for CNAs regions driving many expressions.
Numerous GO terms with low $p$-values, indicative of a high degree of functional similarity, were identified. For example, the targets of 22q13.33 are significantly enriched for GO terms involving the immune response ($p$-value 6E-11), lymphocyte activation (9E-11), and leukocyte activation (5E-10), with the former two pertaining to the adaptive immune response.

Genetic variation at 8q24 has been associated with both breast and prostate cancer \cite{meyer09genetic} and this region is frequently amplified in breast cancer. One of the top ranking hot spots, 8q24.13, has downstream targets enriched for GO terms such as the positive regulation of ubiquitin-protein ligase activity during mitotic cell cycle, and the proteasomal protein catabolic process. Promoter analysis revealed a significantly enriched unknown motif CGGCGACATA (hypergeometric test $p$-value 4.7E-19, FDR on 50 random samples $p$-value 0.0044) in the upstream (-2000 bp) region.  Another hot spot, 9q22.2, exhibits {\it trans}-effects on {\it TOP2A} (17q21) and the MCM family, and has downstream targets enriched for DNA replication and mitotic cell cycle functions. Once again, these findings are in-line with existing knowledge on cancer biology amassed using various approaches.

\subsection{Individual example of copy-number driven expression}

\begin{figure}[t]
\centering
 \sffamily
\textbf{A.}  Correlation analysis of {\it MYC} expression with genome-wide copy number data

\centerline{\includegraphics[width=.5\textwidth]{MYCqval}}

\smallskip
\textbf{B.} Lasso analysis of  {\it MYC} expression with genome-wide copy number data

\centerline{ \includegraphics[width=.5\textwidth]{MYClasso}}
\caption{Example1:  {\it MYC}. The upper plot shows  the $q$-value (a multiple testing corrected $p$-value) of correlation coefficients between copy number and gene expression (red for positive correlations, blue for negative correlations). The lower plot shows the much sparser results from our framework, which illuminates key alterations by removing noise in the data.}
    \label{myc}
\end{figure}

We also examined the impact of copy number on individual target expressions. For the purpose of comparison, correlation of gene expression and genome-wide copy number was performed, and the results after FDR (False Discovery Rate) correction under the Benjamini-Hochberg procedure \cite{benjamini95controlling} is plotted in Fig.\ref{myc}.

Fig.\ref{myc}(a) shows that there is very high correlation between the expression of the oncogene {\it MYC} and copy number 3q26.2, 5q12-5q34. In fact, these regions exhibit nearly as high a correlation as does 8q24 (which denotes the location of  {\it MYC}). Here, the red histogram indicates positive correlation, while the blue histogram indicates negative correlation. Fig.\ref{myc}(b) illustrates that the dependence of  {\it MYC} expression on copy number in 5q21-5q34, when conditioned on its own copy number (8q24), is much smaller than that of copy number in 8q24, and the dependence on 3q26 vanishes. This is due to the fact that co-occurring events, such as co-amplification, have been removed by the conditional dependence assumption in Lasso. Moreover, the effects of copy number in 9q22.2 and 11p11.12 amongst others, shown in Fig.\ref{myc}(a) are completely removed when conditioning on the copy number at 8q24. In contrast, an apparently strong negative regulatory effect from 17q21.31, which harbours {\it WNK4} and the breast cancer susceptibility gene {\it BRCA1} persists, suggesting an interaction between these key genes.

\section{Discussion}

This paper presents an efficient framework for integrating genome-wide copy number and expression data. Both {\it cis}- and {\it trans}- regulatory effects from the genomic level to the transcriptomic level can be quantitatively measured using conditional probabilities. In particular, $L_1$-constrained regression is able to select relevant predictor variables to be included into the statistical model.  In summary, our proposed framework not only allows for a genome-wide description of the effect of gene dosage on gene expression, but also can partition copy-number dependent and independent transcriptional changes. While the former are useful to identify copy number and expression events co-ordinatedly deregulated during cancer progression, the latter are interesting as a starting point for validating putative oncogenes. Thus, our framework not only yields a global view of the interactions and suggests downstream targets of CNAs ripe for validation, it can also be used as a starting point for further focused analyses of the genomic basis of cancer. Our findings in the Chin dataset are currently being validated in another much larger breast cancer dataset, which has not yet been published.

\bibliographystyle{IEEEtran}


\section{Introduction}

It is thought that the Hypercompact or Ultracompact HII (HCHII or
UCHII) regions trace early stages of massive star formation after the
hot molecular core phase. However, the physical processes involved in
this transition are still poorly understood \citep[see for
  reviews:][]{stan2005,liz2008}. The appearance of one HCHII begins
when the Lyman continuum output from a massive young star becomes
sufficient high to ionize its surroundings.  These objects show a
rising continuum spectra, $S_\nu \propto \nu^\alpha$, with a slope
$\alpha$ $\sim$ 1, intermediate between the optically thick and thin
limits, small sizes ($\sim$ 0.01 pc), high temperatures (10$^4$ K),
high densities ($\geq$ 10$^6$ cm$^{-3}$), and very broad radio
recombination lines ($\Delta\nu_{FWHM}$ $\geq$ 40 km s$^{-1}$)
indicating both pressure broadening and the presence of bulk motions
of dense gas \citep[e.g.][]{Sewi2004}. On the other hand, the hot
molecular cores internally heated by massive protostars are small
($\sim$ 0.05 pc) and compact regions of dense molecular gas ($\geq$
10$^6$ cm$^{-3}$) which show very rich millimeter and submillimeter
spectra \citep[see][for reviews on this
  topic]{stan2000,van2004,Cesa2005a}. The hot molecular core phase
begins at much colder temperatures ($\sim$ 100 K) than the HC/UCHII
stage.

There are many cases where the UCHII or even the HCHII regions are
observed in association with hot cores, some of them are: G5.89-0.39
\citep{su2009}; G24 A1 \citep{Bel2006}; G10.47, and G31.41+0.31
\citep{ces2010,Rol2009}; W51e2, W51e8, NGC 7538 IRS1, G28.2 and G10.6
\citep{kla2009}; IRAS 17233-3606 \citep{leu2008}.  However, there are
some others cases where the hot cores do not show the presence of a
strong HCHII \citep[e.g. AFGL490, NGC 7538S, and IRAS
  20126;][]{naka1991,sch2002,cesa2005,san2010} possibly because the
infalling molecular gas partially quenches its formation
\citep[e.g.][]{Wal1995,oso1999} or maybe because the massive central
protostar(s) are/is in an early phase with not high enough
temperatures to ionize their/its surroundings.

W51 North is one of the youngest massive stars within the luminous
cluster W51-IRS2 that is still on the process of formation. This
object is associated with a strong mm and submm source, a hot
molecular core at an excitation temperature of $\sim$ 200 K
\citep{Zhangetal1998}, no centimeter free-free emission at all ($\leq$
5 mJy at 1.3 cm) that could be related with an HC/UCHII or a thermal
jet \citep{Gaumeetal1993, Zhangetal1998, Eisneretal2002}. From this
object emanates a powerful molecular outflow observed at very small
scales by masers spots of SiO and H$_2$O, and at large scales in
thermal SiO($J=5-4$) \citep{Schnepsetal1981, Eisneretal2002,
  Imaietal2002, zapataetal2009}. Recently, \citet{zapataetal2009}
using high angular 1.3 and 0.7 mm continuum, and SO$_2$ line
observations made with the Submillimeter Array (SMA) and the Very
Large Array (VLA) resolved the large dusty and molecular hot core in a
rotating Keplerian ring with a size of $\sim$ 8000 AU \citep[at a
  distance of $\sim$ 6 kpc:] []{Imaietal2002,barba2008,xu2009}. This
structure surrounds a compact central dusty circumstellar disk with a
smaller size of approximately 3000 AU, and where the molecular
northwest(redshifted)-southeast(blueshifted) bipolar outflow
emerges. However, the relationship between both circumstellar
structures was not clear.

In this paper, we present sensitive high angular resolution SMA
submillimeter and millimeter line observations of the W51 North region
that were made in an attempt to understand the nature of the
circumstellar molecular and dusty structures associated with this
young massive protostar.  In Section 2 discuss the observations
undertook in this study. In Section 3, we present and discuss our SMA
millimeter and submillimeter data. Finally, in Section 4, we give the
main conclusions of the observations presented here.

\begin{table*}[ht]
\begin{minipage}[t]{\columnwidth}
\scriptsize \renewcommand{\footnoterule}{}
\caption{Observational and physical parameters of the submillimeter
  and millimeter lines}
\begin{center}
\begin{tabular}{lccccccc}
\hline \hline & Rest frequency & E$_{lower}$ & Range of Velocities &
Linewidth & LSR Velocity\footnote{The linewidth and LSR velocity were
  obtained fitting a Gaussian to the spectra.} & Peak Flux \\ Lines &
[GHz] & [K] & [km s$^{-1}$] & [km s$^{-1}$] & [km s$^{-1}$] & [Jy
  Beam$^{-1}$] \\ \hline \hline C$_2$H$_5$OH(37$_{8,29}$-36$_{9,28}$)
& 216.8121...& 715 & $+$54,$+$67 & 7 & 59 &
0.1\\ HCOCH$_2$OH(62$_{13,49}$-62$_{12,50}$) & 216.7585...& 1476 &
$+$54,$+$66 & 8 & 59 & 0.4\\ CO($J=3-2$)\footnote{Spectra with
  multiple peaks and in absorption about a LSR velocity of 60 km
  s$^{-1}$.} & 345.7959... & 16 & $+$20,$+$120 & 50 & 60 & 3.0
\\ SO(8$_{9}$-7$_8$)\footnote{Spectra with two peaks.}  &
346.5284... & 21 & $+$30,$+$90 & 35 & 59 &
1.0\\ SO$_2$(19$_{1,19}$-18$_{0,18}$) & 346.6521... & 151 &
$+$41,$+$72 & 15 & 58 & 1.0\\ HC$_3$N(38-37) ($\nu_t$=0) &
345.6090... & 306 & $+$51,$+$72 & 9 & 60 & 2.0
\\ HCOCH$_2$OH(68$_{19,49}$-67$_{20,48}$) & 337.3988... & 1512 &
$+$48,$+$73 & 10 & 62 & 2.0 \\ \hline \hline
\end{tabular}
\end{center}
\end{minipage}
\end{table*}

\section{Observations}

\subsection{Millimeter}

The observations were obtained with the SMA\footnote{The Submillimeter
  Array (SMA) is a joint project between the Smithsonian Astrophysical
  Observatory and the Academia Sinica Institute of Astronomy and
  Astrophysics, and is funded by the Smithsonian Institution and the
  Academia Sinica.}  during 2008 January 17. The SMA was in its very
extended configuration, which included 28 independent baselines
ranging in projected length from 30 to 385 k$\lambda$.  The phase
reference center of the observations was R.A. = 19h23m40.05s, decl.=
14$^\circ$31$'$05.0$''$ (J2000.0).  The frequency was centered at
217.1049 GHz in the Lower Sideband (LSB), while the Upper Sideband
(USB) was centered at 228.1049 GHz. The primary beam of the SMA at 230
GHz has a FWHM of about 60$''$.

We detected the lines C$_2$H$_5$OH(37$_{8,29}$-36$_{9,28}$) and
HCOCH$_2$OH(62$_{13,49}$-62$_{12,50}$) in the LSB. See Table 1 for
their rest frequencies and rotational temperatures above of the ground
state. The full bandwidth of the SMA correlator was 4 GHz (2 GHz in
each band). The SMA digital correlator was configured in 24 spectral
windows (``chunks'') of 104 MHz each, with 256 channels distributed
over each spectral window, providing a resolution of 0.40 MHz (0.56 km
s$^{-1}$) per channel. However, in this study we smoothed our spectral
resolution to about 1 km s$^{-1}$.

The zenith opacity ($\tau_{230 GHz}$), measured with the NRAO tipping
radiometer located at the Caltech Submillimeter Observatory, was
$\sim$ 0.15, indicating good weather conditions during the
observations.  Observations of Uranus provided the absolute scale for
the flux density calibration.  Phase and amplitude calibrators were
the quasars 1925+211 and 2035+109.  The uncertainty in the flux scale
is estimated to be 15--20$\%$, based on the SMA monitoring of quasars.
Observations of Uranus provided the absolute scale for the flux
density calibration.  Further technical descriptions of the SMA and
its calibration schemes can found in \citet{Hoetal2004}.

The data were calibrated using the IDL superset MIR, originally
developed for the Owens Valley Radio Observatory
\citep{Scovilleetal1993} and adapted for the SMA.\footnote{The MIR-IDL
  cookbook by C. Qi can be found at
  http://cfa-www.harvard.edu/$\sim$cqi/mircook.html} The calibrated
data were imaged and analyzed in the standard manner using the MIRIAD
and AIPS packages.  We used the ROBUST parameter of the INVERT task of
MIRIAD set to 2 to obtain a slightly better sensitivity while losing
some angular resolution. The resulting image rms noise of line images
was around 30 mJy beam$^{-1}$ for each velocity channel (with a
smoothed size of 1 km s$^{-1}$) at an angular resolution of
$0\rlap.{''}57$ $\times$ $0\rlap.{''}42$ with a P.A. = $57.6^\circ$.

\begin{figure}[ht]
\begin{center}
\includegraphics[scale=0.4]{f1.eps}
\caption{\scriptsize Overlay of the SO$_2$(22$_{2,20}$-22$_{1,21}$)
  moment zero emission (color image), the 7 mm continuum emission from
  the compact dusty disk (brown contours), and the CO($J=3-2$) moment
  zero emission from the bipolar outflow (blue and red contours). The
  SO$_2$ and the 7 mm images were taken from
  \citet{zapataetal2009}. The integrated velocity range for the
  CO($J=3-2$) is from +20 to +50 km s$^{-1}$ for the blueshifted gas,
  and from +70 to +120 km s$^{-1}$ for the redshifted gas. The brown
  contours are the 10\%, 13\%, and 16\% of the peak of the 7 mm
  continuum emission, the peak is 100 mJy beam$^{-1}$.  The
  synthesized beam of the CO($J=3-2$) contour image is 0.71$''$
  $\times$ 0.57$''$ with a P.A. of 87.6$^\circ$ and is shown in the
  bottom right corner. The green triangle marks the position of the
  possible companion reported in this paper. The white dashed line
  shows how the blueshifted side of the outflow seems to be more
  likely energized from the central massive protostar W51 North rather
  than the possible companion.}
\label{fig1}
\end{center}
\end{figure}

\begin{figure*}[ht]
\begin{center}
\includegraphics[scale=0.31]{f2.eps}
\caption{\scriptsize Spectra from the selected millimeter and
  submillimeter lines observed toward W51 North. The red line
  represents in all panels the synthetic spectrum obtained for the
  best fitting solution with the parameters shown in each panel.  The
  spectra was obtained from the average of the total emission in each
  line using the task ``imspect'' of MIRIAD.}
\label{fig2}
\end{center}
\end{figure*}

\begin{figure*}[ht]
\begin{center}
\includegraphics[scale=0.20]{f3.eps}
\caption{\scriptsize Composite images of the molecular emission from
  the circumstellar large disk around the massive protostar W51 North.
  Left: The blue color scale and the contours are representing the
  integrated emission of the SO(8$_{9}$-7$_8$), green ones the
  HC$_3$N(38-37) ($\nu_t$=0), and the red ones the
  HCOCH$_2$OH(62$_{13,49}$-62$_{12,50}$).  The blue contours are from
  5\% to 90\% with steps of 10\% of the peak of the line emission, the
  peak is 53 Jy beam$^{-1}$ km s$^{-1}$.  The green contours are from
  30\% to 90\% with steps of 10\% of the peak of the line emission,
  the peak is 26 Jy beam$^{-1}$ km s$^{-1}$.  The red contours are
  from 30\% to 90\% with steps of 10\% of the peak of the line
  emission, the peak is 3.5 Jy beam$^{-1}$ km s$^{-1}$.  The
  synthesized beams of the millimeter and submillimeter observations
  are shown on the right bottom corner. Right: The light blue color
  scale and the contours are representing the integrated emission of
  the SO$_2$(19$_{1,19}$-18$_{0,18}$), green ones the
  C$_2$H$_5$OH(37$_{8,29}$-36$_{9,28}$), and the orange ones the
  HCOCH$_2$OH(68$_{19,49}$-67$_{20,48}$).  The light blue contours are
  from 30\% to 90\% with steps of 10\% of the peak of the line
  emission, the peak is 38 Jy beam$^{-1}$ km s$^{-1}$.  The green
  contours are from 30\% to 90\% with steps of 10\% of the peak of the
  line emission, the peak is 2 Jy beam$^{-1}$ km s$^{-1}$.  The orange
  contours are from 30\% to 90\% with steps of 10\% of the peak of the
  line emission, the peak is 27 Jy beam$^{-1}$ km s$^{-1}$.  The
  synthesized beam of the three measurements is the same and is shown
  on the right bottom corner.  The triangles represent the position of
  the two possible embedded massive protostars as suggested by the hot
  molecular emission of the HCOCH$_2$OH(62$_{13,49}$-62$_{12,50}$) and
  HC$_3$N(38-37) ($\nu_t$=0). The white dashed line in both panels
  represents the position and orientation where the PV-diagrams shown
  in Figures 4 and 5 were computed. }
\label{fig3}
\end{center}
\end{figure*}

\begin{table*}[ht]
\begin{minipage}[t]{\columnwidth}
\scriptsize \renewcommand{\footnoterule}{}
\caption{Observational parameters of the circumstellar multilayer
  disk}
\begin{center}
\begin{tabular}{lccccccc}
\hline \hline & \multicolumn{2}{c}{Central Position} & Flux Density &
\multicolumn{2}{c}{Deconvolved Size\footnote{The Gaussian fitting was
    obtained with the task JMFIT of AIPS.}} & Approx. cavity
size\\ \hline & $\alpha(2000)$ & $\delta(2000)$ & & & & \\ Specie & [h
  m s] & [$\circ$ $''$ $'$] & [Jy Beam$^{-1}$ km s$^{-1}$] &
[arcsec$^2$]& [Degrees] & [arcsec$^2$] \\ \hline \hline
SO(8$_{9}$-7$_8$) & 19 23 40.074 & 14 31 05.70 & 341$\pm$20 &
1.78$\pm$0.01 $\times$ 1.39$\pm$0.02 & 171$\pm$2 & $\sim$
0.8\\ SO$_2$(19$_{1,19}$-18$_{0,18}$) & 19 23 40.071 & 14 31 05.54 &
234$\pm$15 & 1.65$\pm$0.04 $\times$1.49$\pm$0.03 & 24$\pm$6& $\sim$
0.7\\ HC$_3$N(38-37) ($\nu_t$=0) & 19 23 40.060 & 14 31 05.54 &
88$\pm$15 & 1.36$\pm$0.04 $\times$0.93$\pm$0.04 &135$\pm$7 & $\sim$
0.6\\ C$_2$H$_5$OH(37$_{8,29}$-36$_{9,28}$) & 19 23 40.041 & 14 31
05.49 & 116$\pm$15 & 1.31$\pm$0.03 $\times$0.89$\pm$0.03 & 107$\pm$7&
$\sim$ 0.8\\ HCOCH$_2$OH(68$_{19,49}$-67$_{20,48}$) & 19 23 40.054 &
14 31 05.44 & 105$\pm$13 & 1.14$\pm$0.03 $\times$0.88$\pm$0.03 &
99$\pm$3& --\\ HCOCH$_2$OH(62$_{13,49}$-62$_{12,50}$) &19 23 40.062 &
14 31 05.59 & 9.1$\pm$2 & 0.81$\pm$0.03 $\times$ 0.40$\pm$0.09 &
75$\pm$3 & --\\ ---Companion & 19 23 40.077 & 14 31 05.63 & 8.5$\pm$1
& 1.27$\pm$0.04 $\times$0.39$\pm$0.03 & 74$\pm$2& --\\ \hline & & &
Dusty circumstellar disk\footnote{Data obtained from
  \citet{zapataetal2009}. } & & & \\ 7 mm & 19 23 40.057 & 14 31 05.67
& -- & 0.58$\pm$0.02 $\times$0.27$\pm$0.02 & 70$\pm$6&
--\\ \hline\hline
\end{tabular}
\end{center}
\end{minipage}
\end{table*}

\subsection{Submillimeter}

The observations were obtained with the SMA on 2008 July 13.  At the
time of these observations the SMA had seven antennas in its extended
configuration with baselines ranging in projected length from 30 to
255 k$\lambda$.  The primary beam of the SMA at 340 GHz has a FWHM of
37$''$. The molecular emission from the hot core was found well inside
of the primary beam.

The receivers were tuned to a frequency of 345.796 GHz in the upper
sideband (USB), while the lower sideband (LSB) was centered on 335.796
GHz.  The LSB contained line HCOCH$_2$OH(68$_{19,49}$-67$_{20,48}$)
while the USB the lines SO(8$_{9}$-7$_8$),
SO$_2$(19$_{1,19}$-18$_{0,18}$), CO($J=3-2$), and HC$_3$N(38-37)
($\nu_t$=0). See Table 1 for their rest frequencies and lower level
energies above the ground energy state.  The SMA digital correlator
was configured in 24 spectral windows (``chunks'') of 104 MHz each,
with 128 channels distributed over each spectral window, thus
providing a spectral resolution of 0.81 MHz (0.70 km s$^{-1}$) per
channel.  However, we smoothed our spectral resolution to 1.0 km
s$^{-1}$ per spectral channel.

The zenith opacity ($\tau_{230 GHz}$) was $\sim$ 0.035 -- 0.04,
indicating excellent weather conditions.  Observations of titan
provided the absolute scale for the flux density calibration.  Phase
and amplitude calibrators were the quasars 1751$+$096, and 1925$+$211.
The uncertainty in the flux scale is also estimated to be 15 --
20$\%$, also based on the SMA monitoring of the quasars.

The calibrated data were imaged and analyzed in standard manner using
the MIRIAD, and AIPS packages. We also used the ROBUST parameter set
to 2.  The line image rms noise was around 100 mJy beam$^{-1}$ for
each velocity channel (with a smoothed size of 1 km s$^{-1}$) at an
angular resolution of $0\rlap.{''}78$ $\times$ $0\rlap.{''}58$ with a
P.A. = $84.7^\circ$.

\section{Results}

In Figure 1, we present the disk/ring/outflow system found towards W51
North using observations of the SiO($J=5-4$),
SO$_2$(22$_{2,20}$-22$_{1,21}$), and millimeter continuum by
\citet{zapataetal2009}. However, in this image, we have exchanged the
SiO emission by the emission of the CO($J=3-2$) presented here for the
first time. The CO($J=3-2$) confirms the presence of a single powerful
and collimated (15$^\circ$) bipolar outflow (with the blueshifted side
towards the southeast while the redshifted one to the north) emanating
from the disk/ring structure at similar scales. However, the
CO($J=3-2$) reveals that the outflow extents for a few arcseconds more
on both sides in comparison with the emission traced by SiO($J=5-4$),
and shows a slightly larger radial velocity range ($+$30 to $+$120 km
s$^{-1}$).  In this image, it is more clear to see the deviation of
the blueshifted side of the outflow with respect to the redshifted one
than in \citet{zapataetal2009}.

Figure 2 shows the spectra of the multiple lines reported here with
exception of the CO($J=3-2$).  Additionally, we present synthetic LTE
models for the spectra using XCLASS program
\citep[][]{Comitoetal2005}.  The data were fitted well with the
synthetic model using rotational temperatures, column densities, and
sizes of 250 K, 2 $\times$ 10$^{16}$ cm$^2$, and 1$''$, respectively,
for the lines SO(8$_{9}$-7$_8$), SO$_2$(19$_{1,19}$-18$_{0,18}$),
HC$_3$N(38-37) ($\nu_t$=0), and C$_2$H$_5$OH(37$_{8,29}$-36$_{9,28}$)
and 800 K, 5 $\times$ 10$^{18}$ cm$^2$, and 0.5$''$ for transitions
with high excitation temperatures in the lower energy states,
HCOCH$_2$OH(68$_{19,49}$-67$_{20,48}$) and
HCOCH$_2$OH(62$_{13,49}$-62$_{12,50}$).  This modeling serves as for
avoiding the line confusion, and suggest that hot gas is only found in
the innermost parts of the hot core.  See Table 1 for the rest
frequencies and the physical parameters of the lines. The physical
parameters of all lines were obtained fitting a Gaussian to the line
profile.

In order to be more confident about the detection of the
glycol-aldehyde (HCOCH$_2$OH) toward W51 North, we have fitted
simultaneously many other lines that fall in our 2 GHz bands with our
synthetic LTE model using the same parameters as above. We obtained a
reasonable agreement between the synthetic model and the detected
lines.  A full astrochemical analysis of these lines (and other
species) is beyond the scope of the present paper in which we
concentrate on the spatial distributions of the lines presented in
Table 1.

The spectra of the SO(8$_{9}$-7$_8$) additionally shows
self-absorption possibly due to the molecular gas is falling into the
central object \citep[][]{Zapataetal2008}.

\begin{figure*}[ht]
\begin{center}
\includegraphics[scale=0.37, angle=0]{f4.eps}
\caption{\scriptsize Position-velocity diagrams of the HC$_3$N(38-37)
  (green scale color), and the HCOCH$_2$OH(68$_{19,49}$-67$_{20,48}$)
  (pink scale color) line emission from the circumstellar disk. The PV
  diagrams were computed at a P.A. = 30$^\circ$ (see Figure 3). These
  diagrams are additionally overlaid with the position-velocity
  diagram at same P.A. of our Keplerian disk model (contours). The
  contours are from 20\% to 90\% with steps of 10\% of the peak of the
  line emission of our model. The units of the horizontal axis are in
  arcseconds. The systemic LSR radial velocity of the ambient
  molecular cloud is about 60 km s$^{-1}$, see Table 1. The angular
  and spectral resolutions are shown in each panel on the bottom left
  corner.  The synthesized beam of both images is 0.71$''$ $\times$
  0.57$''$ with a P.A. of 87.6$^\circ$. The spectral resolution was
  smoothed to 1 km s$^{-1}$.  The origin in the horizontal axis is at
  R.A. = 19h23m40.05s, decl.= 14$^\circ$31$'$05.0$''$ (J2000.0).}
\label{fig4}
\end{center}
\end{figure*}

\begin{figure*}[ht]
\begin{center}
\includegraphics[scale=0.4]{f5.eps}\\
\caption{\scriptsize Position-velocity diagrams for the
  SO(8$_{9}$-7$_8$) and the SO$_2$(19$_{1,19}$-18$_{0,18}$) from the
  large circumstellar disk and the outflow.  They were computed at a
  P.A. of 160$^\circ$ (see Figure 3). The contours in both panels are
  from 10\% to 90\% with steps of 10\% of the peak of the line
  emission. The units of the horizontal axis are in arcseconds. The
  systemic LSR radial velocity of the ambient molecular cloud is about
  60 km s$^{-1}$, see Table 1. The angular and spectral resolutions
  are shown in each panel on the bottom left corner.  The synthesized
  beam of both images is 0.71$''$ $\times$ 0.57$''$ with a P.A. of
  87.6$^\circ$. The spectral resolution was smoothed to 1 km
  s$^{-1}$. The origin in the horizontal axis is at R.A. =
  19h23m40.05s, decl.=14$^\circ$31$'$05.0$''$ (J2000.0).}
\end{center}
\label{fig5}
\end{figure*}

The integrated intensity or moment zero maps of the lines associated
with the molecular disk are presented in Figure 3. We have divided
these lines into three groups showing the transitions with energy
levels with moderate (20-160 K), high (300-800 K), and very high
(1400-1600 K) energies and presented them in two different panels.
From these images, it is easy to see how the the molecular emission is
found in layers with the transitions characterized by high excitation
temperatures in their lower energy states (up to 1512 K), and/or
critical densities, being concentrated closer to the central high-mass
protostar W51 North.  This phenomenon confirms clearly that the
heating in W51 North is internal.  Additionally, molecules with
transitions characterized by low and moderate excitation temperatures
{\it i.e.} SO(8$_{9}$-7$_8$), SO$_2$(19$_{1,19}$-18$_{0,18}$),
HC$_3$N(38-37) ($\nu_t$=0), and C$_2$H$_5$OH(37$_{8,29}$-36$_{9,28}$)
show a central cavity of approximately an angular size of 0.7$''$, see
Table 2. These lines trace thus the molecular Keplerian ring (or
toroid) around W51 North reported in SO$_2$ by
\citet[][]{zapataetal2009}, see Figure 1.  The transition of the
molecule SO$_2$ utilized by these authors has a similar excitation
temperature in the lower energy state as those mentioned above. Note
the similarity between the angular size of the inner cavity and the
size of the compact dusty circumstellar disk reported by
\citet[][]{zapataetal2009}, see Table 2.

The molecular emission from the transition with a high excitation
temperature in the lower energy states ({\it i.e.}
HCOCH$_2$OH(68$_{19,49}$-67$_{20,48}$)) is found peaking at the
position of the cavity and showing similar angular sizes to this, see
Table 2. This molecule is thus more intimately related with the
circumstellar compact dusty disk.

The observational parameters of all lines are shown in Table 2.  We
noted from Table 1 and 2, there is a correlation between the
deconvolved sizes of the molecular emission and the excitation
temperatures in lower energy states of each molecular specie. The
molecular emission from transitions characterized by high excitation
temperatures show to be very compact.  This correlation could be
obviously obtained for lines that are optically thick with T$_b$
$\sim$ $\eta$ T$_{ex}$, where the T$_b$ is the brightness temperature,
$\eta$ is the filling factor, and T$_{ex}$ is the excitation
temperature.

In the left panel of Figure 3 the emission of the lines
HCOCH$_2$OH(62$_{13,49}$-62$_{12,50}$) and HC$_3$N(38-37) ($\nu_t$=0)
reveals the presence of a possible warm ``companion'' located to the
northeast of the disk. We have marked the position of this putative
companion with a yellow triangle. This ``companion'' is also observed
in SO$_2$ at the same position, see Figure 3 of
\citet[][]{zapataetal2009}, and is associated with a group of water
maser spots \citep{Imaietal2002,Eisneretal2002}. The position of the
``companion'' was found by fitting a Guassian to its
HCOCH$_2$OH(62$_{13,49}$-62$_{12,50}$)) emission. However, it still
not clear if the ``companion'' is real protostar or if this could be
the result of the interaction of the outflow with a high density zone
of the molecular cloud.

We do not find any clear evidence of outflowing gas activity
associated with this possible ``companion'' \citep[see
  also][]{zapataetal2009}. Although the blueshifted side of the
CO($J=3-2$) outflow has not the same position angle as the redshifted
one, this seems not be ejected from the companion.  We show how the
blueshifted side appears to be more likely ejected from W51 North in
Figure 1. We drew a line that crosses this side of the outflow and
points directly to the dusty compact disk.  The SO$_2$ and SiO
emission indeed show how the outflow is deviated to where the
CO($J=3-2$) is located \citep[][]{zapataetal2009}.  However, more
observations in some other molecular outflow tracers are thus
necessity to firmly discard the existence of a second outflow in W51
North. It is interesting to note that the possible companion is not
observed in the hotter molecular gas tracer
HCOCH$_2$OH(68$_{19,49}$-67$_{20,48}$) suggesting that it may not as
warm as the central massive star.

Figures 4 and 5 show the position-velocity diagrams (PV-diagrams) of
different molecules which trace distinct scales of the disk. In Figure
3, we have marked the orientation and position of the PV-cuts. In the
left panel of Figure 3 we show a white line with a P.A.=30$^\circ$
that corresponds to the PV-cuts shown in Figure 4. On the other hand,
the white line with a P.A.=160$^\circ$ in the right panel corresponds
to the the PV-cuts shown in Figure 5.

Figure 4 shows the molecular emission from HC$_3$N(38-37) and
HCOCH$_2$OH(68$_{19,49}$-67$_{20,48}$) located in the innermost parts
of the disk, while the PV-diagrams of molecules as SO(8$_{9}$-7$_8$),
and SO$_2$(19$_{1,19}$-18$_{0,18}$), which trace the outermost parts,
are presented in Figure 5. The PV-diagrams in Figure 4 reveal that the
hot molecular gas closer to the massive protostar is Keplerian.  In
addition, in this figure we have overlaid the PV-diagram of the LTE
Keplerian disk modeled in \citet[][]{zapataetal2009}, but without an
inner cavity and a smaller size. Both structures shown a very good
correspondence.  The SO(8$_{9}$-7$_8$) and
SO$_2$(19$_{1,19}$-18$_{0,18}$) presented in Figure 5 trace much
larger structures similar to the ring reported in
\citet[][]{zapataetal2009}. These molecules in addition show clearly
two northwest and southeast high velocity extensions excited by the
bipolar outflow mapped in CO($J=3-2$) and SiO($J=5-4$).  The
PV-diagram of the molecule C$_2$H$_5$OH(37$_{8,29}$-36$_{9,28}$) shows
also Keplerian motions, while that obtained from the molecule
HCOCH$_2$OH(62$_{13,49}$-62$_{12,50}$) shows a more compact structure
without a clear morphology, these diagrams are not presented in this
study.

Finally, in Figure 6 we show the spectral energy distribution (SED)
from the centimeter, millimeter, and submillimeter wavelengths of W51
North with data obtained from \citet[][]{zapataetal2009} and Tang et
al. (in prep.). The new value at submillimeter wavelengths presented
here shows also to be good fitted with the $\alpha$ = 2.8 value (where
the flux density goes as S$_\nu$=$\nu^\alpha$) found by
\citet[][]{zapataetal2009}.

\begin{figure}[h!]
\begin{center}
\includegraphics[scale=0.6]{f6.eps}\\
\caption{\scriptsize SED for the source W51 North combining 3.6 cm,
  1.3 cm, 7 mm, 1.3 mm and 0.8 mm continuum data. The millimeter data
  were obtained from \citet{zapataetal2009} and the submillimeter data
  was obtained from Tang et al. (in prep.). The respective error bars
  were smaller than the squares and are not presented here. The line
  is a least-squares power-law fit (of the form S$_\nu$ $\propto$
  $\nu^\alpha$) to the spectrum. Note that the value of $\alpha$=2.8
  obtained by \citet{zapataetal2009} fit very well all millimeter and
  submillimeter measurements. The 1.3 cm point is a upper limit
  (5$\sigma$=5 mJy) obtained from Figure 1 of
  \citet{Eisneretal2002}. The 3.6 cm point is a upper limit
  (5$\sigma$=4 mJy) obtained from Figure 5 of \citet{meh1994}.}
\end{center}
\label{fig6}
\end{figure}

\begin{figure*}[ht]
\begin{center}
\includegraphics[scale=0.75, angle=0]{f7.eps}\\
\caption{\scriptsize A very simple diagram showing the different
  evolutionary phases of a massive protostar. In this diagram, we do
  not included multiplicity for simplicity. }
\end{center}
\label{fig7}
\end{figure*}

\section{Discussion}

\subsection{A multilayer infalling and keplerian disk around W51 North:
Connecting the molecular ring and dusty compact disk.}

The combined millimeter and submillimeter line observations presented
here from the ring and compact dusty disk structures reported in
\citet{zapataetal2009} suggest that we are seeing a single hot
multilayer accreting disk around W51 North (Figure 3). This large
molecular disk connects both circumstellar structures, the dusty
compact disk, observed at 7 and 1 mm, and the molecular ring mapped in
SO$_2$ in a single structure.  Some molecules are found in association
with the dusty compact disk while some others with the ring.

The molecular emission associated with the multilayer disk reveals
that the molecular emission with transitions with a low excitation
temperature in the lower energy states delineates the outer disk and
shows a cavity, while such transitions with a high excitation
temperature in the lower energy state are found only in the innermost
parts of the disk and closer to the central massive protostar(s).  One
would expect to this physical phenomenon happen in disks or cores
around high- and low-mass protostar(s) due to the large optical depths
and higher temperatures only found close to the young star. Up to date
detailed models of the ensuing chemical evolution at these much higher
temperatures do not at present exist.

It still is not clear why we do see an inner cavity on the disk in
some molecules with transitions characterized by low and moderate
excitation temperatures.  This cavity maybe is not a real hole on the
disk rather this could be due to photo-dissociation or opacity effects
of these molecules toward the innermost parts of the disk.  However,
we do not discard the possibility that the cavity could be real and
formed possibly by the tidal effects of the young multiple massive
protostars in the middle of the ``ring'' as in the case of Ori 139-409
located in Orion South where the protostars have lower masses
\citep{zap2010}.  The molecular emission with transitions
characterized by high excitation temperatures in the lower energy
states may arise from the possible compact circumstellar disks inside
of the cavity that with the actual angular resolution are not
resolved.  The presence of a binary system at the center of W51 North
has been suggested by the precession of the molecular outflow at very
small scales found by \citet{Eisneretal2002}. This precession
phenomenon may also explain the difference between the position angles
of the blueshifted and redshifted components of the CO($J=3-2$)
outflow (Figure 1).

The accreting disk hypothesis is supported by the detection of the
powerful outflow that emanates with an almost perpendicular
orientation to the object.  The large disk shows to be Keplerian and
with a contracting velocity of about a few km s$^{-1}$ as revealed by
\citet[][]{zapataetal2009}.  The models presented in Figure 4 were
obtained using a slightly different P.A. = 30$^\circ$, inclination
angle of {\it i}=25$^\circ$, a smaller size (3000 AU), and without an
inner cavity than in \citet[][]{zapataetal2009}.

In \citet[][]{zapataetal2009} and here, we have modeled a contracting
flattened disk in Keplerian rotation with a central hole, using the
disk parametrization from Guilloteau, Dutrey \& Simon (1999).  The
contraction is assumed to have the functional form of free-fall
(i.e. V$_{inf}$ $\propto$ $\frac{1}{\sqrt r}$), with a reference
velocity at the reference radius. However, our resolution is
insufficient to determine the exact functional form.  The model is for
the same molecules and transitions.  A better fit than all other
trials in our recurrence was found until we obtained similar
structures in our model to those imaged (Figure 4).  Most of the
physical parameters in the model were constrained in the process. The
model fits the observations reasonably well.

In particular, the molecules HCOCH$_2$OH(62$_{13,49}$-62$_{12,50}$)
and HC$_3$N(38-37) ($\nu_t$=0) revealed the possible presence of a
warm companion located to the northeast of the disk (Figure 2).  We
suggest that the companion could be a consequence of disk
fragmentation due gravitational instabilities because the disk is
extremely large, see for a reference of this phenomenon
\citet{kra2006,kru2009}.

\subsection{A possible sequence on the evolution of nascent massive stars?}

The multilayer disk imaged in this paper reveals a possible link
between the Hyper/Ultra-compact HII regions and the hot cores.  These
observations suggest a possible sequence on the formation of massive
stars. We present this simple sequence in Figure 7, and discuss it as
follows:

 \begin{itemize}

\item
{\it Phase I:} In the first phase, a large and massive pseudo-disk is
formed together with a bipolar outflow from a large and dense core.
The pseudo-disk is surrounded by an infalling envelope with accretion
rates on the order of 10$^{-3}$ M$_\odot$ yr$^{-1}$.  The molecular
emission arising from the circumstellar pseudo-disk could be observed
as a single structure without a well defined temperature gradient
across it. The pseudo-disks observed toward these objects are not
classical cicumstellar disks associated with low-mass stars which are
centrifugally supported, rather these might be contracting
pseudo-disks \citep[see for example the large and massive contracting
  disks associated with the high-mass protostars AFGL490 and NGC
  7538S:][] {naka1991,san2010}. The pseudo-disk could be also likely
circumbinary or circum-multiple because of the large multiplicity
presented on the massive stars.

\item
{\it Phase II:} The large and massive pseudo-disk formed earlier is
observed now in molecular shells or layers. This object is still
surrounded by an infalling envelope.  A clear temperature gradient is
observed across the disk.  The inner cavity could be formed by
photodissociation due to the high temperatures close to massive stars
or maybe to opacity effects.  In this phase, no free-free emission
from a HCHII region is detected probably because the infalling
molecular gas partially quenches its formation or maybe because the
massive central protostar(s) are/is in an early phase with not high
enough temperatures to ionize their/its surroundings. W51 North is
found here. Another possibility is that maybe there is cluster of late
B-type in middle able of forming a strong multilayer molecular disk,
but not for developing HCHII regions.

\item
{\it Phase III:} A HCHII region appears in the middle of the large
multilayer pseudo-disk. At this stage the central massive protostar
has the sufficient high temperatures ($>$ 10$^4$ K) to ionize its
surroundings, see for example the cases of NGC7538 IRS1
\citep{fra2004,San2009} and MWC349 \citep{tafo2004}. These sources are
maybe variable on time because of the strong accretion from the large
disk continues toward the massive central star, see also G24.78+0.08
A1 \citep{gal2008}. There are some sources that might be situated at
this phase, see for example, the molecular rotating toroids observed
around the UCHII regions G24.78+0.08 A1 and G5.89-0.39
\citep{su2009,Bel2006}. Moreover, \citet{kla2009} found in a group of
HCHII and UCHII regions rotating molecular structures surrounding
them. Therefore, the molecular toroids or rings around an UCHII are
likely the same structures as in Phase II, but just hotter.  The holes
observed in these toroids might be due to the photodissociation of
some molecules, or maybe to the formation of cavities by the presence
of a multiple system of massive stars as mentioned before. In this
phase the masive young star(s) might be still accreting ionized
material as proposed and described by \citet{Ke2006}.

\item
{\it Phase IV:} In this phase the HCHII region expands and ionizes the
remnant molecular disk or toroid, and creates an extended H II region
surrounding one or maybe multiple young massive stars. The sequence
from UCHII to extended HII regions has been discussed by
\citet{gar2004,fra2007}. If the central star(s) is/are still accreting
material, probably ionized material, an outflow show also be present
even at this stage.

 \end{itemize}

\section{Summary}

We have observed the massive and young protostar W51 North with the
SMA using a high angular resolution and sensitivity at millimeter and
submillimeter wavelengths. We give a summary of the results as
follows:

\begin{itemize}
\item
We report the presence of a large and single molecular disk around the
object W51 north, and confirm the existence of a single powerful
outflow emanating from the disk;

\item
The molecular emission from the large disk is observed in layers with
transitions characterized by high excitation temperatures in the lower
states (up to 1512 K) being concentrated closer to the central massive
protostar. The molecular emission at low or moderate excitation
temperatures additionally exhibits a central cavity with an angular
size around 0.7$''$.  This multilayer disk connects the molecular ring
and dusty compact disk reported recently toward this object
\citep[][]{zapataetal2009}. The LTE modeling presented in this study
also confirms that the hot gas is only found in the innermost disk.

\item
A detail study of the kinematics of the molecular gas together with a
model of a circumstellar disk in Local Thermodynamic Equilibrium shows
that the disk is Keplerian, and is accreting fresh gas to the
protostar possibly quenching the formation of a HCHII region;

\item
The thermal emission of the HCOCH$_2$OH reveals the possible presence
of a warm companion located to the northeast of the disk, perhaps
produced by disk fragmentation or maybe is the interaction of the
outflow a high density zone.  The molecular emission of the SO and
SO$_2$ is observed in the circumstellar disk as well as in the
outflow;

\item
We modeled all lines with a LTE synthetic spectra.  For an assumed
source size of 0.5 - 1 arcsec, our modeling yields a column density of
2 - 6 $\times$ 10$^{16-18}$ cm$^{-2}$, a temperature of 250 - 800 K,
and a linewidth of $\sim$ 5 - 10 km s$^{-1}$.

\end{itemize}

We thank the anonymous referee for many valuable suggestions.  This
research has made extensive use of the SIMBAD database, operated at
CDS, Strasbourg, France, and NASA's Astrophysics Data System.  This
research made use of the myXCLASS program
(https://www.astro.uni-koeln.de/projects/schilke/XCLASS), which
accesses the CDMS (http://www.cdms.de) and JPL
(http://spec.jpl.nasa.gov) molecular data bases.

\bibliographystyle{aa}