\section{Introduction}

Quark-Gluon-Plasma (QGP) is a state of matter in the Strong
Interaction system which is actively sought through theoretical study
and heavy ion collision at high energy (or high temperature) and at
low energy (or high baryon density) and must have played an important
role in the evolution of early universe
\cite{Yagi:2005yb}. Unfortunately at physical quark masses, chiral
symmetry restoration/deconfinement transition in QCD from normal
hadronic matter to QGP matter is a cross-over transition and lacks
characteristic behaviors of order parameter in a first order or a
second order thermodynamic phase transition across the thermodynamic
transition. Thus, finding and understanding physical phenomena which
can serve as signatures for QGP matter formation is quite crucial.

Physics of heavy quarkonium in medium may be one of such candidate
phenomena \cite{Matsui:1986dk}. In contrast to large uncertainty in
our understanding of light hadron production and decays, inclusive
productions and decays of heavy quarkonium in proton--proton or
proton--anti-proton collision, which are required for the comparison
with those in heavy ion collisions, is well understood:
Non-Relativistic QCD (NRQCD) factorization in which the scale above
heavy quark (like charm quark or bottom quark) mass is ``integrated
away'' and the smallness of heavy quark velocity at the quarkonium
rest frame offers ``small parameter'' is applicable in heavy
quarkonium \cite{Bodwin:1994jh}. Short distance effects can be
calculated via perturbative expansion in $\alpha_s$ and long distance
non-perturbative matrix elements can be computed either using lattice
NRQCD at zero
temperature\cite{Bodwin:1996tg,Bodwin:2001mk,Bodwin:2005gg} or can be
fitted using experimental data (see e.g.,
\cite{Butenschoen:2010rq}). Based on our understanding of quarkonium
physics in p--p collision, medium effect for quarkonium in heavy ion
collision may be better disentangled. In particular, bottomonium
system offers a better probe for QGP formation since there are clearer
seperations between the long distance physics and the short distance
physics, compared to those in charmonium system due to smaller quark
velocity ($v_b^2 \sim 0.1$ compared to $v_c^2 \sim 0.3$ for charmonium
system) at the quarkonium rest frame.

In recent years, we studied the behavior of S-wave and P-wave
bottomonium states for a temperature range $0.4 T_c \le T \le
2.1 T_c$ using NRQCD in non-zero temperature on anisotropic lattices
($12^3 \times N_t$) with two light quark flavors
\cite{Aarts:2010ek,Aarts:2011sm}. Details on how dynamical gauge field
is generated and what the simulation parameters are can be found in
\cite{Aarts:2010ek,Aarts:2011sm,Morrin:2006tf}. Our calculation method
is explained in \cite{Aarts:2010ek,Aarts:2011sm}. We expect that
in-medium effect such as temperature ($T$) effect or baryon chemical
potential ($\mu$) effect can be accommodated by NRQCD as long as
$\frac{T}{M} \ll 1$ or $\frac{\mu}{M}\ll 1$ where $M$ is the heavy
quark mass (see, e.g. \cite{Hands:2012yy} for baryon chemical
potential efffect in NRQC$_2$D). Indeed our study
\cite{Aarts:2010ek,Aarts:2011sm} meets such expectation. We found that
P-wave bottomonium states ($\chi_{bJ}$) melt immediately above $T_c$,
and bound state signals for S-wave channel ($\eta_b, \Upsilon$)
survive above $T_c$ \cite{Aarts:2010ek}. Closer inspection of the
spectral functions for S-wave correlators ($\chi_b, \Upsilon$), which
is obtained by the Maximum Entropy Method (MEM) with NRQCD kernel,
shows that $(1S)$ peaks appears to survive and the excited state peaks
are suppressed \cite{Aarts:2011sm}. This observation is consistent
with recent CMS experiment which shows sequential suppression of
$(3S), (2S)$ $\Upsilon$ states and survival of $(1S) \Upsilon$
\cite{Chatrchyan:2011pe,Chatrchyan:2012fr}.

Here, we discuss further lattice NRQCD study of ours on bottomonium at
non-zero temperature. First, bottomonium with non-zero momentum is
studied (for full report, see \cite{Aarts:2012ka}). Quarkonium
production in hadron collision proceeds: parton evolution from hadron,
then partonic scattering process (gluon fusion, heavy quark
recombination, gluon Compton scattering, and light quark--light quark
scattering). Among these partonic processes, quarkonium production is
mostly dominated by gluon fusion. Gluon subsequently fragments into
quarkonium. In this case, fragmenting gluon will have large energy and
momentum. If quarkonium is produced in thermal bath, then it will be
moving relative to thermal bath. In contrast, if quarkonium is
produced through recombination of thermalized heavy quark and heavy
anti-quark (albeit less probable), produced quarkonium is at rest in
thermal bath. Schematic view on these two processes is shown in
Fig. \ref{fig:1}. Since a thermal bath defines a preferred rest frame,
``moving'' effect will induce $\frac{1}{2} M v^2$ shift in the mass
compared to the mass of quarkonium produced at rest. This shift may be
larger than experimental resolution for muon detection and may be
observable.

Secondly, we start to investigate systematic errors associated with
our computation in \cite{Aarts:2010ek,Aarts:2011sm}. As a first step,
we investigate how $10 \sim 20 \%$ change in bottom quark mass affects
previous observations made in \cite{Aarts:2010ek,Aarts:2011sm} by
calculating bottomonium correlators with $M a_s = 3.7, 4.14$ and $5.4$
and looking at the correlator ratios.

\begin{figure}[t]
\begin{center}
\epsfig{figure=schem.eps,width=0.3\textwidth}
\end{center}
\caption{Quarkonium production processes: (a)fragmentation from
  gluon (top), (b) recombination of thermalized heavy quark--heavy
  anti-quark (bottom)}
 \label{fig:1}
\end{figure}

\section{S-wave bottomonium moving in thermal bath}

We consider $\eta_b$ and $\Upsilon$ which are moving with regard to
the thermal bath, i.e., S-wave bottomonium states with
momentum. Non-zero momentum quarkonium state is obtained by injecting
momentum to the source for NRQCD quark propagator computation and
projecting out the injected momentum component from the quarkonium
propagator which is obtained by combining the zero-momentum NRQCD
quark propagator with the momentum injected NRQCD quark
propagator. Momenta, $\vec{k} = \frac{2\pi}{L} (n_1,n_2,n_3)$ where
$\vec{n} = (0,0,0), (1,0,0), (1,1,0), (1,1,1), (2,0,0), (2,1,0)$, and
$(2,2,0)$ with $L=12$, have been considered. Since $M a_s = 4.5$ is
used for the bottom quark mass in NRQCD propagator computation and
$M_{\eta(1S)} = 9.391$ GeV and $M_{\Upsilon (1S)} = 9.460$ GeV
\cite{Beringer:1900zz}, these quarkonium momentum range amount to
quarkonium velocity, $0.07 \stackrel{<}{\sim} \;\; v \;\;
\stackrel{<}{\sim} 0.18$ with $\vec{p}^2 a_s^2 = 4 \sum_{i=1,3} \sin^2
\frac{k_i}{2}$.  From the correlators, spectral functions defined as
\begin{equation}
G(\tau, \vec{p}) = \int \frac{d \omega}{2\pi} K(\tau, \omega) \rho
(\omega, \vec{p}), \;\;\;\;\; K(\tau, \omega) = e^{-\omega \tau}
\end{equation}
are computed using MEM. Fig. \ref{fig:2} shows a typical spectral
function as a function of $|\vec{p}|$. For each momentum, the spectral
function exhibits a distinct shape. From this, the peak and width can
be extracted and Fig. \ref{fig:3} show extracted peak values and width
values which corresponds the $\Upsilon (1S)$ state. To a leading order
in $\vec{p}^2$ and $T$, we see that the $1S$ peak roughly follows
\begin{equation}
\Delta E (\vec{p}, T) \simeq \Delta E (\vec{p} = 0, T = 0) +
\frac{\vec{p}^2}{2M} + c (\alpha_s) T \label{eqDeltaE}
\end{equation}
In general, continuum NRQCD dispersion is
\begin{equation}
E = M_0 + \frac{\vec{p}^2}{2M_2} - \frac{(\vec{p}^2)}{8 M_4^3} + \cdots
\end{equation}
Since we use ${\cal O} (v^4)$ NRQCD lagrangian to compute bottomonium
propagator , $M_0$ has ${\cal O} (v^4)$ error, $M_2$ has ${\cal O}
(v^2)$ error, and $M_4$ ${\cal O} (v^0)$ respectively. Using the form
of $M_0 + 4\sum_i( \sin (k_i/2)^2)/2M_2 -
(4\sum_i(\sin(k_i/2)^2)^2)/8M_2^3$ to fit $N_t = 80$ $\Upsilon$
correlator, we get 11.05 GeV for the kinetic mass ($M_2$) of $\Upsilon
(1S)$ state with $M a_s = 4.5$. The fitting form $M_0 + 4\sum_i(\sin
(k_i/2)^2)/2M_2$ gives 11.12 GeV for $M_2$. Thus Eq. (\ref{eqDeltaE})
is suitable to describe our lattice result. On the other hand, the
$1S$ width is independent of $\vec{p}^2$ and is just proportional to
$T$.
\begin{figure}[t]
\begin{center}
\epsfig{figure=rho-vs-w-upsilon-24.eps,width=0.5\textwidth}
\end{center}
\caption{$\Upsilon$ spectral functions on $12^3 \times 24$ for various
  momenta}
 \label{fig:2}
\end{figure}

\begin{figure}[t]
\begin{center}
\epsfig{figure=upsilon_mass_vs_p_2_all_Nt.eps,width=0.45\textwidth}
\epsfig{figure=upsilon_width_vs_p_2_all_Nt.eps,width=0.5\textwidth}
\end{center}
\caption{$\Upsilon (1S)$ energy vs. various momenta from the first
  peak of the spectral functions at each $T/T_c = 2.09, 1.86, 1.68,
  1.40, 1.20, 1.05,$ and $0.42$ for $N_t = 16, 18, 20, 24, 28, 32$,
  $80$ respectively (left). $\Upsilon (1S)$ width vs. various momenta
  from the width of the first peak in spectral functions
  (right)} \label{fig:3}
\end{figure}

This lattice NRQCD finding can be compared with effective theory
prediction. NRQCD bound state at rest in thermal bath below the
melting temperature has been considered in \cite{Brambilla:2010vq} and
non-relativistic bound state which is moving with regard to thermal
bath has been considered in QED (in the rest frame of the bound state,
i.e., bound state in moving thermal bath)
\cite{Escobedo:2011ie}. Finding in \cite{Escobedo:2011ie} for the
width can be summarized as
\begin{equation}
\frac{\Gamma_v}{\Gamma_0} = \frac{\sqrt{1-v^2}}{2v} \log
\left(\frac{1+v}{1-v} \right) \rightarrow 1 - \frac{2}{3} v^2 + {\cal
  O} (v^4)
