# Detection of distant AGN by MAGIC: the transparency of the Universe to high-
energy photons

Marco Roncadelli1, Alessandro De Angelis2, Oriana Mansutti3
1INFN  Sezione di Pavia  via A. Bassi 6  I – 27100 Pavia  Italy
2Dipartimento di Fisica  Università di Udine  Via delle Scienze 208  I – 33100
Udine  and INAF and INFN  Sezioni di Trieste  Italy
3Dipartimento di Fisica  Università di Udine  Via delle Scienze 208  I – 33100
Udine  and INFN  Sezione di Trieste  Italy

###### Abstract

The recent detection of blazar 3C279 by MAGIC has confirmed previous
indications by H.E.S.S. that the Universe is more transparent to very-high-
energy gamma rays than previously thought. We show that this fact can be
reconciled with standard blazar emission models provided photon oscillations
into a very light Axion-Like Particle occur in extragalactic magnetic fields.
A quantitative estimate of this effect explains the observed spectrum of
3C279. Our prediction can be tested in the near future by the satellite-borne
GLAST detector as well as by the ground-based Imaging Atmospheric Cherenkov
Telescopes H.E.S.S., MAGIC, CANGAROO III, VERITAS and by the Extensive Air
Shower arrays ARGO-YBJ and MILAGRO.

##  1 Introduction

As is well known, in the very-high-energy (VHE) band above $100\,{\rm GeV}$
the horizon of the observable Universe rapidly shrinks as the energy further
increases. This comes about because photons from distant sources scatter off
background photons permeating the Universe, thereby disappearing into
electron-positron pairs [1]. The corresponding cross section
$\sigma(\gamma\gamma\to e^{+}e^{-})$ peaks where the VHE photon energy $E$ and
the background photon energy $\epsilon$ are related by
$\epsilon\simeq(500\,{\rm GeV}/E)\,{\rm eV}$. Therefore, for observations
performed by Imaging Atmospheric Cherenkov Telescopes (IACTs) – which probe
the energy interval $100\,{\rm GeV}-100\,{\rm TeV}$ – the resulting cosmic
opacity is dominated by the interaction with ultraviolet/optical/infrared
diffuse background photons (frequency band $1.2\cdot 10^{3}\,{\rm
GHz}-1.2\cdot 10^{6}\,{\rm GHz}$, corresponding to the wavelength range
$0.25\,\mu{\rm m}-250\,\mu{\rm m}$), usually called Extragalactic Background
Light (EBL), which is produced by galaxies during the whole history of the
Universe. Neglecting evolutionary effects for simplicity, photon propagation
is controlled by the photon mean free path ${\lambda}_{\gamma}(E)$ for
$\gamma\gamma\to e^{+}e^{-}$, and so the observed photon spectrum $\Phi_{\rm
obs}(E,D)$ is related to the emitted one $\Phi_{\rm em}(E)$ by

$\Phi_{\rm obs}(E,D)=e^{-D/{\lambda}_{\gamma}(E)}\ \Phi_{\rm em}(E)~{}.$ (1)

Within the energy range in question, ${\lambda}_{\gamma}(E)$ decreases like a
power law from the Hubble radius $4.2\,{\rm Gpc}$ around $100\,{\rm GeV}$ to
$1\,{\rm Mpc}$ around $100\,{\rm TeV}$ [2]. Thus, Eq. (1) entails that the
observed flux is exponentially suppressed both at high energy and at large
distances, so that sufficiently far-away sources become hardly visible in the
VHE range and their observed spectrum should anyway be much steeper than the
emitted one.

Yet, observations have not detected the behaviour predicted by Eq. (1). A
first indication in this direction was reported by the H.E.S.S. collaboration
in connection with the discovery of the two blazars H2356-309 ($z=0.165$) and
1ES1101-232 ($z=0.186$) at $E\sim 1\,{\rm TeV}$ [3]. Stronger evidence comes
from the observation of blazar 3C279 ($z=0.536$) at $E\sim 0.5\,{\rm TeV}$ by
the MAGIC collaboration [4]. In particular, the signal from 3C279 collected by
MAGIC in the region $E<220$ GeV has more or less the same statistical
significance as the one in the range 220 GeV $<E<$ 600 GeV ($6.1\sigma$ in the
former case, $5.1\sigma$ in the latter).

A suggested way out of this difficulty relies upon the modification of the
standard Synchro-Self-Compton (SSC) emission mechanism. One option invokes
strong relativistic shocks [5]. Another rests upon photon absorption inside
the blazar [6]. While successful at substantially hardening the emission
spectrum, these attempts fail to explain why only for the most distant blazars
does such a drastic departure from the SSC emission spectrum show up.

Our proposal – usually referred to as the DARMA scenario – is quite different
[7]. Implicit in previous considerations is the hypothesis that photons
propagate in the standard way throughout cosmological distances. We suppose
instead that photons can oscillate into a new very light spin-zero particle –
named Axion-Like Parlicle (ALP) – and vice-versa in the presence of cosmic
magnetic fields, whose existence has definitely been proved by AUGER
observations [8]. Once ALPs are produced close enough to the source, they
travel unimpeded throughout the Universe and can convert back to photons
before reaching the Earth. Since ALPs do not undergo EBL absorption, the
effective photon mean free path ${\lambda}_{\gamma,{\rm eff}}(E)$ gets
increased so that the observed photons cross a distance in excess of
${\lambda}_{\gamma}(E)$. Correspondingly, Eq. (1) becomes

$\Phi_{\rm obs}(E,D)=e^{-D/{\lambda}_{\gamma,{\rm eff}}(E)}\ \Phi_{\rm
em}(E)~{},$ (2)

from which we see that even a slight increase of ${\lambda}_{\gamma,{\rm
eff}}(E)$ gives rise to a huge enhancement of the observed flux. It turns out
that the DARMA mechanism makes ${\lambda}_{\gamma,{\rm eff}}(E)$ shallower
than ${\lambda}_{\gamma}(E)$ although it remains a decreasing function of $E$.
So, the resulting observed spectrum is much harder than the one predicted by
Eq. (1), thereby ensuring agreement with observations even for a standard SSC
emission spectrum. As a bonus, we get a natural explanation for the fact that
only the most distant blazars would demand $\Phi_{\rm em}(E)$ to substantially
depart from the emission spectrum predicted by the SSC mechanism.

Our aim is to review the main features of our proposal as well as its
application to blazar 3C279.

##  2 DARMA scenario

Phenomenological as well as conceptual arguments lead to view the Standard
Model of particle physics as the low-energy manifestation of some more
fundamental and richer theory of all elementary-particle interactions
including gravity. Therefore, the lagrangian of the Standard Model is expected
to be modified by small terms describing interactions among known and new
particles. Many extensions of the Standard Model which have attracted
considerable interest over the last few years indeed predict the existence of
ALPs. They are spin-zero light bosons defined by the low-energy effective
lagrangian

${\cal L}_{\rm ALP}\ =\
\frac{1}{2}\,\partial^{\mu}\,a\,\partial_{\mu}\,a-\frac{1}{2}\,m^{2}\,a^{2}-\frac{1}{4M}\,F^{\mu\nu}\,\tilde{F}_{\mu\nu}\,a~{},$
(3)

where $F^{\mu\nu}$ is the electromagnetic field strength, $\tilde{F}_{\mu\nu}$
is its dual, $a$ denotes the ALP field whereas $m$ stands for the ALP mass.
According to the above view, it is assumed $M\gg G_{F}^{-1/2}\simeq 250\,{\rm
GeV}$. On the other hand, it is supposed that $m\ll G_{F}^{-1/2}\simeq
250\,{\rm GeV}$. The standard Axion [9] is the most well known example of ALP.
As far as generic ALPs are concerned, the parameters $M$ and $m$ are to be
regarded as independent.

So, what really characterizes ALPs is the trilinear $\gamma$-$\gamma$-$a$
vertex described by the last term in ${\cal L}_{\rm ALP}$, whereby one ALP
couples to two photons. Owing to this vertex, ALPs can be emitted by
astronomical objects of various kinds, and the present situation can be
summarized as follows. The negative result of the CAST experiment designed to
detect ALPs emitted by the Sun yields the bound $M>0.86\cdot 10^{10}\,{\rm
GeV}$ for $m<0.02\,{\rm eV}$ [10]. Moreover, theoretical considerations
concerning star cooling via ALP emission provide the generic bound
$M>10^{10}\,{\rm GeV}$, which for $m<10^{-10}\,{\rm eV}$ gets replaced by the
stronger one $M>10^{11}\,{\rm GeV}$ even if with a large uncertainty [11]. The
same $\gamma$-$\gamma$-$a$ vertex produces an off-diagonal element in the mass
matrix for the photon-ALP system in the presence of an external magnetic field
${\bf B}$. Therefore, the interaction eigenstates differ from the propagation
eigenstates and photon-ALP oscillations show up [12].

We imagine that a sizeable fraction of photons emitted by a blazar soon
convert into ALPs. They propagate unaffected by the EBL and we suppose that
before reaching the Earth a substantial fraction of ALPs is back converted
into photons. We further assume that this photon-ALP oscillation process is
triggered by cosmic magnetic fields (CMFs), whose existence has been
demonstrated very recently by AUGER observations [8]. Owing to the notorious
lack of information about their morphology, one usually supposes that CMFs
have a domain-like structure [13]. That is, ${\bf B}$ ought to be constant
over a domain of size $L_{\rm dom}$ equal to its coherence length, with ${\bf
B}$ randomly changing its direction from one domain to another but keeping
approximately the same strength. As explained elsewhere [14], it looks
plausible to assume the coherence length in the range $1-10\,{\rm Mpc}$.
Correspondingly, the inferred strength lies in the range $0.3-1.0\,{\rm nG}$
[14].

##  3 Predicted energy spectrum

Our ultimate goal consists in the evaluation of the probability
$P_{\gamma\to\gamma}(E,D)$ that a photon remains a photon after propagation
from the source to us when allowance is made for photon-ALP oscillations as
well as for photon absorption from the EBL. As a consequence, Eq. (2) gets
replaced by

$\Phi_{\rm obs}(E,D)=P_{\gamma\to\gamma}(E,D)\,\Phi_{\rm em}(E)~{}.$ (4)

We proceed as follows. We first solve exactly the beam propagation equation
arising from ${\cal L}_{\rm ALP}$ over a single domain, assuming that the EBL
is described by the “best-fit model” of Kneiske et al. [15]. Starting with an
unpolarized photon beam, we next propagate it by iterating the single-domain
solution as many times as the number of domains crossed by the beam, taking
each time a random value for the angle between ${\bf B}$ and a fixed overall
fiducial direction. We repeat such a procedure $10^{.}000$ times and finally
we average over all these realizations of the propagation process.

We find that about 13% of the photons arrive to the Earth for $E=500\,{\rm
GeV}$, representing an enhancement by a factor of about 20 with respect to the
expected flux without DARMA mechanism (the comparison is made with the above
“best-fit model”). The same calculation gives a fraction of 76% for
$E=100\,{\rm GeV}$ (to be compared to 67% without DARMA mechanism) and a
fraction of 3.4% for $E=1\,{\rm TeV}$ (to be compared to 0.0045% without DARMA
mechanism). The resulting spectrum is exhibited in Fig. 1. The solid line
represents the prediction of the DARMA scenario for $B\simeq 1\,{\rm nG}$ and
$L_{\rm dom}\simeq 1\,{\rm Mpc}$ and the gray band is the envelope of the
results obtained by independently varying ${\bf B}$ and $L_{\rm dom}$ within a
factor of 10 about such values. These conclusions hold for $m\ll
10^{-10}\,{\rm eV}$ and we have taken for definiteness $M\simeq 4\cdot
10^{11}\,{\rm GeV}$ but we have cheked that practically nothing changes for
$10^{11}\,{\rm GeV}<M<10^{13}\,{\rm GeV}$.

Our prediction can be tested in the near future by the satellite-borne GLAST
detector as well as by the ground-based IACTs H.E.S.S., MAGIC, CANGAROO III,
VERITAS and by the Extensive Air Shower arrays ARGO-YBJ and MILAGRO.

Figure 1: The two lowest lines give the fraction of photons surviving from
3C279 without the DARMA mechanism within the “best-fit model” of EBL (dashed
line) and for the minimum EBL density compatible with cosmology (dashed-dotted
line) [15]. The solid line represents the prediction of the DARMA mechanism as
explained in the text.

## References

  * [1] G. G. Fazio and F. W. Stecker, Nature 226 135 (1970).
  * [2] P. Coppi and F. Aharonian, Astrophys. J. 487 L9 (1997).
  * [3] F. Aharonian et al. (H.E.S.S. Collaboration), Nature 440 1018 (2006).
  * [4] J. Albert et al. (MAGIC Collaboration), Science 320 1752 (2008).
  * [5] F. W. Stecker, M. G. Baring and E. J. Summerlin, Astrophys. J. 667 L29 (2007). F. W. Stecker and S. T. Scully, Astron. Astrophys 478 L1 (2008).
  * [6] F. Aharonian, D. Khangulyan and L. Costamante, arXiv:0801.3198 (2008).
  * [7] A. De Angelis, M. Roncadelli and O. Mansutti, Phys. Rev. D76 121301 (2007).
  * [8] J. Abraham et al. [Pierre Auger Collaboration], Science 318 939 (2007).
  * [9] J. H. Kim, Phys. Rep. 150 1 (1987). H. Y. Cheng, Phys. Rep. 158 1 (1988).
  * [10] K. Zioutas et al., Phys. Rev. Lett. 94 121301 (2005). S. Andriamoje et al., JCAP 0704010 (2007).
  * [11] G. G. Raffelt, Stars as Laboratories for Fundamental Physics (University of Chicago Press, Chicago, 1996).
  * [12] P. Sikivie, Phys. Rev. Lett. 51 1415 (1983); (E) ibid. 52 695 (1984). L. Maiani, R. Petronzio and E. Zavattini, Phys. Lett. B175 359 (1986). G. G. Raffelt and L. Stodolsky, Phys. Rev. D37 1237 (1988).
  * [13] P. P. Kronberg, Rept. Prog. Phys. 57 325 (1994). D. Grasso and H. Rubinstein, Phys. Rep. 348 163 (2001).
  * [14] A. De Angelis, M. Persic and M. Roncadelli, Mod. Phys. Lett. A23 315 (2008).
  * [15] T. M. Kneiske, T. Bretz, K. Mannheim and D. H. Hartmann, Astron. Astrophys. 413 807 (2004).


