#  QCD effects in non-QCD theories

Vladimir Dzhunushaliev<EMAIL_ADDRESS>Department of Theoretical
and Nuclear Physics, Al-Farabi Kazakh National University, Almaty 050040,
Kazakhstan  Institute of Experimental and Theoretical Physics, Al-Farabi
Kazakh National University, Almaty 050040, Kazakhstan  Academician J. Jeenbaev
Institute of Physics of the NAS of the Kyrgyz Republic, 265 a, Chui Street,
Bishkek 720071, Kyrgyzstan  Vladimir Folomeev<EMAIL_ADDRESS>Institute of
Experimental and Theoretical Physics, Al-Farabi Kazakh National University,
Almaty 050040, Kazakhstan  Academician J. Jeenbaev Institute of Physics of the
NAS of the Kyrgyz Republic, 265 a, Chui Street, Bishkek 720071, Kyrgyzstan

###### Abstract

It is shown that, in some non-QCD theories, there are effects shared by QCD:
(i) in SU(2) Yang-Mills theory containing a nonlinear spinor field, there is a
mass gap; (ii) in SU(3) Proca-Higgs theory, there are flux tube solutions with
a longitudinal electric field required for producing a force binding quarks;
(iii) in non-Abelian Proca-Higgs theories, there exist flux tube solutions
with a momentum directed along the tube axis and particlelike solutions with a
nonvanishing total angular momentum created by crossed color electric and
magnetic fields; in QCD, such configurations may contribute to the proton
spin. We discuss the conjecture that such non-QCD theories might be a
consequence of approximate solution of an infinite set of Dyson-Schwinger
equations describing the procedure of nonperturbative quantization. The
phenomenon of dimensional transmutation for nonperturbative quantization and
the analogy between nonperturbative quantization and turbulence modeling are
also discussed.

mass gap, Proca theory, flux tube and particlelike solutions, angular
momentum, nonperturbative quantization, closure problem

###### pacs:

12.90.+b, 11.90.+t

##  I Introduction

In quantum chromodynamics (QCD), there are many effects caused by the strong
nonlinearity of QCD. This means that such effects must be described using
nonperturbative quantization of strongly nonlinear non-Abelian Yang-Mills
fields. These effects are, for instance, (i) the appearance of a mass gap,
(ii) the existence of tubes connecting quarks and filled with a longitudinal
electric field, (iii) the presence of a contribution to the proton spin coming
from gluon fields, etc. The presence of the mass gap indicates that, for
quantized Yang-Mills fields, there exists a lower limitation on the mass of
particles. The existence of a mass gap in QCD is not an exceptional
phenomenon. For example, in the BCS theory of superconductivity, there is an
energy gap whose existence is caused by the energy profit for a pair of
electrons creating a Cooper pair (evaluating the energy gap for
superconductors can be found in textbooks, see, e.g., Ref. [1]). In the Nambu-
Jona-Lazinio model, there appears the so-called equation for a mass gap
according to which the mass spectrum of mesons starts not from zero but from
some nonzero value which is called a mass gap (see, e.g., the review [2]). In
turn, the existence of the aforementioned tubes is related to the confinement
phenomenon (for a detail review of the confinement problem, see, e.g., Ref.
[3]): such a tube provides so strong attraction between quarks that they are
never found in isolation. It may be noted here that, in the MIT bag model,
there is a string approximation for a flux tube connecting quarks and
containing a longitudinal color electric field (for a review of this model,
see Ref. [4]).

For its part, it was experimentally observed that there exist contributions to
the proton spin coming from quark spins, their orbital angular momenta, as
well as gluon fields [5, 6]. Nearly all these phenomena are described using
lattice calculations (for a detail review on lattice calculations see Ref.
[7]). However, in a certain sense, lattice calculations are similar to an
experiment: they give a certain result that one would like to explain using
some analytic calculations.

In the present paper we wish to discuss the possibility that, in some non-QCD
theories, there may exist effects shared by QCD (here, by QCD, we mean SU(3)
gauge theory with quarks). The appearance of such QCD effects in non-QCD
theories enables one to assume that such non-QCD theories may have some
relevance to QCD. For example, it is possible that such non-QCD theories may
serve as some approximate way to describe nonperturbative effects in QCD. This
would imply that, in nonperturbatively quantized QCD, some approximations
leading to non-QCD theories are possible.

As concrete examples, we will consider here, firstly, the possibility of the
appearance of a mass gap in SU(2) Yang-Mills theory where color fields are
sourced by a nonlinear spinor field (Sec. II). Second, in Sec. III, we will
demonstrate the possibility of obtaining various flux tube solutions (within
SU(3) Proca-Higgs theory) and particlelike solutions with a nonzero total
angular momentum (within SU(2) Proca-Higgs theory). Such configurations may be
used in modeling the confinement phenomenon and contribute to the proton spin.

##  II Mass gap

In this section we demonstrate that the energy spectrum of a “monopole” in
SU(2) Yang-Mills theory with a source in the form of a nonlinear spinor field
has a global minimum corresponding to a mass gap. We have enclosed the word
“monopole” in quotation marks because the asymptotic behavior of the
corresponding solution (sourced by a nonlinear spinor field) differs from that
of the ’t Hooft-Polyakov monopole solution (sourced by a triplet of Higgs
scalar fields). (Detailed description of the properties of magnetic monopoles
can be found, for example, in Ref. [8], and for a review on topological
methods in QCD and of the literature in the field see Ref. [9]). Also, we will
discuss the possibility of the appearance in the theory of a nonlinear spinor
field as a consequence of some approximate description of nonperturbative
effects in non-Abelian quantum field theory.

###  II.1 Energy spectrum of the “monopole”

Here, we closely follow Refs. [10, 11]. The Lagrangian describing a system
consisting of a non-Abelian SU(2) field $A^{a}_{\mu}$ interacting with
nonlinear spinor field $\psi$ can be taken in the form

$\mathcal{L}=-\frac{1}{4}F^{a}_{\mu\nu}F^{a\mu\nu}+i\hbar
c\bar{\psi}\gamma^{\mu}D_{\mu}\psi-
m_{f}c^{2}\bar{\psi}\psi+\frac{\Lambda}{2}g\hbar
c\left(\bar{\psi}\psi\right)^{2}.$ (1)

Here $m_{f}$ is the mass of the spinor field;
$D_{\mu}=\partial_{\mu}-i\frac{g}{2}\sigma^{a}A^{a}_{\mu}$ is the gauge-
covariant derivative, where $g$ is the coupling constant and $\sigma^{a}$ are
the SU(2) generators (the Pauli matrices);
$F^{a}_{\mu\nu}=\partial_{\mu}A^{a}_{\nu}-\partial_{\nu}A^{a}_{\mu}+g\epsilon_{abc}A^{b}_{\mu}A^{c}_{\nu}$
is the field strength tensor for the SU(2) field, where $\epsilon_{abc}$ (the
completely antisymmetric Levi-Civita symbol) are the SU(2) structure
constants; $\Lambda$ is a constant; $\gamma^{\mu}$ are the Dirac matrices in
the standard representation; $a,b,c=1,2,3$ are color indices and
$\mu,\nu=0,1,2,3$ are spacetime indices.

Using the Lagrangian (1), one can find the corresponding field equations

$\displaystyle D_{\nu}F^{a\mu\nu}=$ $\displaystyle\frac{g\hbar
c}{2}\bar{\psi}\gamma^{\mu}\sigma^{a}\psi,$ (2) $\displaystyle
i\hbar\gamma^{\mu}D_{\mu}\psi-m_{f}c\psi+\Lambda
g\hbar\psi\left(\bar{\psi}\psi\right)=$ $\displaystyle 0.$ (3)

In order to find “monopole” solutions and their energy spectrum, we choose the
following Ansätze for the Yang-Mills and spinor fields:

$\displaystyle A^{a}_{t}$ $\displaystyle=$ $\displaystyle 0,\quad
A^{a}_{i}=\frac{1}{g}\left(1-f\right)\begin{pmatrix}0&\phantom{-}\sin\varphi&\sin\theta\cos\theta\cos\varphi\\\
0&-\cos\varphi&\sin\theta\cos\theta\sin\varphi\\\
0&0&-\sin^{2}\theta\end{pmatrix}\,\text{with}\,\,i=r,\theta,\varphi\text{ (in
polar coordinates)},$ $\displaystyle\psi^{T}$ $\displaystyle=$
$\displaystyle\frac{e^{-i\frac{Et}{\hbar}}}{gr\sqrt{2}}\begin{Bmatrix}\begin{pmatrix}0\\\
-u\\\ \end{pmatrix},\begin{pmatrix}u\\\ 0\\\
\end{pmatrix},\begin{pmatrix}iv\sin\theta e^{-i\varphi}\\\ -iv\cos\theta\\\
\end{pmatrix},\begin{pmatrix}-iv\cos\theta\\\ -iv\sin\theta e^{i\varphi}\\\
\end{pmatrix}\end{Bmatrix},$

where $E/\hbar$ is the spinor frequency and the functions $f,u$, and $v$
depend on the radial coordinate $r$ only. Substituting these Ansätze in Eqs.
(2) and (3), one can find the corresponding equations for the functions $f,u$,
and $v$ (for details see Ref. [10]).

In turn, the total dimensionless energy density of the “monopole” under
consideration is

$\tilde{\epsilon}=\frac{1}{\tilde{g}^{2}_{\Lambda}}\left[\frac{{f^{\prime}}^{2}}{x^{2}}+\frac{\left(f^{2}-1\right)^{2}}{2x^{4}}\right]+\left[\tilde{E}\frac{\tilde{u}^{2}+\tilde{v}^{2}}{x^{2}}+\frac{\left(\tilde{u}^{2}-\tilde{v}^{2}\right)^{2}}{2x^{4}}\right],$
(4)

written in terms of the dimensionless variables $x=r/\lambda_{c}$,
$\tilde{u}=u\sqrt{\Lambda/\lambda_{c}g},\tilde{v}=v\sqrt{\Lambda/\lambda_{c}g},\tilde{E}=\lambda_{c}E/(\hbar
c),\tilde{g}^{2}_{\Lambda}=g\hbar c\lambda_{c}^{2}/\Lambda$, where
$\lambda_{c}=\hbar/(m_{f}c)$ is the Compton wavelength. The prime denotes
differentiation with respect to $x$. Correspondingly, the total dimensionless
energy of the “monopole” can be calculated using the formula

$\tilde{W}_{t}\equiv\frac{\lambda_{c}g^{2}}{\tilde{g}^{2}_{\Lambda}}W_{t}=4\pi\int\limits_{0}^{\infty}x^{2}\tilde{\epsilon}dx,$
(5)

where the energy density $\tilde{\epsilon}$ is taken from Eq. (4). The
corresponding profile of the total energy is given in Fig. 1.

Figure 1: A sketch of the energy spectrum of the total energy from Eq. (5) as
a function of the spinor frequency $\tilde{E}$ (for details see Ref. [10]).

###  II.2 Outline of nonperturbative quantization and possible appearance of a
nonlinear spinor field

As we see from Fig. 1, the energy of the “monopole” supported by the nonlinear
spinor field has a minimum which corresponds to a mass gap. In this connection
the natural question arises: whether there exists a relation between the mass
gap obtained and a mass gap in QCD? In other words, whether the appearance of
the mass gap for such a “monopole” solution is coincidental or this is related
somehow to nonperturbative quantization of essentially nonlinear non-Abelian
Yang-Mills fields? If one starts from the conjecture that such a relation does
really exist, it may consist in the fact that the nonlinear spinor field
occurs as a consequence of some approximate description of quantum
nonpertubative effects in QCD.

This may happen as follows. According to W. Heisenberg [12], the procedure of
nonperturbative quantization consists in writing out quantum averages of the
Yang-Mills and Dirac operator equations. But since these equations contain
Green’s functions of different orders, to close the system, one has to write
out equations for higher-order Green’s functions, and so on, ad infinitum. As
a result, this leads to an infinite set of Dyson-Schwinger-type equations for
all Green’s functions. This has been discussed in detail in Ref. [13], where a
technique based on truncating the Dyson-Schwinger equations for quantum
$\phi^{4}$ field theories is presented.

In a more rigorous mathematical language, this looks like