\end{equation}
where $\Gamma_0$ is the width of the bound state not moving with
regard to thermal bath and $\Gamma_v$ is the width of the bound state
moving with the velocity $v$ in thermal bath. With $v^2 < 0.04$,
``moving effect'' in the width is beyond the precision level of our
computation, which explains negligible momentum dependence of the
width in the right figure of Fig. \ref{fig:3}.

\section{Mass dependence of bottomonium correlators}

Our previous bottomonium study was done with bottom quark mass, $M a_s
= 4.5$. The same NRQCD correlator calculation is repeated with $M a_s
= 3.7 (-20 \%), 4.14 (-10 \%)$ and $5.4 (+ 20\%)$ in order to
investigate how small change in bottom quark mass affects bottomonium
spectrum and how important this effect is compared to non-zero
temperature effect. Since lattice NRQCD is an effective theory and is
consistent when $M a_\tau \sim 1$, all these choices ($M a_\tau =
0.62, 0.75, 0.9$ for $M a_s = 3,7, 4.5, 5.4$ respectively) are
expected to be valid.

Lattice NRQCD quarkonium correlator behaves as $A e^{-\Delta E \tau}$
when a bound state is formed and fitting $\Upsilon$ correlator
calculated on $N_t = 80$ configurations gives $\Delta E_\Upsilon (M,T
\sim 0)$. Using experimental value for $M_\Upsilon$, we determine the
unknown constant for NRQCD spectrum by
\begin{equation}
M_\Upsilon^{\rm exp} = E_0 (M) + \Delta E_\Upsilon (M,T = 0)
\end{equation}
If non-zero temperature quarkonium states has a pole like
$\Upsilon (1S)$ state \cite{Aarts:2011sm},
\begin{equation}
M_\Upsilon (T) = E_0 + \Delta E_\Upsilon (M,T) = M_\Upsilon^{\; \rm exp}
+ \Delta E_\Upsilon (M,T) - \Delta E_\Upsilon (M,T = 0).
\end{equation}
So, the ratio of the non-zero temperature quarkonium correlator to the
$N_t = 80$ correlator can reveal the heavy quark mass dependence
$\Delta E (M,T)$ although the true ground state behaviors for $N_t =
80$ correlators and for $N_t = 16$ or $20$ correlators will not set in
at large $\tau$. Fig. \ref{fig:5} shows such ratios for $T = 458 (N_t
= 16)$ (MeV) and $T = 408 (N_t = 18)$ (MeV).

\begin{figure}[t]
\begin{center}
\epsfig{figure=3S1ratioNt16-20.eps,width=0.5\textwidth}
\end{center}
\caption{$\Upsilon$ correlator ratio of $G(\tau;N_t)$ to
  $G(\tau;{N_t = 80})$ for $N_t = 16$ (filled symbol)
  and $18$ (open symbol)}
 \label{fig:5}
\end{figure}

In Fig. \ref{fig:5}, the comparison between the filled symbol and the
empty symbol shows the temperature dependence for the same $M$, and
the comparison between the different symbols shows $M-$dependence for
the same $T$. The correlator ratio change due to 10\% change in the
two different $T$ is larger than due to $10 \sim 20 \%$ change in the
4 different $M$. Since the small change in heavy quark mass is less
important than that due to the temperature change, lattice NRQCD is
valid as an effective theory and systematic error due to heavy quark
mass mis-tuning is small.

\section{Conclusion}

Using NRQCD formalism in non-zero temperature, we computed momentum
dependent $(\eta_b, \Upsilon)$ correlators and calculated spectral
functions using MEM via NRQCD kernel. From the peaks and the width of
momentum dependent spectral function, we find that there are
observable $v^2$ effect on the energy of S-wave state moving in
thermal bath but no observable effect on the width of S-wave state
moving in thermal bath for $v^2_{\rm upsilon} \stackrel{<}{\sim} 0.04$.

Among various sources of systematic errors involved in our result,
bottomonium correlators calculation is repeated with $20 \%$ variation
in bottom quark mass. By studying correlator ratios, We find that the
temperature effect is more important than the heavy quark mass effect
in S-wave bottomonium at the temperature around a few $T_c$, which
suggests lattice NRQCD as an effective theory for bottomonium in
non-zero temperature is a consistent theory.

Improved tuning for bottom quark mass and detecting thermal deviations
from the standard dispersion relation is in progress for $N_f = 2 + 1$
flavor system with smaller spatial lattice spacing and a larger
spatial extent.

\section{Acknowledgements}

We acknowledge the support and infrastructure provided by the Trinity
Centre for High Performance Computing and the IITAC project funded by
the HEA under the Program for Research in Third Level Institutes
(PRTLI) co-funded by the Irish Government and the European Union. The
work of CA and GA is carried as part of the UKQCD collaboration and
the DiRAC Facility jointly funded by STFC, the Large Facilities
Capital Fund of BIS and Swansea University. GA and CA are supported
by STFC. SK is grateful to STFC for a Visiting Researcher Grant and
INFN for the visit to Frascati, and is supported by the National
Research Foundation of Korea grant funded by the Korea government
(MEST) No.\ 2011-0026688. SR is supported by the Research Executive
Agency (REA) of the European Union under Grant Agreement number
PITN-GA-2009-238353 (ITN STRONGnet) and the Science Foundation
Ireland, grant no.\ 11-RFP.1-PHY-3201. DKS is supported in part by US
Department of Energy contract DE-AC02-06CH11357. JIS has been
supported by Science Foundation Ireland grant 08-RFP-PHY1462 and
11-RFP.1-PHY-3193-STTF11.


\section{Introduction}
\label{sec:intro}
Graphs are a common topological representation across a variety of
scientific fields. They are used when relations between a large number
of entities have to be specified in a succinct manner. Chemical
reactions, molecules, social networks (both physical and online) and
the electric grid are some common examples. In many cases, we may have
partial information about the graph, requiring generation of many
graphical realizations consistent with the (partial) characterization.
Such situations arise in case of large online social networks (of
which only a small part can be sampled tractably) or where the data
simply cannot be collected e.g., the web of human sexual relations
(which only allows the estimation of the degree distribution). In
other cases, privacy concerns can prevent the distribution of a graph
for experimentation and or study, e.g., networks of email
communications, critical infrastructure nets, etc. This gives rise to
the problem of constructing ``sanitized'' proxies, that preserve only
some properties of the original graph. Thus, being able to generate
independent graphs conditioned on an incomplete set of graphical
measurements is essential for many applications.

Many graph generation methods aim to preserve the salient features of
graphs~\cite{12sk3a, 12pf4a,55hv1a,62hs1a,96wp2a, kronecker,
Aiello00,Chung2, 12gk3a}.  For statistical analysis, we need
algorithms that can generate uniformly random instances from the space
of graphs with a specified feature.  There has been significant work
on random graph generation with a given degree distribution (DD),
which specifies the number of vertices with a given degree.
In~\cite{ktv}, the problem of generating a graph with a given degree
distribution was reduced to a prefect matching problem, which can be
used to generate random instances by employing results in sampling
perfect matchings~\cite{JeSiVi04,Br86}.  Alternatively, sequential
sampling methods were investigated in~\cite{BlDi11,BaKiSa10}.  These
methods can be compute intensive, and in practice,
Markov chains (MC) are widely used due to their simplicity and
flexibility.  The MC is started using a graph that honors the
specified graphical characteristics, and uses a procedure that can
``rewire'' a graph, to create ensembles of graphs with the same degree
distribution.  Taylor~\cite{81tr1a} showed that edge-swaps could
modify a graph while preserving its DD.
Kannan et al~\cite{ktv} analyzed the mixing time of such a MC, whereas
Gkantsidis~\cite{GkantsidisMMZ03} devised a MC scheme that avoids self
loops. In~\cite{StPi12}, Stanton and Pinar used an MC to generate an
ensemble of graphs using a rewiring scheme that preserved the joint
degree distribution (JDD), which specifies the number of edges between
vertices of specified degrees.  Stanton and Pinar also empirically
analyzed the mixing of the MC using the binary ``time-series'' traced
out by the appearance/disappearance of edges between two labeled
nodes, as the MC executed its random walk in the space of graphs. They
showed that the autocorrelation of the time-series decayed to zero, a
necessary condition for the MC to converge to its stationary
distribution~\cite{Sokal}.

Graph rewiring schemes that preserve DD or JDD are simple (and fast)
and MC chains driving them are easy to construct. However, successive
graphs generated by the MC are only slightly different and the MC has
to proceed for a large number of steps $N$ before the initial graph is
``forgotten.'' The mixing time estimates in~\cite{ktv} take the form
of upper bounds and of $O(|E|^6)$, where $|E|$ is the number of edges
in a graph. Even for small graphs with 1000 edges, the bound is
intractable. In practice, $N$ is chosen arbitrarily. Failure to mix
completely leads to the generation of correlated samples and any
results derived from them are erroneous.  While the method
in~\cite{StPi12} demonstrated the use of autocorrelation to establish
mixing, it is an \emph{a posteriori} test which provides no \emph{a
priori} guidance regarding the length of the MC chain. Further, the
method requires a small, user-specified threshold, below which the
time-series autocorrelation is deemed zero.

In this paper, we construct analytical models that estimate $N$. These
models track the evolution of the binary ``time-series'' formed by the
edges of labeled graphs generated by the MC.  We test the model under
conditions where the DD or the JDD is held constant. In
Sec.~\ref{sec:models}, the model is used to derive an expression for
the mixing time based on when the binary time-series begins to
resemble independent draws from a distribution. The models predict
that the mixing time is proportional to $|E|$, the number of edges in
the graph. The model holds true for a ``representative'' edge and
exploits the constancy of DD (or JDD) in arriving at its
prediction. The model is tested with real graphs in
Sec.~\ref{sec:tests}. In Sec.~\ref{sec:verif}, we develop a
data-driven method, that assumes neither a constant DD (or JDD) nor a
representative edge, to investigate the independence of the edge
time-series. We use this test to verify the assumptions made when
developing mixing time expressions in Sec.~\ref{sec:models}. We
determine the fraction of the edges for which the model prediction of
$N$ is an underestimate, and the consequences of the lack of
stationarity on the ensemble of graphs sampled in this manner.

\section{Theoretical analysis}
\label{sec:models}

The goal of this paper is to provide a mathematically principled
argument for running a MC ${\rm O(|E|)}$ steps to generate independent
graphs with $|E|$ edges. The constant hidden in the big-Oh depends on
the desired accuracy. The graphs so generated may have a prescribed DD
or a JDD; we find that a MC on graphs with a prescribed DD mixes
slightly easier than those where the JDD is preserved. Our empirical results
show that ${\rm 5|E|}$ -- ${\rm 30|E|}$ steps are sufficient for
mixing these MCs. We provide a mathematical justification
for this observation.