Sampling in Combinatorial Spaces with
SurVAE Flow Augmented MCMC

Priyank Jaini Didrik Nielsen Max Welling

Bosch-Delta Lab University of Amsterdam Technical University of Denmark Bosch-
Delta Lab University of Amsterdam

###### Abstract

Hybrid Monte Carlo is a powerful Markov Chain Monte Carlo method for sampling
from complex continuous distributions. However, a major limitation of HMC is
its inability to be applied to discrete domains due to the lack of gradient
signal. In this work, we introduce a new approach based on augmenting Monte
Carlo methods with SurVAE Flows to sample from discrete distributions using a
combination of neural transport methods like normalizing flows and variational
dequantization, and the Metropolis-Hastings rule. Our method first learns a
continuous embedding of the discrete space using a surjective map and
subsequently learns a bijective transformation from the continuous space to an
approximately Gaussian distributed latent variable. Sampling proceeds by
simulating MCMC chains in the latent space and mapping these samples to the
target discrete space via the learned transformations. We demonstrate the
efficacy of our algorithm on a range of examples from statistics,
computational physics and machine learning, and observe improvements compared
to alternative algorithms.

##  1 Introduction

The ability to draw samples from a known distribution is a fundamental
computational challenge. It has applications in diverse fields like
statistics, probability, and stochastic modeling where these methods are
useful for both estimation and inference. These are further useful within the
frequentist inference framework to form confidence intervals for a point
estimate. Sampling procedures are also standard in the Bayesian setting for
exploring posterior distributions, obtaining credible intervals, and solving
inverse problems. The workhorse algorithms in these settings are simulation-
based, amongst which the Markov Chain Monte Carlo Brooks et al., (2011) method
is the most broadly used method. Impressive advances have been made – both in
increasing efficiency and reducing computation costs – for sampling using
Monte Carlo methods over the past century. However, problems in discrete
domains still lack an efficient general-purpose sampler.

In this paper, we consider the problem of sampling from a known discrete
distribution. Inspired by the recent success of deep generative models –
particularly neural transport methods like normalizing flows (Tabak and
Vanden-Eijnden,, 2010; Tabak and Turner,, 2013; Rezende and Mohamed,, 2015)–
in unsupervised learning, we propose a new approach to design Monte Carlo
methods to sample from a discrete distribution based on augmenting MCMC with a
_transport map_. Informally, let $\pi(\boldsymbol{\theta})$ be a discrete
target distribution and $p(\mathbf{z})$ be a simple source density from which
it is easy to generate independent samples e.g. a Gaussian distribution. Then,
a transport map $\mathbf{T}$ from $p(\mathbf{z})$ to
$\pi(\boldsymbol{\theta})$ is such that if $\mathbf{z}_{i}\sim p(\mathbf{z})$
then $\mathbf{T}(\mathbf{z}_{i})\sim\pi(\boldsymbol{\theta})$.

The significance of having such a transport map $\mathbf{T}$ is particularly
consequential: firstly, given such a map $\mathbf{T}$, we can generate samples
from the target $\pi(\boldsymbol{\theta})$. Secondly, these samples can be
generated cheaply irrespective of the cost of evaluating
$\pi(\boldsymbol{\theta})$. Importantly, $\mathbf{T}$ affords us the ability
to sample from the marginal and conditional distributions of
$\pi(\boldsymbol{\theta})$ using $p(\mathbf{z})$ given an appropriate
structure. Indeed, this idea of using neural transport maps based on
normalizing flows have been explored by Parno and Marzouk, (2018) for
continuous densities by learning a diffeomorphic transformation $\mathbf{T}$
from the source density $p$ to the target density $\pi(\boldsymbol{\theta})$
where $\boldsymbol{\theta}\in\mathds{R}^{d}$.

In this paper, we extend this to discrete domains using the recently proposed
SurVAE Flow framework (Nielsen et al.,, 2020). We first learn a transport map
from the discrete space to a continuous space using a surjective
transformation $\mathbf{Q}$. However, such a continuous space is often highly
multi-modal with unfavorable geometry for fast mixing of MCMC algorithms.
Thus, we learn an additional normalizing flow $\mathbf{T}_{\phi}$ that
transforms this _complex_ continuous space to a simple latent space with
density $p(\mathbf{z})$ where sampling is easy. We finally sample from the
desired target distribution by generating samples from the latent space and
using the learned transformations $\mathbf{Q}$ and $\mathbf{T}$ to push-
forward these samples to the target space. Our complete implementation allows
parallelization over multiple GPUs, improving efficiency and reducing
computation time.

The main contributions of this paper are:

  * •

We present a SurVAE augmented general purpose MCMC solver for combinatorial
spaces.

  * •

We propose a new learning objective compared to previous methods in
normalizing flows to train a transport map that “Gaussianizes” the discrete
target distribution.

The rest of the manuscript is organized as follows: We begin in §2 by
presenting a brief background on normalizing flows and setting up our main
problem. In §3, we provide details of our method which consists of two parts:
learning the transport map and, generating samples from the target
distribution. Subsequently, in §4, we put our present work in perspective to
other approaches to sampling in discrete spaces. Finally, we perform empirical
evaluation on a diverse suite of problems in §5 to demonstrate the efficacy of
our method. Our implementation is available at
https://github.com/priyankjaini/discFlowMH

##  2 Preliminaries and Setup

In this section, we setup our main problem, provide key definitions and
notations, and formulate the idea for an MCMC sampler augmented with
normalizing flows for sampling in discrete spaces.

Let $p,q$ be two probability density functions (w.r.t. the Lebesgue measure)
over the source domain $\mathsf{Z}\subseteq\mathds{R}^{d}$ and the target
domain $\mathsf{X}\subseteq\mathds{R}^{d}$, respectively. Normalizing flows
learn a _diffeomorphic_ transformation $\mathbf{T}:\mathsf{Z}\to\mathsf{X}$
that allows to represent $q$ using $p$ via the change of variables formula
(Rudin,, 1987):

$\displaystyle q(\mathbf{x})$
$\displaystyle=p(\mathbf{z})/|\mathbf{T}^{\prime}(\mathbf{z})|$
$\displaystyle=p(\mathbf{T}^{-1}\mathbf{x})/|\mathbf{T}^{\prime}(\mathbf{T}^{-1}\mathbf{x})|,$

where $|\mathbf{T}^{\prime}(\mathbf{z})|$ is the (absolute value) of the
Jacobian (determinant of the derivative) of $\mathbf{T}$. In other words, we
can obtain a new random variable $\mathbf{x}\sim q$ by pushing the source
random variable $\mathbf{z}\sim p$ through the map $\mathbf{T}$. When we only
have access to an i.i.d. sample
$\lbag\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\rbag\sim q$, we can learn
$\mathbf{T}$ and thus $q$ through maximum likelihood estimation:

$\displaystyle\max_{\mathbf{T}\in\mathcal{F}}~{}~{}\frac{1}{n}\sum_{i=1}^{n}\Big{[}-\log|\mathbf{T}^{\prime}(\mathbf{T}^{-1}\mathbf{x}_{i})|+\log
p(\mathbf{T}^{-1}\mathbf{x}_{i})\Big{]}.$

where $\mathcal{F}$ is a class of diffeomorphic mappings. Conveniently, we can
choose _any_ source density $p$ to facilitate estimation e.g. standard normal
density on $\mathsf{Z}=\mathds{R}^{d}$ (with zero mean and identity
covariance) or uniform density over the cube $\mathsf{Z}=[0,1]^{d}$.

This “push-forward” idea has played an important role in optimal transport
theory (Villani,, 2008) and has been used successfully for Monte carlo
simulations. For example Marzouk et al., (2016); Parno and Marzouk, (2018);
Peherstorfer and Marzouk, (2018); Hoffman et al., (2019); Albergo et al.,
(2019) have used normalizing flows for continuous random variables to address
the limitations of HMC which suffers from the chain to mix slowly between
distant states when the geometry of the target density is unfavourable.

Specifically, Parno and Marzouk, (2018) addressed this problem of sampling in
a space with difficult geometry by learning a diffeomorphic transport map that
transforms the original random variable to another random variable with a
simple distribution. Concretely, let our interest be to sample from
$\pi(\boldsymbol{\theta})$ where
$\boldsymbol{\theta}\in\Theta\subseteq\mathds{R}^{d}$. We can proceed by
learning a diffeomorphic map $\mathbf{T}:\mathsf{Z}\to\Theta$ such that
$\tilde{p}(\mathbf{z})=\pi(\boldsymbol{\theta})\cdot|\mathbf{T}^{\prime}(\mathbf{z})|$
where $\mathbf{z}=\mathbf{T}^{-1}(\boldsymbol{\theta})$ such that
$p(\mathbf{z})$ has a simple geometry amenable to efficient MCMC sampling.
Thus, samples can be generated from $\pi(\boldsymbol{\theta})$ by running MCMC
chain in the z-space and pushing these samples onto the $\Theta$-space using
$\mathbf{T}$. The transformation $\mathbf{T}$ can be learned by minimizing the
KL-divergence between a fixed distribution with simple geometry in the z-space
e.g. a standard Gaussian and $\tilde{p}(\mathbf{z})$ above. The learning phase
attempts to ensure that the distribution $\tilde{p}(\mathbf{z})$ is
approximately close to the fixed distribution with easy geometry so that MCMC
sampling is efficient. A limitation of these works that use diffeomorphic
transformations augmented samplers is that they are restricted to continuous
random variables. This is primarily because the flow models used in these
works can only learn density functions over continuous variables.

Uria et al., (2013) introduced the concept of _dequantization_ to extend
normalizing flows to discrete random variables. They consider the problem of
estimating the discrete distribution $\pi(\boldsymbol{\theta})$ given samples
$\lbag\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2},\cdots,\boldsymbol{\theta}_{n}\rbag\sim\pi$
by “lifting” the discrete space $\Theta$ to a continuous one $\mathsf{X}$ by
filling the gaps in the discrete space with a uniform distribution i.e.
$\mathbf{x}\in\mathsf{X}$ is such that
$\mathbf{x}:=\boldsymbol{\theta}+\mathbf{u}|\boldsymbol{\theta}$ where
$\boldsymbol{\theta}\in\Theta$ and
$\mathbf{u}|\boldsymbol{\theta}\sim\mathsf{Uniform}(0,1)$. Subsequently, they
learn the continuous distribution $p_{\mathsf{model}}(\mathbf{x})$ over
$\mathsf{X}$ using a normalizing flow by maximizing the log-likelihood of the
continuous model $p_{\mathsf{model}}(\mathbf{x})$. Theis et al., (2015) showed
that maximizing likelihood of the continuous model is equivalent to maximizing
a lower bound on the log-likelihood for a certain discrete model
$\pi_{\mathsf{model}}:=\int_{\mathbf{u}}p_{\mathsf{model}}(\boldsymbol{\theta}+\mathbf{u})\mathrm{d}\mathbf{u}$
on the original discrete data. Thus, this learning procedure cannot lead to
the continuous model degenerately collapsing onto the discrete data, because
its objective is bounded above by the log-likelihood of a discrete model. Ho
et al., (2019) extended the idea of uniform dequantization to propose
_variational dequantization_ where in $\mathbf{u}|\boldsymbol{\theta}\sim
q(\mathbf{u}|\boldsymbol{\theta})$ instead of the uniform distribution. They
learn the dequantizing distribution, $q(\mathbf{u}|\boldsymbol{\theta})$ from
data by treating it as a variational distribution. They subsequently estimate
$\pi(\boldsymbol{\theta})$ by optimizing the following lower bound:

$\displaystyle\mathbb{E}_{\boldsymbol{\theta}\sim\pi}[\log\pi_{\mathsf{model}}(\boldsymbol{\theta})]\geq\mathbb{E}_{\begin{subarray}{c}\boldsymbol{\theta}\sim\pi\\\
\mathbf{u}\sim
q(\mathbf{u}|\boldsymbol{\theta})\end{subarray}}\Bigg{[}\log\frac{p_{\mathsf{model}}(\mathbf{x})}{q(\mathbf{u}|\boldsymbol{\theta})}\Bigg{]}$
(1)

