\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.