Consider a MC on the space of graphs and two labeled vertices $u$ and
$v$. Under certain circumstances, the existence of an edge $(u,v)$ can
also be described by a Markov chain. Thus we model the behavior of the
MC on the space of graphs in terms of a set of coupled 2-state Markov
chains, each representing edge behavior. We develop the transition
matrix for the (edge) MC, which is exact for DD-preserving
experiments, but heuristically derived for MCs with prescribed
JDDs. We prove an ${\rm O(|E|)}$ mixing time for these 2-state Markov
chains.  This ensures that that after ${\rm O(|E|)}$ MC steps, each
edge behaves as if it comes out of independent draws. While this is
not sufficient for complete mixing, it is a \emph{necessary}
condition.

These mathematical theorems will be empirically validated in the next section,
but we give a quick preview of the method.
Starting from a real
graph, we run a MC (driving a DD- or JDD-preserving rewiring scheme)
for $\ell$ steps, to create a new, nominally independent graph. This
is repeated $M$ times to create an ensemble of graphs. We plot the
distribution of some common graphical quantity (e.g., the number of
triangles in each graph) constructed from the ensemble. This is
repeated for different values of $\ell$ till the plots converge. We
find that $\ell \sim {\rm 10 |E|}$. When $\ell \sim {\rm |E|}$, the MC
is far from convergence.

Below, we review the rewiring schemes used in the paper. This is
followed by derivations of the models of mixing under DD and JDD
preservation.

\subsection{A review of rewiring schemes}
\label{sec:rewire}

Consider an undirected graph $G=(V,E)$, where $|V| = n$, the number of
vertices in the graph, and $|E| = m$, the number of edges.  The degree
distribution of the graph is given by the vector $\mathbf{f}$, where
$\mathbf{f}(d)$ is the number of vertices of degree $d$.  The \emph{joint
degree distribution} is an $n \times n$ matrix $\mathbf{J}$, where
$\mathbf{J}(i,j)$ is the number of edges incident between vertices of degree
$i$ and degree $j$. We use $d_v$ to denote the degree of vertex $v$.

\begin{figure}
  \includegraphics[scale=0.5]{ddswap.png}
  \caption{The DD-preserving swapping operation described in
  Sec.~\ref{sec:rewire}.}
  \label{fig:ddswap}
\end{figure}

In this paper, we use the greedy, JDD-preserving rewiring procedure by
Stanton and Pinar in~\cite{StPi12}. This is analogous to the degree
distribution preserving chain~\cite{ktv,GkantsidisMMZ03}, illustrated in Fig~\ref{fig:ddswap}.
The JDD-preserving rewiring
procedure is shown in Fig~\ref{fig:jddswap}. The procedure does not guard
against parallel edges or self-loops, but the MC itself can be
slightly modified to ensure that these artifacts do not
occur~\cite{StPi12,GkantsidisMMZ03}. We also maintain lists of nodes
and edges indexed by their degree, so that for a given degree $d$, we
can locate a uniform random edge that is incident on a degree-$d$
node. The steps in the rewiring scheme are:

\begin{asparaitem}
  \item Pick an \emph{endpoint} uniformly at random. This is done by
  picking a random edge and choosing one of its endpoints with
  equal probability. Let this vertex be $u_1$. We then choose
  a uniform random neighbor $v$ of $u_1$, to get edge
  $(u_1,v)$.

  \item Pick uniformly at random another vertex $u_2$ such that $d_{u_1} =
  d_{u_2}$, and choose a uniform random neighbor $w$ incident to
  $u_2$. This gives the second edge.

  \item Swap the edges $(u_1,v)$ and $(u_2,w)$ to create $(u_1,w)$ and
  $(u_2,v)$.
\end{asparaitem}
Details of the rewiring scheme and discussions of the ergodicity,
mixing of a MC driving such a rewiring scheme can be found
in~\cite{StPi12}. A pictorial description is in Fig.~\ref{fig:jddswap}.

\begin{figure}
  \includegraphics[scale=0.5]{jdd.png}
  \caption{The JDD-preserving swapping operation described in
  Sec.~\ref{sec:rewire}.}
  \label{fig:jddswap}
\end{figure}

\subsection{Modeling mixing when degree distribution is preserved}
\label{sec:DD}

Consider a MC ${\cal M}$ that generates graphs with a prescribed degree
distribution. Consider two labeled vertices $u$ and $v$, with degrees
$d_u$ and $d_v$. We will show that the appearance/disappearance of $(u,v)$ in  ${\cal M}$ can also be
described by a Markov chain.

\begin{claim}
\label{clm:T}
  The transition matrix of the Markov chain for the edge $(u,v)$ is given by
  \begin{eqnarray}
    \M{T}(d_u, d_v)  & = &  \left( \begin{array}{cc}
                                  1 - \alpha(d_u,d_v) &   \alpha(d_u,d_v) \\
                                  \beta (d_u,d_v)     &  1 - \beta(d_u,d_v)
                               \end{array} \right)  \nonumber

  \label{eqn:T}
  \end{eqnarray}
  where  $\alpha(d_u,d_v)=\frac{d_ud_v}{2m^2}$,   $\beta (d_u,d_v)=1-(1-1/m)^2$, and   the first index in $T$ is state $0$ (no edge).
\end{claim}

\begin{proof}
  Assume that the edge $(u,v)$ does not initially exist. We analyze
  the probability of it appearing after one step. Let $e$ and
  $e^{\prime}$ be the two edges chosen by the rewiring algorithm. $(u,
  v)$ will appear if $u$ is the first endpoints of $e$ and $v$
  is the last endpoint of $e'$ (or vice versa). This leads to four
  possible cases. Consider the case when both the vertices are the
  first endpoints of $e$ and $e^{\prime}$. The probability of this
  event is $d_u/2m \times d_v/2m = d_u d_v / 4m^2$. The other case
  is symmetric and so $\alpha(d_u, d_v) = d_u d_v / 2m^2$.

  Suppose that $(u, v)$ does exist but is removed during the Markov
  transition. This occurs if $(u, v)$ is either $e$ or
  $e^{\prime}$. The probability of this event is $\beta_(d_u, d_v) = 1 - (1-1/m)^2$. Note
  that the values of $\alpha$ and $\beta$ so derived are completely
  independent of the structure of the rest of the graph, and so the
  transition probabilities in Eq.~\ref{eqn:T} satisfy the Markov property.
\end{proof}

\subsection{Modeling mixing when joint degree distribution is preserved}
\label{sec:JDD}

We consider a pair of labeled vertices $(u, v)$ and model the
appearance/disappearance of the edge during a Markov transition.
For convenience, we will assume that $\mathbf{f}(d_u)$ and $\mathbf{f}(d_v)$
are both strictly greater than $1$. The proofs can be extended to these
cases as well. (Observe that when $\mathbf{f}(d_u) = \mathbf{f}(d_v)$, the
edge $(u,v)$ is always present because the JDD is preserved.)

\begin{claim}
\label{clm:beta}
 Assume edge $(u, v)$ is present and is removed using the transition. The probability of  this event is
 \begin{equation}
   \beta(d_u, d_v) = \frac{1}{m} + \frac{\mathbf{f}(d_u) -  1}{2m\mathbf{f}(d_u)} + \frac{\mathbf{f}(d_v) -  1}{2m\mathbf{f}(d_v)}
 \end{equation}
\end{claim}
\begin{proof}
  We follow the JDD-preserving swapping procedure described in
  Sec.~\ref{sec:rewire}. Let $e$ and $e^{\prime}$ be the two edges
  chosen for swapping. If   $(u, v)$ is $e$, it will definitely be
  swapped. The probability of the chosen edge being $(u, v)$ is $1/m$.

  The other alternative is that  $(u, v)$ is $e^{\prime}$.
  This can happen only if the first endpoint chosen has degree $d_u$ (or $d_v$)
  and then $(u,v)$ is chosen. Consider the first case.
  The total number
  of endpoints on degree $d_u$  vertices (excluding $u$) is
  $(\mathbf{f}(d_u)-1)d_u$. The probability that $e$ is chosen
  is the probability that $u$ is chosen as the second endpoint, and
  then $(u,v)$ is chosen. This is equal to $1/(\mathbf{f}(d_u) \times d_u)$.
  The total probability is
  \[  \frac{(\mathbf{f}(d_u) - 1)d_u}{2m} \times \frac{1}{\mathbf{f}(d_u)d_u} = \frac{(\mathbf{f}(d_u) -  1)d_u}{2m\mathbf{f}(d_u)d_u}  = \frac{\mathbf{f}(d_u) -  1}{2m\mathbf{f}(d_u)}\]

  Symmetrically, the randomly chosen endpoint could have had degree
  $d_v$. So the total probability of choosing $e^{\prime} = (u, v)$ is
$$ \frac{\mathbf{f}(d_u) -  1}{2m\mathbf{f}(d_u)} + \frac{\mathbf{f}(d_v) -  1}{2m\mathbf{f}(d_v)} $$
  and the probability of $(u, v)$ being removed is
$$
     \beta(d_u, d_v) = \frac{1}{m} + \frac{\mathbf{f}(d_u) -  1}{2m\mathbf{f}(d_u)} + \frac{\mathbf{f}(d_v) -  1}{2m\mathbf{f}(d_v)}.
$$
\end{proof}

This expression has certain ramifications. Since the JDD is preserved
(and by implication, $\mathbf{f}(d)$), the probability of edge removal does not
vary as the Markov chain proceeds. This satisfies the Markov property
and $\beta$ depends only on $u$ and $v$ i.e., it is independent of the
rest of the graph.

Consider the probability $\alpha$ of $(u, v)$ appearing. This is not
Markovian, as it depends on the structure of the graph. However, this
dependence is weak, and we can obtain a Markovian
approximation using a simple heuristic.

Consider the number of edges $\M{J}(d, d_v)$ between vertices of
degree $d_v$ and an arbitrary $d$. The \emph{expected} number of such
edges incident on a degree $d_v$ vertex is $\M{J}(d,d_v)/\mathbf{f}(d_v)$. The
number of edges between a vertex $v$ and degree $d$ vertices depends
on the graphical structure; however, we will approximate it with the
expected value $\M{J}(d,d_v)/\mathbf{f}(d_v)$.

\begin{claim}
\label{clm:alpha}
  If $(u, v)$ is not present at a stage in the Markov chain, the
  probability $\alpha$ that it appears in the next step is given by
  \begin{equation}
    \alpha(d_u, d_v) = \frac{2\M{J}(d_u,d_v)}{m\mathbf{f}(d_u)\mathbf{f}(d_v)}
  \label{eqn:alpha}
  \end{equation}
\end{claim}