$\displaystyle
i\gamma^{\mu}\left\langle\hat{\psi}_{;\mu}-m\hat{\psi}\right\rangle=$
$\displaystyle 0,$ (6)
$\displaystyle\left\langle\hat{F}^{a\mu\nu}\right\rangle_{;\nu}=$
$\displaystyle-4\pi\left\langle\hat{j}^{a\mu}\right\rangle,$ (7)
$\displaystyle\left\langle\hat{\bar{\psi}}\left(i\gamma^{\mu}\hat{\psi}_{;\mu}-m\hat{\psi}\right)\right\rangle=$
$\displaystyle 0,$ (8)
$\displaystyle\left\langle\hat{A}_{\nu}\left(i\gamma^{\mu}\hat{\psi}_{;\mu}-m\hat{\psi}\right)\right\rangle=$
$\displaystyle 0,$ (9)
$\displaystyle\left\langle\hat{F}^{a\mu\nu}\hat{\psi}\right\rangle_{;\nu}=$
$\displaystyle-4\pi\left\langle\hat{j}^{a\mu}\hat{\psi}\right\rangle,$ (10)
$\displaystyle\ldots$ (11)

Here, the two first equations are the quantum-averaged Dirac and Yang-Mills
operator equations, and all subsequent ones are gotten by multiplying by the
corresponding operators with subsequent quantum averaging. It must be
mentioned that these equations are well known in perturbative quantum field
theory, but there they are always written down using Feynman diagrams, and
this leads to a fundamental difference between the Dyson-Schwinger equations
in perturbative and nonperturbative regimes. The same equations have been
obtained in Ref. [14], using only another notations for Green’s functions.
Notice also that, within this approach, it is of course necessary to study the
problem of unitary evolution for these equations.

It is of great importance to note that these equations greatly simplify in the
stationary case. Apparently, in this case all Green’s functions will be
regular, without involving Dirac delta functions, step functions, etc. This
situation differs in principle from the case of quantum fields changing with
time where one must take into account the finiteness of the velocity of light;
this results in the appearance of singular Green’s functions. Such a situation
occurs, for example, in perturbative quantum field theory where Green’s
functions are singular functions occurring as a result of the fact that in
such a case Green’s functions are polylinear combinations of propagators. In
turn, the propagators are used for a mathematical description of dynamical
processes of a motion and of propagation of quanta.

It is evident that such an infinite set of equations can be solved neither
analytically nor numerically. Therefore, it is necessary to have some
approximate methods of solving Eqs. (6)-(11). One of such methods consists in
using the procedure that enables one to cut off an infinite set of equations
to get a finite one. In this case it is necessary to have a hypothesis about
the behavior of higher-order Green’s function. In turbulence modeling, this
problem is called the closure problem. We will discuss the relation between
nonperturbative quantization and turbulence modeling in more detail below.

Consider a possible scenario of the occurrence of nonlinearity of a spinor
field. One of the equations in the system (6)-(11) contains a nonlinearity
describing the interaction between quarks and the non-Abelian gauge field,

$\mathcal{L}_{i}\sim\left\langle\hat{\bar{\psi}}\hat{A}^{a}_{\mu}\hat{\psi}\right\rangle.$
(12)

If we wish to cut off the infinite set of equations (6)-(11) at just the point
where one may derive separate equations for the gauge field $\left\langle
A^{a}_{\mu}\right\rangle$ and for the spinor field
$\left\langle\psi\right\rangle$, we have to have a hypothesis about an
approximate description of the interaction (12). For this purpose, one can
suggest different approximations for (12), one of which is

$\left\langle\hat{\bar{\psi}}\gamma^{\mu}\hat{A}^{a}_{\mu}\hat{\psi}\right\rangle=\left\langle\hat{\bar{\psi}}\gamma^{\mu}\left(\left\langle\hat{A}^{a}_{\mu}\right\rangle+\widehat{\delta
A^{a}_{\mu}}\right)\hat{\psi}\right\rangle\approx\left\langle\hat{\bar{\psi}}\gamma^{\mu}\hat{\psi}\right\rangle\left\langle
A^{a}_{\mu}\right\rangle+\lambda\begin{cases}\left\langle\hat{\bar{\psi}}\hat{\psi}\right\rangle\left\langle\widehat{\delta
A^{a}_{\mu}}\widehat{\delta A^{a\mu}}\right\rangle,\\\
\left\langle\hat{\bar{\psi}}\hat{\psi}\right\rangle^{2},\\\
\left\langle\hat{\bar{\psi}}\gamma^{\mu}\hat{\psi}\right\rangle\left\langle\hat{\bar{\psi}}\gamma_{\mu}\hat{\psi}\right\rangle,\end{cases}$
(13)

where the last term on the right-hand side represents possible corrections to
the first term needed for an approximate description of the quantum average
given on the left-hand side of this expression; the closure constant $\lambda$
appears here for matching dimensions of all terms. The physical meaning of
such an approximation is that the quantum average of the term describing the
interaction of two fields can be approximately represented as a sum of two
terms. The first term represents the product of two averages of these fields,
and the second one – a nonperturbative quantum correction.

It is of great interest that, in this approximation, there occurs the new
dimensional constant $\lambda$. The presence of such a dimensional constant
is, perhaps, an example of the fact that there is dimensional transmutation in
a nonperturbative regime of quantization, by analogy with the phenomenon of
dimensional transmutation occurring in a perturbative regime when a
dimensional constant appears in a theory where it was initially absent [15].

##  III Flux tubes and particlelike configurations in Proca theories

In this section, we will show that, in Proca theories, there are such QCD
effects:

  * •

A tube filled with a longitudinal non-Abelian electric field. Such a tube is a
necessary ingredient of the confinement phenomenon, since it provides a
constant force between quarks preventing their separation.

  * •

A tube filled with crossed color electric and magnetic fields that give a
nonzero Poynting vector, and hence a momentum directed along the tube axis.
Three such tubes, binding quarks together in a proton, will contribute
correspondingly to the proton spin.

  * •

Particlelike solutions possessing a nonvanishing angular momentum created by
static crossed non-Abelian electric and magnetic fields. It is possible that
this may have a connection with the nature of the angular momentum of gluon
fields contributing to the proton spin.

###  III.1 Tubes with the flux of chromoelectric field

In Ref. [16], within non-Abelian Proca theory with a scalar field, there were
obtained cylindrically symmetric solutions describing infinite tubes with the
flux of a color electric field. Such tubes can be used for a description of
the confinement phenomenon, since they contain a longitudinal electric field
producing a constant attractive force between quarks and preventing their
separation.

The Lagrangian describing such a system consisting of a non-Abelian SU(3)
Proca field $A^{a}_{\mu}$ interacting with nonlinear scalar field $\phi$ can
be taken in the form (in this section we work in units such that $c=\hbar=1$)

$\mathcal{L}=-\frac{1}{4}F^{a}_{\mu\nu}F^{a\mu\nu}-\frac{\left(\mu^{2}\right)^{ab,\mu}_{\phantom{ab,}\nu}}{2}A^{a}_{\mu}A^{b\nu}+\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi+\frac{\lambda}{2}\phi^{2}A^{a}_{\mu}A^{a\mu}-\frac{\Lambda}{4}\left(\phi^{2}-M^{2}\right)^{2}.$
(14)

Here
$F^{a}_{\mu\nu}=\partial_{\mu}A^{a}_{\nu}-\partial_{\nu}A^{a}_{\mu}+gf_{abc}A^{b}_{\mu}A^{c}_{\nu}$
is the field strength tensor for the Proca field, where $f_{abc}$ are the
SU(3) structure constants, $g$ is the coupling constant, $a,b,c=1,2,\dots,8$
are color indices, $\mu,\nu=0,1,2,3$ are spacetime indices. The Lagrangian
(14) also contains the arbitrary constants $M,\lambda,\Lambda$ and the Proca
field “mass tensor” $\left(\mu^{2}\right)^{ab,\mu}_{\phantom{ab,}\nu}$. We
have been using the words “mass tensor” in quotation marks because these
quantities are not actually masses. The quantity
$\left(\mu^{2}\right)^{ab,\mu}_{\phantom{ab,}\nu}$ corresponds to the case of
modified Proca theory where the mass term $m^{2}A^{a}_{\mu}A^{a\mu}$ is
replaced by
$\left(\mu^{2}\right)^{ab,\mu}_{\phantom{ab,}\nu}A^{a}_{\mu}A^{b\nu}$. The
literature in the field offers a few ways for modifying Proca theory (see,
e.g., Ref. [17] where generalized multi-Proca theory with the mass term
$M^{ab}A^{a}_{\mu}A^{b\mu}$ is under consideration). Following this approach,
we have generalized Proca theory by introducing the matrix
$\left(\mu^{2}\right)^{ab,\mu}_{\phantom{ab,}\nu}$ which is a tensor with the
spacetime indices $\mu,\nu$.

Making use of Eq. (14), the corresponding field equations can be written in
the form

$\displaystyle D_{\nu}F^{a\mu\nu}-\lambda\phi^{2}A^{a\mu}=$
$\displaystyle-\left(\mu^{2}\right)^{ab,\mu}_{\phantom{ab,}\nu}A^{b\nu},$ (15)
$\displaystyle\Box\phi=$ $\displaystyle\lambda
A^{a}_{\mu}A^{a\mu}\phi+\Lambda\phi\left(M^{2}-\phi^{2}\right).$ (16)

The required tube is described by the Ansätze

$A^{2}_{t}=\frac{h(\rho)}{g},\;A^{5}_{z}=\frac{v(\rho)}{g},\;A^{7}_{\varphi}=\frac{\rho
w(\rho)}{g},\;\phi=\phi(\rho),$ (17)

where $\rho,z,$ and $\varphi$ are cylindrical coordinates. For such a choice
of the SU(3) Proca field potentials, we have the following electric and
magnetic field intensities:

$\displaystyle E^{2}_{\rho}=-\frac{h^{\prime}}{g},\quad
E^{5}_{\varphi}=-\frac{\rho hw}{2g},\quad E^{7}_{z}=\frac{hv}{2g},$
$\displaystyle H^{2}_{\rho}=-\frac{vw}{2g},\quad H^{5}_{\varphi}=-\frac{\rho
v^{\prime}}{g},\quad
H^{7}_{z}=\frac{1}{g}\left(w^{\prime}+\frac{w}{\rho}\right),$

where the prime denotes differentiation with respect to $\rho$. The most
interesting fact for us is that the tube contains the longitudinal electric
field $E^{7}_{z}\neq 0$ producing the required attraction between quarks
described by the Proca field.

Substituting the potentials (17) in Eqs. (15) and (16), one can get the
corresponding equations for the functions $h$, $v$, $w$, and $\phi$ (for
details see Ref. [16]).

A sketch of the distributions of color electric and magnetic fields is given
in Fig. 2. As we see, there is the longitudinal electric field $E^{7}_{z}$
ensuring the attractive force between quarks and preventing their separation.
Correspondingly, the tube contains a finite flux of this field,

$\Phi_{z}=2\pi\int\limits_{0}^{\infty}\rho E^{7}_{z}d\rho<\infty,$

whose presence is necessary for the creation of a constant attractive force
between quarks.

Figure 2: A sketch of the distributions of color electric fields
$E^{2}_{\rho}$, $E^{5}_{\varphi}$, $E^{7}_{z}$ and of color magnetic fields
$H^{2}_{\rho}$, $H^{5}_{\varphi}$, $H^{7}_{z}$.

###  III.2 Tube with the energy flux/momentum density

In this case we take the Ansätze

$A^{5}_{t}=\frac{f(\rho)}{g},\;A^{5}_{z}=\frac{v(\rho)}{g},\;A^{7}_{\varphi}=\frac{\rho
w(\rho)}{g},\;\phi=\phi(\rho),$ (18)

which give the following components of the electric and magnetic field
intensities:

$\displaystyle E^{2}_{\varphi}=$ $\displaystyle\frac{\rho fw}{2g},\quad
E^{5}_{\rho}=-\frac{f^{\prime}}{g},\quad$ (19) $\displaystyle H^{2}_{\rho}=$
$\displaystyle-\frac{vw}{2g},\quad H^{5}_{\varphi}=-\frac{\rho
v^{\prime}}{g},\quad
H^{7}_{z}=\frac{1}{g}\left(w^{\prime}+\frac{w}{\rho}\right).$ (20)

In this case the Poynting vector $S^{i}$ is already nonzero,

$S^{i}=\frac{\epsilon^{ijk}}{\sqrt{\gamma}}E^{a}_{j}H^{a}_{k}\neq 0,$ (21)

where $\gamma$ is the determinant of the space metric. Substituting the
potentials (18) in Eqs. (15) and (16), one can get the corresponding equations
for the functions $f$, $v$, $w$, and $\phi$ (for details see Ref. [16]). In
turn, substitution of the components of electric and magnetic fields (19) and
(20) in Eq. (21) yields the following expression for the Poynting vector:

$S^{z}=\frac{1}{g^{2}}\left(\frac{df}{d\rho}\frac{dv}{d\rho}+\frac{1}{4}fvw^{2}\right).$
(22)

The graph of the momentum density $S^{z}$ is schematically sketched in the top
panel of Fig. 3. Its nonzero value is caused by the presence of the crossed
electric, $E^{2}_{\varphi}$, and magnetic, $H^{2}_{\rho}$, fields, as well as
the electric, $E^{5}_{\rho}$, and magnetic, $H^{5}_{\varphi}$, fields, whose
typical distributions are depicted in the bottom panels of Fig. 3.

Figure 3: Sketches of the $z$-component of the Poynting vector $S^{z}$ given
by Eq. (22), of the color electric fields $E^{2}_{\varphi}$ and
$E^{5}_{\rho}$, of the color magnetic fields $H^{2}_{\rho},H^{5}_{\varphi}$,
and $H^{7}_{z}$.

###  III.3 Proca balls with angular momentum

In this subsection we follow Ref. [18]. The Lagrangian describing a system
consisting of a non-Abelian SU(2) Proca field $A^{a}_{\mu}$ coupled to a
triplet of real Higgs scalar fields $\phi^{a}$ can be taken in the form

$\mathcal{L}=-\frac{1}{4}F^{a}_{\mu\nu}F^{a\mu\nu}-\frac{1}{2}\left(\mu^{2}\right)^{ab,\mu}_{\phantom{ab,}\nu}A^{a}_{\mu}A^{b\nu}+\frac{1}{2}D_{\mu}\phi^{a}D^{\mu}\phi^{a}-\frac{\Lambda}{4}\left(\phi^{a}\phi^{a}-v^{2}\right)^{2}.$
(23)

Here
$F^{a}_{\mu\nu}=\partial_{\mu}A^{a}_{\nu}-\partial_{\nu}A^{a}_{\mu}+g\epsilon_{abc}A^{b}_{\mu}A^{c}_{\nu}$
is the field strength tensor for the Proca field, where $\epsilon_{abc}$ are
the SU(2) structure constants (the completely antisymmetric Levi-Civita
symbol), $g$ is the coupling constant, $a,b,c=1,2,3$ are SU(2) color indices,
$\mu,\nu=0,1,2,3$ are spacetime indices; $v$ and $\Lambda$ are arbitrary
constants, and $\left(\mu^{2}\right)^{ab,\mu}_{\phantom{ab,}\nu}$ is the Proca
field mass tensor, which we suppose to be symmetric with respect to the color
and spacetime indices. This tensor corresponds to the case of modified Proca
theory in the same sense as described after Eq. (14).

Using Eq. (23), the corresponding field equations can be written in the form

$\displaystyle
D_{\nu}F^{a\mu\nu}+\left(\mu^{2}\right)^{ab,\mu}_{\phantom{ab,}\nu}A^{b\nu}$
$\displaystyle=$ $\displaystyle g\epsilon^{abc}\phi^{b}D^{\mu}\phi^{c},$
$\displaystyle D_{\mu}D^{\mu}\phi^{a}$ $\displaystyle=$
$\displaystyle-\Lambda\phi^{a}\left(\phi^{b}\phi^{b}-v^{2}\right).$

We choose the Ansätze for the Proca and scalar fields in the form

$A^{1}_{t}=\frac{f(\rho,z)}{g},\;A^{1}_{\varphi}=\frac{\rho\,k(\rho,z)}{g},\;A^{3}_{t}=\frac{h(\rho,z)}{g},\;A^{3}_{\varphi}=\frac{\rho\,w(\rho,z)}{g},\;\phi=\left\\{\phi_{1}(\rho,z),0,\phi_{3}(\rho,z)\right\\}$

written in cylindrical coordinates $\\{t,\rho,\varphi,z\\}$. For such a
choice, there are the following nonvanishing color electric and magnetic
fields (physical components):

$\displaystyle E^{1}_{z}=$ $\displaystyle-\frac{f_{,z}}{g},\quad
E^{3}_{z}=-\frac{h_{,z}}{g},\quad E^{1}_{\rho}=-\frac{f_{,\rho}}{g},\quad
E^{3}_{\rho}=-\frac{h_{,\rho}}{g},$ (24) $\displaystyle H^{1}_{z}=$
$\displaystyle-\frac{\rho\,k_{,\rho}+k}{g\rho},\quad
H^{3}_{z}=-\frac{\rho\,w_{,\rho}+w}{g\rho},\quad
H^{1}_{\rho}=\frac{k_{,z}}{g},\quad H^{3}_{\rho}=\frac{w_{,z}}{g},$ (25)

where a comma in lower indices denotes differentiation with respect to the
corresponding coordinate.

The components of the electric field strength given in Eq. (24) can be even or
odd functions of the coordinate $z$. For this reason, we consider here systems
of two types: the cases where the component $E^{1}_{z}$ is an odd/even
function and $E^{1}_{\rho}$ is an even/odd function. We will refer to such
configurations as the type A/B Proca balls (or P-balls), respectively.

For the sake of simplicity, consider the case where

$h=-f,k=-w,\phi_{1}=\phi_{3}=\phi.$ (26)

The corresponding partial differential equations describing such a system can
be found in Ref. [18]. In turn, for the field strengths (24) and (25), the
expression (21) gives the following nonvanishing physical component of the
Poynting vector:

$S_{\varphi}=\frac{2}{g^{2}}\left[f_{,\rho}\left(w_{,\rho}+\frac{w}{\rho}\right)+f_{,z}w_{,z}\right].$

Making use of this expression, one can obtain the expression for a linear
angular momentum density,

$\mathcal{P}_{z}=2\pi\int\limits_{0}^{\infty}S_{\varphi}\rho^{2}d\rho,$ (27)

which describes a distribution of the angular momentum density along the tube
axis. The corresponding distributions are exemplified in Fig. 4. It is seen
from this figure that, for the type A P-balls, the linear angular momentum is
always directed in one direction; this has the result that such
configurations, unlike the type B P-balls, have a nonzero total angular
momentum.

Figure 4: Schematic sketch of the linear angular momentum density
${\mathcal{P}}_{z}$ defined in Eq. (27) [with the component of the Proca field
mass tensor $\mu_{f}^{2}\equiv\left(\mu^{2}\right)^{11,t}_{\phantom{ab,}t}$]
along the tube axis for the P-balls under consideration. The arrows show the
directions of the angular momentum density.

###  III.4 Nonperturbative quantization and the hypothetical appearance of a
mass term

For the above flux tube and particlelike solutions to exist, the presence in
the system of a massive Yang-Mills field (i.e., a Proca field) is absolutely
necessary. In this subsection we would like to discuss the question of
whether, in solving Eqs. (6)-(11) approximately, there can occur additional
terms (in particular, a mass-like term) which might play the role of an
effective mass in Proca theory.

The only way to solve the infinite set of equations (6)-(11) approximately is
the use of the cut-off method for obtaining a finite set of equations.
According to this approach, it is necessary to keep only the first $n$
equations. In this case, the last equation will contain higher-order Green’s
functions ($G_{n+1}$), equations for which are already ruled out from a
consideration. In order that the first $n$ equations form a closed system, it
is necessary to propose a hypothesis about the last Green’s function
$G_{n+1}$. For instance, one can assume that such Green’s function is some
polylinear combination of lower-order Green’s functions. For example, for a
fourth-order Green’s function, this can be represented as

$\begin{split}&\left\langle\hat{A}^{a}_{\alpha}(x_{1})\hat{A}^{b}_{\beta}(x_{2})\hat{A}^{c}_{\gamma}(x_{3})\hat{A}^{d}_{\delta}(x_{4})\right\rangle\approx\left\langle\hat{A}^{a}_{\alpha}(x_{1})\right\rangle\left\langle\hat{A}^{b}_{\beta}(x_{2})\right\rangle\left\langle\hat{A}^{c}_{\gamma}(x_{3})\right\rangle\left\langle\hat{A}^{d}_{\delta}(x_{4})\right\rangle\\\
&+\left(\mu^{2}\right)^{ab}_{\phantom{ab,}\alpha\beta}\left\langle\hat{A}^{c}_{\gamma}(x_{3})\right\rangle\left\langle\hat{A}^{d}_{\delta}(x_{4})\right\rangle+\left(\text{
transmutation }a,b,c,d;\alpha,\beta,\gamma,\delta\right)+C,\end{split}$ (28)

where $\hat{A}^{a}_{\alpha}=\left\langle
A^{a}_{\alpha}\right\rangle+\widehat{\delta A^{a}_{\alpha}}$, the constants
$\left(\mu^{2}\right)^{ab}_{\phantom{ab,}\alpha\beta}$ and $C$ are called
closure constants. Here, the lower line is a nonperturbative quantum
correction to the first term in the right-hand side of the expression (28). An
assumption of this sort enables one to cut off the infinite set of equations
(6)-(11) and to obtain appropriate solutions describing nonperturbative
effects in QCD.

We may notice here a remarkable analogy between such procedure of
nonperturbative quantization and what happens in turbulence modeling when
analysing the Navier-Stokes equation

$\rho\left(\frac{\partial v_{i}}{\partial t}+v_{j}\frac{\partial
v_{i}}{\partial x_{j}}\right)=-\frac{\partial p}{\partial
x_{i}}+\frac{\partial t_{ij}}{\partial x_{j}},$

where $v_{i}$ is the flow velocity, $\rho$ is the fluid density, $p$ is the
pressure, $t_{ij}=2\mu s_{ij}$ is the viscous stress tensor,
$s_{ij}=\frac{1}{2}\left(\frac{\partial v_{i}}{\partial x_{j}}+\frac{\partial
v_{j}}{\partial x_{i}}\right)$, $\mu$ is molecular viscosity, and the velocity
$\vec{v}$ is a fluctuating quantity. Averaging this equation, one can obtain
[19]

$\rho\frac{\partial V_{i}}{\partial t}+\rho V_{j}\frac{\partial
V_{i}}{\partial x_{j}}=-\frac{\partial p}{\partial
x_{i}}+\frac{\partial}{\partial x_{j}}\left(2\mu S_{ji}-\overline{\rho
v_{j}^{\prime}v_{i}^{\prime}}\right),$ (29)

where the mean velocity $V_{i}$ is defined according to the equation
$v_{i}(\vec{x},t)=V_{i}(\vec{x},t)+v^{\prime}_{i}(\vec{x},t)$. Eq. (29) is
usually referred to as the Reynolds-averaged Navier-Stokes equation. The
quantity $\overline{\rho v_{j}^{\prime}v_{i}^{\prime}}$ is known as the
Reynolds-stress tensor and denoted as

$\tau_{ij}=-\overline{\rho v_{i}^{\prime}v_{j}^{\prime}}.$

One can obtain the following equation for the Reynolds stress tensor:

$\begin{split}\frac{\partial\tau_{ij}}{\partial
t}+V_{k}\frac{\partial\tau_{ij}}{\partial x_{k}}=-\tau_{ik}\frac{\partial
V_{j}}{\partial x_{k}}-\tau_{jk}\frac{\partial V_{i}}{\partial
x_{k}}+2\mu\overline{\frac{\partial v^{\prime}_{i}}{\partial
x_{k}}\frac{\partial v^{\prime}_{j}}{\partial
x_{k}}}-\overline{p^{\prime}\left(\frac{\partial v_{i}^{\prime}}{\partial
x_{j}}+\frac{\partial v_{j}^{\prime}}{\partial
x_{i}}\right)}+\frac{\partial}{\partial
x_{k}}\left(\nu\frac{\partial\tau_{ij}}{\partial
x_{k}}+C_{ijk}\right),\end{split}$

where $\nu$ is kinematic viscosity and

$C_{ijk}=\overline{\rho v^{\prime}_{i}\rho v^{\prime}_{j}\rho
v^{\prime}_{k}}+\delta_{jk}\overline{p^{\prime}v^{\prime}_{i}}+\delta_{ik}\overline{p^{\prime}v^{\prime}_{j}}.$

Thus we have six new equations, one for each independent component of the
Reynolds-stress tensor $\tau_{ij}$. However, we have also generated new
unknowns $\overline{\rho v^{\prime}_{i}\rho v^{\prime}_{j}\rho
v^{\prime}_{k}}$, $\overline{\mu\frac{\partial v^{\prime}_{i}}{\partial
x_{k}}\frac{\partial v^{\prime}_{j}}{\partial x_{k}}}$,
$\overline{v^{\prime}_{i}\frac{\partial p^{\prime}}{\partial x_{j}}}$, and
$\overline{v^{\prime}_{j}\frac{\partial p^{\prime}}{\partial x_{i}}}$. In the
same way one can obtain equations for these new unknown cumulants, and so on,
ad infinitum. In order to cut off this infinite set of equations, one uses the
procedure known as the closure problem in turbulence modeling.