Recently, Nielsen et al., (2020) introduced SurVAE Flows that extends the
framework of normalizing flows to include surjective and stochastic
transformations for probabilistic modelling. In the SurVAE flow framework,
dequantization can be seen as a _rounding_ surjective map
$\mathbf{Q}_{\lambda}:\mathsf{X}\to\Theta$ with parameters $\lambda$ such that
$\boldsymbol{\theta}:=\lfloor\mathbf{x}\rfloor$ where the forward
transformation is a discrete surjection
$P(\boldsymbol{\theta}|\mathbf{x})=\mathbb{I}(\mathbf{x}\in\mathcal{B}(\boldsymbol{\theta}))$,
for
$\mathcal{B}(\boldsymbol{\theta})=\\{\boldsymbol{\theta}+\mathbf{u}|\mathbf{u}\in[0,1)^{d}\\}$.
The inverse transformation
$\mathbf{Q}_{\lambda}^{-1}:=q(\mathbf{x}|\boldsymbol{\theta})$ is stochastic
with support in $\mathcal{B}(\boldsymbol{\theta})$. Thus, these ideas of using
dequantization to learn a discrete probability distribution can be viewed as
learning a transformation $\mathbf{Q}_{\lambda}^{-1}$ that transforms a
discrete space $\Theta$ to a continuous space $\mathsf{X}$.

In this paper, we study the problem of sampling in discrete spaces i.e. let
$\Theta\subseteq\\{1,2,\cdots,K\\}^{d}:=[K]^{d}$ be a discrete random variable
in $d$-dimensions with probability mass function (potentially unnormalized)
$\pi(\boldsymbol{\theta}),~{}\boldsymbol{\theta}\in\Theta$. Given access to a
function that can compute
$\pi(\boldsymbol{\theta}),~{}\forall\boldsymbol{\theta}\in\Theta$, we aim to
generate samples
$\lbag\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2},\cdots,\boldsymbol{\theta}_{n}\rbag\sim\pi(\boldsymbol{\theta})$.
The aforementioned works on normalizing flows for discrete data is not
directly applicable in this regime. This is due to the fact that Uria et al.,
(2013); Ho et al., (2019) and Nielsen et al., (2020) estimate
$\pi(\boldsymbol{\theta})$ given samples from the original distribution by
maximizing the log-likelihood of a discrete model. However, we only have
access to the function $\pi(\boldsymbol{\theta})$ and the learning method
given in Equation 1 cannot be used in our setting. In the next section, we
address this by extending the ideas of neural transport MCMC (Parno and
Marzouk,, 2018; Hoffman et al.,, 2019; Albergo et al.,, 2019) and leveraging
SurVAE Flows (Nielsen et al.,, 2020).

##  3 Flow Augmented MCMC

We discussed in Section 2 the utility of normalizing flows augmented Monte
Carlo samplers for overcoming the difficulties posed by unfavourable geometry
in the target continuous distribution, as evidenced by the works of Parno and
Marzouk, (2018); Hoffman et al., (2019) and Albergo et al., (2019). We will
now introduce our method of SurVAE Flow augmented MCMC for sampling in
combinatorial spaces.

Informally, our method proceeds as follows: We first define a rounding
surjective transformation $\mathbf{Q}_{\lambda}:\mathsf{X}\to\Theta$ with
parameters $\lambda$ such that $\mathsf{X}$ is a continuous embedding of the
$\Theta$ space with density $q(\mathbf{x})$. Since the continuous embedding
may be highly multi-modal with potentially unfavourable geometry for efficient
MCMC sampling, we define an additional diffeomorphic transformation,
$\mathbf{T}_{\phi}:\mathsf{Z}\to\mathsf{X}$ with parameters $\phi$, from a
simple latent space $\mathsf{Z}$ to $\mathsf{X}$. Subsequently, we learn these
transformations via maximum likelihood estimation. This concludes the learning
phase of our algorithm. Finally, we generate samples from
$\pi(\boldsymbol{\theta})$ by running MCMC chains to generate samples from the
learned distribution over $\mathsf{Z}$ and pushing-forward these samples to
the $\Theta$ space using $\mathbf{Q}_{\lambda}$ and $\mathbf{T}_{\phi}$. We
elaborate each of these steps below.

$\mathbf{Q}_{\lambda}$ is a surjective transformation that takes as input a
real-valued vector $\mathbf{x}\in\mathsf{X}\subseteq\mathds{R}^{d}$ and
returns the _rounded_ value of $\mathbf{x}$ i.e.
$\boldsymbol{\theta}=\lfloor\mathbf{x}\rfloor$. Thus, the forward
transformation from $\mathsf{X}$ to $\Theta$ is deterministic. The inverse
transformation from $\Theta$ to $\mathsf{X}$ is however stochastic since the
random variable $\mathsf{X}:=\Theta+\mathsf{U}_{|\boldsymbol{\theta}}$ where
$\mathsf{U}_{|\boldsymbol{\theta}}\subseteq[0,1)^{d}$ is given by
$\mathbf{x}:=\boldsymbol{\theta}+\mathbf{u}$ where $\mathbf{u}\sim
q(\mathbf{u}|\boldsymbol{\theta})$. Given access to
$q(\mathbf{u}|\boldsymbol{\theta})$, we can evaluate the density
$q(\mathbf{x})$ at a point $\mathbf{x}\in\mathds{R}^{d}$ exactly as:

$\displaystyle
q(\mathbf{x})=q(\boldsymbol{\theta}+\mathbf{u}):=\pi(\boldsymbol{\theta})\cdot
q(\mathbf{u}|\boldsymbol{\theta})$ (2)

since $\mathbf{x}\in\mathcal{B}(\boldsymbol{\theta})$. Thus, we need to learn
$q(\mathbf{u}|\boldsymbol{\theta})$ in-order to fully specify this surjective
transformation. Since, $q(\mathbf{u}|\boldsymbol{\theta})$ can be any
arbitrary continuous density, we learn $q(\mathbf{u}|\boldsymbol{\theta})$
using a normalizing flow i.e. we learn a diffeomorphic transformation
$\mathbf{T}_{\lambda}(\cdot;\boldsymbol{\theta}):\mathcal{E}\to\mathsf{U}_{|\boldsymbol{\theta}}$
where $\mathcal{E}\subseteq\mathds{R}^{d}$ is standard Gaussian distributed.
Under this setup, we get

$\displaystyle
q(\mathbf{u}|\boldsymbol{\theta})=p(\boldsymbol{\varepsilon})\cdot|\mathbf{T}_{\lambda}^{\prime}(\boldsymbol{\varepsilon};\boldsymbol{\theta})|^{-1}$
(3)

Thus, learning the surjective transformation $\mathbf{Q}$ is equivalent to
learning $q(\mathbf{u}|\boldsymbol{\theta})$ which reduces to learning a flow
$\mathbf{T}_{\lambda}(\cdot;\boldsymbol{\theta})$. Thus, for brevity we use
the informal notation that the transformation
$\mathbf{T}_{\lambda}(\cdot;\boldsymbol{\theta})$ quantizes the continuous
space $\mathsf{X}$ to $\boldsymbol{\theta}$. Finally, using Equation 3, the
density $q(\mathbf{x})$ in Equation 2 for a fixed
$\mathbf{T}_{\lambda}(\cdot;\boldsymbol{\theta})$ can be written as:

$\displaystyle q(\mathbf{x})=\pi(\boldsymbol{\theta})\cdot
p(\boldsymbol{\varepsilon})\cdot|\mathbf{T}_{\lambda}^{\prime}(\boldsymbol{\varepsilon};\boldsymbol{\theta})|^{-1}$
(4)

The density $q(\mathbf{x})$ over the continuous embedding $\mathsf{X}$ of
$\Theta$ learned above can have any arbitrary geometry and may not be
efficient for MCMC sampling. Thus, next we learn a diffeomorphic
transformation $\mathbf{T}_{\phi}$ from a latent space $\mathsf{Z}$ to
$\mathsf{X}$ such that using the change of variables formula we get:

$\displaystyle\tilde{p}_{\phi,\lambda}(\mathbf{z})=q(\mathbf{x})\cdot|\mathbf{T}_{\phi}^{\prime}(\mathbf{T}^{-1}_{\phi}\mathbf{x})|$
(5)

Our complete model, therefore, consists of transformations $\mathbf{T}_{\phi}$
and $\mathbf{T}_{\lambda}(\cdot;\boldsymbol{\theta})$ that transforms
$\mathsf{Z}$ to $\Theta$. The next challenge is to learn $\mathbf{T}_{\phi}$
and $\mathbf{T}_{\lambda}(\cdot;\boldsymbol{\theta})$ such that the induced
density $\tilde{p}_{\phi,\lambda}(\mathbf{z})$ has a simple geometry for
efficient MCMC sampling. We achieve this by forcing
$\tilde{p}_{\phi,\lambda}(\mathbf{z})$ to be close to a standard Gaussian by
minimizing the KL-divergence between $\tilde{p}_{\phi,\lambda}(\mathbf{z})$
and $p(\mathbf{z})$ where
$p(\mathbf{z})=(2\pi)^{-\frac{d}{2}}\mathsf{exp}(-\mathbf{z}^{T}\mathbf{z})$.

Concretely, let
$\mathcal{L}(\lambda,\phi):=\mathsf{KL}\big{(}p(\mathbf{z})~{}||~{}\tilde{p}_{\phi,\lambda}(\mathbf{z})\big{)}$.
$\mathcal{L}(\lambda,\phi)$ can be approximated with the empirical average by
generating i.i.d. samples
$\lbag\mathbf{z}_{1},\mathbf{z}_{2},\cdots,\mathbf{z}_{m}\rbag\sim
p(\mathbf{z})$ giving:

$\displaystyle\mathcal{L}(\lambda,\phi)\approx\frac{1}{m}\sum_{i=1}^{m}\log\frac{p(\mathbf{z}_{i})}{\tilde{p}_{\phi,\lambda}(\mathbf{z}_{i})}$

We thus arrive at the following optimization problem:

$\displaystyle\mathbf{T}_{{\lambda}}^{*},\mathbf{T}_{{\phi}}^{*}$
$\displaystyle:=\arg\min_{\lambda,\phi}~{}\mathcal{L}(\lambda,\phi)\approx\arg\max_{\lambda,\phi}\frac{1}{m}\sum_{i=1}^{m}\log\tilde{p}_{\phi,\lambda}(\mathbf{z}_{i})$

Using Equations (4) and (5), we can rewrite our final objective as:

$\mathbf{T}_{\lambda}^{*},\mathbf{T}_{\phi}^{*}=\arg\max_{\lambda,\phi}\Bigg{(}\frac{1}{m}\sum_{i=1}^{m}\log\pi(\boldsymbol{\theta}_{i})+\log
p(\boldsymbol{\varepsilon}_{i})\\\
-\log|\mathbf{T}_{\lambda}^{\prime}(\boldsymbol{\varepsilon}_{i};\boldsymbol{\theta}_{i})|+|\mathbf{T}_{\phi}^{\prime}(\mathbf{T}_{\phi}^{-1}\mathbf{x}_{i})|\Bigg{)}$
(6)

where $\forall i\in[m],~{}\mathbf{x}_{i}=\mathbf{T}_{\phi}(\mathbf{z}_{i})$,
$\boldsymbol{\theta}_{i}=\lfloor\mathbf{T}_{\phi}(\mathbf{z}_{i})\rfloor$,
$\mathbf{u}_{i}=\mathbf{x}_{i}-\boldsymbol{\theta}_{i}$ and
$\boldsymbol{\varepsilon}_{i}=\mathbf{T}_{\lambda}^{-1}(\mathbf{u}_{i})$.

Learning $\mathbf{T}_{\lambda}^{*},\mathbf{T}_{\phi}^{*}$ results in a density
$\tilde{p}^{*}_{\phi,\lambda}(\mathbf{z})$ that is approximately Gaussian with
a landscape amenable to efficient MCMC sampling. Thus, our sampling phase
consists of running an MCMC sampler of choice with target density
$\tilde{p}(\mathbf{z})$ resulting in samples
$\lbag\mathbf{z}_{1},\mathbf{z}_{2},\cdots,\mathbf{z}_{n}\rbag\sim\tilde{p}^{*}_{\phi,\lambda}(\mathbf{z})$.
We can finally obtain samples
$\lbag\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2},\cdots,\boldsymbol{\theta}_{n}\rbag\sim\pi(\boldsymbol{\theta})$
as
$\boldsymbol{\theta}_{i}=\lfloor\mathbf{T}_{\phi}^{*}(\mathbf{z}_{i})\rfloor$.

We end this section with two important remarks.

######  Remark 1.

In our method, it is possible to by-pass the last step of sampling in the
z-space using MCMC entirely since the trained flow is a generative model that
can be used directly to generate samples from $\pi(\boldsymbol{\theta})$. This
can indeed be done if the learned transformations
$\mathbf{T}_{\lambda}(\cdot;\boldsymbol{\theta})^{*}$ and
$\mathbf{T}_{\phi}^{*}$ result in the density $\tilde{p}(\mathbf{z})$ to be
Gaussian. Thus, there is a natural trade-off here: we can spend only enough
computation to train the flow to learn $\tilde{p}(\mathbf{z})$ that is
suitable for fast-mixing of the MCMC chain and generate samples from
$\tilde{p}(\mathbf{z})$ that is not Gaussian using any sampler of choice or we
can spend a larger amount of compute to learn a flow that perfectly
Gaussianizes the target density $\pi(\boldsymbol{\theta})$. Then, we can
sample directly by sampling from a Gaussian in the z-space and using the
learned transformations to obtain samples from $\pi(\boldsymbol{\theta})$.

######  Remark 2.