\begin{proof} As before, let the two edges chosen by the algorithm
  be $e$ and $e'$, in that order. The edge $(u,v)$ is swapped in
  when either of the four following situations happens.
  Vertex  $u$ is the first (resp. second) endpoint of $e$ and $v$ is the
  second (resp. first) endpoint of $e'$. Or,
  vertex  $u$ is the first (resp. second) endpoint of $e'$ and $v$ is the
  second (resp. first) endpoint of $e$.

  $u$ is the first endpoint of $e$ and $v$ is the second endpoint of $e'$ (or vice versa).

  Consider the following sequence of events. First, $u$ is chosen as an endpoint
  (and we get edge $e$). Then, we choose a vertex $u'$ (who is a neighbor
  of $v$) whose degree is the same as $u$.
  Finally, we select edge $e' = (u',v)$. Swapping will lead to edge $(u,v)$.
  The probability of the first event is $d_u/2m$. The probability of the
  second event depends the number of neighbors of degree $d_u$ incident
  to $v$. Alternately, this is the number of edges connecting $v$
  to degree $d_u$ vertices. Clearly, this depends on the graph structure, but our approximation of this
  quantity is $\M{J}(d_u,d_v)/\mathbf{f}(d_v)$.

  So, the second probability is $\M{J}(d_u,d_v)/(\mathbf{f}(d_v)\mathbf{f}(d_u))$ (because
  there are $\mathbf{f}(d_u)$ neighbors of degree $d_u$, of which $\M{J}(d_u,d_v)/\mathbf{f}(d_v)$
  are neighbors of $v$). The last probability is $1/d_u$. The total probability is
  \[ \frac{d_u}{2m} \times \frac{\M{J}(d_u,d_v)}{\mathbf{f}(d_u)\mathbf{f}(d_v)}
  \times \frac{1}{d_u}= \frac{\M{J}(d_u,d_v)}{2m\mathbf{f}(d_u)\mathbf{f}(d_v)} \]

  Now consider the second case. First, $v'$ (a neighbor of $u$ whose degree is $d_v$) is chosen
  as an endpoint. Then, $u$ is chosen as the neighbor of $v'$, so $e = (u,v')$.
  Then, $v$ is chosen as the first endpoint of the second edge $e'$.
  Swapping will again lead to edge $(u,v)$. The first probability
  is approximated by $\mathbf{J}(d_u,d_v)/(\mathbf{f}(d_u)\mathbf{f}(d_v))$.
  The second probability is $1/d_v$, and the last is $d_v/2m$.
  The probability comes out to be the same as earlier.

  By symmetry, we get the probabilities of the other cases by switching the roles
  of $u$ and $v$, yielding the same expression.
\end{proof}

\subsection{Estimating the mixing time}

Having developed expressions for $\alpha$ and $\beta$ in
Eq.~\ref{eqn:T}, under the assumption of prescribed DD
(Sec.~\ref{sec:DD}) and JDD (Sec.~\ref{sec:JDD}), we use the Markov
transition matrix to estimate the number of steps $N$ to achieve a
stationary distribution, conditional on a tolerance $\epsilon$.

\begin{claim}
\label{clm:markov} Choose $N$ as follows.
Let the final distribution after
running ${\cal M}_{u,v}$ for $N$ steps be $\bp$. Then $\|\bp - \sta_{u,v}\|
< \varepsilon$.
\[ N=\begin{cases}
(m/2)\ln(1/\varepsilon) & \mathrm{DD\;preserving\;\; MCs} \\
m\ln(1/\varepsilon) &   \mathrm{JDD\; preserving\;\; MCs }
\end{cases}\]
\end{claim}

The distribution $\sta_{u,v}$ is the probability of existence/absence
of edge $(u, v)$ once ${\cal M}$ has converged. This implies that after $N
= m\ln(1/\varepsilon)$ we are very close to the converged distribution of
\emph{all} edges, irrespective of the degrees of their terminating
vertices.

\begin{proof}
$\M{T}$ has two eigenvalues, $1$ and $1 - (\alpha(d_u, d_v)
  +\beta(d_u, d_v))$. Let the corresponding unit eigenvectors be
  $\be_1$ and $\be_2$. For notational convenience we will refer to
  $\alpha(d_u, d_v)$ and $\beta(d_u, d_v)$ as $\alpha$ and
  $\beta$. Since $\alpha + \beta > 0$, the eigenvectors form a
  basis. We represent the initial state of the MC as $\bv = c_1 \be_1
  + c_2 \be_2$. After $N$ applications of the transition matrix, the
  distribution of the Markov chain's state is
\begin{eqnarray}
 \bp & = & \M{T}^N \bv = c_1 \M{T}^N \be_1 + c_2 \M{T}^N \be_2 \nonumber \\
     & = & c_1 \be_1 + c_2 \left(1 - (\alpha +\beta) \right)^N \be_2.
\end{eqnarray}

Since $(1 - (\alpha + \beta)) < 1$, the second term decays with $N$,
and the stationary distribution asymptotically approaches $c_1
\be_1$. The 2-norm of the decaying term (the error incurred when we
stop at finite $N$) can be bounded for
$$N = \ln{(1/\epsilon)}/(\alpha+\beta)$$
in the following manner
\begin{eqnarray}
  \|\bp - \sta_{u,v}\|  & = &  \|(1-\gamma)^N c_2 \be_2\|_2
           \leq  (1-\gamma)^{\ln(1/\varepsilon)/\gamma}  c_2 \|\be_2\|_2  \nonumber \\
                        & \leq &  \exp(-\ln(1/\varepsilon)) = \varepsilon \nonumber
\end{eqnarray}
where we have used $\gamma = (\alpha + \beta)$ and
\begin{eqnarray}
  N & = & \frac{m}{2} \ln{\frac{1}{\epsilon}} \geq \frac{1}{\gamma}  \ln{\frac{1}{\epsilon}},
  \mbox{\hspace{3mm} for DD-preserving MC and} \nonumber \\
  N & = & m \ln{\frac{1}{\epsilon}} \geq \frac{1}{\gamma} \ln{\frac{1}{\epsilon}},
  \mbox{\hspace{3mm} for JDD-preserving MC.}
\label{eqn:N}
\end{eqnarray}
The key observation is $\gamma \geq 2/m$ for DD-preserving MC (see
values of $\alpha$ and $\beta$ in Sec.~\ref{sec:DD}) and $\gamma \geq
1/m$ for a JDD-preserving MC (see Sec.~\ref{sec:JDD}).
\end{proof}

\section{Tests with real graphs}
\label{sec:tests}

In this section we perform  empirical studies to find a suitable value of $\epsilon$, i.e., one
which can be used to generate samples of graphs which are
statistically indistinguishable from those generated by a smaller
$\epsilon$. We discriminate graphs using certain graphical metrics.

We run the MC with the DD-preserving rewiring scheme for $N$ steps
before saving the resultant graph. We start the MC using 4 real
graphs - the neural network of \emph{C. Elegans}~\cite{celegans1}
(referred to as ``C. Elegans''), the power grid of the Western states
of US~\cite{celegans1} (called ``Power''), co-authorship graph of
network science researchers~\cite{adjnoun} (referred to as
``Netscience'') and a 75,000 vertex graph of the social network at
Epinions.com~\cite{03ra3a} (``soc-Epinions1''). Their details are in
Table~\ref{tab:graphs}. The first three were obtained
from~\cite{newman_graphs} while the fourth was downloaded
from~\cite{snap}. All the graphs were converted to undirected graphs
by symmetrizing the edges.

We use $N = \{0.5, 2.5, 5, 7.5\}|E|$ corresponding to $\epsilon =
\{0.37, 6.7 \times 10^{-3}, 4.5 \times 10^{-5}, 3.06 \times
10^{-7}\}$. 1000 graphs are generated in this manner, starting with
the real graphs in Table~\ref{tab:graphs}. Thereafter, we calculate
the global clustering coefficient (three times the ratio of the number of
triangles to wedges in a graph)\cite{BaWe00},
the diameter, and the maximum eigenvalue of the Laplacian matrix  of the graph~\cite{Me94}
 for each of members of the ensemble of 1000 graphs. In Fig.~\ref{fig:dd},
we plot these distributions for three of the four graphs in
Table~\ref{tab:graphs}. We clearly see that for $\epsilon = 0.37$
i.e., $N = 0.5|E|$, the distributions for all graphical metrics are
quite different from the distributions obtained with $N > |E|$. The
red and black curves, corresponding to $N = 5|E|$ and $N = 7.5|E|$ are
mostly indistinguishable indicating that the MC has converged to its
stationary distributions. We will henceforth proceed with $\epsilon =
4.5 \times 10^{-5}$, achieved with $N = 5|E|$ when the DD is preserved
during the MC sampling.
\begin{table}[b]
  \caption{Characteristics of the graphs used in this paper. $(|V|,
    |E|)$ are the numbers of vertices and edges in the graph and G-R
    diagnostic is the Gelman-Rubin diagnostic.}
  \label{tab:graphs}
  \begin{ruledtabular}
    \begin{tabular}{lcr}
      \textrm{Graph name} & \textrm{ $(|V|, |E|)$} & \textrm{G-R diagnostic} \\ \colrule
      C. Elegans & (297, 4296)     &  1.05 \\
      Netscience & (1461, 5484)    &  1.02 \\
      Power      & (4941, 13188)   &  1.006 \\
      soc-Epinions1 & (75879, 405740)& 1.06  \\
    \end{tabular}
  \end{ruledtabular}
\end{table}

\begin{figure*}
  \centerline{\includegraphics[width=\textwidth, trim=1cm 7cm 2.cm  7cm,clip=true]
  {Fig1.pdf}}
  \caption{Plots of the distributions of the global clustering
    coefficient, the diameter, and the maximum eigenvalue of the
    graph Laplacian for ``C. Elegans'' (left), ``Netscience'' (middle)
    and ``Power'' (right), evaluated after $0.5|E|, 2.5|E|, 5|E|$ and
    $7.5|E|$ iterations of the Markov chain (green, blue, black and
    red lines respectively). The corresponding values of $\epsilon$
    are in the legend. We see that the distributions converge at
    $\epsilon \sim 1e-5$. In these runs, the DD was held constant
    across all graphical samples.}
  \label{fig:dd}
\end{figure*}

We now turn our attention to the case when the JDD is preserved during
MC sampling. We repeat the same sampling procedure as above, but with
the JDD-preserving rewiring scheme being driven by the Markov
chain. We choose the same sequence of $\epsilon$, corresponding to $N
= \{1, 5, 10, 15\}|E|$, per Eq.~\ref{eqn:N}. The distributions of
global clustering coefficients, graph diameter and maximum eigenvalue
of the graph Laplacian are in Fig.~\ref{fig:jdd}. Again, it is
 clear that the distribution (for any of the graphical
metrics) is far from convergence for $\epsilon = 0.37$ but
indistinguishable for the other values of $\epsilon$. In particular,
the red and black plots lie on top of each other in most of the
subfigures, indicating that $N = 10|E|$ ($\epsilon = 4.5 \times
10^{-5}$) is generally sufficient to drive the MC to its stationary
distribution, as characterized by the aforementioned graphical
parameters.
\begin{figure*}
  \centerline{\includegraphics[width=\textwidth, trim=1cm 7cm 2.cm  7cm,clip=true]
  {Fig2.pdf}}
  \caption{Plots of the distributions of the global clustering
    coefficient, the diameter, and the maximum eigenvalue of the
    graph Laplacian for ``C. Elegans'' (left), ``Netscience'' (middle)
    and ``Power'' (right), evaluated after $|E|, 5|E|, 10|E|$ and
    $15|E|$ iterations of the Markov chain (green, blue, black and
    red lines respectively). The corresponding values of $\epsilon$
    are in the legend. We see that the distributions converge at
    $\epsilon \sim 1e-5$. In these runs, the JDD was held constant
    across all graphical samples.}
  \label{fig:jdd}