Thus, returning to the closure equation (28), we see that the approximate
approach to solving the infinite set of equations of nonperturbative
quantization (6)-(11) leads to the appearance of the closure constants, one of
which may be a mass of a Yang-Mills field. It is quite conceivable that in the
limit $n\rightarrow\infty$ some of such constants will remain. This means that
the nonperturbative quantization leads to the appearance of dimensional
constants which are absent in the initial classical theory; this phenomenon is
called dimensional transmutation [15].

##  IV Discussion and conclusions

In this study, we have shown that some QCD effects can occur in non-QCD
theories which differ from QCD by the presence, for example, a nonlinear
spinor field or a massive Yang-Mills field (a Proca field). For this purpose,
we have discussed hypothetical possibilities of the appearance of such effects
in QCD as a consequence of quantization of essentially nonlinear Yang-Mills
fields. We have found out that when one uses the approximate approach to
solving an infinite set of Dyson-Schwinger equations describing the procedure
of nonperturbative quantization, the question arises as to how one may
describe a higher-order Green’s function in the cut-off set of equations.

In the present paper we have considered two possible ways of solving this
problem related to (i) the use of the nonlinearity describing the interaction
of gauge and spinor fields and (ii) the nonlinearity of the gauge field
itself:

  * •

Green’s function describing the interaction between gauge and spinor fields
(for instance,
$\left\langle\hat{\bar{\psi}}\gamma^{\mu}\hat{A}^{a}_{\mu}\hat{\psi}\right\rangle$)
is approximately described as the sum of two terms. The first term breaks up
the quantum average of the product of two fields as the product of the
averages of these fields, i.e., it is the product of the two-point Green’s
function of the spinor field
$\left\langle\bar{\hat{\psi}}\hat{\psi}\right\rangle^{2}$ by the average of
the gauge field $\left\langle\hat{A}^{a}_{\mu}\right\rangle$. The second term
is the correction, and it consists of the product of two factors: the first
one is the square of the two-point Green’s function
$\left\langle\bar{\hat{\psi}}\hat{\psi}\right\rangle^{2}$, and the second one
is the two-point Green’s function of the Yang-Mills fields
$\left\langle\widehat{\delta A^{a}_{\mu}}\widehat{\delta
A^{a\mu}}\right\rangle$, i.e., it is the dispersion of this field.

  * •

Green’s function of the product of Yang-Mills field potentials is a
combination of lower-order Green’s functions of the gauge field, and the so-
called closure constants appear, among which there may be masses of the gauge
field and $\Lambda_{\text{QCD}}$ (a characteristic parameter coming from QCD).

As concrete examples, we have demonstrated that there exist the following QCD
effects in non-QCD theories:

  * •

In SU(2) Yang-Mills theory with a source in the form of a nonlinear spinor
field, a mass gap is present.

  * •

In non-Abelian Proca theories, there are:

    * –

Infinite tubes containing a longitudinal electric field (see Sec. III.1).

    * –

Infinite tubes with the momentum density created by crossed electric and
magnetic fields and directed along the tube axis (see Sec. III.2).

    * –

Particlelike solutions possessing a nonzero total angular momentum created by
crossed electric and magnetic fields (see Sec. III.3).

    * –

For all the above solutions, the Meissner-like effect is observed: the Proca
fields are expelled by the Higgs fields.

In connection with the conjecture made here, there arise the following
questions requiring more detailed investigations:

  * •

What happens with the closure constants in the limit $n\rightarrow\infty$? If
the closure constants survive in this limit, this may lead to the appearance
of dimensional constants which are absent in the initial classical theory.
This means that there will take place dimensional transmutation for
nonperturbative quantization.

  * •

Whether there exists any generalization of the procedure of renormalization
for the procedure of nonperturbative quantization described in Sec. II.2?

  * •

How deep is the analogy between nonperturbative quantization and turbulence
modeling pointed out in Sec. III.4?

Thus the main result of the present study is that we make the conjecture that
the approximate approach to the procedure of quantization of strongly
interacting fields (the procedure of nonperturbative quantization) may lead to
the appearance of either new nonlinearities or new terms in the corresponding
equations.

## Acknowledgments

The work was supported by the program No. AP14869140 (The study of QCD effects
in non-QCD theories) of the Ministry of Education and Science of the Republic
of Kazakhstan. We are also grateful to the Research Group Linkage Programme of
the Alexander von Humboldt Foundation for the support of this research. We
wish to thank the anonymous referees whose comments helped to improve the
manuscript.

## References

  * [1] M. Tinkham, Introduction to Superconductivity (McGraw-Hill, 1996).
  * [2] S. P. Klevansky, The Nambu-Jona-Lasinio model of quantum chromodynamics, Rev. Mod. Phys. 64, 649 (1992).
  * [3] J. Greensite, An introduction to the confinement problem, Lect. Notes Phys. 821, 1 (2011).
  * [4] P. Hasenfratz and J. Kuti, The Quark Bag Model, Phys. Rept. 40, 75 (1978).
  * [5] J. Ashman et al. [European Muon Collaboration], A Measurement of the Spin Asymmetry and Determination of the Structure Function g(1) in Deep Inelastic Muon-Proton Scattering, Phys. Lett. B 206, 364 (1988) .
  * [6] J. Ashman et al. [European Muon Collaboration], An Investigation of the Spin Structure of the Proton in Deep Inelastic Scattering of Polarized Muons on Polarized Protons, Nucl. Phys. B 328, 1 (1989).
  * [7] C. Ratti, Lattice QCD and heavy ion collisions: a review of recent progress, Rept. Prog. Phys. 81, 084301 (2018).
  * [8] Ya. Shnir, Magnetic Monopoles (Springer, Berlin-Heidelberg-New York, 2005).
  * [9] E. Shuryak, Nonperturbative Topological Phenomena in QCD and Related Theories, Lect. Notes Phys. 977, 1 (2021).
  * [10] V. Dzhunushaliev, V. Folomeev, and A. Serikbolova, Monopole solutions in SU(2) Yang-Mills-plus-massive-nonlinear-spinor-field theory, Phys. Lett. B 806, 135480 (2020).
  * [11] V. Dzhunushaliev, N. Burtebayev, V. N. Folomeev, J. Kunz, A. Serikbolova, and A. Tlemisov, Mass gap for a monopole interacting with a nonlinear spinor field, Phys. Rev. D 104, 056010 (2021).
  * [12] W. Heisenberg, Introduction to the unified field theory of elementary particles (Max-Planck-Institut für Physik und Astrophysik, Interscience Publishers London, New York, Sydney, 1966).
  * [13] C. M. Bender, K. A. Milton, and V. Savage, Solution of Schwinger-Dyson equations for PT symmetric quantum field theory, Phys. Rev. D 62, 085001 (2000).
  * [14] M. Frasca, Differential Dyson–Schwinger equations for quantum chromodynamics, Eur. Phys. J. C 80, 707 (2020).
  * [15] S. R. Coleman and E. J. Weinberg, Radiative Corrections as the Origin of Spontaneous Symmetry Breaking, Phys. Rev. D 7, 1888 (1973).
  * [16] V. Dzhunushaliev and V. Folomeev, Proca tubes with the flux of the longitudinal chromoelectric field and the energy flux/momentum density, Eur. Phys. J. C 80, 1043 (2020).
  * [17] J. Beltran Jimenez and L. Heisenberg, Generalized multi-Proca fields, Phys. Lett. B 770, 16 (2017).
  * [18] V. Dzhunushaliev and V. Folomeev, Proca balls with angular momentum or flux of electric field, Phys. Rev. D 105, 016022 (2022).
  * [19] D. C. Wilcox, Turbulence Modeling for CFD (DCW Industries, Inc. La Canada, California, 1994).


11institutetext:  Institute for Astronomy, University of Hawaii at Manoa, 640
N. Aohoku Place, Hilo, HI 96720, USA
11email<EMAIL_ADDRESS>22institutetext: Tata Institute of Fundamental
Research, Homi Bhabha Road, Mumbai 400005, India 33institutetext: Max Planck
Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

# CO outflows from young stars in the NGC 2023 cluster

G. Sandell 11 B. Mookerjea 22 R. Güsten 33

(Last edit )

###### Abstract

Context. Young early-type HAeBe stars are still embedded in the molecular
clouds in which they formed. They illuminate reflection nebulae, which shape
the surrounding molecular cloud and may trigger star formation. They are
therefore ideal places to search for ongoing star formation activity.

Aims. NGC 2023 is illuminated by the Herbig Be star HD 37903. It is the most
massive member of a small young cluster with about 30 PMS stars, several of
which are Class I objects that still heavily accrete. It might therefore be
expected that they might drive molecular outflows. We examined the whole
region for outflows.

Methods. We analyzed previously published APEX data to search for and
characterize the outflows in the NGC 2023 region. This is the first systematic
search for molecular outflows in this region. Since the outflows were mapped
in several CO transitions, we can determine their properties quite well.

Results. We have discovered four molecular outflows in the vicinity of NGC
2023, three of which are associated with Class I objects. MIR-63, a bright
mid-infrared and submillimeter Class I source, is a binary with a separation
of 2${}^{\prime\prime}\mskip-7.6mu.\,$4 and drives two bipolar outflows
orthogonal to each other. The large southeast–northwest outflow excites the
Herbig-Haro flow HH 247. MIR-73, a Class I object, which is also a far-
infrared source, drives a pole-on outflow. MIR-62 is a Class II object with
strong infrared excess and a luminosity of 7 L⊙. It is not detected in the
far-infrared. The Class I sources have bolometric luminosities of about 20 L⊙
or lower, that is, they are all low-mass stars. One other far-infrared source,
MIR-75, may have powered an outflow in the past because it now illuminates an
egg-shaped cavity.

Conclusions. The four outflows are all powered by young stars and are located
in the immediate vicinity of NGC 2023. They are at a similar evolutionary
stage, which suggests that their formation may have been triggered by the
expanding C ii region.

###### Key Words.:

ISM: jets and outflows – ISM: Herbig-Haro objects – ISM: individual objects:
NGC 2023 – stars: pre-main sequence – stars: formation

##  1 Introduction

The reflection nebula NGC 2023 in the L 1630 dark cloud is a well-studied
reflection nebula that has been the source of several discoveries. The near-
infrared observations by Sellgren (1984) showed that there must be very small
dust grains that are transiently heated by a single photon. These are now
generally known as polyaromatic hydrocarbons, or PAHs. Fullerenes, the
extremely stable C60 and C70 molecules, were also found in NGC 2023 (Sellgren
et al., 2011). More importantly, however, it has served as a testbed for
developing and testing models of photodissociation regions (PDRs) (Black & van
Dishoeck, 1987; Draine & Bertoldi, 1996, 2001; Kaufman, Wolfire & Hollenbach,
2006).

NGC 2023 is illuminated by the young B2 Ve star HD 37903. This Herbig Be star
is the most massive member of a young stellar cluster with $\sim$ 20 – 30 pre-
main-sequence (PMS) stars (Sellgren, 1983; DePoy et al., 1990; Mookerjea et
al., 2009; López-Garcia et al., 2013). The reflection nebula is embedded in
the L 1630 dark cloud at a distance of $\sim$ 400 pc. Most of the PMS stars in
the cluster are at the southern and eastern boundary of NGC 2023, where the
reflection nebula expands into the dense surrounding molecular cloud, creating
a hot PDR (Sandell et al., 2015). Several of the PMS stars are very young
Class I sources (Mookerjea et al., 2009) that are heavily accreting and likely
to drive molecular outflows. However, the outflows in the NGC 2023 region have
not been studied much so far. Malin, Ogura & Walsh (1987) found a chain of HH
objects that is now known as HH 247, in the southwestern boundary of NGC 2023.
They were not able to identify the exciting source of the HH objects, although
they speculated that MIR-62 would be a likely candidate because of its
proximity to the HH objects. Mookerjea et al. (2009) suggested that the HH
objects could be excited by the submillimeter/infrared source MIR-63, but
stated that this was not a clear identification. We have used the extensive
APEX CO data sets from Sandell et al. (2015) to search for and characterize
the molecular outflows from young stars in the NGC 2023 cluster. We find four
molecular outflows and note that even HD 37903 may drive a weak molecular
outflow.

##  2 Observations

###  2.1  APEX and SOFIA observations