As mentioned in Section 2, we can use any density $p(\mathbf{z})$ instead of a
Gaussian for training the transformations
$\mathbf{T}_{\lambda}(\cdot;\boldsymbol{\theta})$ and $\mathbf{T}_{\phi}$. The
main motivation of our method is to “push-forward” $\pi(\boldsymbol{\theta})$
onto a space that is amenable to efficient sampling. An interesting future
work might be to devise learning objectives that explicitly drive the learned
density to have simple geometry that favours “off-the-shelf” samplers.

##  4 Connection to previous works

In Section 3 we introduced a normalizing flow augmented MCMC sampler for
combinatorial spaces. Our method combines surjective transformation for
learning continuous relaxations for discrete spaces and normalizing flows that
map the continuous space to a simpler discrete space for easy MCMC sampling.
In this section, we put both these ideas of continuous relaxations of discrete
spaces and neural transport methods for efficient MCMC sampling in to
perspective with existing work.

### Neural transport samplers for continuous spaces:

As we discussed briefly in Section 2, the neural transport augmented sampler
for continuous variables has been successfully used by Parno and Marzouk,
(2018), Peherstorfer and Marzouk, (2018), Hoffman et al., (2019) and Albergo
et al., (2019). A subtle difference between these works and our work here –
apart from the major difference that the aforementioned works are applicable
only for continuous domains – is the method of training the transport map
itself. These methods train the transport map (or normalizing flow) by
minimizing the Kullback-Liebler divergence between the target density
$\pi(\boldsymbol{\theta})$ and the density $\tilde{\pi}(\boldsymbol{\theta})$
learned via $\mathbf{T}$ by pushing-forward a standard normal distribution.
Albergo et al., (2019) additionally use the reverse KL-divergence for the
target density and the approximate density for training instead. We, on the
other hand, learn by minimizing the KL-divergence in the latent space since we
only have access to a lower-bound of the discrete density
$\tilde{\pi}(\boldsymbol{\theta})$ (cf. Equation 1).

### Continuous relaxations of discrete spaces:

The idea to relax the constraint that random variables of interest only take
discrete values has been used extensively in combinatorial optimization
(Pardalos,, 1996). Such an approach is attractive since the continuous space
affords the function with gradient information, contours, and curvatures that
can better inform optimization algorithms. Surprisingly though, this idea has
not received much attention in the MCMC setting. Perhaps the closest work to
ours in the present manuscript is that of Zhang et al., (2012) who use the
Gaussian Integral Trick (Hubbard,, 1959) to transform discrete variable
undirected models into fully continuous systems where they perform HMC for
inference and the evaluation of the normalization constant. The Gaussian
Integral Trick used by Zhang et al., (2012) can be viewed as specifying a
fixed map from the discrete space to an augmented continuous space. However,
this can result in a continuous space that is highly multi-modal and not
amenable for efficient sampling.

Nishimura et al., (2020) on the other hand map to a continuous space using
_uniform dequantization_ (Uria et al.,, 2013) i.e. filling the space between
points in the discrete space with uniform noise inducing parameters with
piecewise constant densities. 111This corresponds to
$q(\mathbf{u}|\boldsymbol{\theta})=\mathsf{Uniform}(0,1)$ in our method. They
further propose a Laplace distribution for the momentum variables in dealing
with discontinuous targets and argue this to be more effective than the
Gaussian distribution. This work relies on the key theoretical insight of
Pakman and Paninski, (2013) that Hamiltonian dynamics with a discontinuous
potential energy function can be integrated explicitly near the discontinuity
such that it preserves the total energy. Pakman and Paninski, (2013) used this
to propose a sampler for binary distributions which was later extended to
handle more general discontinuities by Afshar and Domke, (2015). Dinh et al.,
(2017) also used this idea on settings where the parameter space involves
phylogenetic trees. A major limitation of Pakman and Paninski, (2013), Afshar
and Domke, (2015) and Dinh et al., (2017) is the fact that these method do not
represent a general-purpose solver since they encounter computational issues
when dealing with complex discontinuities. Similarly, the work of Nishimura et
al., (2020) requires an integrator that works component-wise and is
prohibitively slow for high dimensional problems. Furthermore, it is not clear
if the the Laplace momentum based sampler leads to efficient exploration of
the continuous space which is highly multi-modal due to uniform
dequantization.

### Discrete MCMC samplers:

A main limitation of embedding discrete spaces into continuous ones is that
they can often destroy natural topological properties of the space under
consideration e.g. space of trees, partitions, permutations etc. Titsias and
Yau, (2017) proposed an alternative approach called Hamming ball sampler based
on informed proposals that are obtained by augmenting the discrete space with
auxiliary variables and performing Gibbs sampling in this augmented space.
However, potentially strong correlations between auxiliary variables and the
chain state severely slows down convergence. Zanella, (2020) tried to address
this problem by introducing locally balanced proposals that incorporate local
information about the target distribution. This framework was later called a
_Zanella process_ by Power and Goldman, (2019). They used the insights in
(Zanella,, 2020) to build efficient, continuous-time, non-reversible
algorithms by exploiting the structure of the underlying space through
symmetries and group-theoretic notions. This helps them to build locally
informed proposals for improved exploration of the target space. However,
their method requires explicit knowledge of the underlying properties of the
target space which is encoded in the sampler which can be problematic.

##  5 Experiments

We now present experimental results for our SurVAE Flow augmented MCMC on
problems covering a range of discrete models applied in statistics, physics,
and machine learning. These include a synthetic example of a discretized
Gaussian Mixture Model, Ising model for denoising corrupted MNIST, quantized
logistic regression on four real world datasets, and Bayesian variable
selection for high-dimensional data. The code for our implementation as well
as the experiments is available at https://github.com/priyankjaini/discFlowMH.

We compare our model to two baselines that include a random walk Metropolis-
Hastings algorithm and Gibbs sampling. We further also compare to discrete HMC
(dHMC) (Nishimura et al.,, 2020) although we use the original implementation
released by the authors which is not parallelizable and implements using the
numpy package (Harris et al.,, 2020) and R. In contrast, we implemented our
method and the other baselines in Pytorch (Paszke et al.,, 2019) and is
parallelizable to use multiple GPUs for fast and efficient computation.

For each experiment, we train our model to learn $\mathbf{T}_{\lambda}$ and
$\mathbf{T}_{\phi}$ for $10,000$ iterations with a batch-size of 128, learning
rate of $10^{-3}$ and optimize using the Adam optimizer. We run 128 parallel
chains of Metropolis-Hastings algorithm for $10^{5}$ steps with thinning of 10
(i.e. saving every tenth sample) with no burn-in since our learned flow
already provides a good initialization for sampling. For fair comparison, we
follow a similar setup for the baselines i.e. we run 128 parallel chains for
each method for $10^{5}$ steps with thinning of 10. We use $10^{5}$ burn-in
steps for the baselines in all the experiments. Additionally, for denoising
binarized MNIST using the Ising model, we also experiment with different
number of burn-in steps which are detailed in the results. For dHMC (Nishimura
et al.,, 2020), we used the publicly available code released by the authors.

We compare the efficiencies for all these models by reporting both the mean
effective sample size (ESS)222We chose the mean ESS as a metric because for
our experiments the minimum ESS across features tended to be close to one due
to variables being discrete. across dimensions, mean ESS across dimensions per
minute, and the accuracy or (unnormalized) log-probability of the samples
generated by each sampler for the corresponding downstream task. For our
method (labeled as Flow + MH), we include the training time of the flow i.e.
$(\mathbf{T}_{\lambda},\mathbf{T}_{\phi})$ when evaluating ESS per min to
ensure a fair comparison. We compute the ESS for 16 chains at a time,
resulting in 8 ESS estimates and report mean $\pm$ standard error across these
8 estimates. The major aims of our experimental evaluation presented here are
two-fold: Firstly, our aim is to demonstrate the efficacy of our flow
augmented MCMC over broad applications. In particular, we want to demonstrate
that flow-based generative models present an attractive methodology to sample
from a discrete distribution by modelling the distribution in a continuous
space and sampling in a simpler latent space. Thus, in our experiments here we
have restricted ourselves to using the basic Metropolis-Hastings (MH)
algorithm for sampling to demonstrate the advantages of being able to learn a
transport map from the discrete space to a simple latent space. By using a
more sophisticated sampler like HMC (Duane et al.,, 1987; Neal et al.,, 2011),
we could thus get even better results. Secondly, most samplers (and especially
discrete space samplers) are either inefficient and/or are prohibitively slow
for high-dimensional problems. Thus, through our experiments we want
demonstrate that leveraging the advances of flows to learn high-dimensional
distributions, our sampler is more efficient with significantly better results
than the alternative methods for high-dimensional problems.

### Discretized Gaussian Mixture Model

We first illustrate our method on a 2D toy problem. We define a set of target
distributions by discretizing Gaussian mixture models with 5 components using
a discretization level of 6 bits. In Fig. 1 we compare samples from the target
distributions, samples from the approximation learned by the flow and samples
from the MCMC chain. We observe that the flow gives an initial rough
approximation to the target density – thus greatly simplifying the geometry
for the MCMC method. Next, the samples from the MCMC chain are
indistinguishable from true samples from the target distribution. This also
highlights the trade-off we described in Remark 1 in Section 3.

Target

Flow

Flow+MH

Figure 1: Toy problem. Left: Discretized Gaussian mixture target
distributions. Middle: Approximation learned by flow. Right: Samples from MCMC
chain.

### Ising Model for Denoising Binarized MNIST

We illustrate the application of discrete samplers to un-directed graphs in
this section by considering the removal of noise from a binary $28\times 28$
MNIST image. We take an image from the MNIST dataset and binarize it such that
the true un-corrupted image is described by a vector
$\mathbf{y}\in\\{0,1\\}^{784}$. We obtain a corrupted image $\mathbf{x}$ by
taking this unknown noise-free image $\mathbf{y}$ and randomly flipping the
value of each pixel with probability 0.1. Given this corrupted image, our goal
is to recover the original noise-free image.

Table 1: Results for denoising binary MNIST using the Ising model. We report ESS per $10^{4}$ samples, ESS per minute, and the average log-likelihood (unnormalized) for the generated samples. Sampler | ESS | ESS/min | $\log\pi(\boldsymbol{\theta})$
---|---|---|---
Flow + MH (iters = 1k) | 43.45 $\pm$ 5.07 | 0.15 $\pm$ 0.02 | 4025.67 $\pm$ 7.56
Flow + MH (iters = 2k) | 46.48 $\pm$ 7.02 | 0.15 $\pm$ 0.02 | 4039.92 $\pm$ 6.05
Flow + MH (iters = 5k) | 45.95 $\pm$ 5.34 | 0.13 $\pm$ 0.01 | 4054.25 $\pm$ 5.15
Flow + MH (iters = 10k) | 28.52 $\pm$ 4.17 | 0.07 $\pm$ 0.01 | 4056.82 $\pm$ 5.64
Gibbs (burn-in = 10k) | 7.131 $\pm$ 1.42 | 0.011 $\pm$ 0.002 | 4045 $\pm$ 22.79
Gibbs (burn-in = 20k) | 6.09 $\pm$ 0.82 | 0.009 $\pm$ 0.001 | 4050.08 $\pm$ 16.56
Gibbs (burn-in = 50k) | 7.65 $\pm$ 0.88 | 0.012 $\pm$ 0.001 | 4056.37 $\pm$ 9.94
Gibbs (burn-in = 100k) | 7.82 $\pm$ 1.06 | 0.009 $\pm$ 0.001 | 4060.09 $\pm$ 7.56
discrete MH (burn-in = 10k) | 26.24 $\pm$ 1.84 | 0.07 $\pm$ 0.01 | 3989.26 $\pm$99.77
discrete MH (burn-in = 20k) | 12.01 $\pm$ 2.39 | 0.03 $\pm$ 0.01 | 4019.20 $\pm$ 44.22
discrete MH (burn-in = 50k) | 1.90 $\pm$ 0.17 | 0.004 $\pm$ 0.001 | 4038.96 $\pm$17.95
discrete MH (burn-in = 100k) | 2.40 $\pm$ 0.87 | 0.005 $\pm$ 0.002 | 4044.16 $\pm$ 16.03
dHMC | 17.83 $\pm$ 1.49 | 0.022 $\pm$ 0.001 | 4038.68 $\pm$ 23.55

Following (Bishop,, 2006, Section 8.3), we solve this by formulating it as an
Ising model and using Monte carlo samplers to sample $\boldsymbol{\theta}$
with energy function given by:

$\displaystyle
E(\boldsymbol{\theta},\mathbf{x})=-\beta\cdot\sum_{i,j}\boldsymbol{\theta}_{i}\boldsymbol{\theta}_{j}-\eta\cdot\sum_{i}\boldsymbol{\theta}_{i}\mathbf{x}_{i}$

where $\boldsymbol{\theta}_{i}$ and $\mathbf{x}_{i}$ are the $i^{th}$ pixel in
the sample under consideration $\boldsymbol{\theta}$ and corrupted image
$\mathbf{x}$ respectively. We train the flow for
$\mathsf{iters}=\\{10^{3},2\times 10^{3},5\times 10^{3},10\times 10^{3}\\}$
and run the baselines with $\mathsf{burn-in}=\\{10^{4},~{}2\times
10^{4},~{}5\times 10^{4},~{}10\times 10^{4}\\}$. We report the results in
Table 1 and the samples in Figure 2-3. We The results evidently show that Flow
augmented MCMC significantly outperforms other samplers on both ESS and the
underlying downstream task demonstrating its ability to disentangle complex
correlations present in high-dimensional problems by mapping it to a latent
space whereas other coordinate-wise samplers are not able to handle this
efficiently.

True:

Corrupted:

Flow+MH:

Discrete MCMC:

Gibbs:

Figure 2: Ising model. Left: True MNIST digit and observed corrupted MNIST
digit. Right: Denoised posterior samples from independent MCMC chains.

### Quantized Logistic Regression

Next, we consider the task of logistic regression where the learnable
parameters and biases denoted by $\boldsymbol{\theta}$ taking discrete values.
This problem is particularly applicable for training quantized neural networks
where the weights take discrete values. However, here we restrict ourselves to
a simpler version of quantized logistic regression on four real-world datasets
from the UCI repository that include Iris, Breast Cancer, Wine, and Digits.
For each dataset, we run 5-fold cross-validation – with results averaged
across the 5 folds. In each fold, we train the flow for $10^{4}$ iterations
and run discrete MCMC with a burn-in of $10^{5}$ steps. We consider parameters
quantized to 4 bits and report the results for these in Table 2. We notice
that our method outperforms other methods on ESS but suffers in ESS/min due to
the time it takes to train the flow.

Table 2: Results for quantized logisitic regression on four datasets. We report the ESS per $10^{4}$ samples, ESS per minute, maximum accuracy, and the average log-likelihood for samples generated by each model. Dataset | Parameters | Sampler | ESS | ESS/min | accuracy | $\log\pi(\boldsymbol{\theta})$
---|---|---|---|---|---|---
Iris | 15 | Flow + MH | 923.21 $\pm$ 37.29 | 58.64 $\pm$ 2.37 | 97.4 | -8.55 $\pm$ 1.77
discrete MCMC | 108.53 $\pm$ 1.44 | 26.14 $\pm$ 0.35 | 97.4 | -8.54 $\pm$ 1.77
dHMC | 361.07 $\pm$ 16.97 | 32.85 $\pm$ 0.64 | 97.4 | -8.40 $\pm$ 1.89
Wine | 42 | Flow + MH | 146.25 $\pm$ 3.51 | 6.59 $\pm$ 0.16 | 97.8 | -2.15 $\pm$ 1.63
discrete MCMC | 65.78 $\pm$ 0.46 | 15.05 $\pm$ 0.10 | 97.8 | -2.15 $\pm$ 1.63
dHMC | 89.12 $\pm$ 11.57 | 4.53 $\pm$ 0.98 | 97.7 | -3.54 $\pm$ 0.48
Breast Cancer | 62 | Flow + MH | 15.77 $\pm$ 0.49 | 0.56 $\pm$ 0.02 | 96.8 | -19.85 $\pm$ 3.40
discrete MCMC | 10.46 $\pm$ 0.12 | 2.29 $\pm$ 0.03 | 96.8 | -20.14 $\pm$ 3.46
dHMC | 9.54 $\pm$ 0.28 | 1.38 $\pm$ 0.03 | 96.8 | -21.19 $\pm$ 2.86

### Bayesian Variable Selection

Table 3: Results for Bayesian variable selection on synthetic datasets of 100, 200, and 400 dimensions with only 10, 20, and 40 informative features, respectively. We report ESS per $10^{4}$ samples, ESS per minute, and the average log-likelihood for samples generated by each model. Setting | Sampler | ESS | ESS/min | $\log\pi(\boldsymbol{\theta})$
---|---|---|---|---
$\mathsf{features}=100$, $\mathsf{features}_{\mathsf{informative}}=10$ | Flow + MH | 436.48 $\pm$ 57.82 | 4.78 $\pm$ 0.63 | -1054.97 $\pm$ 15.89
Gibbs | 23.19 $\pm$ 22.20 | 0.12 $\pm$ 0.12 | -1235.02 $\pm$ 13.11
$\mathsf{features}=200$, $\mathsf{features}_{\mathsf{informative}}=20$ | Flow + MH | 213.36 $\pm$ 29.37 | 0.89 $\pm$ 0.12 | -2231.24 $\pm$ 8.28
Gibbs | 8.78 $\pm$ 6.12 | 0.013 $\pm$ 0.009 | -2485.91 $\pm$ 0.14
$\mathsf{features}=400$, $\mathsf{features}_{\mathsf{informative}}=40$ | Flow + MH | 141.45 $\pm$ 17.97 | 0.15 $\pm$ 0.02 | -4627.58 $\pm$ 10.86
Gibbs | 1.46 $\pm$ 0.23 | 0.001 $\pm$ $10^{-4}$ | -5306.61 $\pm$ 0.20

Here, we consider the problem of probabilistic selection of features in
regression where the hierarchical framework allows for complex interactions
making posterior exploration extremely difficult. Following Schäfer and
Chopin, (2013), we consider a hierarchical Bayesian model in which the target
vector $\mathbf{y}\in\mathbb{Z}^{d}$ in linear regression is observed through:

$\displaystyle\mathbf{y}|\beta,\boldsymbol{\theta},\sigma^{2},\mathbf{x}\sim\mathcal{N}(\mathbf{x}\cdot\boldsymbol{\theta}\cdot\beta,\sigma^{2}I)$

with features $\mathbf{x}\in\mathds{R}^{d\times k}$, parameters
$\beta\in\mathds{R}^{k}$ and the binary vector $\boldsymbol{\theta}$ that
indicates the features that are included in the model. To achieve a closed
form for the marginal posterior that is independent of $\beta$ and
$\sigma^{2}$, we make the following choices following Power and Goldman,
(2019):
$\pi(\boldsymbol{\theta}|\mathbf{x},\mathbf{y})=\int\pi(\boldsymbol{\theta}|\mathbf{x},\mathbf{y},\beta,\sigma^{2})\mathrm{d}\beta~{}\mathrm{d}\sigma^{2}$
with conjugates
$p(\beta|\sigma^{2},\boldsymbol{\theta})=\mathcal{N}(0,\nu^{2}\sigma^{2}I_{n}\boldsymbol{\theta})$,
$p(\sigma)=\text{I}\Gamma(\frac{2}{2},\frac{\alpha w}{2})$ and,
$p(\boldsymbol{\theta})=\mathsf{Uniform}(\\{0,1\\}^{k})$ where
$\text{I}\Gamma$ is the inverse-gamma distribution and we set $\nu,w,\alpha$
as in George and McCulloch, (1997). For this setting, we create high-
dimensional synthetic datasets rather than the low-dimensional datasets
usually used to demonstrate the efficacy of our method. We present the results
in Table 3.

##  6 Conclusion

In this paper, we presented a flow based Monte Carlo sampler for sampling in
combinatorial spaces. Our method learns a deterministic transport map from a
discrete space to a simple continuous latent space where it is efficient to
sample. Thereby, we sample from the discrete space by generating samples in
the latent space and using the transport map to obtain samples in the discrete
space. By learning a map to a simple latent space (like standard Gaussian),
our method is particularly suited for high-dimensional domains where
alternative samplers are not efficient due to the presence of more complex
correlations. This is also reflected in our implementation which is faster and
efficient as demonstrated by our results on a suite of experiments. In the
future, it will be interesting to devise learning strategies for the transport
map that explicitly pushes the latent space to have certain desirable
properties for efficient sampling. Another direction could be to extend the
framework of SurVAE Flow layers to incorporate underlying symmetries and
invariance in the target domain.