\end{figure*}

In summary, a tolerance of $\epsilon = 4.5 \times 10^{-5}$ provides
converged/stationary distributions of global clustering coefficient,
graph diameter and maximum eigenvalue of the graph Laplacian. This
corresponds to running the MC on the space of graphs for $N = 5|E|$
steps when DD is preserved and for $N = 10|E|$ steps when JDD is
preserved. These results provide an empirical proof that the models in
Sec.~\ref{sec:DD} and Sec.~\ref{sec:JDD} provide a good \emph{a
priori} estimate of the steps to run a MC to obtain stationary
results.

\section{Verification of model assumptions}
\label{sec:verif}

In this section, we check our assumption that a MC on the space of
graphs can be approximated by coupled 2-state Markov chains. We also
devise a test to ensure that the MC does not converge to a local mode
in the distribution of graphs. This can occur in case of multimodal
distributions.

\subsection{Verifying edge-by-edge convergence}
\label{sec:ebe}

The equation for $N$ (Eq.~\ref{eqn:N}) is based on a heuristic and has
to be verified. Further, the equation is strictly valid for an edge
and it is unlikely that after $N$ steps of the MC, \emph{all} the edges
have decorrelated and converged to their stationary distribution. The
residual number of correlated edges have to be identified and their
impact on graphical properties quantified.

Consider the behavior of an arbitrary edge (between labeled vertices)
during a long walk executed by a MC on the space of graphs. During
this walk, the edge appears and disappears, thus describing a binary
time-series $\{Z_t\}$. If the MC on the space of graphs reaches a
stationary distribution in ${\rm O(|E|)}$ steps, then thinning the
binary time-series by a factor of ${\rm O(|E|)}$ should result in a
time-series that resembles independent draws of zeros and ones. We
develop a nonparametric test of independence for such time-series.
This method exploits Sokal's observation that if an MC reaches a
stationary distribution, the samples it draws are uncorrelated i.e.,
the autocorrelation is small~\cite{Sokal}. It is general and does not
depend on the preservation of DD or JDD in the graphs. While this
test has been in existence for some time~\cite{92rl2a,96ra1a}, it is
not well known in the network generation literature.

Assume that the time-series $\{Z_t\}$ is very long i.e., it takes $K
\gg N$ steps. The time-series will be correlated as observed
in~\cite{StPi12}. Consider too, a $k$-thinned chain $\{Z^k_t\}$,
obtained by retaining every $k^{th}$ member of $\{Z_t\}$.  The thinned
chain will have a far smaller autocorrelation and upon further
thinning will begin to resemble independent draws from a
distribution. If Eq.~\ref{eqn:N} is correct, $k = N$ should yield a
time-series that resembles independent draws \emph{more} than a
first-order Markov process. This forms the basis of our test of
independence. We fit models of two processes, independence and
first-order Markov, to the data and compute the
log-likelihood. The ``goodness'' of model fit is determined by computing
the Bayesian Information Criterion (BIC). This paper is the first
application of this technique to graphs, though it has been used
in other contexts.

Let $x_{ij}$ be the number of $(i, j), i, j \in (0, 1)$ transitions in
the $k-$thinned chain $\{Z^k_t\}$. The $x_{ij}$ are used to populate
$X$, a $2 \times 2$ contingency table. The table entries are
normalized by the length of the $k$-thinned chain ($K/k -1$) to
provide us with the empirical probabilities $p_{ij}$ of an $(i,j)$
transition in $\{Z^k_t\}$. Let $\widehat{p_{ij}}$ and
$\widehat{x_{ij}} = (K/k -1) \widehat{p_{ij}}$ be the predictions
of expected values of the table entries provided by a model. Then the
goodness-of-fit of the model can be provided by a likelihood ratio
statistic (called the $G^2$-statistic; Chapter 4.2 in~\cite{07bfh3a})
and a Bayesian Information Criterion (BIC) score:
\begin{eqnarray}
  G^2 & = & -2 \sum_{i=0}^{i=1} \sum_{i=0}^{i=1} x_{ij} \log\left( \frac{\widehat{x_{ij}}}{x_{ij}}\right),
  \nonumber \\
  BIC & = & G^2 + q \log\left( \frac{K}{k} - 1\right),
\label{eqn:ll}
\end{eqnarray}
where $q$ is the number of parameters in the model used to fit the
table data. Log-linear models are generally used to model tabular data
(Chapter 2.2.3 in~\cite{07bfh3a}). The log-linear models' predictions
for table entries generated by independent sampling and a first-order
Markov process are
\begin{eqnarray}
\log(p_{ij}^{(I)}) & = & u^{(I)} + u^{(I)}_{1,(i)} + u^{(I)}_{2, (j)}
 \mbox{\hspace{2mm} and \hspace{2mm}}  \nonumber \\
\log(p_{ij}^{(M)}) & = & u^{(M)} + u^{(M)}_{1,(i)} + u^{(M)}_{2, (j)} + u^{(M)}_{12,(ij)},
\label{eqn:llm}
\end{eqnarray}
where superscripts $I, M$ indicate an independent and a Markov process
respectively. The maximum likelihood estimates (MLE) of the parameters
($u^{(W)}_{b,(c)}$) are available in closed form in literature (Chapters 2.2.3 and
3.1.2 in~\cite{07bfh3a}) and we reproduce them below:
\begin{equation}
\widehat{x_{ij}^{I}}  =  \frac{(x_{i+}) (x_{+j})}{x_{++}}
\mbox{\hspace{2mm} and \hspace{2mm}}
\widehat{x_{ij}^{M}}  =  x_{ij},
\label{eqn:llmans}
\end{equation}
where $x_{i+}$ and $x_{+j}$ are the sums of the table entries in row
$i$ and column $j$ respectively. $x_{++}$ is the sum of all entries
(i.e., $K/k - 1$, the number of transitions observed in $\{Z^k_t\}$,
or the total number of data points). We compare the fits of the two
models thus:
\begin{eqnarray}
  \Delta BIC & =  &  BIC^{(I)} - BIC^{(M)} \nonumber \\
             & =  & -2 \sum_{i=0}^{i=1} \sum_{i=0}^{i=1} x_{ij}
               \log\left( \frac{\widehat{x_{ij}^{(I)}}}{x_{ij}}\right)
               - \log\left(\frac{K}{k} - 1\right). \nonumber \\
 \label{eqn:bic}
\end{eqnarray}
Above, we have substituted $\widehat{x_{ij}^{(M)}} = x_{ij}$ and the
fact that the log-linear model for a Markov process has one more
parameter than the independent sampler model.  Large BIC values
indicate a bad fit. A negative $\Delta BIC$ indicates that an
independent model fits better than a Markov model.

This test is applied as follows. We construct a thinned binary
time-series $\{Z^k_t\}$ for $k = |E|$ for each of the edges. The
$\Delta BIC$ is computed and edges with negative $\Delta BIC$ are
deemed to have become independent after $k$ steps of the Markov
chain. The test is repeated with a higher $k$ till all the edges are
deemed to have become independent. A C++ implementation of this test,
using a binary time-series as input, is available
at~\cite{papersoftware}.

However, a difficulty appears at this point. We see that the accuracy
of the log-linear models in Eq.~\ref{eqn:llm} and~\ref{eqn:llmans}
depend on the length of the chain being tested i.e., if the chain is
too short, the ratio of likelihoods and the BIC will be wrongly
calculated. Thus, before we populate the contingency table, we ensure
that the edge time-series is long enough.

Consider an edge time-series obtained from the $K$-step MC. One can
empirically compute the edge-mean $\overline{Z_t}$ i.e., the
probability of existence of the edge. An estimate of the same quantity can be
obtained from a $k-$thinned version of the chain i.e.,
$\overline{Z^k_t}$. If the thinned time-series are sufficiently long,
then the estimates of the edge mean, $\overline{Z^k_t}$, computed for
different values of $k$, approximately describe a normal distribution,
with mean $q$ (the edge mean) and a variance $\sigma^2$.

We wish to ensure that a $k$-thinned MC is sufficiently long, i.e.,
the empirical edge-mean $\overline{Z^k_t}$ can be estimated to lie
inside $q \pm r$ with confidence $s$ i.e.$ P( q - r \le
\overline{Z^k_t} < q + r ) = s$, or
\begin{equation}
 \left( \frac{r}{\Phi^{-1}\{ 0.5(1+s) \}} \right)^2 = \sigma^2
\label{eqn:sigma}
\end{equation}
where $\Phi$ is the cumulative distribution function for a standard
normal distribution and $r$ is a tolerance on the empirical edge-mean
calculated using a $k$-thinned MC.

Consider that thinning the $K$-step time-series by a factor $k^\prime$
renders it a first-order MC (see Appendix~\ref{ap:som} on how this may
be done). The contingency table entries (Eq.~\ref{eqn:llmans}) then
indicate the number (or proportions) of 0-1 and 1-0 transitions of the
2-state Markov chain model of the time-series. It is then trivial to
compute the entries $\alpha$ and $\beta$ of the Markov transition
matrix. If $n^\prime$ is the length of the $k^\prime$-thinned time-series, then
the variance of the $\overline{Z^k_t}$, $\sigma^2$ can be written in
terms of $\alpha, \beta$,
\[ \sigma^2 = \frac{\alpha \beta (2 - \alpha - \beta)}{n^\prime (\alpha +
  \beta)^3} \]
and so, using Eq.~\ref{eqn:sigma},
\begin{equation}
  n^\prime = \frac{\frac{\alpha \beta (2 - \alpha - \beta)}{(\alpha + \beta)^3}}
           {\left( \frac{r}{\Phi^{-1}\{ 0.5(1+s) \}} \right)^2},
\end{equation}
i.e., the $k^\prime$-thinned time-series must be at least $n^\prime$ in
length to provide an estimate $\overline{Z^k_t}$ of $q$ (computed from
the $K$-step time-series) within tolerance $r$ with the required
confidence $s$. In this paper, we use $r = 0.01$ and $s = 0.95$.