All the Atacama Pathfinder EXperiment, APEX111APEX, the Atacama Pathfinder
Experiment is a collaboration between the Max-Planck-Institut für
Radioastronomie, Onsala Space Observatory (OSO), and the European Southern
Observatory (ESO).(Güsten et al., 2006), and the SOFIA (Stratospheric
Observatory for Infrared Astronomy) CO(11–10) and [C ii] observations were
discussed in detail in Sandell et al. (2015). With APEX, we observed CO(3 –2),
CO(4–3), CO(6–5), CO(7–6), and 13CO(3–2), which are crucial for this paper. We
additionally examined the SOFIA observations of [O i] 63 $\mu$m, which were
presented in Mookerjea et al. (2023). [C ii] is generally not an outflow
tracer and is not discussed further. Although the [O i] fine structure line at
63 $\mu$m can be strongly enhanced in shocks and is often seen in jet-like
molecular outflows, it was not detected in any of the outflows either. In
Table 1 we list the basic information of the spectral lines that we analyzed
for molecular outflow activity. Further details of the observations can be
found in Sandell et al. (2015) and Mookerjea et al. (2023). All the spectral
lines discussed in this paper were calibrated in main-beam brightness
temperature. In order to compare the CO maps more easily, we smoothed the
CO(6–5) and CO(7–6) data by regridding them to the same spatial resolution as
the CO(4–3) map, that is, 14${}^{\prime\prime}\mskip-7.6mu.\,$4\. These
smoothed maps are used in all the figures, tables, and analyses presented in
this paper.

Table 1: Observing setup. The quoted rms sensitivities are in Tmb for a velocity resolution of 0.5 km s-1. $\theta_{\rm FWHM}$ is the half-power beam width, and $\eta_{\rm mb}$ is the main-beam coupling efficiency.  Receiver | Molecular transition | Frequency | $\theta_{\rm FWHM}$ | $\eta_{\rm mb}$ | rms
---|---|---|---|---|---
|  | [GHz] |  |  | [K]
FLASH+ |  13CO(3–2) | 330.587965 | 18.5 | 0.68 | 0.3
FLASH+ | CO(3–2) | 345.795990 | 17.7 | 0.68 | 0.2
FLASH+ | CO(4–3) | 461.040768 | 13.3 | 0.58 | 0.4
CHAMP+ | CO(6–5) | 691.473076 | 9.1 | 0.49 | 1.3
CHAMP+ | CO(7–6) | 806.651806 | 7.7 | 0.48 | 2.5
GREAT L1 | CO(11–10) | 1267.014486 | 23.0 | 0.67 | 1.1
upGREAT HFA |  [O i] 63 $\mu$m | 4744.777490 | 6.3 | 0.63 | 1.6
Table 2: Coordinates and UKIDSS magnitudes of the outflow sources in NGC 2023.
For saturated sources (marked with a), we quote 2MASS magnitudes.

Source | $\alpha_{2000}$ | $\delta_{2000}$ | $Z$ | $Y$ | $J$ | $H$ | $K$
---|---|---|---|---|---|---|---
| [h m s] |  [∘ ′ ″] | [mag] | [mag] | [mag] | [mag] | [mag]
MIR-62 | 5:41:36.390 | -2:16:46.1 |  13.663 $\pm$ 0.002 |  12.818 $\pm$ 0.001 |  11.567 $\pm$ 0.001 |  10.049 $\pm$ 0.023a |  18.83 $\pm$ 0.02a
MIR-63a | 5:41:37.178 | -2:17:17.3 | … | … | … |  14.822 $\pm$ 0.010 |  11.24 $\pm$ 0.01a
MIR-63b | 5:41:37.028 | -2:17:18.0 | … | … |  18.70 $\pm$ 0.18 |  14.127 $\pm$ 0.006 |  11.39 $\pm$ 0.01a
MIR-73 | 5:41:44.601 | -2:16:06.6 |  19.38 $\pm$ 0.13 |  17.92 $\pm$ 0.04 |  15.832 $\pm$ 0.012 |  12.212 $\pm$ 0.001 |  10.38 $\pm$ 0.03a
MIR-75 | 5:41:44.75 | -2:15:54.8 |  17.39 $\pm$ 0.03 |  15.80 $\pm$ 0.01 |  14.63 $\pm$ 0.05 |  111.44 $\pm$ 0.05a |  99.99 $\pm$ 0.03a

a magnitude from 2MASS

###  2.2 UKIDDs images

We retrieved the UKIRT Infrared Deep Sky Surveys (UKIDDS)222The UKIDSS project
is defined in Lawrence et al. (2007). UKIDDS uses the UKIRT Wide Field Camera
(WFCAM; Casali et al., 2007) and a photometric system described in Hewett et
al. (2006). calibrated images as well as pipeline-reduced photometry of NGC
2023 from the WFCAM science archive333http://was.roe.ac.uk. NGC 2023 was
observed as part of the Galactic Clusters Survey on 2005 October 14, in all
UKIDDS filters. The near-infrared counterpart to MIR-63, a binary with a
separation of 2${}^{\prime\prime}\mskip-7.6mu.\,$4, did not have pipeline-
reduced photometry. We manually performed photometry using GAIA444GAIA is a
derivative of the Skycat catalogue and image display tool, developed as part
of the VLT project at ESO. Skycat and GAIA are free software under the terms
of the GNU copyright. and APT555http://www.aperturephotometry.org on both
components using 2″ apertures. The photometry for the outflow sources
discussed in this paper is given in Table 2.

###  2.3 Herschel PACS archive data

NGC 2023 was observed twice with the Herschel Space Observatory666Herschel was
an ESA space observatory with science instruments provided by European-led
Principal Investigator consortia and with important participation from NASA.
in the guaranteed time key program, the Herschel Gould Belt Survey (HGBS)
(André et al., 2010; Bontemps et al., 2010). On operational day (OD) 668 (2011
March 13) NGC 2023 was observed with fast scanning (60″/sec) in parallel mode
with PACS (70/160 $\mu$m; AOR 1342215984) and SPIRE. We retrieved the PACS 70
$\mu$m level 2.5 JScanam image from the Herschel data archive. The field was
additionally observed with medium scan speed (20″/sec) and an orthogonal cross
scan on OD 513 (2010 October 8) with PACS (100/160 $\mu$m; AORs 1342206080 and
1342206081). These data sets have a much better image quality. Level 2.5
JScanam images were retrieved for both 100 $\mu$m and 160 $\mu$m from the
Herschel data archive. All these data sets were discussed in more detail by
Stutz et al. (2013). Although PACS photometry of NGC 2023 sources was
published by several groups (Furlan et al., 2016; Könuves et al., 2020;
Fischer et al., 2020), the deduced flux densities vary strongly, and we found
one source, MIR-75, which is visible both at 70 $\mu$m and 100 $\mu$m, but has
no entry in any of the above mentioned papers, nor is it found in the PACS
point source catalog at IRSA777https://irsa.ipac.caltech.edu. We therefore
performed photometry using a combination of point-spread function fitting on
background-subtracted images using GAIA and MIRIAD (Sault et al., 1995). The
PACS photometry is listed in Table 3.

Table 3: Positions and flux densities of Herschel PACS sources in the vicinity of NGC 2023. Positions are derived from the 100 $\mu$m image.  Name |  $\alpha$(2000.0) |  $\delta$(2000.0) |  S(70 $\mu$m) |  S(100 $\mu$m) |  S(160 $\mu$m)
---|---|---|---|---|---
| [h m s] |  [∘ ′ ″] | [Jy] | [Jy] | [Jy]
MIR-63 | 05:41:37.165 | -2:17:16.3 |  36.5 $\pm$ 3.0 |  35.8 $\pm$ 1.0 |  33.2 $\pm$ 3.8
MIR-73 | 05:41:44.580 | -2:16:05.2 |  13.4 $\pm$ 2.0 |  16.0 $\pm$ 0.7 |  18.2 $\pm$ 1.4
MIR-75 | 05:41:44.760 | -2:15:55.3 |  14.3 $\pm$ 0.4 |  15.4 $\pm$ 1.1 | …

##  3 Molecular outflows

In the radio regime, CO is by far the best tracer of molecular outflows. CO is
one of the most abundant molecules after H2, which has no dipole moment and
which therefore has no observable transitions at radio wavelengths. CO is
easily excited in both diffuse and dense gas. The most accurate mass estimates
are obtained by observing at least two rotational transitions and one line of
13CO (Bachiller & Gómez-González, 1992). We have observed five rotational
transitions of 12CO and one 13CO transition. Molecular outflows can also be
studied in a variety of other molecular lines, especially in strong shocks
(see e.g. Bachiller & Pérez Gutiér, 1997), but even in these regions, the CO
emission is quite strong. For the purpose of this study, there was no need to
observe any other molecular tracer, especially since the nondetection of [O i]
63 $\mu$m, a good shock tracer, demonstrates that there are no strong shocks
in any of the outflows.

Figure 1: Contours of high-velocity CO emission plotted in 2 km s-1 wide
intervals from a center velocity of 10.5 km s-1 and overlaid on the 8 $\mu$m
IRAC image in color. HD 37903 is at an offset 0″, 0″. The left panels show
high-velocity emission in CO(3–2), and the right panels show this in CO(4–3).
The blue- and redshifted velocity ranges for the CO emission are indicated in
the top right corner in each plot. The positions of the young stars MIR-63,
MIR-62, and MIR-73 are labeled in the top left panel, and the Herbig Haro
objects are plotted with black squares. The hatched lines in the CO(3–2) panel
for intermediate velocities are drawn to outline the bipolar outflow from
MIR-63. The topmost panels show three linearly spaced contours from 0.4 K to
1.1 K in CO(3–2), and from 0.65 K to 1.68 K and 1.125 K in red- and
blueshifted CO(4–3), respectively. At high velocities, $\pm$ 7.5 km s-1, we
plot four contours the from 0.6 K for CO(3–2) and 0.7 K for CO(4–3), with the
highest red and blue contours at 2.01 K and 2.41 K for CO(3–2) and 2.12 K and
2.46 K for CO(4–3). For low and intermediate velocities, we plot six contours
starting at 0.8 K for both CO(3–2) and CO(4–3) at intermediate velocities,
$\pm$ 5.5 km s-1, and at low velocities, $\pm$ 3.5 km s-1. The contours start
at 3 K and 2.5 K to avoid contamination from the molecular cloud. For
intermediate velocities, the highest contours are at 6.1 K and 6.0 K and 4.9 K
and 6.1 K for CO(3–2) and CO(4–3), and for low velocities, they are 11.6 K and
16.3 K and 12.0 K and 15.9 K for CO(3–2) and CO(4–3), respectively.

###  3.1 Overview

In order to identify potential outflows more easily, we created channel maps
from the CO(3–2) and CO(4–3) images and overlaid the blue- and redshifted
channel maps on top of an 8 $\mu$m IRAC image in color (see Figure 1). This
enabled us to see how the blue- and red-shifted high velocity gas are related
to each other and allowed us to identify potential bipolar outflows. In Figure
1, we show blue- and redshifted gas in the velocity range 2 – 8 km s-1and 13 –
21 km s-1. We omitted the emission from 8 – 13 km s-1, which is dominated by
cloud emission. These velocity ranges capture almost all of the high-velocity
gas. The emission was averaged over 2 km s-1 wide velocity intervals spaced
symmetrically around a center velocity of 10.5 km s-1. We only plot the part
of the images with high-velocity gas. After we identified potential outflow
sources, we examined the spectra in the vicinity of that source and explored
position-velocity diagrams, which is a good way to obtain information about
the velocity structure of the outflows.

The dominant outflow source is MIR-63. It is also a submillimeter (MM 3) and a
far-infrared source. In the near-infrared, it is a heavily reddened binary. At
low redshifted velocities, 13 – 15 km s-1, a wide-angle redshifted outflow
lobe is visible, centered on MIR-63 at a PA888Position Angle (PA) is measured
counterclockwise from north of $\sim$ 130°, with a more highly collimated
blueshifted (6 – 8 km s-1) lobe at a PA of $\sim$ -10°, which appears to form
a misaligned bipolar outflow. Some redshifted emission also surrounds the
Class I objects MIR-73 and MIR-75 at $\sim$ 120″ east of HD 37903, where there
are no young stars. At these velocities, no gas is observed to be associated
with the Herbig Haro flow HH 247999The Herbig Haro flow HH 247 was originally
called HH 1, HH 4 and HH 5 by Malin, Ogura & Walsh (1987), where blueshifted
CO emission would be expected.

