\documentclass{article}
\usepackage[affil-it]{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
\usepackage{url}
\usepackage{hyperref}
\hypersetup{colorlinks=false,pdfborder={0 0 0}}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[english]{babel}

\newcommand{\truncateit}[1]{\truncate{0.8\textwidth}{#1}}
\newcommand{\scititle}[1]{\title[\truncateit{#1}]{#1}}



\begin{document}

\title{Effects of an Embedded B-star Wind on the Properties of the Near-by Cloud: Ophiuchus}


\author{Hope How-Huan Chen}
\affil{Affiliation not available}

  
\author{Alyssa Goodman}
\affil{Affiliation not available}
  



\date{\today}

\maketitle

\selectlanguage{english}
\begin{abstract}
\textbf{Abstract}. $\rho$ Ophiuchii is a group of five B-stars, embedded in a nearby molecular cloud:  Ophiuchus, at a distance of $\sim$ 119 pc. A ``bubble''-like structure is found in dust thermal emission around $\rho$ Oph. The circular structure on the H$\alpha$ map further indicates that this bubble is physically connected to the source at the center. The goal of this paper is to estimate the impact of feedback from these embedded B-stars on the molecular cloud, by comparing the energy associated with the material entrained in the bubble to the total turbulent energy of the cloud. In this paper, we combine data from the COMPLETE Survey, which includes $^{12}$CO (1-0) and $^{13}$CO (1-0) molecular line emission from FCRAO, an extinction map derived from 2MASS near-infrared data using the NICER algorithm, and far-infrared data from IRIS (60/100 $\mu$m) with data from the Herschel Science Archive (PACS 100/160 $\mu$m and SPIRE 250/350/500 $\mu$m). With the wealth of data tracing different components of the cloud, we try to determine the best strategy to derive physical properties and to estimate the energy budget in the shell and in the cloud. We also experiment with the hierarchical Bayesian-fitting technique introduced by \citet{uffmann_Goodman_Launhardt_2012} in an effort to eliminate the bias in the derived column densities and/or temperatures induced by noise in the far-IR data. We find that the energy entrained in the bubble is $\sim$ 12 \% of the total turbulent energy of the Ophiuchus molecular cloud. This fraction is similar to the number\citet{Arce_2011} give for the Perseus molecular cloud, and it suggests the non-negligible role of B-stars in driving the turbulence in clouds. We expect that a complete survey of ``bubbles'' in the Ophiuchus cloud will reveal the importance of B-star winds in molecular clouds.%
\end{abstract}%



\section{Introduction} 
Stars go through many periods of mass loss throughout their lives. Young embedded protostars drive powerful collimated outflows, and as they evolve, the mass loss rate decreases and the outflows become less collimated \citep{Bontemps_1996,Arce_2006}. Wide-angle winds, sometimes coexistent with a more collimated, jet-like component, are observed around evolved pre-main sequence stars (T-Tauri and Herbig Ae/Be stars). It was realized early on in the study of winds that both the outflows driven by protostars and winds from T-Tauri stars have the potential to significantly affect the dynamics and structure of the parent cloud \citep{Norman_1980}, and the impact has been observed and analysed by Arce et al (2007) and \citet{Arce_2010}. However, the original bubble model proposed by \citet{Norman_1980} lost its appeal when T-Tauri stars were more often observed to have an effect on their surroundings through what are now well-known bipolar outflows.

However, parsec-scale circular cavities (typically referred to as ``shells'' or ``bubbles'') are regularly found in regions of high-mass star formation, and are likely created by spherical winds from high-mass stars \citep{Churchwell_2006,Churchwell_2007,Beaumont_2010}. High-mass stars evolve faster than low-mass stars, and they reach the main sequence while they are still accreting material and embedded in their parent cloud. Mass loss rates of $\sim 10^{-6}$ M$_{\odot}$ yr$^{-1}$ that can create these bubbles have been observed around high-mass stars during their main sequence phase. Previously, people thought that the same process did not happen in regions of intermediate- and/or low-mass star formation. However, \citet{Arce_2011} find that in the Perseus cloud it is possible, and frequent, for B or later types of main sequence stars and evolved pre-main sequence stars to drive spherical winds that can have a significant impact on the natal cloud. In their estimates of the impact, comparisons between the energy entrained in these shells and the turbulent energy in the cloud show that the embedded stellar winds have enough energy to disturb the cloud and contribute to the turbulence therein.

Reliable estimates of density and  dynamics are critical to measuring the significance of wind-cloud interactions. \citet{Goodman_2009} find that mass estimates can vary significantly due to biases inherent to various techniques for measuring column density in molecular clouds, including near-infrared extinction, mid/far-infrared thermal emission and molecular line emission.  These uncertainties and biases in mass determination can lead to similarly uncertain and/or biased estimates of momentum and energy.  \citet{Goodman_2009} find that not only do different types of observations trace different components (e.g. warm/cold gas/dust) in molecular clouds, but they also have sensitivity limits that are different from one another. These limits are dependent on the physical environment, and extreme care should be taken extreme when using (any of) these tracers to estimate the density and other physical properties.

The recent release of data taken by the Herschel Space Observatory has inspired many works examining the validity and caveats of using far-infrared thermal emission in estimates of column density. One of the most discussed issues is the artificial anti-correlation between the temperature and the emissivity spectral index that can be caused by (even small amounts of) noise in real data, when certain fitting procedures are used. The boost of computational power and new algorithms in the past decade has allowed researchers to use parameter-estimation techniques that are less likely to cause false correlations amongst fitted parameters.  \citet{Kelly_2012} have proposed a hierarchical Bayesian-fitting method and tested it on  Herschel data, and they claim that the Bayesian procedure (as opposed to least-squares fitting) reduces (false) temperature-opacity anti-correlation when fitting models for thermal dust emission at multiple wavelenghts.

In this paper we focus on an embedded stellar wind in one of the most near-by molecular clouds: Ophiuchus. A shell structure is found around a group of B-stars ($\rho$ Oph; \S2). Multi-wavelength observations are used in the analysis of the shell. These include near-infrared data from the Two-Micron All-Sky Survey (2MASS), mid/far-infrared maps released by the Improved Reprocessing of the IRAS Survey (IRIS), the Herschel Space Observatory and $^{12}$CO (1-0) and $^{13}$CO (1-0) line data taken at Five College Radio Astronomy Observatory (FCRAO) from the COoridnated Molecular Probe Line Emission Thermal Emission \citep[COMPLETE;]{Ridge_2006} project (\S3). These data sets allow us to reproduce a similar analysis of the limits on using them to trace the column density in molecular clouds as in \citet{Goodman_2009} (\S4.1). In this paper, we also experiment with the hierarchical Bayesian-fitting method \citep{Kelly_2012} on the Herschel maps at PACS 70/100 $\mu$m and SPIRE 250/350/500 $\mu$m (\S4.2). This is the first time the algorithm has been applied on an extended source. Then we conclude with an estimate of the impact caused by the embedded stellar wind (\S5). The result reconfirms the idea proposed by \citet{Arce_2011} that stellar winds driven by B or even later types of stars can contribute to the turbulence in molecular clouds.


\section{Target}
The Ophiuchus cloud is one of the closest molecular clouds, sitting at a distance of $119 \pm 6$ pc \citep{Lombardi_2008}. It hosts a few star forming regions, including the well-studied $\rho$ Oph cluster. Fig. 1 shows the entire Ophiuchus cloud at mid- and far-infrared wavelengths, as well as the $^{12}$CO line emission and H$\alpha$ emission. A shell-like structure can be identified at infrared wavelengths and in the molecular line observations. The H$\alpha$ emission appears in a filled circular shape that fits in the interior of the shell. The center of both the shell and the H$\alpha$ emission coincides with the position of $\rho$ Ophiuchii, which is a group of B-type stars. The shell-like structure, the H$\alpha$ emission and its proximity to B-stars all indicate that the shell is driven by an embedded stellar wind. 

$\rho$ Ophiuchii \citep[16h 25m 35.12s, -23d 26m 49.82s, J2000,]{van_Leeuwen_2007}, not to be confused with the YSO cluster of the same name lying more than a degree to the south, is a neighborhood of five B-type stars \citep{Cordiner_2013}. It includes two pairs of close binaries ($\rho$ Oph AB and $\rho$ Oph DE, at 111$^{+12}_{-10}$ pc and 135$^{+12}_{-10}$ pc respectively) and a single star at the northern edge of the neighborhood ($\rho$ Oph C, at $\sim$ 125$^{+14}_{-11}$ pc). The distances of $\rho$ Oph and the Ophiuchus molecular cloud indicate that at least $\rho$ Oph C is embedded in the cloud, if not all five B-stars considering the uncertainties of the distance measurements.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/all1/all1}
\caption{{The Ophiuchus cloud seen in emission at different wavelengths/from different lines.%
}}
\end{center}
\end{figure}

\section{Data}

\subsection{Near-Infrared Extinction Map from 2MASS}
We use the extinction map presented in \citet{Ridge_2006} and \citet{Alves_2006}, made by applying the Near-Infrared Color Excess Revisited technique (NICER; \citet{Lombardi_2001}) to the Two Mircon All-Sky Survey (2MASS) point source catalog. The resulting K-band extinction is then converted to V-band extinction by the relation $A_K = 0.112 A_V$ (\citet{Rieke_1985}). The derived extinction map has a fixed resolution of 3$'$, with a Nyquist-sampled pixel size of 1.5$'$. The map is made public on the COordinated Molecular Probe Line Extinction Thermal Emission (COMPLETE) project website. 

The 2MASS/NICER extinction map is accurate up to $A_V \simeq$ 10 mag, limited by the number of background stars that can be used to derive the extinction. The confusion of more evolved YSOs embedded in the cloud with the background stars (with their J, H, and K band colors in the same range) limits the accuracy of the extinction on the low end. See \S3.4 for more details on how we deall with these limits and the related data editing.

\subsection{Far-Infrared Dust Thermal Emission}
Mid-infrared and far-infrared emission in molecular clouds are mainly produced by thermal emission from dust. To create spectral-energy distributions (SEDs) suitable for the model-fitting described below,  we use flux measurements from the Improved Reprocessing of the IRAS Survey (IRIS) and data extracted from the Herschel Science Archive.

\subsubsection{IRAS/IRIS}
IRIS flux density maps at 60 and 100 $\mu$m, in contrast with earlier releases of the IRAS all-sky maps, have been corrected for the effects of zodiacal dust and striping. The gain and offset have also been calibrated \citep{Miville_Deschenes_2005}. The IRIS maps of Ophiuchus have a resolution of approximately 2.5', close to the 3' resolution of the 2MASS/NICER extinction map. The noise in the IRIS maps is approximately 0.03 MJy sr$^{-1}$ and 0.06 MJy sr$^{-1}$ at 60 and 100 $\mu$m, respectively. This is $\leq$ 1 \% of the overall flux density in the Ophiuchus region.

\subsubsection{Herschel PACS/SPIRE}
The Herschel Gould Belt Survey (Herschel key program, Andr\'{e} et al 2007) has completed observations of 13 molecular clouds. The data have been made public in the Herschel Science Archive. For the Ophiuchus cloud, the observations were taken in the PACS/SPIRE parallel mode, granting us wavelength coverage over PACS 70/100 $\mu$m and SPIRE 250/350/500 $\mu$m. We processed the archival data in HIPE 10.0.1, using the standard pipeline for extended sources. The calibration of the offset level is done by looking for the faintest region on the maps and subtracting the mean of the lowest flux densities. We check the calibration by comparing the resulting maps to the Planck maps at 857 and 545 GHz. Since the original observation of Ophiuchus is divided into three fields (containing the cluster, the Eastern streamer and the Northern streamer), the three maps are mosaiced together using the Montage package provided by the Infrared Processing and Analysis Center (IPAC).

\subsection{Molecular Line Emission from FCRAO}
We use the $^{12}$CO (1-0) and $^{13}$CO (1-0) line maps of the COMPLETE survey \citep{Ridge_2006} to measure column density. Both lines were observed simultaneously using the Five College Radio Astronomy Observatory (FCRAO) telescope. The resolution is 46" and the pixels in the maps are Nyquist sampled with a pixel size of 23". The main beam efficiency is assumed to be 0.45 and 0.49, respectively for $^{12}$CO and $^{13}$CO. The flux calibration uncertainty is assumed to be 15 \% (\citet{Pineda_2008}). The channel maps of these two molecular lines are also used to assess the kinematics in the cloud. The observation provides a velocity resolution of 0.12 km s$^{-1}$ for both lines.

\subsection{Data Editing}
In this work we analyze two versions of our full dataset, pruned differently. 

The first version includes every pixel that is covered in every data set that we use. No other editing is done for the first version. Since the coverage of the molecular line maps is the smallest among the observations described in the previous section, the constraint is mainly set by the boundary of the molecular line maps.

In the second version, we follow Pineda et al (2008) and edit out the pixels that are less reliable. We first require that the integrated $^{12}$CO (1-0) and $^{13}$CO (1-0) emission is positive, and the main beam temperature at the peak of each line is larger than 10 and 5 times the rms noise of the $^{12}$CO and $^{13}$CO line profile, respectively. We then assume that the $^{12}$CO line is broader kinematically than the $^{13}$CO line and exclude pixels with a line width of $^{12}$CO (1-0) smaller than 0.8 times the linewidth of $^{13}$CO (1-0). Lastly, to address the problem that evolved embedded YSOs can be confused with background stars, we exclude pixels with a stellar density larger than 10. Following the second set of criteria, we obtain a map that has 10825 pixels. This is 4276 pixels fewer than the map we get from the first set of data editing.

To make maps that can be compared to each other, we consistently convolve all maps to a resolution of 3', which is the resolution of the 2MASS/NICER extinction. The convolution is done using a two dimensional Gaussian beam with FWHM = 3'. Then the maps are regridded using the Montage package (by IPAC) to match the 1.5'-pixels on the 2MASS/NICER extinction map. 

\section{Analysis}

\subsection{Tracing Column Density with Different Tracers}
To estimate the turbulent energy of the Ophiuchus cloud, we first estimate its physical properties, including density and temperature. Various tracers have been used in deriving the column density of the molecular gas. For star forming regions like the Ophiuchus, people often use three tracers: near-infrared extinction, mid/far-infrared dust thermal emission, and molecular line emission as proxies for the molecular gas. However, as noted in \citet{Goodman_2009}, since none of these three tracers are the result of physical processes directly related to the dominant component of the molecular gas, which is H$_2$, one has to be extremely careful when choosing a tracer and when comparing measurements from one with another. Following the analysis done in \citet{Goodman_2009}, we examine the distributions of column densities derived from these tracers for the Ophiuchus cloud.

As described in \S 3, we use the 2MASS-based extinction map, the IRIS data, the Herschel data and the FCRAO observation of $^{12}$CO/$^{13}$CO line emission in our analysis. Fig. 2 shows the column density maps based on these tracers, given in units of equivalent V-band magnitude.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/av2/av2}
\caption{{The resulting column density map in units of equivalent V-band extinction (mag), using algorithms described in \S3.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/temp4/temp4}
\caption{{The temperature map based on Herschel using $\chi^2$-based SED fitting, and the same map based on IRIS 60/100 $\mu$m.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/av-scatter1/av-scatter1}
\caption{{Scatter plots of column density estimates based on different tracers used in this paper. The dashed line in the lower left plot is the cutoff (lowest) value of column density, in units of equivalent V-band extinction, which is measurable with the molecular line emission. The red solid line in each plot shows shows the 1-1 relation.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.98\columnwidth]{figures/hist/hist}
\caption{{Distribution of column density estimates in units of equivalent V-band extinction, based on tracers used in this paper. The blue curve in each plot is the normal profile fitted to the distribution of 2MASS/NICER extinction, shown in the first plot. Curves in subsequent plots, coded in red and green, are the normal profiles fitted to distributions shown in these plots. The standard deviation $\sigma$ for the normal profiles are represented by the horizontal bars in the first three plots. In the last plot, the green histogram shows the distribution of the column density estimate based on N($^{13}$CO), and the dark gray filled histogram shows the distribution based on the integrated $^{13}$CO (1-0) line emission. See \S3.3 for details about the fitting methods.%
}}
\end{center}
\end{figure}

\subsubsection{2MASS Extinction}
As discussed in \citet{Goodman_2009}, although it is expected that in a region as large and complicated as the Ophiuchus the gas-to-dust ratio can vary by a significant amount, the 2MASS/NICER-based extinction measurement might still be the most reliable tracer of column density amongst all the tracers presented in this paper. However, we have to treat it with extreme care, especially towards the very diffuse and very dense regions. In both cases, the uncertainty comes from the number of background stars that can be used to derive the extinction. This being said, the 2MASS extinction is used as the ``true'' extinction measurement and helps us determine the conversion from various types of density measurement (for example, the optical depth and the mass column density) to the equivalent V-band extinction.

Fig. 5 shows the distribution of the 2MASS extinction. It has a clear log-normal-like profile. \citet{Padoan_1997} point out that this log-normal-like profile is equivalent to the $\sigma(A_V)-A_V$ relation first discovered by \citet{Lada_1994}. Thoraval et al. (1997) find that this is the result of small structures below the resolution limit. 

\subsubsection{IRIS Two-point Fit of Dust Thermal Emission}
The basic method we use to derive the dust temperature and column density from the IRIS 60 and 100 $\mu$m flux densities follows \citet{Arce_1999} and \citet{Schnee_2005}. The dust temperature is determined by the ratio of flux densities at 60 and 100 $\mu$m. We then calculate the column density by estimating the opacity along each line of sight. This calculation depends on the values of two parameters: the emissivity spectral index and the conversion factor between the optical depth at 100 $\mu$m and the extinction.

The modified blackbody equation is used to model the dust thermal emission at mid/far-infrared wavelengths. It assumes that each line of sight has a uniform temperature and a uniform volume density, and that the material along the line of sight is in local thermal equilibrium (LTE). After accounting for radiative transfer along the line of sight, the modified blackbody equation can be written as the blackbody equation modified by a power-law emissivity spectral index. At the dust temperature $T_d$, we have the modified blackbody equation:

\begin{equation}
F_\lambda = B_\lambda(T_d)\,(1-\exp{(-\tau)})
\end{equation}

The optical depth $\tau$ can be written in terms of a power-law emissivity, which depends on the wavelength: $\tau = N_d\,\alpha\,\lambda^{-\beta}$ (\citet{Schnee_2005}), where $N_d$ represents the column density of dust grains, $\alpha$ is the constant that relates the flux to the optical depth of the dust, and $\beta$ is the emissivity spectral index. In the optically thin limit, we can calculate the flux from:

\begin{equation}
F_\lambda = B_\lambda(T_d)\,N_d\,\alpha\,\lambda^{-\beta}
\end{equation}

Here we have measurements in IRAS band 3 (60 $\mu$m) and band 4 (100 $\mu$m). To derive the column density, we first calculate the dust temperature from the modified blackbody equation. We can write the ratio of the flux densities at 60 $\mu$m and 100 $\mu$m as:

\begin{equation}
R = \frac{F_{60}}{F_{100}} = 0.6^{-(3+\beta)}\,\frac{\exp{(hc/\lambda_{60}kT)}-1}{\exp{(hc/\lambda_{100}kT)}-1}
\end{equation}

Once a reasonable value of $\beta$ is given, we can calculate the dust temperature. From there, we calculate the optical depth at 100 $\mu$m, assuming again that the lines of sight are optically thin.

\begin{equation}
\tau_{100} = \frac{F_{100}}{B_{100}(T_d)}
\end{equation}

We then convert the 100 $\mu$m optical depth to the V-band extinction using

\begin{equation}
A_V = X\,\tau_{100}
\end{equation}

The V-band extinction provides an estimate of the column density. Here we follow the discussion in \citet{Goodman_2009} and use the 2MASS extinction map as the ``true'' extinction measurement (with limits discussed in \S 4.1.1) to fit for a conversion factor $X$. The equivalent V-band extinction can be further converted to the column density.

Previous works have shown that the dust thermal emission at mid/far-infrared wavelengths is dominantly contributed by Big Grains (BGs; \citet{Stepnik_2003}). However, at wavelengths less than 100 $\mu$m, the observed flux density contains a significant amount of emission from Very Small Grains (VSGs), which we can no longer describe with a simple modified blackbody equation. \citet{Schnee_2008} examine the contribution of VSGs and finds that not only does a large amount of emission at mid-infrared wavelengths come from the VSGs, but also the contribution of VSGs varies systematically with $T_d$ and $N_d$. In the Perseus cloud, the contribution of VSGs at 60 $\mu$m ranges from 77 to 83 \%. We adopt a modification factor of 0.2 to account for the ``true'' thermal emission from the BGs at 60 $\mu$m, which we can model with the modified blackbody equation.

In addition to the uncertainty in the VSG contribution of emission at 60 $\mu$m, there is another uncertainty caused by the emissivity spectral index $\beta$. First we note that the modified blackbody equation assumes a single temperature and volume density for each line of sight. \citet{Shetty_2009} point out that the variation in temperature along the line of sight can affect the shape of the observed SED, even when each dust grain is in local thermal equilibrium and emits according to the blackbody equation. The result is that the modified blackbody equation can no longer recover the shape of SEDs, and that the emissivity spectral index in the empirical power law varies from one line of sight to another. \citet{Kelly_2012} further address that making $\beta$ a free parameter in the calculation does not solve this problem. Any errors in the data from observations and/or calibration will cause the derived $T_d$ and $\beta$ from the modified blackbody equation to show an artificial anti-correlation. This is because some degree of redundancy exists between increasing $T_d$ and decreasing $\beta$ (and vice versa) when we cannot measure the exact SED without any uncertainty. This false anti-correlation has led many works to assume the anti-correlation \emph{a priori} in their derivation (including \citet{Schnee_2005}, but not \citet{Schnee_2008}). For us to derive the dust temperature and column density from the IRIS 60 and 100 $\mu$m, the same trend between $T_d$ and $\beta$ forces us to use a single value for $\beta$ throughout the entire region, due to the lack of data points to fit $T_d$, $N_d$, $\beta$, and the conversion factor $X$ simultaneously. We assume $\beta = 2$ empirically from observations of cold regions (T \< 35 K; \citet{Dupac_2001}) and leave the discussion of temperature variation and the effects of noise to \S 4.2, where we enjoy the wealth of data from Herschel observations.

Fig. 2 and 3 show the equivalent V-band extinction and temperature from the calculation described above. Fig. 4 shows the distribution of the IRIS-based V-band extinction against that based on 2MASS/NICER after the data editing described in \S 3.4. We see from Fig. 4 that there are sub-populations within this distribution. Linking the data with the pixel positions on the sky makes it clear that these sub-populations correspond to different parts of the cloud. Fig. 6 shows this linked data view using the entire dataset without any editing. The dominance of data points along the 1-1 line is from the main cloud, including most of L1729, L1712 to the east and L1688 around the cluster. Points sitting slightly above the 1-1 line are mostly from around the ``shell'' region. Another population that corresponds to the ``northern streamer'' (L1755 and L1765) can be identified when we include the entire data set. The existence of ``sub-populations'' on the plot of $A_V$-$A_V$ hints at the fact that the conversion factor $X$ is different for different regions. Similar results have been reported in \citet{Pineda_2008} for the CO-based extinction. This result was not seen in \citet{Schnee_2005} due to their cropping and smoothing. We adopt the smoothing of a 3' beam, instead of the 5' beam used in \citet{Schnee_2005} and include the very east end of the Ophiuchus cloud. Both the smoothing and cropping adopted in this paper increase the possibility of including more pixels with divergent physical properties.

Fig. 5 presents the distribution of the IRIS-based extinction, which seems askew from the log-normal profile of the 2MASS extinction. This is the result of uncertainties in both the VSG emission and the emissivity spectral index $\beta$. The uncertainty in the amount of VSG emission at 60 $\mu$m affects the determination of $T_d$ using the ratio of flux densities at 60 and 100 $\mu$m. On the other hand, fixing $\beta$ = 2 affects both $T_d$ and $\tau_{100}$, and thus the column density. It is difficult to disentangle these two effects, which again were not seen in \citet{Schnee_2005} (and subsequently \citet{Goodman_2009}) due to the cropping, the smoothing and their \emph{a priori} assumption about the anti-correlation bewteen $\beta$ and $T_d$.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/glue-IRIS/glue-IRIS}
\caption{{{\bf Left} The scatter plot of column density estimates in units of equivalent V-band extinction, based on the 2MASS/NICER extinction and the IRIS 60/100 $\mu$m data. {\bf Right} Regions corresponding to the three ``sub-populations'' identified on the plot of A$_V$-A$_V$. The regions are overlaid on the 2MASS extinction map in gray scale. The plots are made with Glue (Beaumont et al 2012; see \url{http://glueviz.org} for information about Glue).%
}}
\end{center}
\end{figure}

\subsubsection{Herschel SED Fit of Dust Thermal Emission}
The Herschel data for Ophiuchus covers two of the PACS bands (70 and 160 $\mu$m) and all three of the SPIRE bands (250, 350, and 500 $\mu$m). This wealth of data enables us to perform SED fitting to the modified blackbody equation. We follow \citet{Arzoumanian_2011} and use the following form of the modified blackbody equation, for the convenience of comparison between this work and previous works done by them and other Herschel teams.

\begin{equation}
F_\lambda = B_\lambda(T_d)\,0.1\,\left(\frac{\nu}{\text{1000 GHz}}\right)^\beta\,\Sigma_d
\end{equation}

With this equation, we adopt the emissivity anchor in \citet{Arzoumanian_2011} at 1000 GHz and derive the mass column density $\Sigma_d$ instead of the number column density. These do not change the conclusions from the discussion in the previous section (\S 4.1.2). We deal with the uncertainty due to the VSG emission by abandoning the 70 $\mu$m data, which also suffers from the large uncertainty of its own calibration, especially the background subtraction. Here we also try to set $\beta$ as a free parameter in the fit but find that freeing $\beta$ makes uniquely determing $T_d$ and $N_d$ difficult (see Fig. 7 for the temperature and column density maps). This is again due to the anti-correlation discussed in the previous section, which exists whenever there is any amount of uncertainty in the data. Since our map of the Ophiuchus region is mosaiced from three different observations, it is difficult to track the effect of uncertainties. We see at least three ``sub-populations'' of data points showing up on the plot of $A_V$-$A_V$ against 2MASS extinction (Fig. 8), when $\beta$ is allowed to vary in the fitting. This is a natural result of the false anti-correlation between $\beta$ and $T_d$ and the variation of physical properties across different regions. To solve the complexity caused by a free $\beta$ and make the result comparable to that based on other tracers, we fix $\beta$ to 2 (same as the value used to fit the IRIS data) in our $\chi^2$-based fitting. We also seek alternative fitting schemes, including the hierarchical Bayesian method (Kelly et al 2012), the result of which is presented in \S 4.2.

Fig. 2 and Fig. 3 show the result of the fitting where $\beta = 2$. Fig. 4 presents a similar scatter plot of the Herschel-based extinction against the 2MASS extinction. Compared to Fig. 8, it is clear that fixing $\beta$ causes the trend to converge well to the 2MASS extinction, despite the different physical properties across the cloud. The effect of a fixed $\beta$ is also obvious when comparing the distribution of Herschel-based to then IRIS-based extinction (Fig. 5, compared to Fig. 9). The distribution of Herschel-based extinction is less askew and agrees better with the 2MASS extinction.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/bad-her1/bad-her1}
\caption{{{\bf Left} The column density map in units of equivalent V-band extinction, based on the Herschel data using a $\chi^2$-based fitting algorithm with a varying emissivity spectral index $\beta$. {\bf Right} The temperature map using the same algorithm.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/Glue-her/Glue-her}
\caption{{{\bf Left} The scatter plot of column density estimates in units of equivalent V-band extinction, based on the 2MASS/NICER extinction and the Herschel data, with a varying spectral index $\beta$. {\bf Right} Regions corresponding to the three ``sub-populations'' identified on the plot of A$_V$-A$_V$. The regions are overlaid on the 2MASS extinction map in gray scale. The plots are made with Glue (Beaumont et al 2012; see \url{http://glueviz.org} for information about Glue).%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.98\columnwidth]{figures/hist-varyB/hist-varyB}
\caption{{The distribution of column density in units of equivalent extinction, based on the Herschel data using a $\chi^2$-based algorithm where the emissivity spectral index $\beta$ is set as a free parameter.%
}}
\end{center}
\end{figure}

\subsubsection{FCRAO Molecular Line Emission}
We follow the analysis in \citet{Pineda_2008} to derive the equivalent V-band extinction from the $^{12}$CO (1-0) and $^{13}$CO (1-0) line observations. We start by assuming that the $^{12}CO$ is optically thick. We then calculate the excitation temperature using the main beam temperature at the peak of $^{12}CO$:

\begin{equation}
T_{\text{ex}} = \frac{5.5\,\text{[K]}}{\ln{\left(1+5.5\,\text{[K]}/\left(T_{\text{max}}(^{12}\text{CO})+0.82\,\text{[K]}\right)\right)}}\;\;\;\text{,}
\end{equation}

where $T_{\text{max}}(^{12}\text{CO})$ is the brightness temperature at the peak. If the excitation temperature of the $^{13}CO$ (1-0) line is the same as the $^{12}CO$ (1-0), we can calculate the optical depth of the $^{13}CO$ (1-0) line:

\begin{equation}
\tau(^{13}\text{CO}) = - \ln\left[1 - \frac{T_{\text{max}}(^{13}\text{CO})/5.3\,\text{[K]}}{1/\left(e^{5.3\,\text{[K]}/T_{\text{ex}}}-1\right) - 0.16}\right]
\end{equation}

Similarly, here $T_{\text{max}}(^{12}\text{CO})$ is the main beam brightness temperature at the peak of $^{13}CO$. Using the definition of column density in Rohlfs \& Wilson (1996), we can derive the $^{13}CO$ column density from the two equations above:

\begin{equation}
N(^{13}\text{CO}) = 3.0\times10^{14}\;\left[\frac{\tau(^{13}\text{CO})}{1-e^{-\tau(^{13}\text{CO})}}\right]\,\frac{W(^{13}\text{CO})}{1-e^{-5.3\,\text{[K]}/T_{\text{ex}}}}\;\text{[cm}^{-2}\text{]}\;\;\;\text{,}
\end{equation}

where $W(^{13}\text{CO})$ is the integrated intensity along the line of sight in units of [K km/s]. Similarly, we then fit the conversion from both $N(^{13}\text{CO})$ and $W(^{13}\text{CO})$ to the 2MASS extinction.

\begin{equation}
\begin{aligned}
A_V &= X_{N(^{13}\text{CO})}\,N(^{13}\text{CO}) \\
A_V &= X_{W(^{13}\text{CO})}\,W(^{13}\text{CO})
\end{aligned}
\end{equation}

Fig. 4 shows the scatter of the derived extinction based on $N(^{13}\text{CO})$. The dashed line shows the cutoff (lowest) value of column density measurable with $^{13}CO$. This cutoff value is 4.3 mag for the Ophiuchus cloud, compared to 1.67 mag for the Perseus. It is clear that no ``sub-population'' is seen here in the molecular line-based extinction. This suggests that the variation of the conversion factors inherent to the ``sub-populations'' seen in the thermal emission-based extinctions is the result of different dust-to-gas ratios in different parts of the cloud. This is consistent with the limit of using dust tracers for the column density, as presented in \citet{Goodman_2009}. For a region as complicated as the Ophiuchus (where multiple components with very different physical properties present), the variation in the dust-to-gas ratio affects the conversion from dust to the total column density.

Fig. 5 shows the distribution of the $CO$-based extinction. Although the profile does not seem as log-normal as other dust tracers, the mean and the variance agree with the 2MASS extinction measurements.

\subsection{Deriving Physical Properties from Herschel Data}
The Herschel data provides wavelength coverage from 70 to 500 $\mu$m. This allows us to experiment with different schemes to derive the temperature and column density. As mentioned in the previous section (\S 4.1), the modified blackbody equation widely used to derive physical properties from mid/far-infrared observations assumes a single temperature and volume density along each line of sight but treats each pixel independently from others. This assumption does not agree with what actually happens in molecular clouds and systematically affects the results of the fitting. Both the variation of physical properties along the line of sight (especially the variation of temperature) and the noise in the data can seriously bias the distribution of the derived parameters. Recently, multiple authors (Loredo 2004; Hogg et al 2010) have proposed various Bayesian fitting schemes to solve this problem. Here we examine the effects of varying temperature and column density along the line of sight with a simple model, and then experiment with the hierarchical Bayesian-fitting method(Kelly et al 2012) on our Herschel data.

First we examine the effects of varying temperature and column density along the line of sight by building a ``toy'' column of 2 pc. We divide this column into 100 slabs, each with a thickness of 0.02 pc along the length of the toy column. Each of these slabs has its own temperature, volume density and the emissivity spectral index $\beta$. Temperatures and \emph{volume} densities in these slabs are drawn independently from two uncorrelated log-normal distributions to create an artificial anti-correlation distribution in temperature-density space (Fig. 10). This roughly simulates the empirical distribution of the temperature and column density in the interstellar medium. The emissivity spectral index $\beta$ is drawn again independently from a normal distribution. These parameters are then randomly placed in the slabs, creating our ``toy'' line of sight. Fig. 10 also shows the temperature, column density and $\beta$ along this line of sight. We then place an observer in front of the toy column and synthesize an observation by calculating the radiative transfer, assuming that each slab is in local thermal equilibrium and emits according to the modified blackbody equation (that is, the assumptions for the modified blackbody equation are still true within the slab but not across the entire line of sight). We sample the resulting SED at Herschel wavelengths (at PACS 70/100/160 $\mu$m and SPIRE 250/350/500 $\mu$m) and fit a modified blackbody equation to these data points. The result is presented in Fig. 11, showing that the fitted temperature and column density are nowhere near the mean or the mode of the simulated parameters. One can also see in Fig. 11 that fitting with data at wavelengths further away from the peak gives us numbers closer to the mean values, although the mean values are no longer representative of the physical environment in a column of materials when there are multiple components along the line of sight. This result shows that with varying parameters along the line of sight, the resulting SED can not be fully described by the modified blackbody equation. This agrees with \citet{Shetty_2009}, which also suggests that fitting with data at wavelengths further away from the peak of the SED is more likely to give the ``true'' values of these parameters. Notice that we have not added any noise yet, and the result purely comes from the variation of the temperature, the column density and the emissivity spectral index. 

To deal with the effects of noise in the data, we experiment with the hierarchical Bayesian-fitting method, proposed in \citet{Kelly_2012}, on our Herschel data. In this Bayesian-based fitting scheme, Kelly et al. try to deal with 1) how the SED parameters for individual pixels are generated from distributions of these parameters and 2) how the measured fluxes are generated from these SED parameters for individual pixels, at the same time. They model the distributions of $\log{(N_d)}$, $\log{(T_d)}$ and $\beta$ with a multivariate Student-T distribution and then model the distributions of the logarithmic multiplicative errors (representing calibration uncertainties) and the additive errors (representing observational errors) with another two Student-T distributions. This is basically equivalent to modeling each of these parameters with a single population in parameter space and in the distribution of observed fluxes. This fitting scheme is called hierarchical and differs from other traditional Bayesian fitting schemes because it deals with the distribution in parameter space and the real observation at the same time. In \citet{Kelly_2012}, the hierarchical Bayesian fitting method is tested with a simulated cube and a real Herschel observation of a protostellar core. In both cases it solves the biased T-$\beta$ anti-correlation and successfully models the effects of noise in the data.

Using the code proposed by \citet{Kelly_2012}, we fit the temperature, the column density and $\beta$ for the Ophiuchus cloud. Although the code has been successfully implemented in a few single source cases, this is the first time that the hierarchical Bayesian-fitting method has been applied to a very extended source. Fig. 12 shows the distribution of the resulting temperature and $\beta$ for a sample region with 128 $\times$ 128 pixels, compared to that of traditional $\chi^2$-based fitting. We see that the false anti-correlation caused by the existence of the noise is resolved, and the result of the hierarchical Bayesian method shows a much milder anti-correlation.

Since running the code is computationally expensive and the code is meant for a single population of fluxes, we propose to run the code in parallel on smaller overlapping regions (``boxes'') and mosaic the result back together. We experiment with this procedure by dividing the region close to the center of the shell into 42 65 $\times$ 65 pixels boxes. The mosaicing is then done using the Montage package provided by the Infrared Processing and Analysis Center. Fig. 13 shows the mosaiced map of the temperature and column density. The same T-$\beta$ plot is also presented, and it shows that the mosaicing spreads out the relatively constrained T-$\beta$ correlations in each box over a larger extent on the T-$\beta$ plot. However, compared to the $\chi^2$-based result (Fig. 12; also see \S4.1.3), the result of the hierarchical Bayesian fitting is still more constrained (due to the assumption of the distributions) and has a relatively milder anti-correlation.

We conclude that the hierarchical Bayesian-fitting method solves the problem of bias due to the noise, but the assumption that $\log{(T)}$, $\log{(N)}$ and $\beta$ are from one multivariate distribution does not apply to an extended and complicated region like the Ophiuchus. The procedure of using smaller boxes makes it perform better within each box, but the mosaicing process erases the original design of the fitting method. The result is no longer modeled by a single multivariate Student-T distribution after the mosaicing process tries to adjust the parameters in each box to fit a global plane. Embedded sources that are often seen in molecular clouds affect the fitting result by dragging the assumed distribution as well. Notice also that the hierarchical Bayesian-fitting method (or any fitting scheme) does not/cannot solve the problem of variations along the line of sight without knowledge about the distribution along that direction. In the future, a more careful modeling of the distributions of the physical parameters in molecular clouds will have to be done before we can apply the hierarchical Bayesian-fitting scheme to extended sources.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/toy/toy}
\caption{{{\bf Upper} The 2 dimensional histogram showing the probability distribution of the distributions that we draw the temperatures and volume densities from. The two histograms show distributions of temperatures and volume densities that are drawn to build up the toy line of sight. {\bf Lower} The toy line of sight. The observer are placed at the left-hand side facing right, thus observing the entire column of material.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/SEDs/SEDs}
\caption{{The synthetic SED and the SEDs calculated from the fitted parameters. The synthetic SED is represented by the blue solid line, with blue dots showing fluxes sampled at the Herschel wavelengths. The red solid lines show the SEDs calculated based on parameters when fitting to fluxes at all five wavelengths (in lightest red), the longest four wavelengths and the longest three wavelengths (in darkest red). The fitted temperature and the spectral index $\beta$ are written in corresponding colors. Mean values of the temperatures and $\beta$s in the toy line of sight are written in blue.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Bay-exp/Bay-exp}
\caption{{{\bf Left} The plot of temperature-spectral index $\beta$ showing the results from $\chi^2$-based fitting and from hierarchical Bayesian-fitting, for the 128$\times$128 box. The red dots are randomly drawn from the Monte Carlo chains, with each chain corresponding to one pixel on the map. The distributions of the Monte Carlo chains provide a natural estimate of the uncertainty of the fit. {\bf Right} The same plot for the 42 65$\times$65 boxes, after mosaiced together. Notice that in this plot, since we use the mean values of the Monte Carlo chains to create the maps, the red dots here are showing mean values shifted by the mosaicing process. This is why the scatter in the data seems smaller.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/bay-chi2/bay-chi2}
\caption{{The column density maps and the temperature maps using {\bf Left)} the hierarchical Bayesian-fitting method and {\bf Right)} $\chi^2$-based fitting.%
}}
\end{center}
\end{figure}

\section{Results}

\subsection{Energy/Momentum Estimation}
To estimate the contribution of the embedded B-star wind to the turbulence in the Ophiuchus, we adopt the near-infrared extinction for the mass calculation (reasons discussed above in \S4.1) and the molecular line emission for the velocity calculation. The molecular line emission provides the resolution as well as the dispersion along the velocity axis. Thus, we base our calculation on the 2MASS extinction map and the FCRAO $^{12}$CO/$^{13}$CO maps. In this section, no data editing is applied, so the estimates of mass, energy and momentum in the molecular cloud are the upper limits. The data editing does not affect the calculation for the shell significantly.

\subsubsection{The Cloud}
The mass within the cloud and the shell is estimated by calculating the mass within projected regions using the column density map. As discussed in \S4.1, we adopt the 2MASS extinction map in our mass calculation. We then calculate the mass column density from the equation (\citet{Arce_2001}):

\begin{equation}
\Sigma = \mu\,N_{H_2} = \mu\,(9.4\times10^{20})A_V\;\text{g/cm}^2
\end{equation}

$\mu$ is the average mass per particle and is $\sim$ 2.72 times mass of a hydrogen atom (a.m.u.) for molecular clouds in our Milky Way (ref?). We then obtain the total mass from $M = \Sigma\times \text{Area}$.

In the case of the Ophiuchus, we first calculate the threshold extinction value used to determine the edge of the cloud by calculating the average extinction in places with no obvious far-infrared emission. This threshold value is $\sim$ 0.93 mag. ($\sim$ $8.7 \times 10^{20}$ cm$^{-2}$) on our 2MASS extinction map. We then select the projected region of our cloud to include places where the extinction is larger than 0.93 mag. This threshold value is then subtracted from the extinction within the selected region. The total mass of the Ophiuchus cloud is estimated to be $\sim$ 6000 M$_{\sun}$.

The turbulent velocity of the cloud is estimated using the $^{12}$CO (1-0) line emission. Since the line profile in a molecular cloud at this scale varies from one part of the cloud to another, it is unsuitable to fit the profile and to calculate the velocity dispersion from the fit. A natural way to calculate the velocity dispersion is the second central moment:

\begin{equation}
\sigma_v = \sqrt{\left(\frac{\int{I(\mathbf{p}, v)\,(v(\mathbf{p}, v)-v_{cent}(\mathbf{p}))^2\,\text{d}(\mathbf{p}, v)}}{\int{I(\mathbf{p}, v)}\,\text{d}(\mathbf{p}, v)}\right)}\;\;\;\text{,}
\end{equation}

where $I(\mathbf{p}, v)$ and $v(\mathbf{p}, v)$ are the intensity and the velocity at each position $\mathbf{p}$ and velocity $v$, respectively. $v_{cent}(\mathbf{p})$ is the velocity centroid and is calculated from:

\begin{equation}
v_{cent} = \frac{\int{I(\mathbf{p}, v)\,v(\mathbf{p}, v)\,\text{d}(\mathbf{p}, v)}}{\int{I(\mathbf{p}, v)\,\text{d}(\mathbf{p}, v)}}
\end{equation}

This is also the first moment, which can be understood as the mean velocity weighted by the intensity at each velocity. In the data cube obtained at the FCRAO, the velocity resolution is 0.21 km/s, smaller than the typical thermal linewidth in molecular clouds ($\sim$ 1 km/s). Thus, we can use the discrete sum to calculate the integrals above without losing any information.

The average of the second central moment in the Ophiuchus cloud is $\sim$ 0.97 km/s. The rms velocity of the turbulence is $\sqrt{3}$ times the velocity dispersion and is thus $\sim$ 1.68 km/s on average over the cloud. To calculate the total momentum and energy in the turbulence, we match the maps of the second moment and the 2MASS/NICER extinction. We then calculate the momentum/energy within each pixel by using the extinction and the second moment at the pixel. For regions outside the coverage of the FCRAO observations, we adopt the mean second moment (0.97 km/s) for the turbulent velocity. Summing up the energy $E_i = 0.5\,M_i\,v_{turb,\,i}^2$ over the entire cloud gives us an estimate of the total momentum and energy in the turbulence. These are 9900 M$_{\sun}$ km/s and $8.82\times10^{46}$ erg, respectively.


\subsubsection{The Shell}
The ``bubble'' model (Silk 1985) of an embedded stellar wind predicts that we should observe a hot and relatively diffuse interior and a warm dense shell. In the case of the embedded stellar wind, the ``thickness'' of the shell is often much smaller than the radius of the shell (\citet{Churchwell_2007}; \citet{Arce_2011}), thus making the opacity of the shell smaller along lines of sight through the center of the shell. As a result, the projection of the shell on the sky often looks like a ``ring.'' Due to the variation of the density and the irregularity of the shape of the cloud, the ring is likely incomplete and asymmetric. Estimating mass from the column density within a projected area is thus a lower bound, and serves as a good proxy for the total mass of the shell when the shell is optically thin and symmetric. To estimate the mass of the shell, we first overlay the temperature and density maps to determine an elliptical region with a finite width in the warm and dense part of the map. Notice that we completely exclude the cluster region to avoid confusion between the gas in the shell and within the cluster (Fig. 14, on the 2MASS extinction map). After this projected region of the shell is determined, we use the 2MASS extinction map and convert the extinction magnitude into the mass column density using Eq. 11. Summing up the mass column density within the shell mass region then gives us the shell mass, estimated to be 462 $M_{\sun}$. This is a lower bound, since the cluster region we exclude when choosing the shell mass region likely contains gas within the shell.

The expanding velocity of the shell is determined from the average spectrum of the shell. The symmetry of the shell predicts that the velocity we observe with the molecular line emission within the ``ring'' is symmetric around some system velocity. That is, for each ``particle'' in the shell, the velocity is determined by ($v_{system}$ + $v_{expansion}\,\sin{\left(\cos^{-1}{\left(R/R_0\right)}\right)}$ + $v_{intrinsic}$), where $v_{intrinsic}$ is drawn from a velocity distribution characteristic of the molecule's intrinsic (microscopic) motion. Assuming that the line emission from $^{12}$CO (1-0) traces the gas component of the shell and that the intrinsic velocity distribution can be approximated by a normal distribution, we find $v_{expansion}$ from the residual of fitting the $^{12}$CO line with $N(v_{system}, \sigma)$. Fig. 15 shows the residual as well as the fitted Gaussian line profile, with an asymmetry of the residual due the opacity. $v_{residual}$ is $\sim$ 1.2 km/s, which is $v_{expansion}\,\sin{\left(\cos^{-1}{\left(R/R_0\right)}\right)}$ in the equation above. With the shell mass region selected to have an inner radius $\sim$ 0.98 pc and an outer radius $\sim$ 1.61 pc, we find $v_{expansion}$ to be $\sim$ 1.52 km/s. Notice that $v_{residual}$ is a lower bound since the residual calculated from the $^{12}$CO (1-0) line profile is likely probing the gas component at a larger radius than the inner radius due to the opacity effect.

From these estimates, we obtain the momentum and the kinetic energy of the shell. These are given by $P_{shell} = M_{shell}\,v_{expansion}$ and $E_{shell} = 0.5\,M_{shell}\,v^2_{expansion}$, respectively. For this shell, the momentum is $\sim$ 702 $M_{\sun}\,km/s$, and the energy is $\sim$ $1.06\times10^{46}\,erg$. This is $\sim$ 12 \% of the total turbulent energy of the cloud (compared to $\sim$ 15 \% for the largest shell in Perseus; Fig. 16). Notice again that since the estimates for the mass and the expansion velocity are both the lower limits of possible values, the momentum and energy of the shell presented here are likely smaller than the real values as well.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/shell-2mass1/shell-2mass1}
\caption{{The region selected as the ``shell'' in the estimation of the shell mass and energy, overlaid on the 2MASS/NICER extinction map.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/v-exp/v-exp}
\caption{{{\bf Left} The average $^{12}$CO line profile in green and the Gaussian fit to the center of the line in orange. The velocities between the dashed lines are the velocities (and their corresponding flux densities) used in the Gaussian fit. {\bf Right} The residual of the $^{12}CO$ line over the Gaussian fit in green, and the same fit in orange. The dashed lines show peaks of the residual (green) and the fit (orange).%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/bar-energy/bar-energy}
\caption{{Energy in clouds and in shells, in units of 10$^{46}$ erg.%
}}
\end{center}
\end{figure}

\subsection{Impact of the Shell}
To estimate the impact of the shell on the cloud, we first compare the energy in the shell with the gravitational binding energy. With an effective radius of 5 pc (from the approximated geometrical mean of the cloud extent), we calculate the gravitational binding energy to be $3.97\times10^{47}$ erg, which is 30 times larger than the kinetic energy in the shell. Clearly, the shell does not have the energy to unbind the entire cloud.

One way to assess whether or not the embedded stellar wind is powerful enough to drive the turbulence is by comparing the energy injection rate ($\dot E_W$) of the shell to the turbulence dissipation rate ($L_{turb}$). The intuitive way to estimate the energy injection rate of the shell is to assume a timescale during which the wind from the embedded B-stars has been accumulating the kinetic energy that we observe within the shell. Following Arce et al. (2011), we assume $\tau_W$ $\sim$ $1\times10^6$ yr. We then have $\dot E_W = E_{shell}/\tau_W \sim 3.4\times10^{33}$ erg/s. Another way to assess this is following Eq. 3.7 in \citet{McKee_1989}:

\begin{equation}
\dot E_W = 0.5\,(\dot M_W\,v_W)\,v_{rms}\;\;\;\text{,}
\end{equation}

where $v_{rms}$ is the rms velocity of the turbulence in the cloud. $\dot M_W$ and $v_W$ are the mass loss rate and the wind velocity respectively. The mass loss rate of the shell can be estimated from the shell momentum we calculated in \S5.1 by assuming $\tau_W$ and a wind velocity $v_W$.

\begin{equation}
\dot M_W = \frac{P_{shell}}{\tau_W\,v_W}
\end{equation}

Assuming $\tau_W$ $\sim$ $6\times10^6$ yr and $v_W$ $\sim$ 200 km/s, we get a mass loss rate $\dot M_W$ $\sim$ $2.3\times10^{-5}$ M$_{\sun}$/yr. With the rms velocity of the turbulence $\sim$ 1.68 km/s, we get an energy injection rate of $\sim$ $10^{33}$ erg/s.

The turbulence dissipation rate depends on the dissipation time, which is related to the free-fall time and is often written as $t_{diss} = \eta\,t_{ff}$. In \citet{McKee_1989} and McLow (1999), numerical simulations give $\eta$ in the range of 1 to 10. Assuming $\eta$ = 5 and a volume density of $10^3$ cm$^{-3}$, we have $t_{diss}$ $\sim$ $5\times10^6$ yr. We can then calculate the turbulence dissipation rate:

\begin{equation}
L_{turb} = \frac{E_{turb}}{t_{diss}}
\end{equation}

This gives $L_{turb}$ $\sim$ $10^{33}$ erg/s, which is on the same order of magnitude as the energy injection rate of the shell. These estimates indicate that the spherical wind from B-stars ($\rho$ Oph) has the potential for driving the turbulence in the Ophiuchus cloud. This result agrees with \citet{Arce_2011}.

\subsection{Comparing to Perseus}
Table 1 shows the physical properties of the shell discussed here and the largest shell (CPS 5) in the Perseus (\citet{Arce_2011}). We find that while the momenta of the two shells are comparable, the shell in the Ophiuchus cloud is more massive with a lower expansion velocity. This results in its smaller kinetic energy. However, since the total energy in the turbulence is smaller in the Ophiuchus cloud, the ratios of the shell energy to the turbulence energy in the cloud are comparable to each other. In the comparison of the energy injection rate and the turbulence dissipation rate, the shell in the Ophiuchus is powerful enough to drive the turbulence by itself, while in the Perseus, all shells with confidence scores larger than 4 (\citet{Arce_2011}) are required to get an energy injection rate at the same order of magnitude as the turbulence dissipation rate. We think that this hints at the fundamental difference between the two clouds and the shell(s) therein. The Ophiuchus cloud is slightly smaller than the Perseus cloud in terms of its mass and size, and the larger mass/smaller velocity of the shell in the Ophiuchus cloud probably indicates that the shell has been expanding for a longer period of time.\selectlanguage{english}
\begin{table}
\begin{tabular}{lcc}
\hline
& \textbf{Shell in Ophiuchus} & \textbf{CPS5 in Perseus} \\   
Radius (IR) [pc] & 1.30 & 2.65 \\ 
Mass [M_{\sun}] & 462 & 53 \\ 
Expansion Velocity [km/s] & 1.2 & 3 \\ 
Momentum [M_{\sun} km/s]& 702 & 315 \\ 
Kinetic Energy [erg]& 1.06$\times10^{46}$ & 1.88$\times10^{46}$ \\ 
... of total turbulent energy & (12.0\%) & (11.8\%) \\ 
Energy Injection Rate [erg/s] & 2.46$\times10^{33}$ & 1.49$\times10^{32}$ \\
\hline
\end{tabular}
\caption{{Comparison between the shell in Ophiuchus and the shell with the largest radius, CPS5, in Perseus (Arce et al 2011)}}
\end{table}

\section{Conclusions}
In this work, multi-wavelength observations enable us to analyse the derivation of physical properties in the interstellar medium. By doing this, we provide a reliable estimate of the energy and momentum in the cloud and the shell-like structure of an embedded B-star wind. We find that, although this one shell is not powerful enough to sustain the turbulence in the entire cloud, it does inject a substantial amount of energy into the cloud. Here we list our conclusions:

\begin{itemize}
\item Three of the most common ways to trace the column density in molecular clouds include the near-infrared extinction, mid/far-infrared thermal emission and molecular line emission. Comparing the results using data from 2MASS, IRIS, Herschel and FCRAO, we find that the 2MASS/NICER-based extinction measurement serves as the most reliable tracer of column densities in the cloud. However, each of these different tracers, including the 2MASS/NICER-based extinction, suffers from its own limits. This is because none of these tracers directly trace the most dominant component of mass in molecular clouds: the molecular hydrogen. The variation of dust-to-gas ratio and the dust properties affect the effectiveness of both dust tracers, while the excitation criterion (for example, the critical density and temperature) makes it difficult to assess the conversion from the line emission to the gas column density. In dense regions, the molecular line emission can be optically thick or depleted, making it unsuitable for tracing the column density. Following \citet{Goodman_2009}, we propose that when choosing tracers of the column density, one should always consider the physical environment of the target region. 
\item Two of the complexities of using the mid/far-infrared emission as a density tracer are the variation of physical properties along the line of sight and the noise in the data. The commonly used modified blackbody equation assumes a single temperature, column density and the emissivity spectral index $\beta$ throughout the entire line of sight. This makes it difficult to describe the observed SED with the modified blackbody equation, especially when we free both the temperature and $\beta$ in the fitting. By building a toy line of sight, we find that the temperature and $\beta$ derived by fitting the modified blackbody equation are not representative of a line of sight with a varying temperature and $\beta$. On the other hand, we experiment with the hierarchical Bayesian-fitting method in the hope of decreasing the bias caused by noise in the data. We find that while it decreases the bias we see in previous works, including the anti-correlation between the temperature and $\beta$, the original code used in \citet{Kelly_2012} is not readily applicable to extended sources. To build a hierarchical Bayesian-fitting method that treats the distributions in the parameter space and of the observed fluxes at the same time, we need a more sophisticated model for these distributions than the multivariate Student-T distribution used in \citet{Kelly_2012}.
\item With a better understanding of various density tracers, we calculate the total turbulent energy in the Ophiuchus cloud and the energy entrained in the embedded B-star wind. By comparing the two, we find that the energy in the shell strucuture caused by the B-star wind amounts to $\sim$ 12 \% of the total turbulent energy in the cloud. The wind's energy injection rate is also comparable to the dissipation rate of the turbulent energy in the entire cloud. Comparing this result to shells in the Perseus cloud (Arce et al. 2011), we find that the shell in the Ophiuchus has energy comparable to the largest few shells in the Perseus and its energy injection rate is comparable to the \emph{sum} of all shells in the Perseus. We conclude that the embedded stellar winds can contribute to the turbulent energy in molecular clouds. In future simulations, embedded winds driven by B or even later types of stars should be included as important sources of turbulence in molecular clouds.
\end{itemize}

\selectlanguage{english}
\FloatBarrier
\bibliographystyle{plain}
\bibliography{bibliography/converted_to_latex.bib%
}

\end{CJK}\end{document}