We now apply this test to a Markov chain on graphs initiated with
those mentioned in Table~\ref{tab:graphs}. In Fig.~\ref{fig:ind} we
see the thinning factors required to render the $k$-thinned
time-series $\{Z^k_t\}$ resemble independent draws from a
distribution. In Fig.~\ref{fig:ind} (left) we see that thinning by
$4|E|$ is sufficient to render $\{Z^k_t\}$ independent, for all three
graphs tested, when the DD is preserved by the MC. In
Fig.~\ref{fig:ind} (right) we see results for the JDD-preserving
counterpart. It is clear that higher thinning factors are required for
independence, but a thinning factor of $10|E|$ seems to be
sufficient. Thus empirically, the factor-of-two difference in $N$, as
predicted by the models in Eq.~\ref{eqn:N}, and as seen empirically in
Fig.~\ref{fig:dd} and Fig.~\ref{fig:jdd}, are corroborated by a
separate test that does not assume DD or JDD preservation.
\begin{figure*}[t]
  \centerline{
    \includegraphics[width=0.5\textwidth, trim=1cm 6cm 2.cm 6cm,clip=true] {Fig3a.pdf}
    \includegraphics[width=0.5\textwidth, trim=1cm 6cm 2.cm 6cm,clip=true] {Fig3b.pdf}}
  \caption{Fraction of edges testing independent, for ``C. Elegans'',
  ``Netscience'' and ``Power'' for various values of the thinning
  factor $k$. Left: Results for the case when DD is preserved. We see
  that with a thinning factor of $k = 4|E|$ all the edges have
  converged to their stationary distribution. Right: Results for the
  case when JDD was preserved. We see that $k = 10|E|$ ensures that at
  least 95\% of the edges become independent.}
  \label{fig:ind}
\end{figure*}

Convergence of distributions of global clustering coefficients,
diameter and maximum eigenvalues do not automatically indicate that
the time-series of \emph{all} the edges in a graph resemble
independent draws. This is clearly seen in Fig.~\ref{fig:jdd} for the
graph ``Power'' (bottom row). In Fig.~\ref{fig:jdd}, $N \sim 10|E|$
($\epsilon = 4.5\times10^{-5}$) was sufficient to construct
distributions that do not change appreciably for greater $N$. However,
from Fig.~\ref{fig:ind}, we see that about 2\% of the edges are still
not independent, but presumably close to being so, since about $15|E|$
steps of the MC result in all edges becoming independent.

In order to check this issue, we performed a test with the
soc-Epinions1 graph. We generate an ensemble of 1000 graphs using $N =
30|E|$ and the method in Sec.~\ref{sec:JDD} i.e., where the JDD of the
initial graph was preserved as samples were generated. We also run a
long $K$-step MC ($K = 21.6\times10^6|E|$) and compute the thinning
factor $k$ required to render each edge's time-series similar to
independent draws. Since the graph has 405,740 edges, the independence
test was performed for only 10\% of the edges, chosen randomly from
the initial graph. In Fig.~\ref{fig:soc_ind} we plot a histogram of
$k/|E|$. We clearly see that while by $k = 30|E|$ about 90\% of the
edges have become independent, there are quite a few ``pathological''
edges that are very far from becoming independent (one such edge
becomes independent at $k = 720|E|$).
In Fig.~\ref{fig:soc} we plot the
distributions for the diameter and the maximum eigenvalue for $N =
30|E|, 150|E|, 270|E|$ and $390|E|$. We see that the two
distributions, while not identical, do not change much. Thus
distributions of the chosen metrics are robust (insensitive) to an
ensemble which might contain graphs with a few ``hard-to-converge''
edges. Thus, if one is interested in generating proxy graphs where
anonymity is critical (e.g. proxies of email traffic in an
organization), an exhaustive, edge-by-edge checking of independence,
using the method described above, may be necessary. On the other hand,
when one only desires a collection of roughly independent graphs with
a prescribed DD or JDD, the far simpler approach of Sec.~\ref{sec:DD}
or~\ref{sec:JDD} may suffice.

\begin{figure}
  \includegraphics[width=0.5\textwidth, trim= 1cm 6cm 1cm 6cm, clip=true]{Fig4.pdf}
  \caption{A histogram of $k/|E|$, where $k$ is the number of steps
    the MC has to take to render the time-series of a particular edge
    resemble independent draws, calculated for the soc-Epinions1
    graph. We see that while $30|E|$ steps render about 90\% of the
    edges independent, there are a few edges which are still very far
    from becoming independent.These results are for a case where the
    JDD was preserved.}
\label{fig:soc_ind}
\end{figure}

\begin{figure*}
  \centerline{
    \includegraphics[width=0.5\textwidth, trim= 1cm 6cm 1cm 6cm, clip=true]{Fig5a.pdf}
    \includegraphics[width=0.5\textwidth, trim= 1cm 6cm 1cm 6cm, clip=true]{Fig5b.pdf}}
  \caption{Plots of the distribution of graph diameter (left) and
  maximum eigenvalue of the graph's Laplacian (right) for $N = 30|E|,
  150|E|, 270|E|$ and $390|E|$, for the soc-Epinions1 graph (green,
  blue, black and red lines respectively). We see that distributions
  are very similar, even though, from Fig.~\ref{fig:soc_ind}, it is
  evident that a few edges are far from being independent. Thus these
  two graphical metrics are not very sensitive graphs that are almost,
  but not quite, independent. These results are for a case where the
  JDD was preserved.}
  \label{fig:soc}
\end{figure*}
\subsection{The Gelman-Rubin test for convergence of Markov chains}

Finally we address the question whether the MC generating the ensemble
of graphs is sampling from a particular mode of the distribution of
graphs (in case the distribution is multimodal). This can, in
principle, occur since the MC is always started from the same starting
graph. Convergence to the stationary distribution (and not just to one
mode) can be tested by starting the MC from an overdispersed set of
points and checking whether the same distributions are realized,
irrespective of the starting point. In our case, we will choose three
points and perform this test only for the JDD-preserving case.

In order to generate a set of overdispersed starting graphs, we run a
Markov chain for $N = 10,000|E|$ steps, after initializing it with a
real network. This serves as the second over-dispersed starting point
for the MC and the initial condition for the third point (realized
after another 10,000$|E|$ steps). Three concurrent MC are initialized
with these graphs, and we calculate the Gelman-Rubin (G-R)
diagnostic~\cite{92gr1a} using the binary edge time-series. Values of
the diagnostic between 1 and 1.1 indicate that the states of the
concurrent Markov chain are not dependent on the starting location. We
performed this test for all 4 graphs; the corresponding G-R
diagnostics are tabulated in Table~\ref{tab:graphs}.

\section{Conclusions}
\label{sec:concl}

We have developed a method for generating independent realizations of
graphs which share a particular graphical property. In this paper, we
have demonstrated the method in the case where (1) the degree
distribution was preserved over all realizations and (2) where the
joint degree distribution was held constant. The graphs are generated
using a MC approach, which drives a graph ``rewiring'' scheme. The
rewiring scheme is responsible for preserving the shared graphical
properties when it generates a new graphical realization, from an old
one. Our method involves running the Markov chain for $N$ steps before
extracting a graph; the chain is run repeatedly to generate an
ensemble of realizations, given an initial graph. We have developed
models and closed form expressions for $N$ that allows a coupled
2-state Markov chain \emph{of an edge} to converge to its stationary
distribution. This is a necessary, but not sufficient, condition for
the Markov chain \emph{on the space of graphs} to ``forget'' the
starting graph and generate an independent realization. We find that a
Markov chain requires $5|E| - 30|E|$ steps to mix fully, i.e., to
provide converged distributions of graphical properties like the global
clustering coefficient, graph diameter and maximum eigenvalues of the
graph Laplacian. We find that the Markov chains that preserve the
degree distribution are easier to mix than those where the joint
degree distribution is held constant.

We verified our model (for $N$) by performing a test for the
independence of the Markov chain described by any edge between two
labeled vertices in the graph, as the Markov chain on the space of
graphs executes its random walk. This test is not dependent on any
heuristic or graphical properties. The test uses the time-series of
occurrence/non-occurrence of an edge between labeled vertices, thins
them by a factor of $N$, and fits it with a first-order Markov and an
independent sampling model. The goodness of fit of the two models,
determined by their individual BICs, is used to select between
them. The Markov chain is considered mixed for any $N$ for which the
independence sampling model is selected over the Markov model. We find
that even for $N$ that provide converged distributions of the
graphical properties above, a small fraction of the edges may still be
significantly correlated (as determined by the independence test). The
robustness of the graphical properties indicate that if fully
independent graphical realizations are desired, the laborious process
of testing each edge may be unavoidable.

Finally, we repeated our tests by starting parallel Markov chains from
widely dispersed starting points and computing the Gelman-Rubin
diagnostic for them. We find that the parallel chains converge to the
same stationary distribution. This check was performed to ensure that
our distributions of graphical properties was not being driven by the
starting point of the Markov chain.

Apart from the theoretical results summarized above, we also provide
an open-source C++ implementation of the non-parametric test, along
with sample problems to illustrate its use~\cite{papersoftware}.

\section*{Acknowledgments}
This work was funded by the Applied Mathematics program of the United States Department of Energy, Office
of Science and by an Early Career Award from the Laboratory Directed
Research \& Development (LDRD) program at Sandia National
Laboratories. Sandia National Laboratories is a multiprogram
laboratory managed and operated by Sandia Corporation, a wholly owned
subsidiary of Lockheed Martin Corporation, for the United States
Department of Energy's National Nuclear Security Administration under
contract DE-AC04-94AL85000.


\section{Introduction}

In the formation of low-mass stars, disks and jets play a major role in controlling mass  accretion
 and removing excess angular momentum.
Specifically, low-mass protostellar jets are expected to be launched and collimated via magneto-centrifugal processes at distances $<<$10 AU from Young Stellar Objects or YSOs \citep[e. g.][]{Konigl:2000p607}.
  However, the balance of forces in proximity to high-mass YSOs has
  been difficult to infer because high extinctions, clustering, and
  large distances hinder attempts to resolve the gas structure and
  dynamics at small radii ($<$100 AU) where outflows are launched and
  collimated from accretion disks.