This picture changes at higher velocities. Faint blueshifted emission centered
on the HH 1 cluster is visible there, which was most likely hidden in the
relatively strong emission at near cloud velocities from the molecular cloud
surrounding NGC 2023. Blue- and redshifted CO emission is also visible toward
MIR-73 and MIR-75 at intermediate velocities, $\pm$ 5.5 km s-1, but not at
higher velocities. A closer examination shows that MIR-73 powers this low-
velocity outflow. The blueshifted outflow lobe from MIR-63, which extended all
the way to HD 37903 at low velocities, is far more compact at higher
velocities and appears to turn by 90° at the position of MIR-62 (see, e.g.,
CO(4–3) in the velocity range 4 – 6 km s-1). At higher velocities, it seems to
align with a jet-like feature that is apparently associated with MIR-62,
suggesting that it is a separate outflow driven by MIR-62.

At higher velocities, the redshifted outflow from MIR-63 appears to be more
collimated, with two prominent peaks east of the star. The northern peak,
which dominated at low velocities (13 – 15 km s-1), is still strong at 15 – 17
km s-1, but we also see a southern peak that completely dominates the emission
at higher velocities. A line from the northern peak through MIR-63, that is,
at a PA of 84°, that is extended westward passes through HH suggesting that
the northern redshifted peak forms a bipolar outflow with HH 1, which is part
of the blueshifted counterflow. Likewise, a line from the southern peak
through MIR-63, that is, PA = 118 °, that is extended westward passes through
HH 4. The northern and southern redshifted peaks appear to be part of a single
outflow, whose blueshifted counterflow is the Herbig Haro flow HH 247. This is
illustrated in Fig. 1 in the CO(3–2) panel for intermediate velocities, that
is, the velocity range 4 – 6 km s-1 in blue- and 15 – 17 km s-1 in redshifted
CO. We outlined the outflow by the hatched black lines drawn through the two
eastern redshifted CO peaks and continuing through MIR-63 to the blueshifted
side of the outflow, coinciding with HH 4 and HH 1.

The blueshifted outflow lobe north of MIR-63 must be a separate outflow
powered by MIR-63. It would therefore be expected to have a redshifted
counterflow to the south, which is indeed the case. At low velocities, a
tongue of emission protrudes from MIR-63 at a PA of $\sim$170°, which agrees
well with the northern blueshifted outflow, which is at a PA of $\sim$ -10°.
MIR-63 therefore drives two bipolar outflows, a southeast-northwest outflow,
and a north-south one. This is expected because we know that it is a binary.

Before we start to derive some basic parameters of these outflows, we review
what is known about the outflow sources and what we have learned from the
near- and far-infrared photometry presented in this paper.

###  3.2 Outflow sources

Figure 2: UKIDDS color image (JHK) of the MIR-63 and the MIR-62 region in the
southern part of NGC 2023. Both stars are labeled. The very red binary system
in the middle of the image is MIR-63. The western component of the binary
coincides with the visible star. MIR-62 is $\sim$ 65″ north of MIR-63. We have
labeled the bar Sheffer et al. (2011) called SR, the southern ridge, which
extends east of the star at a PA $\sim$ 110°. There is a dark lane (cavity) on
the opposite side. We outlined the extent of the blue outflow lobe centered on
the southern ridge with a hatched white ellipse. The HH object that Malin,
Ogura & Walsh labeled HH 1 is also labeled. In this image, it appears as a
faint extended red feature. The size of the image is 183.6″ $\times$ 101.6″.
Figure 3: Panel of CO spectra (CO(3–2), 13CO(3–2), CO(4–3), and CO(6–5))
toward MIR-63, the southern and northern peak in the red outflow, and the
blueshifted emission peak near HH 1. MIR-63 is at an offset
-18${}^{\prime\prime}\mskip-7.6mu.\,$3,
-104${}^{\prime\prime}\mskip-7.6mu.\,$7, the southern peak at
41${}^{\prime\prime}\mskip-7.6mu.\,$3,
-133${}^{\prime\prime}\mskip-7.6mu.\,$0, the northern peak at +30″,-97″, and
HH 1 at -97″, -110″. All offsets are relative to HD 37903. The spectra toward
the northern peak, where the CO emission is very strong, are scaled by a
factor of 0.6. The gray vertical line marks the systemic velocity ($\sim$ 10.0
km s-1). The spectra for the different CO transitions are offset by 20 K.

####  3.2.1 MIR-63

MIR-63 is a bright mid-infrared source that Mookerjea et al. (2009) found to
coincide with the submillimeter source MM 3. MM 3 was first detected by
Wyrowski et al. (2000) at 3 mm with BIMA and also at 850 $\mu$m and 450 $\mu$m
with SCUBA. At the time, there was no known optical or near-infrared
counterpart, and therefore Wyrowski et al. assumed that it was most likely a
cold starless core. Mookerjea et al. (2009) clearly identified it with MIR-63,
however, and they also associated Sellgren’s star D with the same mid-infrared
source. They classified MIR-63 as a Class I object and suggested that it might
be exciting the Herbig Haro flow HH 247, originally named HH 5, HH 4, and HH 1
by Malin, Ogura & Walsh (1987). We now know that MIR-63 is a heavily reddened
binary in the near-infrared (Fig. 2). The eastern component, MIR-63a, is
associated with a nebulous arc and is only seen in the H and K bands. This
star dominates the emission in the mid- and far-infrared. The western
component, MIR-63b, which similarly is heavily reddened, is also detected at J
(Table 2). MIR-63b coincides with the optically visible star. It clearly
contributes to the emission at mid-infrared wavelengths and appears to have
some emission even at 70 $\mu$m. In the far-infrared, MIR-63a is quite strong
and was detected by PACS in all bands (Table 3).

Although the source appeared to be extended in SCUBA observations, SCUBA did
not have the angular resolution to resolve the source. Recently, MIR-63 was
observed with ALMA and the VLA as part of the VANDAM survey of Orion
protostars (Tobin et al., 2020). It was not detected by the VLA, but the
continuum emission at 870 $\mu$m was resolved by ALMA into two sources,
separated by 2$\aas@@fstack{\prime\prime}$4 (Tobin et al., 2020). The brighter
eastern component coincides with our far-infrared source MIR-63a, while the
fainter component coincides with MIR-63b. They marginally resolve MIR-63a as
0$\aas@@fstack{\prime\prime}$19 $\times$ 0$\aas@@fstack{\prime\prime}$12; PA =
3.1°, while b is unresolved. MIR-63 was also included in the ALMA Survey of
Orion Planck Galactic Cold Clumps (ALMASOP), which observed 72 young dense
cores at 1.3 mm in three different array configurations in continuum and
molecular lines (Dutta et al., 2020). They found MIR-63 to be extended and
identified four embedded sources in the cloud core. Two of them are compact
and coincide with MIR-63a and MIR-63b. These are the same sources that Tobin
et al. (2020) detected at 870 $\mu$m. The other two were extended, $\sim$1″ to
$>$3″, and appear to be starless cores. They are within 500 au of the two
young heavily accreting Class I protostars, both of which drive outflows.
Their temperatures are therefore likely to be around 25 to 30 K. A recent
paper Luo et al. (2023) reanalyzed the Dutta et al. ALMA data for MIR-63 and
claimed that it drives a complicated molecular outflow that is viewed pole-on.
Their analysis of the starless cores showed that they are transient structures
for any reasonable assumption of gas temperature. As we already saw in the
previous section, MIR-63 drives two molecular outflows, a large southeast-
northwest (SE-NW) outflow, and a more compact north-south (N-S) outflow.
Neither of them are close to pole-on. The SE-NW component has an inclination
of $\sim$ 70°(see below). The N-S outflow probably has a similar inclination
because the red- and blueshifted outflow lobes are spatially well separated.

We derive a bolometric luminosity of 33.0 $\pm$ 1.2 L⊙ using the photometry in
Tables 2 & 3 supplemented with Spitzer IRAC and MIPS photometry from Mookerjea
et al. (2009). Additionally, we used the tabulated Spitzer IRS photometry
points from Furlan et al. (2016). The latter showed that MM 3, or MIR-63a, has
a fairly strong silicate absorption feature, indicating that the disk is seen
somewhat edge-on. The bolometric temperature is around 200 K for an assumed
foreground extinction of 10 mag. Furlan et al. (2016) derived a bolometric
luminosity of 28.2 L⊙. This is slightly lower than what we derive, most likely
because of differences in the IRAC/MIPS flux densities. The modeling by Furlan
et al. (2016) suggests that the outflow has an inclination of 70°, which
sounds very plausible.

In the overview section, we interpreted the outflow emission from MIR-63 seen
in Fig. 1 as being composed of two outflows. The large redshifted outflow lobe
east of the source is one bipolar outflow, and the blueshifted lobe is faint
but still visible in the Herbig-Haro flow HH 247. This outflow is most likely
driven by MIR-63a, which was resolved with ALMA by Tobin et al. (2020) to have
a roughly north-south oriented disk (PA $\sim$3°). The second outflow goes
roughly north-south, with the blueshifted outflow lobe to the north and the
redshifted lobe to south, partly overlapping with the large eastern redshifted
outflow lobe. The position-velocity diagram that goes through the northeastern
redshifted peak and the HH objects to the west, that is, PA 84°, (Fig. 8),
shows that to the east, the outflow lobe terminates at $\sim$ 80″ from MIR-63,
while the blueshifted outflow lobe extends to $\sim$100″. The position-
velocity diagram through the southeastern redshifted peak (Fig. 9), that is,
at a PA of 118°, shows that the redshifted emission is still quite strong at
the edge of our CO maps, at 100″ from MIR-63, indicating that we do not fully
cover the southeastern outflow lobe. In the following, we assume that the
redshifted outflow has a total extent of 110″. There is some blueshifted
emission to the west, but except for the high-velocity gas associated with HH
1 the blueshifted outflow is largely invisible.

At the position of MIR-63, CO(3–2), CO(4–3), and CO(6–5) show strong blue- and
redshifted high-velocity wings extending from about -5 to +28 km s-1 (Fig. 3).
The systemic velocity, determined from 13CO(3–2), is $\sim$ 10 km s-1. Only
faint high-velocity emission at near cloud velocities is seen in CO(7–6), and
CO(11–10) shows no sign of high-velocity emission. The 13CO(3–2) emission
shows clear redshifted high velocity wings not only on MM 3, but also on the
northern and southern redshifted peak (Fig. 3), indicating that the redshifted
high-velocity gas is somewhat optically thick. At the blueshifted high-
velocity peak, which roughly coincides with the HH 1 cluster, the blueshifted
emission is stronger in CO(6–5) and CO(4–3) than in CO(3–2), indicating that
the blueshifted high-velocity gas in this region is warmer than in the rest of
the outflow. Even CO(7–6) shows a blueshifted emission wing extending to a
velocity of -8 km s-1. This is the only region where high-velocity emission is
seen in CO(7–6).

The N-S outflow, PA $\sim$ -10° from MIR-63, is less extended. The blueshifted
outflow lobe, which extends to $\sim$ 60″ from MIR-63, shows high-velocity gas
to $\sim$ -4 km s-1 (total velocity 14 km s-1), while the redshifted gas
extends to +28 km s-1(total velocity 18 km s-1; see Fig. 10). However, since
the redshifted outflow lobe overlaps with the redshifted southeastern outflow
lobe, some of the emission could also be from the southeastern outflow. There
is still high-velocity gas at the edge of our map, indicating that the outflow
extends farther than 45″ from MIR-63. In our analysis, we assumed that it
terminates 55″ from MIR-63. As we already mentioned, the N-S outflow is most
likely driven by MIR-63b, the optically visible star in the MIR-63 core.

Figure 4: CO spectra toward MIR-62 at an offset
-30${}^{\prime\prime}\mskip-7.6mu.\,$0,
-74${}^{\prime\prime}\mskip-7.6mu.\,$6, and at the tip of the blue outflow
lobe at an offset +4${}^{\prime\prime}\mskip-7.6mu.\,$0,
-81${}^{\prime\prime}\mskip-7.6mu.\,$1\. The spectra at the tip of the outflow
are scaled by a factor of 0.6, as indicated in the figure. The CO transitions
are offset by 20 K for clarity. The 13CO(3–2) spectra are plotted in red. The
vertical gray line marks the systemic velocity ($\sim$ 10 km s-1) on MIR-62.

####  3.2.2 MIR-62