## References

  * Afshar and Domke, (2015) Afshar, H. M. and Domke, J. (2015).  Reflection, Refraction, and Hamiltonian Monte Carlo.  In Advances in neural information processing systems, pages 3007–3015.
  * Albergo et al., (2019) Albergo, M., Kanwar, G., and Shanahan, P. (2019).  Flow-based Generative Models for Markov Chain Monte Carlo in Lattice Field Theory.  Physical Review D, 100(3):034515.
  * Bishop, (2006) Bishop, C. M. (2006).  Pattern Recognition and Machine Learning.  springer.
  * Brooks et al., (2011) Brooks, S., Gelman, A., Jones, G., and Meng, X.-L. (2011).  Handbook of Markov Chain Monte Carlo.  CRC press.
  * Dinh et al., (2017) Dinh, V., Bilge, A., Zhang, C., and Matsen IV, F. A. (2017).  Probabilistic Path Hamiltonian Monte Carlo.  arXiv preprint arXiv:1702.07814.
  * Duane et al., (1987) Duane, S., Kennedy, A. D., Pendleton, B. J., and Roweth, D. (1987).  Hybrid Monte Carlo.  Physics letters B, 195(2):216–222.
  * George and McCulloch, (1997) George, E. I. and McCulloch, R. E. (1997).  Approaches for Bayesian Variable Selection.  Statistica sinica, pages 339–373.
  * Harris et al., (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., et al. (2020).  Array programming with NumPy.  Nature, 585(7825):357–362.
  * Ho et al., (2019) Ho, J., Chen, X., Srinivas, A., Duan, Y., and Abbeel, P. (2019).  Flow++: Improving Flow-based Generative Models with Variational Dequantization and Architecture Design.  ICML, 2019.
  * Hoffman et al., (2019) Hoffman, M., Sountsov, P., Dillon, J. V., Langmore, I., Tran, D., and Vasudevan, S. (2019).  Neutralizing Bad Geometry in Hamiltonian Monte Carlo using Neural Transport.  ICML, 2019.
  * Hubbard, (1959) Hubbard, J. (1959).  Calculation of Partition Functions.  Physical Review Letters, 3(2):77.
  * Marzouk et al., (2016) Marzouk, Y., Moselhy, T., Parno, M., and Spantini, A. (2016).  Sampling via Measure Transport: An Introduction.  In Ghanem, R., Higdon, D., and Owhadi, H., editors, Handbook of Uncertainty Quantification, pages 1–41. Springer.
  * Neal et al., (2011) Neal, R. M. et al. (2011).  MCMC using Hamiltonian Dynamics.  Handbook of markov chain monte carlo, 2(11):2.
  * Nielsen et al., (2020) Nielsen, D., Jaini, P., Hoogeboom, E., Winther, O., and Welling, M. (2020).  SurVAE Flows: Surjections to Bridge the Gap between VAEs and Flows.  arXiv preprint arXiv:2007.02731.
  * Nishimura et al., (2020) Nishimura, A., Dunson, D. B., and Lu, J. (2020).  Discontinuous Hamiltonian Monte Carlo for Discrete Parameters and Discontinuous Likelihoods.  Biometrika, 107(2):365–380.
  * Pakman and Paninski, (2013) Pakman, A. and Paninski, L. (2013).  Auxiliary-Variable Exact Hamiltonian Monte Carlo Samplers for Binary Distributions.  In Advances in neural information processing systems, pages 2490–2498.
  * Pardalos, (1996) Pardalos, P. M. (1996).  Continuous Approaches to Discrete Optimization Problems.  In Nonlinear optimization and applications, pages 313–325. Springer.
  * Parno and Marzouk, (2018) Parno, M. and Marzouk, Y. (2018).  Transport Map Accelerated Markov Chain Monte Carlo.  SIAM/ASA Journal on Uncertainty Quantification, 6(2):645–682.
  * Paszke et al., (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S. (2019).  Pytorch: An imperative style, high-performance deep learning library.  In Wallach, H., Larochelle, H., Beygelzimer, A., dAlché Buc, F., Fox, E., and Garnett, R., editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc.
  * Peherstorfer and Marzouk, (2018) Peherstorfer, B. and Marzouk, Y. (2018).  A Transport-based Multifidelity Pre-Conditioner for Markov Chain Monte Carlo.
  * Power and Goldman, (2019) Power, S. and Goldman, J. V. (2019).  Accelerated Sampling on Discrete Spaces with Non-Reversible Markov Processes.  arXiv preprint arXiv:1912.04681.
  * Rezende and Mohamed, (2015) Rezende, D. J. and Mohamed, S. (2015).  Variational inference with normalizing flows.  In ICML.
  * Rudin, (1987) Rudin, W. (1987).  Real and Complex Analysis.  McGraw-Hill, 3rd edition.
  * Schäfer and Chopin, (2013) Schäfer, C. and Chopin, N. (2013).  Sequential Monte Carlo on Large Binary Sampling Spaces.  Statistics and Computing, 23(2):163–184.
  * Tabak and Turner, (2013) Tabak, E. G. and Turner, C. V. (2013).  A family of nonparametric density estimation algorithms.  Communications on Pure and Applied Mathematics, 66(2):145–164.
  * Tabak and Vanden-Eijnden, (2010) Tabak, E. G. and Vanden-Eijnden, E. (2010).  Density estimation by dual ascent of the log-likelihood.  Communications in Mathematical Sciences, 8(1):217–233.
  * Theis et al., (2015) Theis, L., Oord, A. v. d., and Bethge, M. (2015).  A Note on the Evaluation of Generative Models.  arXiv preprint arXiv:1511.01844.
  * Titsias and Yau, (2017) Titsias, M. K. and Yau, C. (2017).  The Hamming Ball Sampler.  Journal of the American Statistical Association, 112(520):1598–1611.
  * Uria et al., (2013) Uria, B., Murray, I., and Larochelle, H. (2013).  RNADE: The Real-Valued Neural Autoregressive Density Estimator.  In Advances in Neural Information Processing Systems, pages 2175–2183.
  * Villani, (2008) Villani, C. (2008).  Optimal Transport: Old and New.  Springer.
  * Zanella, (2020) Zanella, G. (2020).  Informed Proposals for Local MCMC in Discrete Spaces.  Journal of the American Statistical Association, 115(530):852–865.
  * Zhang et al., (2012) Zhang, Y., Ghahramani, Z., Storkey, A. J., and Sutton, C. A. (2012).  Continuous Relaxations for Discrete Hamiltonian Monte Carlo.  In Advances in Neural Information Processing Systems, pages 3194–3202.

##  Appendix A Additional Results

We here provide additional results for the Ising model experiment in Sec. 5.
The results are shown in Fig. 3.

(a) True MNIST digit and observed corrupted MNIST digit.

(b) Flow+MH.

(c) Discrete MCMC.

(d) Gibbs.

Figure 3: Original and corrupted MNIST digit together with denoising posterior
samples from different MCMC methods. Samples are from independent MCMC chains.


# Infall and outflow detections in a massive core JCMT 18354-0649S

Tie Liu11affiliation: Department of Astronomy, Peking University, 100871,
Beijing China<EMAIL_ADDRESS><EMAIL_ADDRESS>, Yuefang Wu11affiliation:
Department of Astronomy, Peking University, 100871, Beijing China;
<EMAIL_ADDRESS><EMAIL_ADDRESS>, Qizhou Zhang22affiliation: Harvard-
Smithsonian Center for Astronomy, 60 Garden St., Cambridge, MA 02138, USA ,
Zhiyuan Ren11affiliation: Department of Astronomy, Peking University, 100871,
Beijing China<EMAIL_ADDRESS><EMAIL_ADDRESS>, Xin Guan11affiliation:
Department of Astronomy, Peking University, 100871, Beijing China;
<EMAIL_ADDRESS><EMAIL_ADDRESS>33affiliation: Current affiliation: I.
Physikal. Institut, Universität zu Köln, Zülpicher St.77,D-50937 Köln, Germany
and Ming Zhu44affiliation: National Astronomical Observatories, Chinese
Academy of Sciences, Beijing, 100012

###### Abstract

We present a high-resolution study of a massive dense core JCMT 18354-0649S
with the Submillimeter Array. The core is mapped with continuum emission at
1.3 mm, and molecular lines including CH3OH ($5_{23}$–$4_{13}$) and HCN (3–2).
The dust core detected in the compact configuration has a mass of
$47~{}M_{\odot}$ and a diameter of $2\arcsec$ (0.06 pc), which is further
resolved into three condensations with a total mass of $42~{}M_{\odot}$ under
higher spatial resolution. The HCN (3–2) line exhibits asymmetric profile
consistent with infall signature. The infall rate is estimated to be
$2.0\times 10^{-3}~{}M_{\odot}\cdot$yr-1. The high velocity HCN (3-2) line
wings present an outflow with three lobes. Their total mass is
$12~{}M_{\odot}$ and total momentum is $121~{}M_{\odot}\cdot$km s-1,
respectively. Analysis shows that the N-bearing molecules especially HCN can
trace both inflow and outflow.

Massive core:pre-main sequence-ISM: molecular-ISM: kinematics and dynamics-
ISM: jets and outflows-stars: formation

††slugcomment: Accepted by APJ

##  1 Introduction

Studies of high-mass star formation have received much attention during recent
years. One of the main questions is whether massive stars form through an
accretion-disk-outflow process, similar to low-mass counterparts (Shu, Adams,
& Lizano 1987), or via collision-coalescence (Wolfire & Cassinelli 1987;
Bonnell, Bate, & Zinnecker. 1998). Studying the characteristics of massive
cores at the early stages is critical for understanding their formation
process. High-Mass Protostellar Objects (HMPOs) are precursors of UC Hii
regions, and represent an essential phase in high-mass star formation
(Churchwell 2002). HMPOs often have strong dust emission and high bolometric
luminosity. But their radio emission is weak or non detectable at a level of
approximately 1 mJy (Molinari et al. 1996, 2000; Sridharan et al. 2002;
Beuther, Schilke, & Menten. 2002; Wu et al. 2006).. Their natal clouds have
not been affected significantly by the star forming process. Thus, they
present the information about the early kinematic processes of high mass star
formation.

The dense core JCMT 18354-0649S was first detected in an ammonia survey of
high-mass star forming regions with Max-Planck-Institut für Radioastronomie
(MPIfR) 100 m telescope at Effelsberg (Wu et al. 2006), and was later
confirmed by the observation with the Submillimeter Common-User Bolometric
Array (SCUBA) of James Clerk Maxwell telescope (JCMT) (Wu et al. 2005).
Another SCUBA core which harbors a UC H ii region G35.4NW is located about
$1\arcmin$ north of JCMT 18354-0649S. The kinetic distance of the two SCUBA
cores is 5.7 or 9.6 kpc (Wu et al. 2005), and 5.7 kpc was adopted in this
paper. Core JCMT 18354-0649S has no counterpart in radio continuum. Multiple
lines towards this core including HCN (3–2), H13CO+ (3–2), and C17O (2–1)
reveal typical ”blue profile” (Wu et al. 2005), indicating that the core is
undergoing gravitational collapse (Keto, Ho & Haschick. 1988; Zhou et al.
1993; Zhang, Ho, & Ohashi. 1998; Wu & Evans. 2003; Wu et al. 2005; Fuller,
Williams, & Sridharan. 2005; Wyrowski 2006; Wu et al. 2007; Birkmann et al.
2007; Klaassen & Wilson 2007; Sun & Gao 2008; Velusamy et al. 2008). The core
is also associated with a near-infrared point source, corresponding to a star
of 6-11 $M_{\odot}$ (Zhu et al.2010, submitted). Carolan et al. (2009)
observed sixteen different molecular line transitions including CO, HCN, HCO+
and their isotopes in this region, and modeled the source with a chemically
depleted rotating envelope collapsing onto a central protostellar source which
has evolved sufficiently to generate a molecular outflow. All the evidence
suggests that JCMT 18354-0649S is forming high mass protostellar
object(s)(HMPO). However, single-dish observations with a resolution of
$15\arcsec-40\arcsec$ can not reveal detailed kinematics in the core at a
distance as large as 5.7 kpc. In this paper we report the results of a high
resolution study with the Submillimeter Array (SMA111Submillimeter Array 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.) in order to probe the
details of the core and its kinematic features. The observation and initial
results are presented in sections 2 and 3. Properties of infall and outflow
motions are discussed further in section 4, and a brief summary is given in
section 5.

##  2 Observations

The observations of JCMT 18354-0649S was carried out with the SMA in July 2005
with seven antennas in its compact configuration and in September 2005 with
six antennas in its extended configuration. The 345 GHz receivers were tuned
to 265 GHz for the lower sideband (LSB) and 275 GHz for the upper sideband
(USB). The frequency spacing across the spectral band is 0.8125 MHz or $\sim$1
km s-1 for both configurations. The phase reference center of both
observations was R.A.(J2000) = 18h38m08.10s and DEC.(J2000) = -$6\arcdeg
46\arcmin 52.17\arcsec$.

In the observations with the compact configuration, Jupiter, Uranus and QSO
3c454.3 were observed for antenna-based bandpass correction. An amplitude
offset was found on some baselines and the baseline-based errors in bandpass
were further corrected using the point source QSO 3c454.3. QSOs 1741-038 and
1908-201 were employed for antenna-based gain correction. Uranus was observed
for flux-density calibration. The synthesized beam size is $3.76\arcsec\times
2.72\arcsec$ (PA=-54$\arcdeg$).

For the extended configuration, QSO 3c454.3 was used as bandpass calibrator,
QSOs 1741-038 and 1908-201 as gain calibrators, and Uranus as a flux
calibrator, respectively. The synthesized beam size is about $1\arcsec$.

Miriad222http://carma.astro.umd.edu/miriad was employed for calibration and
imaging. The 1.3 mm continuum data were acquired by averaging all the line-
free channels over both the 2 GHz of upper and lower spectral bands. MIRIAD
task ”selfcal” was employed to perform self-calibration on the continuum data.
Since the dust emission is weak, self-calibration with phase only was
performed. The gain solutions from the self-calibration were applied to the
line data.

The continuum data combined from both configurations yield a synthesized beam
of $1.63\arcsec\times 1.28\arcsec$ (PA=-81.4$\arcdeg$), and 1 $\sigma$ rms of
2.5 mJy in the naturally weighted maps. HCN (3-2) and CH3OH
($5_{23}$–$4_{13}$) were detected in the compact configuration.

The shortest baseline in compact configuration observations is 16.5 m,
corresponding to a spatial scale of $20\arcsec$. Spatial structures more
extended than this limit, such as HCN maps close to the cloud velocity, would
be filtered out. HCN (3–2) data in JCMT archive (Wu et al. 2005; Carolan et
al. 2009) were used to recover the missing flux. The JCMT archive data were
reduced using the KAPPA and GAIA packages in the STARLINK suite. The JCMT beam
size for HCN (3–2) was $18.3\arcsec$, and the main beam efficiency was 0.69.
The combination of the SMA compact and JCMT HCN (3–2) data was done using the
task ”immerge” in MIRIAD.

##  3 Results

###  3.1 Dust core

The 1.3 mm continuum images are shown in Fig.1. The left panel is obtained in
the compact array and the right panel from the combined data of the compact
and extended configurations. An elongated core is revealed with the 1.3 mm
continuum emission observed with the compact array, and is further resolved
into three condensations by the continuum emission using the combined data
from both configurations. The three condensations are named as MM1, MM2 and
MM3. The peak position of MM1 is R.A.(J2000)=$18^{\rm h}38^{\rm m}08.1^{\rm
s}$, DEC.(J2000)=$-6\arcdeg 46\arcmin 52.98\arcsec$. MM2 peaks at
R.A.(J2000)=$18^{\rm h}38^{\rm m}07.9^{\rm s}$, DEC.(J2000)=$-6\arcdeg
46\arcmin 51.36\arcsec$, and MM3 peaks at R.A.(J2000)=$18^{\rm h}38^{\rm
m}08.05^{\rm s}$, DEC.(J2000)=$-6\arcdeg 46\arcmin 51.36\arcsec$.

From the fit of an elliptical Gaussian, the core revealed by the compact array
is found to be elongated from south-east to north-west. It has an average FWHM
diameter of 0.06 pc ($\sim 2\arcsec$) at a distance of 5.7 kpc, smaller than
the beam size of the compact configuration. The total integrated flux is 0.47
Jy. The total dust and gas mass can be obtained with the formula
$M=S_{\nu}D^{2}/\kappa_{\nu}B_{\nu}(T_{d})$, where $S_{\nu}$ is the flux at
1.3 mm, D is the distance, and $B_{\nu}(T_{d})$ is the Planck function. We
adopt a dust opacity $k_{1330}$=$1.4\times 10^{-2}$ cm2g-1 at 1.3 mm
calculated from Ossenkopf & Henning (1994) with a dust opacity index
$\beta=2$. Here the ratio of gas to dust is taken as 100. Molecular line CH3OH
($5_{23}$–$4_{13}$) is detected, and its emission peak coincides with the dust
core very well (see Sec.3.2). The upper energy level of the CH3OH
($5_{23}$–$4_{13}$) line is 57 K above the ground, indicating a relatively
warm conditions. Assuming $T_{d}$=57 K, a total dust and gas mass of
$47~{}M_{\odot}$ is derived. A beam-average gas/dust density amounts to
$2.0\times 10^{6}$ cm-3, which is larger than $1.1\times 10^{6}$ cm-3 obtained
from single-dish telescope (Wu et al. 2005).

The three condensations (MM1, MM2 and MM3) have total integrated flux of 0.42
Jy, leading to a total mass of $42~{}M_{\odot}$ (assuming Td = 57 K as above).
MM1 has a diameter of 0.02 pc, and a mass of $30~{}M_{\odot}$. MM1 is
centrally concentrated and compact, while MM2 and MM3 are much more diffuse
and extended.

With UKIRT (United Kingdom Infrared Telescope) Zhu et al. (2010, in
preparation) detected three near-infrared sources IRS1a, IRS1b, IRS1c in H, K
and L bands. Their positions are marked with crosses in Fig.1. IRS1a ($18^{\rm
h}38^{\rm m}08.135^{\rm s}$, $-6\arcdeg 46\arcmin 51.57\arcsec$) lies about
$1.6\arcsec$ north-east of MM1. IRS1b ($18^{\rm h}38^{\rm m}08.026^{\rm s}$,
$-6\arcdeg 46\arcmin 56.24\arcsec$) and IRS1c ($18^{\rm h}38^{\rm
m}07.929^{\rm s}$, $-6\arcdeg 46\arcmin 55.34\arcsec$) are about $3\arcsec$
southwest from MM1 and are much fainter.

###  3.2 Gas core

The molecular line CH3OH ($5_{23}$–$4_{13}$) is detected in the compact
configuration. Fig.2 presents its spectrum at three positions and the
integrated emission overlaid on the 1.3 mm continuum image. The central
velocity of the CH3OH ($5_{23}$–$4_{13}$) spectra is 96.7 km s-1, which is
taken as the systemic velocity of the core. The central velocity of CH3OH
($5_{23}$–$4_{13}$) does not shift at different positions (see Fig.2), which
should exclude rotation at the core. The P-V diagram of CH3OH
($5_{23}$–$4_{13}$) is shown in Fig.3, indicating a compact gas core without
rotation. The emission center of CH3OH ($5_{23}$–$4_{13}$)
(R.A.(J2000)=$18^{\rm h}38^{\rm m}08.092^{\rm s}$, DEC.(J2000)=$-6\arcdeg
46\arcmin 52.318\arcsec$) coincides with MM1 very well. While there are no
CH3OH components corresponding with MM2 and MM3. The deconvolved size of the
gas core revealed by CH3OH ($5_{23}$–$4_{13}$) is $3.78\arcsec\times
2.76\arcsec$ (PA=-29$\arcdeg$), comparable to the synthesized beam size of the
compact array.

The HCN (3–2) (265.886GHz) spectra obtained from the SMA compact configuration
and from the data combined from both compact and extended configurations are
presented in the left panel of Fig.4. Both of the two spectra are averaged
over a region of $5\arcsec\times 5\arcsec$, which show a redshifted absorption
dip and broad wings. The line profiles observed with the SMA and JCMT, as well
as a combination of the two are presented in the right panel of Fig.4. All the
spectra in the right panel of Fig.4 are convolved with the JCMT beam
($18.3\arcsec$) for comparing. One can see that the SMA compact array
observations recover less than 10$\%$ of JCMT flux around the systematic
velocity, but recover more than 30$\%$ flux at the wings. The combination of
the SMA and JCMT data recovers more than 70$\%$ of the JCMT flux at all the
velocity channels.

###  3.3 Kinematic signatures of lines

####  3.3.1 Infall motion

The left panel of Fig.4 shows the most prominent feature (”blue profile”) of
the HCN (3–2) line at the core. The absorption gap is more than 8 km s-1 wide,
ranging from 93 to 101 km s-1. Fig.5 presents the channel maps of the HCN
(3–2) emission from 80 km s-1 to 109 km s-1 constructed from the combined
data, which is convolved with the beam of the SMA compact configuration. The
absorption is obvious in the velocity range (95,99) km s-1. The absorption dip
is also clearly seen in the P-V diagrams (see Fig.6), which is much deeper
than that revealed by the single-dish observation (Wu et al. 2005).

From the left panel of Fig.4, it is clearly to see that the centeral velocity
of the absorption dip (98 km s-1) is redshifted from the systematic velocity
(96.7 km s-1) by 1.3 km s-1. Such a blue asymmetric line profile where the
blue emission peak is at a higher intensity than the red one is a collapse
signature of molecular cores (Zhou et al. 1993). The spectra constructed from
JCMT data and the combined data (the right panel of Fig.4) also show
significant ”blue profile”, confirming the existence of infall motions.

####  3.3.2 Molecular outflow

Besides the absorption dip, the HCN (3–2) line exhibits remarkable broad wings
extending more than 40 km s-1. High-velocity gas also can be easily identified
in P-V diagrams of the HCN (3–2) emission along the direction of
P.A.=15$\arcdeg$ and P.A.=90$\arcdeg$ as shown in Fig.6. The HCN (3–2)
emission obtained from SMA compact configuration is integrated from 80 to 87
km s-1 for the blue lobe and from 103 to 109 km s-1 for the red lobe,
respectively. The contour map of the integrated flux are shown in Fig.7. As in
the channel maps (Fig.5), we can see several clumps in each lobe in the
integrated map. The integrated HCN (3–2) emission seems to comprise an
S-shaped structure from north-east to south. Another jet-like structure
extended more than $10\arcsec$ is also seen at the west of the continuum
emission center. IRS1a seems to be the driving source of the outflow.

The southern redshifted lobe (S-lobe) comprises two clumps named ”Clump1” and
”Clump2”. In the north-east blueshifted lobe (NE-lobe), two clumps are also
found and named ”Clump3” and ”Clump4”. These clumps are distributed along the
direction of the outflow and likely to be outward gas knots. They are probably
not physically related with other stellar sources except the driving source
though Clump2 is close to IRS1b and IRS1c.

##  4 Discussion

###  4.1 Infall motion

Although the HCN emission is extended over a region larger than the compact
configuration beam, the infall region is still difficult to confine due to the
contamination of the outflow. Since the size of the gas core traced by CH3OH
($5_{23}$–$4_{13}$) is comparable to the compact configuration beam size, the
beam size of the compact configuration was taken as the radius ($R_{in}$) of
the infall region (Wu et al. 2009). The kinematic mass infall rate can be
calculated using dM/dt=$4{\pi}n\mu_{G}m_{H_{2}}R_{in}^{2}V_{in}$, where
$V_{in}$, $\mu_{G}=1.36$, $m_{H_{2}}$ and n=$2.0\times 10^{6}$ cm-3 are the
infall velocity, the mean molecular weight, the H2 mass, and the beam-average
gas/dust density, respectively. The infall velocity $V_{in}$ is 1.3 km s-1 by
comparing the systemic velocity (96.7 km s-1) and the velocity of the
redshifted absorbing dip (98 km s-1) in the HCN (3-2) spectrum (Welch et al.
1987), leading to a kinematic mass infall rate of 2.0$\times
10^{-3}~{}M_{\odot}\cdot$yr-1. In core G10.6-0.4, the redshifted NH3 indicates
large infall velocity $5.0\pm 1.7$ km s-1 at about 0.05 pc, and a mass infall
rate as high as 5$\times 10^{-3}~{}M_{\odot}\cdot$yr-1 (Keto, Ho & Haschick.
1987). Also with NH3 inverse lines, Zhang & Ho (1997) obtained a high infall
velocity $\sim 3.5$ km s-1 within a region smaller than 0.02 pc towards core
W51e2. Large infall velocities ($>1.5$ km s-1) and mass infall rates
($>1\times 10^{-3}~{}M_{\odot}\cdot$yr-1) were also detected towards G10.47
and G34.26 with HCO+ (4–3) line (Klaassen & Wilson 2007). In core G19.61+0.23,
an infall velocity of 2.5 km s-1 and a mass infall rate as high as $6.1\times
10^{-3}~{}M_{\odot}\cdot$yr-1 were derived (Wu et al. 2009). It seems high
mass infall rate is required by high-mass star formation. The results of the
core JCMT 18354-0649S are comparable with those of the above sources. For
comparison, the $V_{in}$ from pure free-infall assumption is also derived with
the formula $V_{in}^{2}~{}=~{}2GM/R_{in}$. The pure free-infall velocity
inferred is 2.9 km s-1, larger than the infall velocity obtained from the
spectrum.

Wu et al. (2005) obtained a small infall velocity ($\sim 0.3~{}km~{}s^{-1}$)
at a radius of 4$\arcsec$. The absorption dip of HCN (3–2) line seen by the
SMA is much deeper and broader than that observed by JCMT. However, the
kinematic mass infall rate $\dot{M}_{in}$ ($2.0\times
10^{-3}~{}M_{\odot}\cdot$ yr-1) obtained here is well coincident with that
obtained with JCMT (Wu et al. 2005), $3.4\times 10^{-3}~{}M_{\odot}\cdot$
yr-1.

###  4.2 Properties of HCN (3-2) outflow

The column density of HCN at each velocity channel in each outflow lobe can be
obtained through (Garden et al. 1991):

$N_{HCN}(v)=\frac{3k}{8\pi^{3}B\mu^{2}}\frac{exp[hBJ(J+1)/kT_{ex}]}{(J+1)}\frac{(T_{ex}+hB/3k)}{1-exp(-h\nu/kT_{ex})}\int\tau_{v}dv$
(1)

Where $v$ is the central velocity of the channel relative to the systemic
velocity, the rotational constant B=44.315976 GHz and permanent dipole moment
$\mu=3$ debye for HCN, the velocity channel width is smoothed to be
$1~{}km~{}s^{-1}$. Assuming HCN emission in the line wings to be optically
thin and excitation temperature of Tex=30 K (Wu et al. 2004), the optical
depth $\tau_{v}$ can be derived with the equation:

$\tau_{v}=\frac{kT_{r}(v)}{h\nu}(\frac{1}{exp(h\nu/kT_{ex})-1}-\frac{1}{exp(h\nu/kT_{bg})-1})^{-1}$
(2)

where $T_{r}(v)$ is the excess brightness temperature of HCN(3–2) emission at
$v$. Adopting $X_{HCN}=[HCN]/[H_{2}]=1\times 10^{-10}$ (Carolan et al. 2009),
the mass of each lobe at $v$ can be calculated with:

$M(v)=X_{HCN}^{-1}\mu_{G}m_{H_{2}}D^{2}{\int}N_{HCN}(v)d{\Omega}$ (3)

where D, $\Omega$ are the cloud distance and the solid angle. Thus the total
mass of each lobe is given by $M=\sum$$M(v)$, the total momentum by
$P=\sum$$M(v)v$, and the energy by $E={\frac{1}{2}}\sum$$M(v)v^{2}$. The
dynamical timescale $t_{dyn}$ is estimated as $R/V_{max}$, where R is the
outflow extent, and $V_{max}$ is the maximum velocity of the outflow lobe. The
mechanical luminosity L, and the mass-loss rate $\dot{M}$ are calculated as
L=E/t, $\dot{M}=P/(tV_{w}$), where the wind velocity $V_{w}$ is assumed to be
500 km s-1 (Lamers et al. 1995). The derived parameters are listed in Table.1.
The total mass, momentum, energy of the three lobes are 12 $M_{\odot}$, 121
$M_{\odot}\cdot$ km s-1 and $1.3\times 10^{46}$ erg, respectively. The average
dynamical timescale is about $1.6\times 10^{4}$ yr, and the total mass-loss
rate $1.6\times 10^{-5}~{}M_{\odot}\cdot$ yr-1. The outflow is massive with
parameters similar to that of IRAS 05274+3345E and the other outflows detected
towards five massive star formation regions (Zhang et al. 2007b; Klaassen &
Wilson 2008).

The Position-Velocity diagram at the left panel of Fig.6 shows that at low
velocities, the NE-lobe and S-lobe both have compact morphology near the core
center. At higher velocities the southern lobe becomes further away from the
center. As shown in Fig.7, the outflow axis traced by Clump1 and Clump2
differs from that traced by Clump3 and Clump4. The different outflow
orientations in the large and small scales may be attributed to the precession
of the outflow axis (Su et al. 2007). From the right panel of the P-V diagram,
a high-velocity component (V $<$ 85 km s-1) with velocities decreasing with
distance from the protostar, and a second component tracing the low-velocity
material (V $>$ 85 km s-1) extending about $15\arcsec$ along the axis of the
W-lobe are clearly seen. Such convex spur PV structure was also revealed in a
simulation of a pulsed jet driven outflow (Lee et al. 2001).

###  4.3 HCN — tracer of both inflow and outflow motions

HCN is among the most abundant molecular species with a high critical density
larger than $10^{6}$ cm-3 (for HCN (1–0)) (Carolan et al. 2009), and is
believed to trace dense molecular cores. HCN is detected in both low mass
class 0 and I sources (Park, Kim, & Ming. 1999; Yun et al. 1999), and high-
mass hot cores (Boonman et al. 2001).

HCN is thought as a good tracer of inflow motions (Wu & Evans. 2003). The
infall asymmetry in the HCN spectra is found to be more prevalent, and more
prominent than in any other previously used infall tracers such as CS (2–1),
DCO+ (2–1), and N2H+(1–0) during a survey toward 85 starless cores (Sohn et
al. 2007). Among the small group of pre- and protostellar objects in L1251B,
infall signature was also detected in the HCN emission (Lee et al. 2007). HCN
also traces inflow motions very well in massive star-formation regions. Wu and
Evans. (2003) found 12 sources showing ”blue profile” in the HCN lines during
a spectroscopic survey of 28 massive cores with water maser. Besides HCN,
other nitrogen bearing molecules such as N2H+ are also tracers of inflow
motions (Tsamis et al. 2008; Schnee et al. 2007; Crapsi et al. 2005). Recently
inverse P Cygni profile of CN line in hot cores was found (Zapata et al. 2008;
Wu et al. 2009). These results suggest nitrogen bearing molecular species be
good tracers of inflowing motions in star-formation regions.

Outflows traced by HCN are often detected not only in low-mass star-formation
regions but also in massive star-formation regions (Bachiller, Gutiérrez, &
Pérez. 1997; Choi 2001; Su et al. 2007; Zhang et al. 2007a). HCN outflow of
the core JCMT 18454-0649S is another good sample. Additionally, Zhu et al.
(2010 in preparation) and Cyganowski (2008) found excess emission at 4.5
$\micron$ at the position of source IRS1a, which is close to the center of the
NE-lobe and S-lobe. Such excess emission at the 4.5 $\micron$ band could be
shock-excited. In IRAS 20126+4104, HCN emission is also found to be closely
related to the shock-excited near-IR H2 knots and was identified to be
associated with shock wings (Su et al. 2007). The inner clumps (Clump1 in the
S-lobe and Clump3 in the NE-lobe) of the core JCMT 18354-0649S should also be
coincident with shocks. In fact, models have already demonstrated a dramatic
increase of HCN molecules, during the intense interaction between outflow and
ambient gas, or slow shock front (Mitchell 1984; Nejad, Williams, & Charnley.
1990). In this process sulfur and nitrogen react with hydrocarbons to produce
various compounds, wherein HCN abundance gets higher than the rest of the
products (Nejad, Williams, & Charnley. 1990). Thus HCN may trace the outflow
even better than sulfur containing molecules.

##  5 Summary

Both dust continuum at 1.3 mm and CH3OH emission detected with SMA reveal a
compact core in JCMT 18354-0649S. The core observed with the compact
configuration has a mass of $47~{}M_{\odot}$ and an average density of
$2.0\times 10^{6}$ cm-3. With the combination of the compact and extended
configurations, the core is resolved to three condensations with a total mass
of $42~{}M_{\odot}$.

HCN (3-2) spectra exhibit an infall signature in this region. The red shifted
absorption seen in the SMA observation is deeper and broader than that in the
JCMT observation. The infall rate is $2.0\times 10^{-3}~{}M_{\odot}\cdot$
yr-1. High velocity gas is detected in HCN (3-2) emission. The outflow has
three lobes and their total mass is 12 $M_{\odot}$ and momentum of 121
$M_{\odot}\cdot$ km s-1. The average dynamical timescale and the total mass-
loss rate are about $1.6\times 10^{4}$ yr and $1.6\times
10^{-5}~{}M_{\odot}\cdot$ yr-1, respectively. All the findings indicate a
high-mass protostar is forming via rapid accretion. Our results suggest that
nitrogen bearing molecules especially HCN are good for probing both infall and
outflows.

## Acknowledgment

We are grateful to the SMA staff. We also thank Dr. Shengli Qin. for his help
with data reduction and discussion. This work was funded by Grants of NSFC No
10733030 and 10873019.

## References

  * Arce et al. (2007) Arce, H. G., Shepherd, D., Gueth, F., Lee, C.-F., Bachiller, R., Rosen, A., Beuther, H., 2007, arXiv:astro-ph/0603071
  * Bachiller, Gutiérrez, & Pérez. (1997) Bachiller, R., Gutiérrez, M., & Pérez. 1997, ApJ, 487, L93
  * Beuther, Schilke, & Menten. (2002) Beuther, H., Schilke, P., & Menten, K. M. 2002, ApJ, 566, 945
  * Birkmann et al. (2007) Birkmann, S. M., Krause, O., Hennemann, M., Henning, Th., Steinacker, J., Lemke, D. 2007, A&A, 474, 883
  * Bonnell, Bate, & Zinnecker. (1998) Bonnell, I. A., Bate, M. R., & Zinnecker, H. 1998, MNRAS, 298, 93
  * Boonman et al. (2001) Boonman, A. M. S., Stark, R., Van der tak, F. F. S., Van dishoeck, E. F., Van der wal, P. B., Scháfer F., De lange, G., & Laauwen, W. M. 2001, ApJ, 553, L63
  * Carolan et al. (2009) Carolan, P. B., Khanzadyan, T., Redman, M. P., Thompson, M. A., Jones, P. A., Cunningham, M. R., Loughnane, R. M., Bains, I., Keto, E. 2009, MNRAS, 400, 78
  * Choi (2001) Choi, Minho. 2001, ApJ, 553, 219
  * Churchwell (2002) Churchwell, E. 2002, ARA&A, 40, 27
  * Crapsi et al. (2005) Crapsi, A., Caselli, P., Walmsley, C., M., Myers, P., C., Tafalla, M., Lee, C., W., Bourke, T., L. 2005, ApJ, 619,379
  * Cyganowski (2008) Cyganowski, C. J. et al. 2008, AJ, 136, 2391
  * Fuller, Williams, & Sridharan. (2005) Fuller, G. A., Williams, S. J., & Sridharan, T. K. 2005, A&A, 442, 949
  * Garden et al. (1991) Garden, R. P., Hayashi, M., Gatley, I., Hasegawa, T., & Kaifu, N. 1991, ApJ, 374, 540
  * Keto, Ho & Haschick. (1987) Keto, E. R., Ho, P. T. P., & Haschick, A. D. 1987, apj, 318, 712
  * Keto, Ho & Haschick. (1988) Keto, E. R., Ho, P. T. P., & Haschick, A. D. 1988, ApJ. 324, 920
  * Klaassen & Wilson (2007) Klaassen, P. D. & Wilson, C. D. 2007, ApJ, 663, 1092
  * Klaassen & Wilson (2008) Klaassen, P. D., & Wilson, C. D. 2008, ApJ, 684, 1273
  * Lamers et al. (1995) Lamers, H. J. G. L. M., Snow, T. P., Lindholm, D. M. 1995, ApJ, 455, 269
  * Lee et al. (2001) Lee, C.-F., Stone, J. M., Ostriker, E. C., & Mundy, L. G. 2001, ApJ, 557, 429
  * Lee et al. (2007) Lee, J.-E., Di Francesco, J., Bourke, T. L., Evans, N. J. II & Wu, J. 2007, ApJ, 671, 1748
  * Lester et al. (1985) Lester, D. F., Dinerstein, H. L., Werner, M. W., Harver, P. M., Evans, N. J. II, Brown, R. L. 1985, ApJ, 296, 565L
  * Mitchell (1984) Mitchell, G. F. 1984, ApJ, 287, 665
  * Molinari et al. (1996) Molinari, S., Brand, J., Cesaroni, R., Palla, F. 1996, A&A, 308, 573
  * Molinari et al. (2000) Molinari, S., Brand, J., Cesaroni, R., & Palla, F. 2000, A&A, 355, 617
  * Motoyama & Yoshida (2003) Motoyama, K., & Yoshida T. 2003, MNRAS, 344, 461
  * Myers et al. (1996) Myers, P. C., Mardones, D., Tafalla, M., Williams, J. P. & Wilner, D. J. 1996, ApJ, 465, L133
  * Nejad, Williams, & Charnley. (1990) Nejad, L. A. M., Williams, D. A., & Charnley, S. B. 1990, MNRAS, 246, 183
  * Ossenkopf & Henning (1994) Ossenkopf, V., Henning, Th. 1994, A&A, 291, 943
  * Park, Kim, & Ming. (1999) Park, Y. S., Kim, J., & Ming Y. C. 1999, ApJ, 520, 223
  * Schnee et al. (2007) Schnee, S., Caselli, P., Goodman, A., Arce, H. G., Ballesteros-Paredes, J., Kuchibhotla, K. 2007, ApJ, 671, 1839
  * Shu (1977) Shu, F. H. 1977, ApJ, 214, 488
  * Shu, Adams, & Lizano (1987) Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23
  * Sohn et al. (2007) Sohn, J. J., Lee, C. W., Park, Y.-S., Lee, H. M., Myers, P. C., & Lee Y. 2007, ApJ, 664,928
  * Sridharan et al. (2002) Sridharan, T. K., Beuther, H., Schilke, P., Menten, K. M., Wyrowski, F. 2002, ApJ, 566, 931
  * Su et al. (2007) Su, Y.-N., Liu, S.-Y., Chen, H.-R., Zhang, Q., Cesaroni, R. 2007, ApJ, 671,571
  * Sun & Gao (2008) Sun, Y., & Gao, Y. 2008, MNRAS, 392, 170
  * Tsamis et al. (2008) Tsamis, Y. G., Rawlings, J. M. C., Yates, J. A., Viti, S. 2008, MNRAS, 388, 898
  * Velusamy et al. (2008) Velusamy, T., Peng, R., Li, D., Goldsmith, P. F., Langer, William D. 2008, ApJ, 688, L87
  * Welch et al. (1987) Welch, Wm. J., Dreher, J. W., Jackson, J. M., Terebey, S., Vogel, S. N. 1987, Science, 238, 1550
  * Wolfire & Cassinelli (1987) Wolfire, M. G. & Cassinelli, J. P. 1987, ApJ, 319, 850
  * Wu & Evans. (2003) Wu, J., Evans N. J. II. 2003, ApJ, 592, L79
  * Wu et al. (2004) Wu, Y., Wei, Y., Zhao, M., Shi, Y., Yu, W., Qin, S., Huang, M. 2004, A&A, 426, 503
  * Wu et al. (2005) Wu, Y., Zhu, M., Wei, Y., Xu, D., Zhang, Q., & Fiege, J. D. 2005, ApJ, 628, L57
  * Wu et al. (2006) Wu, Y., Zhang, Q., Yu, W., Miller, M., Mao, R., Sun, K., & Wang, Y. 2006, A&A, 450, 607
  * Wu et al. (2007) Wu, Y., Henkel, C., Xue, R., Guan, X., Miller, M. 2007, ApJ, 669, L37
  * Wu et al. (2009) Wu, Y., Qin, S.-L., Guan, X., Xue, R., Ren, Z., Liu, T., Huang, M., Chen, S. 2009, ApJ, 697, L116
  * Wyrowski (2006) Wyrowski, F., Heyminck, S., G$\ddot{u}$sten, R., Menten, K. M. 2006, A&A, 454, L95
  * Yun et al. (1999) Yun, J. L., Moreira, M. C., Afosso, J. M. & Clemens, D. P. 1999, ApJ, 118, 990
  * Zapata et al. (2008) Zapata, L. A., Palau, A., Ho, P. T. P., Schilke1, P., Garrod, R. T., Rodr guez, L. F., & Menten, K. 2008, A&A, 479, L25
  * Zhang & Ho (1997) Zhang, Q., & Ho, P. T. P. 1997, ApJ, 488,241
  * Zhang, Ho, & Ohashi. (1998) Zhang, Q., Ho, P. T. P., & Ohashi, N. 1998, ApJ, 494, 636
  * Zhang et al. (2007a) Zhang, Q., Sridharan, T. K., Hunter, T. R., Chen, Y.; Beuther, H., Wyrowski, F. 2007a, A&A, 470, 269
  * Zhang et al. (2007b) Zhang, Q., Sridharan, T. K., Hunter, T. R., Chen, T., Beuther, H., & Wyrowski, F. 2007b, ApJ, 658, 1152
  * Zhou et al. (1993) Zhou, S., Evans, N. J. II., Koempe, C., Walmsley, C. M. 1993, ApJ, 404, 232

Table 1: Outflow parameters of each lobe. outflow | Vmax | $t_{dyn}$ | Mass | momentum | Energy | $L$ | $\dot{M}_{out}$
---|---|---|---|---|---|---|---
| km s-1 | (10${}^{4}~{}yr$) | ($M_{\odot}$) | ($M_{\odot}$ km s-1) | (10${}^{45}~{}erg$) | ($L_{\odot}$) | ($10^{-6}M_{\odot}$ yr-1)
southern lobe | 11.7 | 1.4 | 3.3 | 29 | 3.0 | 1.8 | 4.0
north-eastern lobe | 15.3 | 1.0 | 4.2 | 42 | 4.5 | 3.6 | 8.0
western lobe | 15.3 | 2.3 | 4.6 | 50 | 5.5 | 2.1 | 4.0

Figure 1: The 1.3 mm continuum emission towards JCMT 18354-0649S. The left one
is obtained with the compact configuration. The rms level is 3 mJy beam-1 (1
$\sigma$). The contours are at -6, 3, 6, 12, 21, 33, 48, 66, 87$\sigma$. The
right panel gives the contours of the continuum emission combined from both
configurations. The rms level is 2.5 mJy beam-1 (1 $\sigma$) and the contours
are at -6, 3, 6, 12, 21, 33, 48$\sigma$. The three near-infrared sources are
marked with crosses. Figure 2: The lower-left panel is the contours of the
CH3OH integrated intensity overlaid on the 1.3 mm continuum image (grey
scale). The contours start from 30$\%$ in steps of 10$\%$ of the peak emission
(15 Jy beam${}^{-1}\cdot$km s-1). The three near-infrared sources are marked
with crosses. The beam-averaged spectrum of CH3OH at three positions are
presented in the other panels. The gaussian fit towards each spectrum is shown
with solid lines. Figure 3: Position-velocity diagrams of CH3OH along a P.A.
of 0$\arcdeg$. The contour levels are from 15% to 90% in steps of 15% of the
peak intensity in both panels. The intensity at the peak is 3.73 Jy beam-1.

Figure 4: The HCN (3-2) spectra. Left: the solid black line exhibits the
spectrum constructed from SMA compact array and the dashed line shows the
spectrum obtained from combining the compact and extended data together. Both
of the two spectra are integrated over a region of $5\arcsec\times 5\arcsec$.
Right: the solid black line shows the spectrum constructed from the combined
SMA and JCMT data, which is convolved with the JCMT beam ($18.3\arcsec$); the
dashed line shows the spectrum from JCMT only; the dash-dotted gray line shows
the spectrum obtained with the SMA compact array and convolved with the JCMT
beam. The vertical dashed lines in both panels mark the position of the
systematic velocity (96.7 km s-1). Figure 5: Combined JCMT and SMA HCN (3-2)
channel maps from 80 km s-1 to 109 km s-1, which is convolved with the beam of
the SMA compact configuration ($3.76\arcsec\times 2.72\arcsec$,
PA=-54$\arcdeg$). The contours are in steps of 0.5 Jy beam-1 (3$~{}\sigma$)
from 0.5 Jy beam-1 (3$~{}\sigma$). The velocity of each channel is plotted at
the upper-left of each panel, and the beam size at the lower-right. The
positions of MM1, MM2 and MM3 are marked with crosses.

Figure 6: Position-velocity diagrams of HCN (3-2) outflow observed by the SMA
along a P.A. of 15$\arcdeg$ (left panel) and 90$\arcdeg$ (right panel). The
image is smoothed to 2 km s-1 velocity resolution. The contour levels are from
15% to 90% in steps of 15% of the peak intensity in both panels. The peak is
1.84 Jy beam-1 in the left panel and 1.92 Jy beam-1 in the right panel. The
four clumps are labeled by solid lines with arrows. The clumps are distinguish
by the thick dashed lines. Figure 7: The high-velocity HCN (3-2) intensity
contours overlaid on the 1.3 mm continuum image, integrated from 80 to 87 km
s-1 for the blueshifted lobes (solid contours) and from 103 to 109 km s-1 for
the redshifted lobe (dashed contours), with contours from 30$\%$ in steps of
10$\%$ of the peak emission. The peak is 8.55 Jy beam-1$\cdot$km s-1 for the
blueshifted lobes and 7.16 Jy beam-1$\cdot$km s-1 for the redshifted lobe. The
empty and solid ellipses in the lower-right corner represent the synthesized
beams of HCN (3-2) emission and 1.3 mm continuum emission combined from both
configurations, respectively.