Radio source I in Orion BN/KL  is the nearest example of a high-mass YSO.
The close distance and exceptional richness in molecular maser tracers enabled the best detailed picture of  gas dynamics  in a massive YSO  on scales comparable with the Solar System.
In particular, detailed imaging of  vibrationally-excited states SiO
masers with the Very Long Baseline Array (VLBA) enabled resolving the
surface of an edge-on disk as well as a wide-angle outflow in the
collimation and launch region at radii 10-100~AU
(\citealt{Matthews:2010p781}; Figure~\ref{Obs_figure}). Additionally,
at larger radii (100-1000 AU), the wide-angle wind appears to
collimate in a straight bipolar flow, as traced by H$_2$O and $\nu$=0 SiO
maser and thermal emission \citep{Plambeck:2009p10531, Greenhill:2012p10712}. Measurement of
the 3D velocity fields of SiO masers provides a clear evidence of rotation  in both the disk and the wind along the same axis.
\cite{Cunningham:2005p8732}.
explored the importance of stellar winds from a luminous central source and proposed the idea that SiO masers in Source I may arise
from the interaction of a spherical wind with an infalling envelope. However, \cite{Matthews:2010p781} argued that such a model could not
explain in detail the gas dynamics and physical conditions as revealed by SiO masers and alternately proposed
that the SiO masers arise from material
ablated from the surface of an accretion disk and discussed
qualitatively the applicability of various classes of disk-wind models
to Source I.
They find evidence to support that magnetic fields may play a role in driving and/or shaping the disk-wind from Source I, based on details of gas dynamics on small scales, like curved and helical trajectories traced by maser motions.
Assuming a magnetic field is present, a magneto-centrifugally launched wind \citep[e.g.,][]{Blandford:1982p892,
  Konigl:2000p607, Fendt:2006p574} would be a natural candidate for powering a disk-wind from Source I. The latter thus  becomes an ideal laboratory to extend the MHD models from solar-like stars to  massive YSOs.
 Recently, \cite{Vaidya:2011p8992} explored the interplay of
 radiation forces and magnetic fields to show that radiative line driving
 from luminous stars with masses $>30\,{\rm M}_{\sun}$ does modify the dynamics
 of outflows. However, since the luminosity for a $10\,{\rm M}_{\sun}$ binary Orion Source I would be less than a single $20\,{\rm M}_{\sun}$,
we adopt that model in the ideal MHD regime
 (neglecting radiative force) to numerically simulate the dynamics of
 the disk-wind inside 100 AU.

\begin{figure} \centering
\includegraphics[width=0.8\columnwidth]{f1.pdf}
\caption{Velocity field (in km/s) of the SiO $v$=1,2 $J$=1-0 maser emission in Orion Source I as observed
with the VLBA \citep{Matthews:2010p781}. The stellar velocity is $\sim$6 km s$^{-1}$.
The bulk of the SiO emission is located within four
"arms" of an X-shaped pattern, while a $\sim$14 AU thick band with no SiO emission harbors the 7mm continuum source \citep{Reid:2007p9519, Goddi:2011p8577}. The continuum emission traces an ionized disk
while the SiO masers probe material ablated from the surface of the disk in a  wind.}
\label{Obs_figure}
\end{figure}

\section{Model Assumptions}
\label{modelA}
 \begin{figure}
 \centering
 \includegraphics[width=1.0\columnwidth]{f2.pdf}\\
 \caption{Cartoon Model for Source I.}
 \label{SourceI_cartoon}
 \end{figure}
The initial conditions adopted to model the disk-wind from Orion Source I are reported in Table~\ref{ParaSummary} and are derived  from either observations and/or modeling of the system.
Source I is believed to be an eccentric, hard  binary with a total
mass of $2\times10$ M$_{\odot}$,
based on studies of the dynamics of YSOs and molecular gas in Orion BN/KL, as well as N-body simulations \citep{Goddi:2011p8577,Bally:2011p9419}.
In particular, proper motion measurements of radio sources  provided strong evidence that a dynamical interaction occurred about 500~years ago between massive YSOs Source~I and BN  \citep[e. g.,][]{Goddi:2011p8577}.
 N-body numerical simulations showed that the dynamical interaction
 between a binary of two 10\,${\rm M}_{\sun}$\,stars (Source I) and a single 10\,${\rm M}_{\sun}$\,star
  (BN) may lead to simultaneous ejection of both stars  and the BN/KL outflow
as well as binary hardening,  while preserving the original circumbinary disk  around Source I
 \citep{Bally:2011p9419, Moeckel:2011p8574}.
Mechanical energy conservation (and virial theorem) implies an orbit semi-major axis $a_0=
 2-5$ AU, depending on the estimate of the kinetic energy of the
 runaway stars and the molecular gas flow \citep{Goddi:2011p8577,
   Bally:2011p9419}.
For simplicity, we neglect the eccentricity
  and assume that the system is an equal-mass circular
binary, comprised of two $10\,{\rm M}_{\sun}$ stars  orbiting in a
circle of radius $a_0 = 3.5$\,AU.

A high-resolution study of H$\alpha$ and [O I]$\lambda$6300 line profiles in a pre-main-sequence spectroscopic binary has shown evidence of a bipolar jet as being launched by the whole binary  as opposed to each star individually  \citep{Mundt:2010p5}. Similarly, we assume  that the binary members  orbit inside the inner radius  of a truncated circumbinary disk and that  the wind  is launched at this radius.
The cartoon figure representing our model assumption is shown in Fig.~\ref{SourceI_cartoon}.

To study the wind launching from Source I,
we perform axi-symmetric ideal MHD simulations in cylindrical
coordinates (r,$\phi$,z)  using the PLUTO code \citep{Mignone:2007p644}.
We essentially prescribe a hydrostatic corona
(density and pressure) threaded with force-free magnetic field.
An adiabatic equation of state, $P=\frac{\gamma-1}{\gamma}\rho^{\gamma}$, with $\gamma = 5/3$, relates the pressure with the density in the flow.
For the central gravity, we treat the system of equal mass binary members as an equivalent $5\,{\rm M}_{\sun}$ (reduced
mass) object at their barycenter. The accretion disk which is treated as a boundary forms the
launching base of the wind \cite[e.g.,][]{Ouyed:1997p634}
Further, we choose the wind launching point to be at a distance equal to the radius of the binary orbit,
i.e. $r_0$ = $a_0$, where the material is rotating with a
sub-Keplerian velocity $V_{\phi} = \chi v_0$.
Further, the disk is considered to be
denser than the hydrostatic corona and this density contrast is
prescribed by $\delta$. In a realistic (thick and hot)
disk, the central gravity is balanced radially by contributions from thermal
pressure and centrifugal rotation. This can be ensured by
consistently deriving the sub-Keplerianity($\chi$) based on the choice
of density contrast. For example, smaller values of $\delta$ implies a
hotter (denser) disk and a larger contribution
of thermal pressure resulting in the disk to rotate with sub-Keplerian
speeds to maintain radial balance. We chose $\delta \sim 3$ as a
reference value implying the underlying disk to rotate at a speed 0.8
times the Keplerian speed.

The magnetic field strength is prescribed by the value of plasma beta (the ratio of thermal pressure to
the magnetic pressure) at $r_0$, $\beta_0$.
In order to put constraints on the magnetic field
strengths in Source I, we study the wind
properties with varying $\beta_0$   (see Sect.~\ref{magfpara}).
Finally, we assume the number density at $r_0$ to be  $n_0$ = $3 \times 10^{10} \rm cm^{-3}$.
\cite{Reid:2007p9519} argued the  disk in Source I to be ionized
with a typical electron density of $\sim 10^{7}\rm cm^{-3}$ (their "Model C"),  which corresponds to an
ionization fraction of $\sim 10^{-3}$.

\begin{table}
\begin{center}
\caption{Parameters chosen for
  the reference run to study MHD driven outflow from Orion Source I. }
\begin{tabular}{ l l }
\hline\hline
\noalign{\smallskip}
Quantity & Reference Value \\
\hline
\noalign{\smallskip}
Masses of Binary[$m_{1}$, $m_{2}$]   & $10\,{\rm M}_{\sun}, 10\,{\rm M}_{\sun}$\\
\noalign{\smallskip}
Orbital Separation [2$a_0$] & 7\,AU\\
\noalign{\smallskip}
Reduced Stellar Mass [$\mu_{*}$]& $5\,{\rm M}_{\sun}$\\
\noalign{\smallskip}
Wind Launching radius [$\rm{r}_0$] & 3.5\,AU\\
\noalign{\smallskip}
Keplerian Velocity at $\rm{r}_0$ [$\rm{v}_0$] & 35.7 km s$^{-1}$\\
\noalign{\smallskip}
Volume Density at $\rm{r}_0$ [$\rho_0$] & $5.0\times10^{-14}$
g\,cm$^{-3}$\\
\noalign{\smallskip}
Plasma Beta at $\rm{r}_0$ [$\beta_0$] & 1.0 \\
\noalign{\smallskip}
Density Contrast [$\delta$] & 2.85 \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
\end{tabular}\label{ParaSummary}
\end{center}
\end{table}

The flow is freely allowed to leave the outer boundaries (r,z = 50
$a_0$, 150 $a_0$) and symmetric boundary conditions are applied on the r=0 axis.
The poloidal component of the magnetic field at the outer boundaries
is set to avoid artificial collimation by ensuring zero current
($\nabla \times \vec{B} = 0$; \citealt{Ustyugova:1999p655, Porth:2010p1041}).
 The critical boundary at z=0 mimics the inflow of gas and magnetic flux from the disk surface.
Here, we ensure the disk to be in pressure equilibrium with the
medium but with a higher density (i.e. colder than the medium).
The angular speed of the magnetic field line in the disk is equal to the Keplerian
value and the poloidal velocity is updated each time with the value at
the grid cell just above the boundary.
 Further, ideal MHD approximation permits to set the flow velocity along the poloidal magnetic field ($v_{\rm  p} \parallel B_{\rm p}$).
Thus, the mass flux at the disk boundary is not fixed but consistently derived by causal
interaction between the outflowing gas and the boundary
value. Detailed description of these boundary conditions are given in
\citet{Vaidya:2011p8992}.

\section{Results and Discussion}
\subsection{MHD driven wind profiles}

\begin{figure*} \centering
\includegraphics[width=1.9\columnwidth]{f3.pdf}
\caption{Number density (left), poloidal velocity
  (middle) and azimuthal velocity (right) of the reference run in
  steady-state. Temperature contours (black) for T =
  500\,K, 1000\,K, 2000\,K and 5000\,K are overlaid in the left
  panel. The middle panel shows the positions of
  three MHD critical surfaces - magneto-sonic surface (yellow),
  Alfv$\acute{e}$n surface (magenta), magneto-fast surface (black).
  White arrows in the right panel represents the velocity
  vectors which are parallel to magnetic field lines shown by solid
  red lines.}
\label{phyvelref}
\end{figure*}