MIR-62 (Sellgren’s star C; Haro 5-49) is also a binary with a separation of
$\sim$ 0${}^{\prime\prime}\mskip-7.6mu.\,$35 (Sheffer et al., 2011; Kounkel et
al., 2016). It has strong infrared excess and was classified as a Class I/II
object by Mookerjea et al. (2009). Megeath et al. (2012) classified it as a
Class II source and estimated a foreground extinction of $\sim$ 5 mag. It was
not detected by PACS at 70 $\mu$m, and it is not a submillimeter source.
Scarrott, Rolph & Mannion (1989) found that it illuminates a small bipolar
reflection nebula oriented north-south (see the JHK image in Fig. 2.) This
would suggest that it might drive a north-south outflow, but as we described
above, it drives an east-west outflow. However, if the outflow is largely in
the plane of the sky, then light would be reflected from the cavity walls of
the east-west outflow. This is most likely the case. Kounkel et al. (2019)
found it to have a spectral type of M3. We determine a bolometric luminosity
of 6.7 $\pm$ 0.1 L⊙ and a bolometric temperature of $\sim$1900 K for an
assumed foreground extinction of 5 mag.

Fig. 1 showed that the star appears to power a blueshifted outflow lobe to the
east. This outflow lobe is centered on a bright ridge pointing toward MIR-62
(see Fig. 2). Sheffer et al. (2011) argued that this is a PDR feature and
named it the southern ridge (SR) because studies of vibrationally excited H2
emission (McCartney et al., 1999; Burton et al,, 1990) and purely rotational
H2 lines (Sheffer et al., 2011) showed that these lines are dominated by PDR
emission in the southern ridge and are not due to shock excitation. It is a
strange coincidence that the southern ridge is so well aligned with the
blueshifted outflow lobe if it is not associated with the outflow. The eastern
end of the ridge is 35″ from MIR-62, while we estimate the length of the
outflow to be $\sim$ 45″ (0.087 pc) from the position-velocity diagram in Fig.
11.

The outflow velocities are rather modest. At the tip of the outflow, they
extend from $\sim$ -2 to +16 km s-1 (Fig. 4). There is some redshifted
emission at the position of MIR-62, but there is no sign of a redshifted
counterflow (Fig. 11). The UKIDDS JHK color image in Fig. 2 shows a narrow
dark structure west of the star, which is approximately aligned with the
blueshifted outflow lobe, and which opens up to a faintly limb-brightened
cavity farther northwest (like a wine glass). This is most likely created by
the counterflow from MIR-62. It is seen even more clearly in the JHK image of
Sheffer et al. (2011), their Fig. 1.

Table 4: Physical parameters of the outflows. Dynamical timescales and mass-
loss rates are computed using mass weighted velocities. All derived parameters
are corrected for inclination. The formulae for the inclination correction
factors are given in Bachiller et al. (1990).

Parameter | MIR-63a | MIR-63b | MIR-62 | MIR-73
---|---|---|---|---
| red | blue | red | blue | blue | red | blue
Maximum outflow velocity (km s-1) | 18 | 15 | 18 | 14 | 12 | 5.6 | 5.4
Inclination (°) | 70a | 70a | 70 | 70 | 70 | 17.6 | 2.4
Temperature (K) | 30 | 75 | 30 | 30 | 40 | 30 | 40
Dynamical time scale (103 yr) | 5.6 | n.a. | 3.1 | 4.3 | 3.2 | 14.4 | 12.3
Mass (10-3 M⊙) | 130 | 2.4b | 52 | 48 | 45 | 2.4 | 1.5
Momentum (10-2 M⊙ km s-1) | 160 | n.a. | 63 | 47 | 44 | 0.99 | 0.67
Force (10-5 M⊙ km s-1 yr${}^{-1})$ | 10.0 | n.a. | 7.4 | 3.8 | 4.8 | 0.055 | 0.073
Energy (1036 J) | 50.4 | n.a. | 8.6 | 5.0 | 1.7 | 0.014 | 0.031
Mass loss rate (10-6 M⊙ yr-1) | 23 | n.a. | 17 | 11 | 14 | 0.17 | 0.12

a Outflow inclination from Furlan et al. (2016); b Mass of HH 1 condensation
only

####  3.2.3 MIR-73

MIR–73, on the eastern side of the reflection nebula (Fig. 5), was classified
as a Class I object by Mookerjea et al. (2009). It has strong near-infrared
excess and is relatively bright in the far-infrared. It was easily detected by
PACS at all bands (Table 3). It was not detected at 850 $\mu$m by SCUBA
(Johnstone, Matthews & Mitchell, 2006), but Furlan et al. (2016) reported a
marginal detection at 350 $\mu$m. Furlan et al. (2016) derived a bolometric
luminosity of 17.8 L⊙, and a bolometric temperature of 277 K for an assumed
foreground extinction of 16 mag. This agrees reasonably well with what we
derive, a bolometric luminosity of 20.9 $\pm$ 0.1 L⊙ and a bolometric
temperature of about 260 K for a foreground extinction of 10 mag.

Figure 5: UKIDDS color image (JHK) of the MIR-73 and MIR-75 region east of HD
37903. Both are Class I sources with strong infrared excess, and both were
detected by PACS. The stars are labeled. MIR-73, the very red source in the
middle of the image, drives a compact almost pole-on outflow. MIR-75 does not
appear to drive a molecular outflow, but sits at the apex of a bubble of
emission, which we outlined with a dashed red ellipse. The hole west of it is
an image artifact. The size of the image is
41${}^{\prime\prime}\mskip-7.6mu.\,$6 $\times$
37${}^{\prime\prime}\mskip-7.6mu.\,$6\.  Figure 6: CO spectra toward MIR-73.
The left panel shows CO(3–2), and the right panel shows CO(4–3). The 13CO(3–2)
spectrum is also plotted in red in the left panel. The vertical gray line
marks the systemic velocity, $\sim$ 10.4 km s-1. The narrow emission feature
at 5.2 km s-1 is due to an IR dark cloud filament and not related to the
outflow from MIR-73.  Figure 7: High-velocity emission of CO(6–5) overlaid on
a 100 $\mu$m PACS image in color. The blueshifted CO is integrated from +6 to
+ 8 km s-1, and the redshifted CO is integrated from 13–15 km s-1. This
CO(6–5) map was created with the same angular resolution as CO(4–3) or
14${}^{\prime\prime}\mskip-7.6mu.\,$4\. MIR-73 is placed at an offset 0″, 0″ ,
and MIR-75 is $\sim$ 12″ north of it. Both are marked by black stars.

The systemic velocity of the molecular gas surrounding MIR-73 is $\sim$ 10.4
km s-1(see Fig. 6). The blue- and redshifted emission lobes partly overlap.
The separation between the red- and blueshifted peak is
10$\aas@@fstack{\prime\prime}$1 , suggesting that the outflow is seen almost
pole-on (see Fig. 7). The outflow models by Cabrit & Bertout (1986) make an
inclination angle of 20° appear about right. The blueshifted peak is almost
centered on the star, offset by only 1$\aas@@fstack{\prime\prime}$2 to the
west, while the redshifted outflow lobe is off by
8$\aas@@fstack{\prime\prime}$9 to the east. We therefore assign an inclination
of 2.4° to the blueshifted outflow and leave 17.6° for the redshifted outflow.
The outflow velocities are very modest, $\sim$ 5 – 6 km s-1 for both outflow
lobes (see Fig. 6).

####  3.2.4 MIR-75

MIR-75 was classified as a Class I object by Mookerjea et al. (2009). We see
no high-velocity CO toward the star, but it sits at the apex of a limb-
brightened cavity on the western side of the star (see Fig. 5). This egg-
shaped cavity is seen in all UKIDDS bands, suggesting that it is illuminated
by reflected light from MIR-75. It appears likely that a past outflow may have
created the cavity, although there is currently no evidence for outflow
activity. MIR-75 is also detected with PACS at 70 and 100 $\mu$m and is most
likely present at 160 $\mu$m as well, but the emission from the star cannot be
separated from the strong emission from the surrounding cloud. We derive a
bolometric luminosity of 3.6 $\pm$ 0.2 L⊙ and a bolometric temperature of 825
$\pm$ 9 K, assuming that the foreground extinction is 10 mag, which is similar
to the foreground extinction toward MIR-73.

There is no sign of a cavity on the eastern side, possibly because it is
obscured by the IR dark cloud filament just east of the star. Another
possibility is that MIR-75 is near the edge of the NGC 2023 cloud, so that the
outflow expanded more freely into the low-density cloud envelope.

####  3.2.5 HD 37903

Wyrowski et al. (2000) found a weak extended ($\sim$ 40″) continuum source at
8 GHz centered at $\sim$ -10″, -12″ relative to HD 37903. They interpreted
this continuum emission as free-free emission from a density enhancement
ionized by the Lyman continuum flux from HD 37903. However, PDR tracers such
as [C ii] and [O i] do not show any enhancement from this continuum source
(Sandell et al., 2015; Mookerjea et al., 2023). As we showed in Fig. 1, a
blueshifted emission feature lies southwest of the star extending to $\sim$ 4
km s-1, that is, $\sim$ 6.5 km s-1 from the systemic velocity. This
blueshifted emission peaks roughly at -16″, -24″ from HD 37903 and has a size
of $\sim$ 50″. It is in the same direction as the low brightness continuum
emission, suggesting that they may be related. Since we see it in CO(3–2),
CO(4–3) and CO(6–5), we can estimate the temperature assuming LTE and
optically thin emission (see Section 3.3). This gives us an excitation
temperature of $\sim$ 30 K. If the blueshifted CO emission were associated
with the ionized emission feature, the gas would be expected to be warmer. The
cause for this blueshifted CO emission is unclear. Maybe it is just being
pushed by the stellar wind from HD 37903. The total mass of the mass of the
blueshifted gas$\sim$ 0.008 M⊙, which is on the same order as mass in the
MIR-73 outflow (see Table 4).

###  3.3 Properties of the outflows

Since we detect all outflows in CO(3–2), CO(4–3), and CO(6–5) and in most
cases do not see any high-velocity gas in CO(7–6), we know that the outflowing
gas must be warmer than 20 K, but probably not much warmer than 40 K. There is
one exception, however. The blueshifted high-velocity gas associated with the
HH 1 cluster is detected even in CO(7–6), and here, we derive a temperature of
$\sim$ 75 K. The redshifted counterflow is much colder. To estimate the column
densities, we integrated over the outflows in 2 km s-1 wide velocity intervals
from -1 to 22 km s-1, except for the velocity range 8 – 13 km s-1, which is
dominated by strong molecular cloud emission. For regions in which the
outflows overlap, we assumed equal contribution from each outflow. We repeated
this for 13CO(3–2), although we only detected 13CO(3–2) line wings in the near
cloud velocity ranges, that is, from 6 – 8 km s-1 in blueshifted emission and
13 – 15 km s-1 in redshifted emission. We can now solve for the optical depth
from the ratio of CO/ 13CO by assuming that both have the same excitation
temperature and that we know the isotope ratio [C/13C], which we assume is 70
(see, e.g., Sandell et al. (2015)). In this way, we find that all outflows
except MIR-73 are somewhat optically thick at low velocities, with optical
depths in the range 2.5 – 6 and a median of around 3.5. For higher velocities,
our upper limits suggest that the emission is largely optically thin. We can
now determine the excitation temperatures and CO column densities by assuming
that CO is in LTE. We can now use the formulae for calculating column
densities given by Mangum & Shirley (2015) for optically thin emission and
assuming the excitation temperature is the same for all velocities over each
outflow lobe. To do this, we calculated the column density for a grid of
temperatures from 25 – 50 K for each CO transition, and then we determined the
resulting temperature within the errors with the same column density for each
CO transition. For the initial analysis, we gave the near cloud velocities
lower weight because we know that they are optically thick. After we
constrained an excitation temperature (Table 4), we calculated optical depth
corrections for CO(4–3) and CO(6–5) to obtain optical depth-corrected column
densities. For each velocity bin, we took the average of the total column
densities derived for each transition. The masses given in Table 4 assume that
the CO/H2 ratio is 10-4 and that the total atomic weight per hydrogen atom is
1.4 (Kauffmann et al., 2008).

We calculated the outflow parameters largely following the procedures outlined
by Lada (1985), Bachiller et al. (1990), and Choi, Evans & Jaffe (1993). To
derive the mass, momentum, and and energy, we integrated over the same
velocity intervals as we use in Fig. 1. All outflow parameters were corrected
for inclination using the formulas given in Bachiller et al. (1990) (see also
Cabrit & Bertout (1992)). To determine the dynamical timescales, we used the
mass-weighted outflow velocity rather than the maximum velocity. The
difference compared to using the maximum outflow velocity for our outflow
sources is smaller than 50%. The outflow parameters are given in Table 4.

