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