At the onset of our reference run, material at the base of the disk is ejected in vertical
direction along the magnetic field. This material propagates outwards
producing bow shocks in the ambient medium.
As this shock propagates outside of the domain, the flow collimates due to magnetic forces
until a steady equilibrium is attained, beyond which no further
evolution of magnetic field lines takes place.
The steady-state values of the density, poloidal and azimuthal velocities of the wind are shown in
Figure~\ref{phyvelref}, along with the magnetic field lines (in red), the poloidal velocity versors (white
arrows), temperature contours and MHD critical surfaces. The magnetic
field lines ideally mark the morphology of the outflow.
The left panel of Figure~\ref{phyvelref}  shows the steady-state
number density of the fluid, whose values range from 10$^{8}$\,cm$^{-3}$ to 10$^{11}$\,cm$^{-3}$. Overlaid on the density
map are the temperature contours, with values ranging from 5000\,K closer to the central source to 100\,K at large
distances. The high density and temperatures close to the axis reflect the fact that the initial hydrostatic conditions are fixed at the
boundary inside of the wind launching radius (this ensures numerical stability of the flow very close to the axis).
The middle panel of Figure~\ref{phyvelref} shows the  poloidal  wind velocity, whose values range from 8.0 to 32.0 km s$^{-1}$.
The lines over-plotted show  three critical surfaces in a steady-state MHD flow,
viz. the magneto-sonic surface (yellow), the magneto-fast surface (black),
and the Alfv$\acute{e}$n surface (red).
The flow starts sub-magneto-sonically from the disk and is
accelerated to super fast speeds mainly via Lorentz forces.
The Alfv$\acute{e}$n surface marks the locus
of points in the flow where ideally the kinetic energy of the flow is
converted to magnetic energy. This results in bending of the field
lines thereby generating a toroidal field and corresponding hoop stress
which collimates the flow.
The right panel of Figure~\ref{phyvelref} finally shows the azimuthal component of the wind velocity, with values ranging from 10
km\,s$^{-1}$ at the outer edge to 25 km\,s$^{-1}$ close to the inner
wind launching point.
The transfer of angular momentum from the wind to the medium may impart rotation to the latter, as observed in SiO masers (Sec.~\ref{siocomp}).

\subsection{Comparison between MHD wind models and  wind properties observed from SiO masers}
\label{siocomp}
We propose that the SiO masers in Orion Source I
are excited as a moderately rotating disk-wind interacts with the ambient
medium via slow MHD shocks and transfer of angular momentum.
Similar mechanisms have been successfully applied to explain the excitation of
$\rm H_{2}$O  masers \citep{Kaufman:1996p9607}.
Here, we qualitatively compare the steady-state properties of a
single-fluid, rotating MHD wind (launched from an underlying accretion disk)
with the observed physical and dynamical properties of the outflow molecular component probed by SiO masers in  Source I.

Three physical properties of the wind from the MHD modeling are consistent with SiO maser observations:
temperature/density profiles, morphology, and velocity field.

  The conditions for optimum excitation of the SiO maser emission are  n(H$_2$)$\sim 10^{9} \rm cm^{-3}$ and T$\lesssim 2000\,\rm K$, respectively, as derived
 from a radiative transfer analysis by \citet{Goddi:2009p8571}. In our
 simulations, the regions close to the inner launching point (z,r
 $\lesssim$  15\,AU) have values of densities (10$^{10-11}$ cm$^{-3}$) and temperatures (2000-5000 K)  beyond the critical values required for SiO maser excitation
 (Fig.~\ref{phyvelref}, left panel). This is  consistent with the 14 AU thick "dark band" seen in
 Fig.~\ref{Obs_figure} where no masers are observed. The intermediate region (15\,AU $<$ z,r $<$
 60\,AU) has typical densities of $10^{9} \ \rm cm^{-3}$ and $T \sim 1000-2000\,K$ (Fig.~\ref{phyvelref}, left panel). Here, the wind can provide a
 strong IR background source for pumping the SiO masers in the shocks
 produced by its interaction with the ambient molecular medium.

In the steady-state, the magnetic collimation and
 acceleration result in an open funnel shaped
 structure of field lines  (Fig.~\ref{phyvelref}).
 Such a configuration of magnetic flux surfaces
 naturally results in a disk-wind with wide-angle bipolar morphology, which
 resembles the X-shaped wind structure around Source I revealed by SiO maser emission (Fig.~\ref{Obs_figure}).

Regarding the velocity field, the dynamical properties of molecular outflows can be extrapolated from the single-fluid MHD wind by analyzing the momentum transfer \citep{Masson:1993p9661}. In this study, the density of the MHD driven wind  exceeds the mean density of its surrounding medium. Such dense and cold winds carry a large amount of momentum and impart a significant fraction of it to the ambient medium, resulting in the  shocked  material to have velocities of the same order as the  wind itself  \citep{Masson:1993p9661}, i.e. in the range of 10 - 30 km\,s$^{-1}$ (Fig.~\ref{phyvelref}, middle panel).
Consistently, \cite{Matthews:2010p781} measured the 3D velocities of the molecular gas flows around
Source I, to vary in the range from 5 km s$^{-1}$ to
25 km s$^{-1}$, with an average of 14 km s$^{-1}$.
In addition, SiO masers show rotation seen as a velocity gradient in the
radial velocities of the four arms of the "X" tracing the disk-wind,
which is consistent with the trend obtained for $V_{\phi}$ (Figure~\ref{phyvelref}, right panel).

Therefore, the simulated MHD winds can drive molecular flows by interacting via
shocks with the ambient medium.

\begin{figure} \centering
\includegraphics[width=0.9\columnwidth]{f4.pdf}
\caption{Radial profiles of the magnetic field [in mG] at z =
  20\,AU [solid] and z = 60\,AU [dashed] for the reference run, obtained after the MHD wind has reached a steady-state.}
\label{phyBtot}
\end{figure}

\subsection{Magnetic field strength
\label{magfpara}
In our models, the magnetic field strength is governed by $\beta_0$.
We have carried out simulation runs with three different
values of $\beta_0$ viz. 1.0, 3.0 and 10.0.
Figure~\ref{phyBtot}
shows the radial profiles of  magnetic field strength at z = 20\,AU
and z = 60\,AU for different values of $\beta_0$ and $\delta$=
2.85.
We find that the magnetic field strength varies in the range 0.2--0.8 G (depending on $\beta_0$) and decreases
with radius and/or vertical height.
From this study, we observe that the changes in density and
temperature are not significant with varying $\beta_0$ (e.g., $\Delta T \sim 100-200\,K$ when $\beta_0$ increases by ten fold).
On the other hand, poloidal and azimuthal velocities decrease substantially with
increasing values of $\beta_0$. Since lower field strengths imply lower
Lorentz force, the flow is not efficiently accelerated resulting in
lower velocities. The typical velocities obtained from $\beta_0 > 3.0$
are not enough to explain the observed velocities in Source I.
Therefore, we assume as a reference the magnetic field strength corresponding to
$\beta_0 = 1$, $\sim$0.35 G (at r, z =40\,AU, 20\,AU).
Assuming equipartition, \citet{Matthews:2010p781} quote a magnetic field strength of $\sim$0.3\,G for
Source I, based on kinematics of SiO masers. This value is consistent
with a measurement  of OH 1665 MHz masers across a $\sim 10^4$ AU region around Source I, which provided field strengths of 1.8 to 16.3 mG \citep{Cohen:2006p9565}.
Since the OH masers are believed to arise from material with
$n(H_2)\sim 10^7$ cm$^{-3}$ \citep{Gray:1992p9570} and the magnetic
field strength is expected to scale as the square root of the gas
density,  a simple extrapolation from the large scales (probed by OH masers) to the small scales (probed by SiO masers)
gives field strengths of order the value predicted by the equipartition assumption.
These empirical estimates  lie comfortably within the range predicted by our
model, but direct measurements via the Zeeman effect could provide a more reliable comparison  between models and observations.

\subsection{Model Limitations}
\label{limitations}

The main qualitative features of a steady MHD driven disk-wind
are consistent with the physical and dynamical properties of the molecular flow probed by SiO masers.
A more quantitative analysis requires however further theoretical and observational studies. In this section, we list the
limitations of the present model and the scope for future work.

(i) The assumption of treating the disk as a boundary allows us to follow wind dynamics till
   a steady state is reached, but neglects the impact of  disk
  dynamics on wind propagation. A more realistic model should allow
   both disk and   outflow to  evolve simultaneously
 \cite[e.g.,][]{Casse:2002p1142, Sheikhnezami:2012p10641}. Also, we
 neglect any feedback on the disk-wind from a stellar wind (Source I has however  a rather modest $L_{\rm bol} \sim 10^{4}\,L_{\odot}$;
 \citealt{Gezari:1998p10535}).

(ii) In our ideal MHD model, we have assumed a single fluid MHD driven wind and have ignored effects due to diffusion of
  neutral molecules across the magnetic field.
Though the regions
  close to the star are hot and potentially have a high
  degree of ionization, consistent with an ideal approximation,
  a more accurate model should incorporate effects due
  to diffusion between the neutral and ionized medium.
  Future simulations will relax the ideal MHD approximation and
   explore  effects of magnetic diffusivity.

(iii) We have ignored any feedback due to molecular chemistry.
  Ideally, a direct comparison between MHD models and maser observations   would require a   detailed multi-fluid chemical model. Recently, conditions for molecular survival in MHD
disk-winds  were computed using chemical evolution of a self-similar
MHD wind model \citep{Panoglou:2012p10039}.
However, such a chemical model
  with micro-physics of maser excitation is quite complex. As a follow-up, we plan to include non-equilibrium cooling and separately follow the evolution of
neutral and ionized gas fraction in the flow.

(iv) Including  in the modeling the elements discussed above,  requires measurement of relevant physical quantities within the launching region ($\lesssim
  100$\,AU), i.e. ionization degree, molecular abundances, and
  magnetic field strengths.
High resolution and sensitive observations of  Hydrogen radio recombination lines (with JVLA and ALMA) will help us constrain the dynamical and physical properties of ionized gas. Further, radiative transfer analysis of rotational molecular lines (from ALMA) can constrain molecular gas properties and abundances. Finally, polarimetric observations of SiO masers (with VLBI) will enable to measure magnetic field strength and orientation at the disk/jet interface.

\section{Summary}
We have explored the plausibility of MHD driven wind
model to explain the disk-wind properties in Orion Source I as probed by SiO masers. In
particular, we have performed axisymmetric MHD wind launching
simulations treating Source I as an equal mass
binary and choosing initial parameters  derived from observations.
We propose that the SiO masers are excited as the MHD driven wind interacts with the ambient
molecular medium in form of shocks.
Four elements of consistency support this hypothesis:

\begin{enumerate}
\item
Our model consistently reproduces the bipolar morphology of the disk-wind traced by the SiO masers.

\item
The poloidal and azimuthal velocities from the simulations are in the range of 10 -- 30 km s$^{-1}$, which
agree very well with the observed average 3D velocity of 14 km s$^{-1}$ and with the rotational component of SiO masers.

\item
The physical conditions in the modeled wind are consistent with the requirements for
pumping SiO masers.

\item
We predict  field strengths of 0.2-0.8 G within 100 AU, in agreement
with empirical estimates from SiO masers.
\end{enumerate}

We conclude that the main  dynamical features of the SiO masers disk-wind in Orion Source I can be
 explained by the MHD  wind model.
In a broader context,  Source I has provided a unique laboratory
 to extend the standard MHD models from low-mass to high-mass YSOs, indicating that  in both regimes
similar physics may describe the transfer of matter and angular momentum between disks and outflows.

\section*{Acknowledgments}
We thank Christian Fendt and Lynn Matthews for their valuable
comments. B.V thanks the STFC for funding

\bibliographystyle{mn2e}