In the review by Bachiller & Gómez-González (1992), they found that a
comparison of different methods for determining the accuracy of the outflow
parameters by Cabrit and her collaborators (Cabrit & Bertout, 1986; Cabrit,
Goldsmith & Snell, 1988; Cabrit & Bertout, 1992) indicated that the most
accurate methods estimate the mass within a factor of 2, the momentum within a
factor of 10, and the force within a factor of 20. Even with the extensive CO
data set we have, we do not appear to do much better. Our mass estimates are
uncertain because of the overlap of outflows and because we average the
emission over the whole outflow lobe. This results in large uncertainties in
the optical depth correction for gas at low velocities. Furthermore, we are
somewhat overcautious in setting the boundary between outflow emission and the
ambient cloud, which may contain the bulk of the outflowing gas, resulting in
a uncertainty of about a factor of 3. For other outflow parameters, the
uncertainty in inclination is the dominant error source. We assumed an
inclination angle of 70° for the two MIR-63 outflows as well as for MIR-62.
However, the inclination angle for these sources could be as high as 80° or as
low as 50°, resulting in uncertainty factors of 4 for momentum, 14 for energy,
and 18 for the outflow force, just from the uncertainty in inclination. For
MIR-73, the almost pole-on outflow, the errors in momentum and energy due to
inclination are small, while the dynamical timescale and force is uncertain by
a factor of 10 or more.

##  4 Discussion

All the molecular outflows discovered in the NGC 2023 region are in the
immediate periphery of the C ii region powered by HD 37903. Furthermore, they
all lie in the eastern and southern part of nebula, where the C ii region
expands into the dense surrounding molecular cloud. This suggests that their
formation may have been triggered by the C ii region. Elmegreen & Lada (1977)
introduced the concept of triggered star formation when they investigated the
formation of OB associations, which appear to follow a sequential star-forming
mechanism. In this scenario, an H ii region ionizes the surrounding molecular
gas and drives a shock front into the cloud, accumulating a dense layer of gas
between the ionized gas and the shock front. This compressed gas layer
eventually becomes unstable and forms a new OB subgroup. Depending on the
density structure of the surrounding cloud, the compressed gas layer may
trigger star formation in several different ways, ranging from radiation-
driven implosion of pre-existing clouds to a collect-and-compress scenario of
the gas layer (see, e.g., Deharveng et al. (2010)). The latter mechanism would
also work for a C ii region, but instead of triggering OB stars, it would
trigger low-mass stars. A good example of collect-and-collapse triggered star
formation is the bright H ii region RWC 79, which shows a dust emission ring
surrounding the ionized gas, in which Zavagno et al. (2006) found two massive
fragments that were diametrically opposite to each other, each containing
bright near-infrared sources. They identified 12 luminous Class I sources
toward the most massive of these dust fragments. It is often difficult,
however, to judge whether the expanding H ii region has triggered new star
formation or just revealed pre-existing dense cores in the surrounding cloud.
NGC 2023, however, looks very similar to RWC 79 in the sense that all the
young Class I sources are near the periphery of the C ii region. Therefore,
their formation may well have been triggered by the expanding C ii region.

It is clear that when we study the properties of outflows from young stars,
both high spatial resolution and sensitivity are required to see large-scale
structures and small details. With our moderate spatial resolution, $\sim$
15″, we can characterize all the outflows in the vicinity of NGC 2023, but we
cannot say with certainty whether it is MIR-63a or b that drives the large
southeast-northwest outflow. The CO(2–1) images by Luo et al. (2023) are
heavily filtered by the ALMA interferometer and do not give a clear view of
the two outflows. We assigned the outflow to MIR-63a because Tobin et al.
(2020) found it to be elongated north-south in continuum at 0.87 mm. If the
continuum emission traces the accretion disk, it would be roughly orthogonal
the outflow. It is also the more luminous star in the binary system.
Nevertheless, it would be very useful to combine the ALMA data with single-
dish data in order to better determine the morphology and origin of each
outflow in the region in which they are launched.

##  5 Summary and conclusions

The reflection nebula NGC 2023 is illuminated by the Herbig Be star HD 37903.
It is the most massive member of a small young cluster with about 30 PMS
stars, a few of which are Class I objects and still accrete heavily. We
therefore expected that they might power molecular outflows. By analyzing
previously published data, we found four molecular outflows. These four
outflows are all powered by young stars and are all located in the immediate
vicinity of the reflection nebula. All the outflows have roughly similar ages,
suggesting that the stars driving them were formed about the same time, which
supports the hypothesis that they may have been triggered by the expanding C
ii region. There are no Class 0 sources in the vicinity of NGC 2023, although
young protostars like this exist in the large L 1630 molecular cloud, but they
are obviously unrelated to the young NGC 2023 cluster.

###### Acknowledgements.

We would like to thank an anonymous referee for the valuable comments and
suggestions which helped to improve the manuscript. BM acknowledges the
support of the Department of Atomic Energy, Government of India, under Project
Identification No. RTI 4002. We also thank Dr. W. Varricattu for advice on how
to analyze the UKIDDS images. The paper is based in part on observations made
with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA).
SOFIA is jointly operated by the Universities Space Research Association, Inc.
(USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI)
under DLR contract 50 OK 0901 to the University of Stuttgart. This work is
based in part on data obtained as part of the UKIRT Infrared Deep Sky Survey.

## References

  * André et al. (2010) André, Ph., Men‘shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102
  * Bachiller & Pérez Gutiér (1997) Bachiller, R., & Pérez Gutiér, M. 1997, ApJ, 487, L93
  * Bachiller & Gómez-González (1992) Bachiller, R., Gómez-González, J. 1992. A&A Rev., 3, 257
  * Bachiller et al. (1990) Bachiller, R., Cernicharo, J., Martín-Pintado, J., et al. 1990, A&A, 231, 174
  * Black & van Dishoeck (1987) Black, J. H., & van Dishoeck, E. F. 1987, ApJ, 322, 412
  * Bontemps et al. (2010) Bontemps, S., André, Ph., Könyves, V., et al. 2010, A&A, 518, L85
  * Burton et al, (1990) Burton, M. G., Geballe, T. R., Brand, P. W. J. L., et al. 1990, ApJ, 352, 625
  * Cabrit & Bertout (1992) Cabrit, S., & Bertout, C. 1992, A&A, 261, 274
  * Cabrit, Goldsmith & Snell (1988) Cabrit, S., Goldsmith, P. F., & Snell, R. L. 1988, ApJ, 334, 196
  * Cabrit & Bertout (1986) Cabrit, S., & Bertout, C. 1986, ApJ, 307, 313
  * Casali et al. (2007) Casali, M., Adamson, A., Alves de Oliveira, C., et al. 2007, A&A, 467, 777
  * Choi, Evans & Jaffe (1993) Choi, M., Evans II, N. J., & Jaffe, D. T. 1993, ApJ, 417, 624
  * Deharveng et al. (2010) Deharveng, L., Schuller, F., Anderson, L. D., et al. 2010, A&A, 523, A6
  * DePoy et al. (1990) DePoy, D. L., Lada, E. A., Gatley, I., et al. 1990 ApJ, 356, L55
  * Draine & Bertoldi (2001) Draine, B. T., & Bertoldi, F. 2000, in Molecular hydrogen in Space, Cambridge, UK: Cambridge University Press, Eds, F. Combes & G. Pineau des Forr̂ets, p. 131
  * Draine & Bertoldi (1996) Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269
  * Dutta et al. (2020) Dutta, F., Lee, C.-F., Liu, T., et al. 2020, ApJS, 251,20
  * Elmegreen & Lada (1977) Elmegreen, B. G., & Lada, C. J. 1977, ApJ, 214, 725
  * Fischer et al. (2020) Fischer, W. J., Megeath, S. T., Furlan, E., et al. 2020, ApJ, 905, 119
  * Fleming et al. (2010) Fleming, B., France, K., Lupu, R. E., & McCandliss, S. R. 2010, ApJ, 725, 159
  * Furlan et al. (2016) Furlan, E., Fischer, W. J., Ali, B., et al. 2016, ApJS, 224, 5
  * Güsten et al. (2006) Güsten, R., Nyman, L.-Å, Schilke, P., et al. 2006, A&A, 454, L13
  * Hewett et al. (2006) Hewett, P. C., Warren, S. J., Leggett, S. K., et al. 2006, MNRAS, 367, 454
  * Johnstone, Matthews & Mitchell (2006) Johnstone, D., Matthews, H., & Mitchell, G. F. 2006, ApJ, 639, 259
  * Kauffmann et al. (2008) Kauffman, J., Bertoldi, F. Bourke, T. l., et al. 2008, A&A, 487, 993
  * Kaufman, Wolfire & Hollenbach (2006) Kaufman, M. J.., Wolfire, M. G., & Hollenbach, D. J. 2006, ApJ, 644, 283
  * Kounkel et al. (2019) Kounkel, M., Covey, K., Moe, M., et al. 2019, AJ, 157, 196
  * Kounkel et al. (2016) Kounkel, M., Megeath, S. T., Poteet, C. A., et al. 2016, ApJ, 821,52
  * Könuves et al. (2020) Könuves, V., André, Ph., Arzoumanian, D., et al. 2020, A&A, 635, A34
  * Lada (1985) Lada, C.J. 1985, ARA&A, 23, 267
  * Lawrence et al. (2007) Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599
  * López-Garcia et al. (2013) López-Garcia, M. A., López-Santiago, J. Albacete-Colombo, J. F., et al. MNRAS, 429, 775
  * Luo et al. (2023) Luo, Q.-y., Liu, T., Lee, A. T., et al. 2023, ApJ, 952, L2
  * Malin, Ogura & Walsh (1987) Malin, D. F., Ogura, K., & Walsh, J. R. 1987, MNRAS, 227, 361
  * Mangum & Shirley (2015) Mangum, J. G., & Shirley, Y. 2015, PASP, 127, 266
  * McCartney et al. (1999) McCartney, M.S. K., Brand, P. W. J. L., Burton, M. G., et al. 1999, MNRAS, 307, 315
  * Megeath et al. (2012) Megeath, S. T., Gutermuth, R. Muzerolle, J., et al. 2012, AJ, 144, 192
  * Mookerjea et al. (2023) Mookerjea, B., Sandell, G., Güsten, R., et al. 2023, MNRAS, 525, 5468
  * Mookerjea et al. (2009) Mookerjea, B., Sandell, G., Jarrett, T. H., & McMullin, J. P. 2009, A&A, 507, 1485
  * Sandell et al. (2015) Sandell, G., Mookerjea, M., Güsten, R., et al. 2015, A&A, 578, A41
  * Sault et al. (1995) Sault, R. J., Teuben, P. J., & Wright, M. C. H. 1995, in Astronomical Society of the Pacific Conference Series, Vol. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes, 433–+
  * Scarrott, Rolph & Mannion (1989) Scarrott, S. M., Rolph, C. D., & Mannion, M. D. 1989,MNRAS, 237, 1027
  * Sellgren et al. (2011) Sellgren, K., Werner, M. W., Ingalls, J. G., et al. 2011, in EAS Publications Series, Volume 46, p. 209
  * Sellgren (1984) Sellgren, K. 1984, ApJ, 277, 623
  * Sellgren (1983) Sellgren, K. 1983, AJ, 88, 1985
  * Sheffer et al. (2011) Sheffer, Y., Wolfire, M. G., Hollenbach, D. J., et al. 2011, ApJ, 741, 45
  * Stutz et al. (2013) Stutz, A. M., Tobin, J. J., Megeath, S. T., et al. 2013, ApJ, 767, 36
  * Tobin et al. (2020) Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130
  * Wyrowski et al. (2000) Wyrowski, F., Walmsley, C. M., Goss, W. M., et al. 2000, ApJ, 543, 245
  * Zavagno et al. (2006) Zavagno, A., Deharveng, l., Comerón, F., et al. 2006, A&A, 446, 171

##  Appendix A Position-velocity diagrams

We plot position-velocity diagrams of CO(3–2) and CO(4–3) for the large
southeastern outflow from MIR-63, that is, cutting through both the northern
and southern red peaks as well as through the outflow from MIR-62.

Figure 8: Position-velocity diagram plotted with contours and overlaid with
grayscale through the northern redshifted peak of the large outflow powered by
MIR-63a, i.e., at a PA $\sim$84°.  Figure 9: Position-velocity diagram plotted
with contours and overlaid with grayscale through the southern redshifted peak
of the large southeastern outflow powered by MIR-63a, i.e., at a PA
$\sim$118°,  Figure 10: Position-velocity diagram plotted with contours and
overlaid with grayscale of the north-south outflow powered by MIR-63b, i.e.,
at a PA of -10°.  Figure 11: Position-velocity plot with contours and overlaid
with grayscale along the outflow axis of the blueshifted outflow from MIR-62,
i.e., at a PA $\sim$110°.