text:
# 1 Introduction

## 1 Introduction

MACRO was a large area multipurpose underground experiment designed to search for rare events in the penetrating cosmic radiation. These included the study of atmospheric neutrinos and their oscillations, high energy $`(E_\nu \stackrel{>}{}\mathrm{\hspace{0.17em}\hspace{0.17em}1}\mathrm{GeV})`$ muon neutrino astronomy, indirect searches for WIMPs and search for low energy ($`E_\nu \stackrel{>}{}\mathrm{\hspace{0.17em}\hspace{0.17em}7}`$MeV) stellar collapse neutrinos. The detector was placed in the Gran Sasso laboratory, located on the highway Rome-Teramo , 120 km east of Rome. The lab consists of three underground halls, each about 100 m long. It is at an altitude of 963 m above sea level, is well shielded from cosmic rays by a mean rock thickness of $`3700`$ m.w.e.; the minimum is $`3150`$ m.w.e: this defines the minimum muon energy at the surface at $`1.3\mathrm{TeV}`$ in order to reach MACRO. The average residual muon energy and the muon flux are $`320\mathrm{GeV}`$ and $`1\mathrm{m}^2\text{h}^1`$, respectively.

MACRO was composed of three sub-detectors: liquid scintillation counters, limited streamer tubes and nuclear track detectors. Each one of them could be used in “stand-alone” and in “combined” mode. A cross section of the detector is shown in Fig. 1. Notice the division in the lower and in the upper part (this was often referred to as the Attico); the inner part of the Attico was empty and lodged the electronics . The mass of the lower MACRO was $`4200`$t, mainly in the form of boxes filled with crushed Gran Sasso rock. The detector had a modular structure: it was divided into six sections referred to as supermodules. Each active part of one supermodule had a size of $`12.6\times 12\times 9.3\mathrm{m}^3`$ and had a separate mechanical structure and electronics readout. The full detector had global dimensions of $`76.5\times 12\times 9.3\mathrm{m}^3`$ and provided a total acceptance to an isotropic flux of particles of $`10000\mathrm{m}^2\mathrm{sr}.`$ The total mass was $`5300`$t.

Data were taken from March 1989 till april 1994 with the detector under construction and from the middle of 1994 till the end of 2000 with the full detector.

## 2 Atmospheric neutrino oscillations

If neutrinos have non-zero masses, one has to consider the weak flavour eigenstates $`\nu _\mathrm{e},\nu _\mu ,\nu _\tau `$ and the mass eigenstates $`\nu _1,\nu _2,\nu _3`$. The weak flavour eigenstates $`\nu _ł`$ are linear combinations of the mass eigenstates $`\nu _\mathrm{m}`$ through the elements of the mixing matrix $`U_{lm}`$:

$$\nu _ł=\underset{m=1}{\overset{3}{}}U_{lm}\nu _m$$
(1)
In the simple case of only two flavour eigenstate neutrinos ($`\nu _\mu ,\nu _\tau `$) which oscillate with two mass eigenstates ($`\nu _2,\nu _3`$) one has

$$\{\begin{array}{cc}\nu _\mu =\hfill & \nu _2\mathrm{cos}\theta _{23}+\nu _3\mathrm{sin}\theta _{23}\\ \nu _\tau =\hfill & \nu _2\mathrm{sin}\theta _{23}+\nu _3\mathrm{cos}\theta _{23}\end{array}$$
(2)
where $`\theta _{23}`$ is the mixing angle ($`\theta `$ will be used in the following). In this case one may easily compute the following expression for the survival probability of a $`\nu _\mu `$ beam

$$P(\nu _\mu \nu _\mu )=1\mathrm{sin}^22\theta _{23}\mathrm{sin}^2\left(\frac{E_2E_1}{2}t\right)=1\mathrm{sin}^22\theta _{23}\mathrm{sin}^2\left(\frac{1.27\mathrm{\Delta }m^2L}{E_\nu }\right)$$
(3)
where $`\mathrm{\Delta }m^2=\mathrm{m}_3^2\mathrm{m}_2^2`$ and L is the distance travelled by the muon neutrino from production to detection.

High energy primary cosmic rays interact in the upper atmosphere producing pions and Kaons, which by decay yield muons and muon neutrinos, $`\pi \mu \nu _\mu `$; further decays of the muons lead to electron and muon neutrinos, $`\mu \nu _\mu \nu _łł`$. The neutrinos are produced in a spherical shell at about 10- 20 km from the earth surface. Upgoing $`\nu _\mu `$’s may interact in the rock below MACRO or inside its lower part leading to upgoing muons, $`\nu _\mu N\mu ^++\mathrm{}`$.

Upward going muons are identified by the streamer tube system (for tracking) and the scintillator system (for time-of-flight measurement). A rejection factor of at least $`10^7`$ is needed in order to separate upgoing muons from the background due to the downgoing muons. Fig. 1 shows sketches of the different neutrino event topologies analyzed: Upthroughgoing muons, Upsemicontained (also called Internal Upgoing muons, IU), Upgoing Stopping muons (UGS), Internal Downgoing muons (ID). The numbers of events measured and expected for the three topologies are given in Table 1. All data samples deviate from MC expectations; the deviations point to neutrino oscillations.

### 2.1 Upthroughgoing muons

The upthroughgoing muons come from $`\nu _\mu `$ interactions in the rock below the detector; the $`\nu _\mu `$’s have an average energy $`\overline{E}_\nu \mathrm{\hspace{0.17em}\hspace{0.17em}50}\mathrm{GeV}`$ . The data refer to a livetime of 6.16 years (full detector equivalent). The data deviate in absolute value and in shape from Monte Carlo (MC) predictions, see Fig. 2. We studied a large number of possible effects that could affect our measurements: no significant systematic problems exist in the detector or in the data analyses.

The measured data have been compared with MC simulations. For the upthroughgoing muon simulation, the neutrino flux computed by the Bartol group was used . The cross sections for the neutrino interactions were calculated using the deep inelastic parton distributions of ref. . The muon propagation to the detector was done using the energy loss calculation in standard rock. The total systematic uncertainty on the expected muon flux, obtained adding in quadrature the errors from neutrino flux, cross section and muon propagation, is estimated to be 17 %. This uncertainty is mainly a scale error; the error on the shape of the angular distribution is $`5\%`$.

Fig. 2 shows the zenith angle ($`\mathrm{\Theta }`$) distribution of the measured flux of upthroughgoing muons. The Monte Carlo expectation for no oscillations is shown as a dashed line. Assuming $`\nu _\mu \nu _\tau `$ oscillations, the best fit parameters are $`\mathrm{\Delta }m^2=2.510^3\mathrm{eV}^2`$ and $`sin^22\theta =1`$ with a probability of 66%; the result of the fit is the solid line in Fig. 2. The probability for no-oscillations is 0.4 %.

Fig. 3a shows the allowed region for $`\nu _\mu \nu _\tau `$ oscillations in the $`sin^22\theta \mathrm{\Delta }m^2`$ plane, computed according to ref.  for the upthroughgoing muon events;

our region is compared with those obtained by the SuperKamiokande (SK)  and Soudan 2  experiments.

### 2.2 Matter effects. $`\nu _\mu \nu _\tau `$ versus $`\nu _\mu \nu _{sterile}`$

Matter effects due to the difference between the weak interaction effective potential for muon neutrinos with respect to sterile neutrinos would produce a different total number and a different zenith distribution of upthroughgoing muons . In Fig. 4 the measured ratio between the events with $`1<cos\mathrm{\Theta }<0.7`$ and the events with $`0.4<cos\mathrm{\Theta }<0`$ is shown as a black point. In this ratio most of the theoretical uncertainties on neutrino flux and cross section cancel. The remaining theoretical error combined with the experimental error is estimated to be $`7\%`$. We measured 305 events with $`1<cos\mathrm{\Theta }<0.7`$ and 206 with $`0.4<cos\mathrm{\Theta }<0`$; the ratio is R = $`1.48\pm 0.13_{stat}\pm 0.10_{sys}`$. For $`\mathrm{\Delta }m^2=2.510^3\mathrm{eV}^2`$ and maximal mixing, the expected value of the ratio for $`\nu _\mu \nu _\tau `$ is $`R_\tau =1.72`$ while for $`\nu _\mu \nu _s`$ is $`R_{sterile}=2.16`$. The maximum probabilities $`P_{best}`$ to find a value of $`R_\tau `$ and of $`R_{sterile}`$ smaller than the expected ones are 9.4 % and 0.06 % respectively. The ratio of the maximum probabilities is $`P_{best_\tau }/P_{best_{sterile}}=157`$, so that $`\nu _\mu \nu _s`$ oscillations are disfavoured at 99% c.l. compared to the $`\nu _\mu \nu _\tau `$ channel .

### 2.3 $`\nu _\mu `$ energy estimates by MCS of upthroughgoing muons

The oscillation probability is a function of the ratio $`L/E_\nu `$, Eq. 3. Even if a precise measurement of $`E_\nu `$ is not possible, MC simulations have shown a correlation between $`E_\nu `$ and the residual muon energy $`E_\mu `$. $`E_\nu `$ was thus estimated for each event by measuring the muon energy $`E_\mu `$, via Multiple Coulomb Scattering (MCS) of the iduced muon in the absorbers. The r.m.s. of the lateral displacement for a muon crossing the apparatus on the vertical is $`\sigma _{MCS}10\mathrm{cm}/\mathrm{E}_\mu (\mathrm{GeV})`$. The muon energy $`E_\mu `$ can be measured up to a saturation point, occurring when $`\sigma _{MCS}`$ is comparable with the detector space resolution. A first analysis was made studying the deflection of upthroughgoing muons with the streamer tubes in “digital mode”. This method could reach a spatial resolution of $`1\mathrm{cm}`$ which implies a maximum measurable energy of 10 GeV . As the interesting energy region for atmospheric neutrino oscillations spans from $`1\mathrm{GeV}`$ to tens of GeV, it is important to improve the spatial resolution of the detector to push the saturation point as high as possible. For this purpose, a second analysis was performed with the streamer tubes used in “drift mode”. To check the electronics and the feasibility of the analysis, two tests were performed at the CERN PS-T9 and SPS-X7 beams. The space resolution achieved is $`3\mathrm{mm}`$, a factor 3.5 better than in the first analysis. For each muon, seven MCS sensitive variables were given in input to a Neural Network (NN) previously trained to estimate the muon energy with MC events of known energy crossing the detector at different zenith angles. The method allowed to separate the upthroughgoing muons in 4 subsamples corresponding to average neutrino energies of 12, 20, 50 and 102 GeV, respectively. The comparison of the 4 zenith angle distributions with the predictions of the no oscillations MC shows a disagreement at low energies (where there is a deficit of vertical events), while the agreement is restored at the higher neutrino energies. The distribution of the ratios $`R=(Data/MC_{noosc})`$ obtained by this analysis is plotted in Fig. 5 as a function of $`log_{10}(L/E_\nu )`$  . The black points with error bars are the data; the vertical extent of the shaded areas represents the uncertainties on the MC predictions for $`\nu _\mu \nu _\tau `$ oscillations with maximal mixing and $`\mathrm{\Delta }m^2=2.510^3\mathrm{eV}^2`$. The horizontal dashed line is the expectation without oscillations. The black square point was obtained from the low energy IU sample.

### 2.4 Low energy data.

The Internal Upgoing (IU) muons come from $`\nu _\mu `$ interactions in the lower apparatus; for these events two scintillation counters are intercepted; the T.o.F. is applied to identify upward going muons. The average parent neutrino energy for these events is 4.2 GeV. For neutrino oscillations one expects a reduction by about a factor of two in the flux, without any distortion in the shape of the angular distribution. This is what is observed in Fig. 6a .

The upstopping muons (UGS) are due to external $`\nu _\mu `$ interactions yielding upgoing muons stopping in the detector. The semicontained downgoing muons (ID) are due to $`\nu _\mu `$-induced downgoing muon tracks with vertex in the lower MACRO (Fig. 1). The two types of events are identified by topological criteria; the lack of time information prevents to distinguish the two sub-samples. An almost equal number of UGS and ID events is expected. In case of oscillations, the flux of the UGS muons should be reduced by 50%; no reduction is instead expected for the semicontained downgoing events (coming from neutrinos with path lengths of $`20\mathrm{km}`$); therefore one expects a global reduction of 25%.

The number of events and the angular distributions are compared with MC predictions in Table 1 and Figs. 6a,b. The data show a uniform deficit for the whole angular distribution with respect to predictions, $`50\%`$ for IU, 75% for ID + UGS; there is good agreement with the predictions based on neutrino oscillations with the parameters obtained from the upthroughgoing muons. The average value of the double ratio $`R=(Data/MC)_{IU}/(Data/MC)_{ID+UGS}`$ over the measured zenith angle range is $`R0.77\pm 0.07`$; the error includes statistical and theoretical uncertainties; $`R=1`$ is expected in case of no oscillations.

## 3 Search for Astrophysical Sources of High Energy Muon Neutrinos

High energy $`\nu _\mu `$’s are expected to come from several galactic and extragalactic sources. Neutrino production requires astrophysical accelerators of charged particles and some kind of astrophysical beam dumps.

The excellent angular resolution of our detector allowed a sensitive search for upgoing muons produced by neutrinos coming from celestial point sources, with a negligible atmospheric neutrino background. An excess of events was searched for around the positions of known sources in $`3^{}`$ (half width) angular bins. This value was chosen so as to take into account the angular smearing due to multiple muon scattering in the rock below the detector and by the energy-integrated angular distribution of the scattered muon, with respect to the neutrino direction.

The pointing accuracy of the detector was checked with the observation of the shadow of primary cosmic rays (and thus on downgoing muons) by the Moon and the Sun  . It was also checked with the study of the correlations between muons detected by MACRO and extensive air showers detected by the EASTOP experiment located above the Gran Sasso massif .

In a total livetime of 6.16 y we obtained a total of 1356 events, see Fig. 7a . The 90% c.l. upper limits on the muon fluxes from specific celestial sources lay in the range $`10^{15}10^{14}\mathrm{cm}^2\mathrm{s}^1`$, see Fig. 7b. The solid line is our sensitivity vs declination. Notice that we have two cases, GX339-4 ($`\alpha =255.71^o`$, $`\delta =48.79^o`$) and Cir X-1 ($`\alpha =230.17^o`$, $`\delta =57.17^o`$), with 7 events: in Fig. 7b they are considered as background, therefore the upper flux limits are higher; but they could also be indications of signals .

We searched for time coincidences of our upgoing muons with sources of $`\gamma `$-ray bursts as given in the BATSE 3B and 4B catalogues, for the period from April 91 to December 2000 . No statistically significant time correlation was found.

We have also searched for a diffuse astrophysical neutrino flux, for which we establish a flux upper limit at the level of $`1.510^{14}\mathrm{cm}^2\mathrm{s}^1`$ .

## 4 Indirect Searches for WIMPs

Weakly Interacting Massive Particles (WIMPs) could be part of the galactic dark matter; they could be intercepted by celestial bodies, slowed down and trapped in their centers, where WIMPs and anti-WIMPs could annihilate and yield $`\nu _\mu `$ upthroughgoing muons. The annihilations in these celestial bodies would yield neutrinos of GeV- TeV energy, in small angular windows from their centers. For the Earth we have chosen a $`15^o`$ cone around the vertical: we find 863 events. The MC expectation for atmospheric $`\nu _\mu `$ without oscillations gives a larger number of events. We set a conservative flux upper limit assuming that the measured number of events equals the expected ones. We obtain the $`90\%`$ c.l. limits for the upgoing muon flux shown in Fig. 8a (it varies from 0.8 to 0.5 $`10^{14}\mathrm{cm}^2\mathrm{s}^1`$) . If the WIMPs are identified with the smallest mass neutralino, our limits may be used to constrain the neutralino mass, following the model of ref. , Fig. 8a.

A similar procedure was used to search for muon neutrinos from the Sun, using 10 search cones from $`3^o`$ to $`30^o`$. In the absence of statistically significant excesses the muon flux upper limits are at the level of $`1.5210^{14}\mathrm{cm}^2\mathrm{s}^1`$. The limits are shown in Fig. 8b as a function of the WIMP (neutralino) mass.

## 5 Neutrinos from Stellar Gravitational Collapses

A stellar gravitational collapse (GC) of the core of a massive star is expected to produce a large burst of all types of neutrinos and antineutrinos with energies of $`730`$ MeV and with a duration of $`<10\mathrm{s}`$. The $`\overline{\nu }_e`$’s can be detected via the process $`\overline{\nu }_e+pn+e^+`$ in the liquid scintillator. About $`100÷150`$ $`\overline{\nu }_e`$ events should be detected in our 580 t scintillator for a stellar collapse at the center of our Galaxy.

We used two electronic systems to search for bursts of $`\overline{\nu }_e`$’s from stellar gravitational collapses. The first system was based on the dedicated PHRASE trigger, the second one was based on the ERP general trigger. Both had energy thresholds of $`7\mathrm{MeV}`$ and recorded pulse shape, charge and timing informations. Immediately after a trigger, the PHRASE system lowered its threshold to 1 MeV, for a duration of $`800\mu \mathrm{s}`$, in order to detect (with a $`25\%`$ efficiency) the $`2.2\mathrm{MeV}`$ $`\gamma `$ released in the reaction $`n+pd+\gamma _{2.2\mathrm{MeV}}`$ induced by the neutron produced in the primary process.

A redundant supernovae alarm system was in operation, alerting immediately the physicists on shift. We defined a general procedure to alert the physics and astrophysics communities in case of an interesting alarm. Finally, a procedure to link the various supernovae observatories around the world was set up .

The effective MACRO active mass was $`580`$t; the live-time fraction was $`97.5\%`$; we kept always ruming at least 1-2 supermodules. No stellar gravitational collapses were observed in our Galaxy from the beginning of $`1989`$ to the end of 2000 .

## 6 Conclusions

The MACRO detector took data from 1989 to the end of year 2000 and it obtained important results in all the items listed in the proposal, in particular:

* Atmospheric neutrino oscillations. Analyses of different event topologies at different energies, the exploitation of the muon Coulomb multiple scattering in the detector give strong support to the hypothesis of $`\nu _\mu \nu _\tau `$ oscillations with $`\mathrm{\Delta }m^2=0.0025eV^2`$ and maximal mixing. The statistical and systematic significance of our global data is well above at the 5 standard deviation level. We are at present performing some global analyses using also the new data on atmospheric muon flux measurements.
* High energy muon neutrino astronomy. MACRO has been highly competitive with other underground experiments thanks to its large acceptance and good angular accuracy. It has been limited by its livetime and the size of the detector.
* Search for bursts of $`\overline{\nu }_e`$ from stellar gravitational collapses. MACRO was sensitive to supernovae events in our Galaxy, it started the SN WATCH system, and for a certain time it was the only detector in operation.
* Sensitive indirect searches have been carried out for possible Dark Matter candidates like the neutralinos, looking for upgoing muons from the center of the Earth and from the Sun.

It may be pointed out that MACRO also obtained the best flux upper limits for GUT Magnetic Monopoles over the widest $`\beta `$ range .

### Acknowledgments

I thank the members of the MACRO Collaboration for their cooperation; in particular I acknowledge the discussions and the help from the members of the Bologna group, in particular Drs. Y. Becherini and M. Giorgini. I thank Ms Anastasia Casoni for typing and correcting the manuscript.


format: text

text:
# Velocity statistics from spectral line data: effects of density-velocity correlations, magnetic field, and shear

## 1 Introduction

The very large Reynolds numbers (defined as the ratio of the inertial force to the viscous force acting on a parcel of gas) clearly suggest that the interstellar medium is turbulent. Understanding of this turbulence is crucial to the correct and complete description of many physical processes that take place in the ISM: molecular cloud dynamics, star formation, heat transfer, magnetic reconnection, accretion disks, cosmic ray propagation and diffusion, just to mention a few.

To adequately describe turbulence we need to use statistical methods, to extract the underlying regularities and reject incidental details. So far, studies of turbulence statistics have been most successful using interstellar scintillations, and have provided important information of the density statistics on scales $`10^810^5cm`$, (see Narayan & Goodman 1989; Spangler & Gwinn 1990). These studies are done with the ionised media, and are restricted to the density fluctuations. We expect to observe power-law power spectrum of density in a turbulent medium, however density fluctuations alone fail to distinguish between fossil and active turbulence. The distribution of sizes of sand grains in a beach follows a power-law but it is not turbulence. Therefore would be better to study a more direct and dynamically important statistic: velocity fluctuations. Studies of the neutral media have included velocity information, e.g. using spectral line widths, centroids of velocity (Miesch & Bally, 1994; Scalo, 1987); or more fancy statistics (Heyer & Schloerb, 1997; Rosolowsky et al., 1999; Brunt & Heyer, 2002a, b).

Turbulence statistics is a very important link between theory and observations. Recent breakthroughs in turbulence studies have been possible with the use of scaling laws (see reviews by Cho & Lazarian 2003; Cho, Lazarian & Vishniac 2003, and references therein). Scaling laws predictions can be confronted with observations.

A discussion of various approaches to the study of turbulence with spectral line data can be found in Lazarian (1999). In particular the problem of the contribution of velocity and density fluctuations to the emissivity statistics was addressed in Lazarian & Pogosyan (2000, hereafter LP00), where analytical description of the power spectrum of the emissivity in velocity channels was obtained. This study provides a technique to extract statistics of density and velocity from spectral line data cubes. H I has been chosen as a test case because it was possible to disregard self-absorption, and consider the emissivity linearly depending on the density (see review in Lazarian, Pogosyan & Esquivel 2002). The technique, termed as velocity-channel analysis (or simply, VCA), is based on the variations of the power spectra in velocity channels with the change of velocity resolution. Predictions in LP00 have been confirmed through observations (Stanimirović & Lazarian, 2001). The first tests using numerical data in Lazarian et al. (2001) also confirmed the predictions in LP00. However they did not address the problems of numerical noise, effects of magnetic field and galactic shear. Moreover the study of velocity-density correlations was performed over a very limited dynamical range.

In this work we subject the predictions in LP00 to further scrutiny. In particular, we systematically analyse the effects of velocity-density correlations, and unlike Lazarian et al. (2001) we include a study of the effect of shear. In addition we address the issue of an anisotropic turbulent cascade, with their implications on the VCA. We also analyse the sources of numerical noise.

In §$`2`$ we review the basic problem, we then (§$`3`$) test the extraction of the spectral indices using the VCA, and obtain the density-velocity correlations. In §$`4`$ we study the effect of shear in the VCA. The anisotropic cascade, and a possible new method to extract information of the magnetic field direction can be found in §$`5`$. And we finally provide a summary of our results (§$`6`$).

## 2 Statistics of spectral line data and VCA

Spectral line observations contain important information of density and velocity. The problem is that both velocity fluctuations and density fluctuations contribute to the intensity of a spectral line at a given velocity, and their separation is far from trivial.

### 2.1 3D turbulence statistics

Turbulence, due to its stochastic nature is studied using statistical tools, such as 3D correlation and structure functions (see Monin & Yaglom 1975). For instance, the density correlation function can be expressed as

$$\xi (𝒓)=\rho (𝒙)\rho (𝒙+𝒓),$$
(1)
where $`𝒙`$ is the spatial position ($`xyz`$ in Cartesian coordinates), $`𝒓`$ is a spatial separation (or ‘lag’), and the angular brackets denote an average to be performed over all $`𝒙`$ space. The power spectrum provides an alternative description and is related to the correlation function (CF) as

$$P(𝒌)=d𝒓\mathrm{e}^{i𝒌\mathbf{}𝒓}\xi (𝒓),$$
(2)
where $`𝒌`$ is the wave-number, related to the scale under study like $`k=2\pi /r`$, and the integration is performed over all 3D space. For power-law statistics, the $`N`$ dimensional power spectrum ($`P_Nk^n`$) and the correlation ($`\xi _Nr^m`$) have indices related as $`n=Nm`$, in other words111Notice that the spectral indices contain the sign. This notation is unfortunate (contrary to a conventional spectral index of the form $`P_Nk^n`$, with $`n`$ positive). However, it was introduced in LP00, and for consistency we adhere to it.

$$(spectralindex)=(spatialdimensionsCFindex).$$
(3)
To describe velocity fields is customary to use structure functions:

$$SF(𝒓)=[𝒗(𝒙)𝒗(𝒙+𝒓)]^2.$$
(4)
They provide a more formal treatment of velocity fields (Monin & Yaglom 1975). In the case of isotropic velocity fields it differs from the correlation function only by a constant. Therefore for power-law statistics it has the same spectral index $`m`$ as the correlation function (i. e. $`SFr^m`$). For instance, Kolmgorov turbulence scales as $`vr^{1/3}`$, which corresponds to $`m=2/3`$.

### 2.2 Statistics in position-position-velocity (PPV) space

In spectroscopic observations we don’t observe the gas distribution in real space coordinates $`𝒙(x,y,z)`$. Instead the intensity of the emission in a given spectral line is measured towards some direction in the sky, given by $`𝑿(x,y)`$, and at a given line-of-sight (LOS) velocity $`v`$. Thus, the coordinates of observational data are $`(𝑿,v)`$. Therefore, observational data are contained in, so called, position-position-velocity (or simply ‘PPV’) cubes. If we identify the $`z`$ coordinate with the LOS, then the relation between the real space and PPV descriptions is that of a map $`(𝑿,z)(𝑿,v)`$.

At a given position $`𝒙`$ the LOS component of the velocity $`v`$ can be decomposed as a regular flow velocity $`u_{reg}(𝒙)`$, plus a thermal component $`v_{thermal}`$ and a turbulent component $`u(𝒙)`$. This way the observed Doppler shifted atoms follow a Maxwellian distribution of the form

$$\varphi _v(𝒙)\mathrm{d}v=\frac{1}{(2\pi \beta )^{1/2}}\mathrm{exp}\left\{\frac{[vv_{reg}(𝒙)u(𝒙)]^2}{2\beta }\right\}\mathrm{d}v,$$
(5)
where $`\beta =\kappa _BT/m`$, $`m`$ is the atomic mass, $`T`$ the gas temperature, and $`\kappa _B`$ the Maxwell-Boltzmann constant.

The density of emitters in PPV space $`\rho _s(𝑿,v)`$ can be obtained by integration along the LOS:

$$\rho _s(𝑿,v)\mathrm{d}𝑿\mathrm{d}v=\left[dz\rho (𝒙)\varphi (𝒙)\right]\mathrm{d}𝑿\mathrm{d}v,$$
(6)
where $`\rho (𝒙)`$ is the mass density of the gas in spatial coordinates. This expression just counts the number of atoms along the LOS, at a given position $`𝑿`$ with a $`z`$ component of velocity in the interval $`[v,v+\mathrm{d}v]`$, and the limits of integration are defined by the extent of spatial distribution of emitting gas.

The statistics available in PPV are the correlation functions, $`\xi _s\rho _s(𝑿_\mathrm{𝟏},v_1)\rho _s(𝑿_\mathrm{𝟐},v_2)`$, and the power spectra (this can be 1D, 2D or 3D) of emissivity fluctuations. Relations between this statistics and the underlying velocity and density spectra were established in LP00, where the 2D power spectrum in velocity channels was used. This spectrum can be directly obtained with radio-interferometric observations (see Lazarian 1995).

### 2.3 Velocity-channel analysis (VCA)

This technique was developed in LP00 to extract the velocity and density spectral indices of turbulence from spectral line data. LP00 provide analytical predictions for the emissivity power spectral indices in velocity channels, as a function of the velocity resolution employed. The relative contribution to the intensity in a velocity channel (also referred as ‘slice’) to the total intensity fluctuations changes in a regular fashion with the width of the velocity slice. This is easy to understand qualitatively since more velocity fluctuations within a channel are averaged out as we increase the thickness of the channel, decreasing the relative contribution of velocity.

The notation that we will use is such that: $`n`$ is the 3D density spectral index ($`P_\rho k^n`$); $`\mu `$ the 3D velocity spectral index ($`P_vk^\mu `$); $`\gamma `$ corresponds to the spectral index in velocity slices ($`P_nk^\gamma `$); and $`m`$ is the structure function index of the velocity. From eq. (3) is clear that

$$\mu =3m.$$
(7)
For a given 3D density power spectrum two distinct regimes are present when: a) $`n>3`$ (‘shallow’ density case), and (b) $`n<3`$ (‘steep’ density case). A summary of the predictions in LP00 for the intensity spectral index in 2D velocity channels is given in Table 1.

In both cases (shallow or steep) the power-law index ‘gradually steepens’ with the increase of the velocity slice thickness. In the thickest velocity channels all the velocity information is averaged out and we naturally get the density spectral index $`n`$. On the other hand, the velocity fluctuations dominate in thin channels. For instance, in the case of a steep density ($`n<3`$) the structure seen in the velocity channels is determined only by velocity fluctuations (is not explicitly dependant on the index $`n`$). This immediately sends a warning against identifying structures in PPV as actual density enhancements (‘clouds’). The spectral velocity index can be obtained if we combine the density velocity index from the very thick channels and the structure function index $`m`$ from thin channels (using Table 1 and eq. 7). Notice that the notion of thin and thick slices depends on the turbulence scale under study, and the same slice can be thick for small-scale turbulent fluctuations and thin for large scale ones. The formal criterion for the slice to be thick is that the dispersion of turbulent velocities on the scale studied should be less than the velocity slice thickness. Otherwise the slice is thin.

From Table 1 one may observe that for any given density spectral index, the spectrum of fluctuations becomes shallower as the velocity field gets steeper. This is somewhat counterintuitive. But can be understood as follows: a steeper velocity index corresponds to relative more power on the largest scales. A steeper velocity power-spectrum means a larger velocity difference across scales, and correspondingly to a larger velocity dispersion. In PPV space, large eddies, because their velocity dispersion, will be spread over a larger velocity range than smaller eddies. In individual channels, the number of emitters associated with large eddies will decrease relative to that of small eddies as the velocity contrast increases. Therefore, in velocity slices, the power at small $`k`$ numbers (large scales, from large eddies) decreases for a steeper velocity, rendering in a shallower power-spectrum.

In following sections we perform tests of the VCA using numerical data, and take a step further exploring the sources of uncertainty of and the effects of shear in the VCA.

## 3 Testing the VCA

Numerical simulations provide good means of theory testing. Using numerical data, we can measure directly the velocity and density statistics, and we can produce synthetic spectra to analyse them as an observer would do. However, we need to be aware of the differences of synthetic and real data. For instance, only a limited number of points is available with numerics. This, as we show later, results in the most important source of uncertainty.

### 3.1 The data

Along the paper we make use of three types of data to construct the PPV cubes and analyse the power spectrum:

1. MHD simulations. We use a 3rd-order hybrid essentially-non-oscillatory (ENO) code to simulate fully developed turbulence. It was performed on a $`216^3`$ Cartesian grid, it is isothermal, compressible, and with a mean magnetic field in the $`x`$ direction. To reduce spurious oscillations near shocks, we combine two ENO schemes. When variables are sufficiently smooth, we use the 3rd-order weighted ENO scheme (Jiang & Hu, 1999) without characteristic mode decomposition. When opposite is true, we use the 3rd-order Convex ENO scheme (Liu & Osher, 1998). We use a three-stage Runge-Kutta method for time integration. We solve the ideal MHD equations in a periodic box:

$`{\displaystyle \frac{\rho }{t}}+\mathbf{}(\rho 𝒗)=0,`$
$`{\displaystyle \frac{𝒗}{t}}+𝒗\mathbf{}𝒗+\rho ^1(c_s^2\rho )(\mathbf{\times }𝑩)\mathbf{\times }𝑩/4\pi \rho =𝒇,`$ (8)
$`{\displaystyle \frac{𝑩}{t}}\mathbf{\times }(𝒗\mathbf{\times }𝑩)=0,`$
with $`\mathbf{}𝑩=0`$ and an isothermal equation of state. Here $`𝒇`$ is a random large-scale driving force, $`\rho `$ the density, $`c_s`$ the sound speed, $`𝒗`$ the velocity, and $`𝑩`$ the magnetic field. The rms velocity $`\delta V`$ is maintained to be approximately unity, so that $`𝒗`$ can be viewed as the velocity measured in units of the rms velocity of the system, and $`𝑩/\sqrt{4\pi \rho }`$ as the Alfvén speed in the same units. The time $`t`$ is roughly in units of the large eddy turnover time ($`L/\delta V`$) and the length in units of $`L`$, the scale of the energy injection. The magnetic field consists of the uniform background field and a fluctuating field: $`𝑩=𝑩_0+𝒃`$. The Alfvén velocity of the mean field is roughly the same as the rms velocity. The average Mach number is $`2.5`$.

The outcomes are the density and velocity fields. The PPV cubes constructed with this data have a correlation between the velocity and the density fields which is self-consistent with MHD evolution. These and are used to test the applicability of the VCA where the basic assumption was to disregard such correlation.
2. ‘Modified’ data. Due to numerical limitations, the power spectrum of the simulations mentioned above doesn’t show exact power-law behaviour over the desired inertial range, making difficult to test the VCA, which deals with power-law spectra. For that reason we apply the same procedure as in Lazarian et al. (2001) to modify the power spectrum to follow a power-law preserving most of the phase correlations. The procedure consists in replacing the amplitudes of the Fourier transform of the data so they follow a power-law while keeping the phases intact (were most of the spatial information is). The set of simulations used here show a larger well-developed inertial range than those in Lazarian et al. (2001); therefore the ‘correction’ introduced is smaller, especially for the $`11/3`$.
3. Gaussian fields. We use a standard procedure to generate Gaussian fields with prescribed power-law spectrum (see Bond, Kofman & Pogosyan 1996). We generate a Fourier representation of the field

$$F(k)=ka_kP(k)^{1/2}\mathrm{e}^{ikr},$$
(9)
where $`a_k`$ are independent Gaussian variables with dispersion $`1`$, and a mean of $`0`$ for velocity fields or $`1`$ for density fields; in our notation the wave-number $`k=2\pi L/\lambda `$, where $`\lambda `$ is the wavelength associated with $`k`$, and $`L`$ is the box size. To obtain the field in configuration space a Fourier transform is applied. Velocity and density fields generated this way are entirely uncorrelated. However the facility of generating as many cubes as we want in a short amount of time is useful when we want to improve the statistics of our sample.

Comparing (ii) and (iii) allows to test the importance of density-velocity correlations.

With the data cubes mentioned above we produced synthetic spectra (i. e. PPV ‘cubes’555Strictly speaking they are parallelepipeds as the sizes need not be the same in all directions.) to analyse them using VCA. The procedure to construct the PPV cubes is as follows: We generate a matrix of the same dimensions in $`𝑿`$ as the $`𝒙`$ data and a given number of velocity bins (velocity channels) as the third dimension. The velocity bins are equally spaced ranging from the minimum value of the velocity ($`v_{min}`$) to the maximum value of the velocity ($`v_{max}`$) in the whole velocity data cube. We then sum all the elements of the density cube in a position $`𝑿`$ that correspond to each velocity bin. This way we simulate observational data, neglecting self absorption, and assuming that the thermal width of the emitting gas is always smaller than the thickness of the velocity channels. An introduction of an intrinsic width (equivalent to thermal) to the emitters is discussed later.

### 3.2 Density and velocity correlations

LP00 considered density and velocity either uncorrelated or correlated to maximal degree. What would happen to intermediate cases? We computed the correlation of the density and the velocity fields as

$$C(𝒓)=\frac{\rho (𝒙)𝒖(𝒙+𝒓)\mathbf{}𝒓/|𝒓|}{\sigma _\rho \sigma _𝒖}.$$
(10)
Where $`𝒓`$ is the spatial separation (‘lag’), $`𝒙`$ is the spatial position, $`\rho `$ the density field, $`𝒖`$ is the velocity field, $`\sigma _\rho `$ and $`\sigma _𝒖`$ are the standard deviations of the density and velocity fields respectively. The average is to be performed over all $`𝒙`$ space and angles of $`𝒓`$. For numerical economy Lazarian et al. (2001) used Fourier transform techniques to obtain similar information (although restricted to the LOS velocity component). Here this is done directly using equation (10), which is applied to the original output from the MHD simulations. We found moderate correlations, the maximum correlation (in absolute value) was of $`0.22`$. This correlation is shown in Fig. 1, where we plotted the average projection of the cross-correlation on the $`x`$, $`y`$, $`z`$ axis and the average correlation as a function of the separation $`r`$.

Clearly the velocity and density correlations are very different along and perpendicular to the magnetic field direction. Some implications of this anisotropy are discussed in §5.

### 3.3 Applying the VCA

In this subsection we will start the test of the VCA generating several data cubes with known (‘modified’) spectral indices, and arbitrary velocity resolutions that will be shown to range from very thick, to excessively thin. This is to check the general trends predicted in LP00. Then with a more careful choice of the velocity thickness of the channel, we will narrow down the measured spectral indices to check the numerical value with the analytical expectation. A more detailed analysis of the sources of error introduced with numerics, with an example case of Kolmogorov indices, follows in the next subsection.

We produced density data cubes with spectral indices modified to $`2.5`$ and $`11/3`$, LOS velocity cubes with indices of $`3.3`$, $`11/3`$ (Kolmogorov), and $`4.0`$. With each velocity and the density fields we constructed PPV cubes with velocity resolutions, ranging from $`150`$ to $`1`$ number of channels. The 2D power spectra in velocity channels (averaged over all velocity channels for each PPV cube) is presented in Fig. 2.

We should stress, that whether a slice in velocity space can be considered thin or thick depends not only on the slice width $`\delta v=(v_{max}v_{min})/N=\mathrm{\Delta }v/N`$ (where $`N`$ is the number of channels) but on the scale as well. For power-law statistics the squared dispersion, $`\sigma _L^2`$, over the box size, $`L`$, equals $`CL^m`$ (here $`m`$ is the velocity structure function index related to the spectral index as described in eq 3); similarly the velocity dispersion squared, $`\sigma _r^2`$, over the scale $`r`$ is $`Cr^m`$. Therefore the criteria for a channel to be considered thin ($`\delta v<\sigma _r`$) from the largest scale $`L`$ up to the scale $`r`$ translates into $`N>\frac{\mathrm{\Delta }v}{\sigma _L}\left(\frac{L}{r}\right)^{m/2}`$. Thus, we do not expect a pure power-law result for a slice of a given width. Indeed, the largest scales, on the order of the size of the cube, are almost always in a thin regime; unless, $`N<\frac{\mathrm{\Delta }v}{\sigma _L}56`$. Figure 2 is consistent with this behaviour.

In Figure 2 we can see that the power spectrum has an additional deviation from the power-law at small scales. For that reason, to measure numerical values of the spectral indices, we have restricted ourselves to the wave-number range from $`\mathrm{log}k1.3`$ to $`2.2`$, where the power spectra in Fig. 2 is close to a power-law. The slopes (spectral indices) of the power spectra obtained with a linear least-squares method are shown in Table 2.

From Fig. 2 (and Table 2) we can see, indeed, a steepening of the power-spectrum as we increase the thickness of the velocity slices. However, the spectral-indices in thin channels do not appear to reach an asymptotic value, instead they keep getting shallower as the velocity resolution increases. This behaviour arise due to the limited number of emitters along the LOS: with only $`216`$ emitters along each LOS, to distribute in a large number of velocity channels, there will be many channels empty. This produces artificial sharp edges, and introduce high frequency components. Rendering in a shallower spectrum at large wave-numbers. To avoid this problem, in Lazarian et al. (2001) it was suggested as optimal to search for thin slice asymptotics in the channels that are just thin enough,

$$N\frac{\mathrm{\Delta }v}{\sigma _L}\left(\frac{L}{r}\right)^{m/2}$$
(11)
at the smallest scale $`r`$ of interest. The expression in eq. (11) is somewhat tricky, because we need to know the value of $`m`$, to get the number of channels to be used to best determine $`m`$! But, it allows to check with emulated observations if the predictions hold. Our results, in Table 2, are in agreement with the prediction in Lazarian et al. (2001) that more channels are required when steeper velocity fields are studied. For real observations the problem is alleviated because in that case we are not limited by the number of emitters and we can use the maximum velocity resolution without having empty channels. In Lazarian et al. (2001) the factor of $`\mathrm{\Delta }v/\sigma _L`$ in eq. (11) was estimated to be $`2.5`$, in the assumption of Gaussian statistics. In practice we can use the standard deviation of the velocity in the LOS velocity cube for $`\sigma _L`$, and $`\mathrm{\Delta }v`$ is the range of velocity covered in our PPV cubes. This is important to keep in mind because a side effect of the modification of the spectral index is an appreciable change in the velocity range of the velocity cubes. To get the factor ($`L/r`$) we need to take the maximum wave-number used to obtain the power-spectra. For the results presented in Table 2 it was $`k10^{2.2}`$, so $`L/r=k/2\pi 25`$. This way we generated PPV cubes with the velocity resolution according to eq. (11), and measure the power spectrum in velocity channels, in Table 3 we display the results.

From Table 3 we can see that the measured spectral indices within $`15\%`$ error of the predictions, but always shallower than expected. This is an indication, that the noise arising from the limited number of emitters is present even in channels barely thin, although not near as large as that present for PPV cubes constructed with more channels (see Table 2).

Could deviations above be due to correlations between the density and the velocity fields? The answer is no: to illustrate this point, we can compare the power spectrum of the MHD simulations to that of an ensemble of synthetic PPV cubes with the same spectral indices but no velocity-density correlations whatsoever. For this comparison to work we must use PPV cubes constructed from data with the same power-law inertial range. For that reason, we modified the spectral index of the density field from the simulations to $`11/3`$ (close to the original index, so the modification on the velocity-density correlations is minimal). With this density field, we then constructed PPV cubes using the modified velocity output from the simulations and $`60`$ synthetic velocity fields. All the velocity cubes used had the same spectral index (also $`11/3`$). The power spectra in velocity channels is shown in Fig. 3.

In this figure we can appreciate that apart of some ‘wiggles’ at small values of $`k`$, due to lack of enough statistics (see discussion below) the power-spectra is practically identical whether the density field is correlated with the velocity field or not.

### 3.4 More on the numerical limitations

As we saw in the previous section, one important practical problem with the determination of the spectral indices through VCA is to measure with accuracy the slopes (precisely the spectral indices) in the plots. We can notice from each plot in Fig. 2, two departures from well behaved power-laws. At the left end of the spectrum (low $`k`$ numbers, large spatial scales) we have some random wiggling, while at the right end of the curves we see a systematic flattening of the power spectrum. The latest is increasingly evident for PPV cubes formed with larger number of channels (with thinner velocity channels), and was referred as ‘shot-noise’ in Lazarian et al. (2001).

#### 3.4.1 Departures from power-law at low wave-numbers

In the spirit of avoiding other sources of uncertainty, we will use simple Gaussian fields with prescribed power-law spectrum as velocity fields, mapped with a constant density field to produce PPV cubes. We generated an ensemble of $`60`$ Gaussian cubes with a prescribed spectral index of $`11/3`$. Then we applied to them the VCA, that is, we obtained PPV cubes with different number of channels (different velocity resolution) and compute the 2D power spectra in velocity channels. The results are shown in Fig. 4.

From the figure is evident that these departures from strict power-laws at low $`k`$ numbers are due to the lack of sufficient statistics at the largest scales. Power spectra in velocity channels for individual velocity cubes differ very little from intermediate to large $`k`$ numbers, but differ significantly at small $`k`$, showing some ‘wiggles’. However, the ensemble average converges to a power-law that extends from the largest spatial scales to intermediate spatial scales, until a systematic flattening of the power spectrum at large $`k`$ becomes evident. This problem may ocurr in actual observations, where the low wave-numbers (large scales) are affected by the general size and shape of the object under study.

#### 3.4.2 Departures from power-law at high wave-numbers

One of the most important limitations that we to construct spectral lines with numerics, is the finite (and very small) number of emitters in a particular direction; in our simulations we have $`216`$ density elements in each LOS to distribute in velocity channels. If we impose a high velocity resolution in the PPV cubes (i. e. large number of channels), the spectral lines instead of being smooth will have sharp edges and empty channels. This artificial jumps introduce high frequency components in the spectrum and therefore an increase of the power spectrum at large wave-numbers. In Lazarian et al. (2001) this was regarded as ‘shot-noise’. This is true in the sense that the flattening becomes larger as we increase the velocity resolution. But is more complicated that the usual shot-noise, that can be easily described in terms of Poisson statistics. This flattening of the spectrum is particularly problematic because it increases as we go into thinner channels, and the predictions in LP00 to extract the velocity spectral index were done precisely in an asymptotic regime for thin channels.

A natural way to avoid the problem of this flatening of the power spectrum is simply to increase the number of emitters along the LOS. This can be accomplished increasing the size of the data cubes, which translates into more computing power required to analyse the power-spectrum (and more importantly more computing power to generate larger MHD cubes). With the resources at our hand, we increased the number of emitters along each LOS by creating an ensemble of $`50`$ 2D Gaussian velocity fields of size $`1536^2`$. All with spectral index of $`8/3`$, Kolmogorov in 2D. Analogously to the 3D ensemble we simulated observations and analysed the power spectrum in velocity slices. In this case instead of PPV cubes we have simply ‘PV arrays’, and the power spectrum computed in slices is one-dimensional. In Fig. 5 we compare the ensemble average power-spectrum in velocity channels for both, 3D and 2D cases.

It’s clear the increase in wave-number range acquired with the 2D data, in which we have about half decade more. The relevant result of this comparison is to show that the flattening of the power spectrum for the thinnest channels is present in both cases, at large $`k`$. But for the 3D case this flattening is present earlier than in the 2D case.

It’s important to stress that this is NOT a real dynamical effect, the cubes used to construct the PPV cubes are just Gaussian velocity and constant density.

Another way to alleviate the problem of limited number of emitters is to include an intrinsic width to each emitting element. We acomplished this task by adding some Gaussian smoothing on each spectral line (i. e. convolving a Gaussian profile along each LOS). This is somewhat equivalent to introduce a thermal width to the emitters666We used this width as free parameter, therefore we rather use the term ‘intrinsic’ instead of ‘thermal’ width, wich stricltly speaking, is determined by physical mechanisms beyond the scope of our simulations.. Some examples of the spectral lines in our PPV cubes, with and without an intrinsic width are shown in Figure 6.

Introducing an intrinsic width to the PPV cubes makes the spectral lines to look more realistic because it smooths sharp edges due to the limited numerical resolution. However, it also limits the resolution that can be used to look for thin channel asymptotics. The criterion in LP00 for thin or thick channels, was based on an ‘effective width’ of the velocity channels, which included not only the thickness of the channel but also a thermal width $`v_T`$ as:

$$\delta v_{eff}\left(\delta v^2+2v_T^2\right)^{1/2}.$$
(12)
The additional intrinsic width limits the resolution attainable in two ways. For a fixed maximum $`k`$ scale the velocity channels need to be narrower to remain in the thin regime, with a minimum number of channels given by:

$$N\left[\frac{\sigma _L^2}{\mathrm{\Delta }v^2}\left(\frac{r}{L}\right)^m2\frac{v_T^2}{\mathrm{\Delta }v^2}\right]^{1/2}.$$
(13)
And for a fixed velocity resolution it reduces the wave-number up to which the channel remains thin to

$$k2\pi \left[\frac{1}{\sigma _L^2}\left(\frac{\mathrm{\Delta }v^2}{N^2}+2v_T^2\right)\right]^{1/m}.$$
(14)
The criteria in eq.(14) is more restrictive, beyond that scale the velocity channels are entering in a transition to the thick regime. The condition in eq.(13) is less restrictive, in the sense that any slice with that minimum number of channels is to be considered thin. However, there is no reason to slice narrower than that, beacuse the thermal width smears velocity fluctuations on smaller scales. If we go to yet thinner channels we reduce the amplitude of the signal, but do not see any further change of spectral index. In Figure 7 we show the power spectra in velocity channels for PPV cubes where we included an intrinsic width similar to the width of the barely thin channels (thin down to the scale corresponding to the limit we chose previously to measure the power spectra, $`\mathrm{log}k2.2`$).

From the figure we can notice immediately that the rising tail at high wave-numbers dramatically decreases. In the cases where the velocity fluctuations at small scale structure are relativelly small (in Fig. 7, $`\mu =4`$) the transition from thin to thick channels as we go into smaller scales (large $`k`$) is evident. This corresponds to predictions in LP00. We computed the spectral indices in velocity channels, for our previous inertial range with slices of thickness according to eq.(13), and also for the conservative inertial range from eq.(14) for the same velocity resolutions as before. In both cases we observe that the measured spectral indices doesn’t change significantly from the earlier estimates with barely thin channels (presented in Table 3). This is understandable because the intrinsic width we included alleviates the problem of small number of emitters only on scales smaller than such width.

This shot-like noise presents an obstacle if we want to obtain the velocity index of data cubes generated synthetically. However, we must stress that the problem is not expected to be present in real observations, where the number of emitters along the LOS is much larger. This was confirmed through observations of the Small Magellanic Cloud (Stanimirović & Lazarian, 2001). The spectral indices there, measured in thinner and thinner channels do tend to an asymptotic value, not being affected by ‘shot-like noise’ allowed to measure accurately the power-spectrum in velocity channels.

## 4 Effect of the shear on the VCA

The VCA is sensitive to density and velocity fluctuations at a particular scale. On the other hand, shear as any organized motion, like galactic rotation, would correspond to the smaller wave-numbers (the largest scales, sometimes much larger that the largest scale in a PPV cube). If we choose the right range in $`k`$ to perform the VCA, we should avoid any interference of this ordered motions with the analysis. To explore numerically the effect of a linear shear on the VCA we used modified data cubes of spectral indices of $`11/3`$ for both velocity and density fields. We chose the LOS to be along the $`z`$ axis and added a component to the velocity cube (the $`v_z`$ cube) following

$$v_z(x,y,z)=v_{0z}(x,y,z)+C\sigma _{v_z}\frac{x}{L}\widehat{𝒛},$$
(15)
where $`v_z`$ is the new velocity cube used to generate the PPV cubes, $`v_{0z}`$ is the original velocity cube before introducing the shear, $`\sigma _{v_z}`$ is the velocity dispersion of $`v_{0z}`$, $`L`$ is the cube size in the spatial coordinates ($`x`$, $`y`$), and $`C`$ is a parameter we introduce to vary the magnitude of the shear introduced. The shear introduced is the same on all scales, and equal in magnitude, to $`C\sigma _{v_z}/L`$. We generated PPV cubes adding a linear shear with magnitude up to $`1`$, $`3`$, $`5`$, $`10`$, and $`20`$ times maximum velocity dispersion over the largest scale (i. e. with a velocity component with $`C=1`$, $`3`$, $`5`$, $`10`$, and $`20`$ in eq. 15). As before, we measure the spectral indices in the range less affected by numerical noise, $`\mathrm{log}k`$ from $`1.3`$ to $`2.2`$. The range used is somewhat arbitrary, but the purpose of the study at this point is just comparison with the case of no shear. The measured slopes are summarized in Table 4.

The introduced shear is the same on all scales, while the turbulent shear scales as $`k^{2/(m2)}`$. The magnitude of the external shear is the same as the magnitude of the turbulent shear when $`k2\pi C^{3/2}`$. Therefore the external shear contributed initially to the power spectrum only at the largest scales, but eventually became important at intermediate to small scales as its magnitude increased (see Fig. 8).

Our results show that the effect of an external shear becomes important only when it is comparable to the local turbulent shear. In our example that was when the maximum shear reached about $`5`$ times the velocity dispersion over the largest scale. The relative differences of the spectral indices measured compared with the case of no shear are below $`10\%`$, for the cases of less shear ($`C5`$). It is noticeable that the larger distortion happens at lower $`k`$, even when its magnitude is comparable to the turbulent shear at smaller scales. This is expected, given that the shear is regular, and represents a large scale motion. And holds if we measure the slope of the power spectrum in intermediate spatial scales, which are available using synthetic data. With actual observations, we can go further into smaller scales and avoid the effects of a stronger shear. The shear due to turbulence dominates in the galaxy, and its the effect to the VCA is marginal.

## 5 Anisotropic turbulent cascade

Magnetic field plays a crucial role in the dynamics of interstellar turbulence, and makes the turbulent cascade anisotropic (Montgomery, 1982; Higdon, 1984). In a turbulent magnetized media, the energy on the large scale motions is larger than that on the small scale structures; however the local magnetic field strength remains almost the same. Thus, is easier to bent magnetic field lines on the large scale by turbulent motions, but on the small scale becomes more difficult. Another way to think about this phenomenon is noting that hydrodynamic motions can easily mix magnetic field lines on the small scale, rather than bent them, in the direction perpendicular to the mean magnetic field (see a discussion in Lazarian & Vishniac 1999). As a result we get elongated eddies, relative to the magnetic field that become even more elongated as we go further into smaller scales. A model for incompressible anisotropic turbulence is that of Goldreich & Sridhar (1995). This model has been supported by numerical simulations (e.g. Cho & Vishniac 2000; Maron & Goldreich 2001; Cho, Lazarian & Vishniac 2002). Compressible anisotropic MHD turbulence has been studied in Lithwick & Goldreich (2001) and Cho & Lazarian (2002, 2003).

### 5.1 Magnetic Field Direction: Anisotropy in Correlation Functions

In an isotropic turbulent cascade correlations depend only on the separation between points. Contours of equal correlation are circular in that case. If on the other hand the turbulence is anisotropic, two point statistics are also anisotropic. In this case contours of iso-correlation become elongated, with symmetry axis given by the magnetic field direction. The technique presented here is similar to one proposed to study magnetic field direction using synchrotron maps (Lazarian, 1992). Iso-contours of correlation functions can allow us to determine the direction of the magnetic field. We calculate two point correlation for different 2D maps generated with our data. The correlation function between the two scalar functions $`f(𝒙)`$ and $`g(𝒙)`$ can be computed from

$$C(𝒓)=\frac{(f(𝒙)f(𝒙))\mathbf{}(g(𝒙+𝒓)g(𝒙))}{\sigma _f\sigma _g},$$
(16)
where $`𝒓`$ is the separation between two points or ‘lag’, $`𝒙`$ is the spatial separation, $`\sigma _f`$ and $`\sigma _g`$ are the standard deviations of $`f`$ and $`g`$ respectively. Here the average is performed on all the radial components and azimuthal angles.

With the original cubes from simulations (no artificial modification of the spectral indices), we calculated the auto-correlation (correlation in eq. (16) with $`f(𝑿)=g(𝑿)`$) on various 2D maps: integrated intensity (integrated along the LOS, equivalent to a 2D map of column density); emissivity at the centroid of the lines, obtained as:

$$j_c(𝑿)=\frac{j(𝑿)vdv}{j(𝑿)dv},$$
(17)
where the emissivity $`j(𝑿)`$ in our case is proportional to the first power of the density; and maps of emissivity in individual velocity channels. Contour maps of iso-correlations calculated of these maps are given in Fig. 9.

We can see that anisotropy is present in the contours plotted. It is evident in the case of the emissivity at the centroids of the lines, as well as in the integrated emissivity. In individual velocity slices, we found that the majority of the velocity channels show this anisotropy in the expected direction (elongation of the contours parallel to the mean B field), but in some slices we see the elongation of correlation contours tilted from the mean magnetic field direction. For Fig. 9 we choose a representative slice to illustrate this misleading behaviour.

The main purpose of this subsection was to show an alternative method to obtain the direction of the magnetic field, and also to illustrate that anisotropic cascade is present in our simulations.

The anisotropy based technique could be a unique tool for magnetic field studies, especially when other techniques, e.g. based on dust alignment fail (see Lazarian 2000 for a review on grain alignment).

### 5.2 Anisotropic cascade & VCA.

It’s important to test the VCA, which was developed considering isotropic turbulence, in a more realistic case where the magnetic field determines a preferential direction of the turbulent motions. With that purpose we constructed PPV cubes choosing the LOS to be parallel to the mean magnetic field (along the $`x`$ axis) and compare the results with the case where the magnetic field is perpendicular to the LOS ($`z`$ axis). For this exercise we use a data cubes corrected to Kolmogorov spectral indices ($`11/3`$ in 3D). The power spectrum in velocity channels is shown in Fig. 10.

In the figure we can see that the slopes of the spectrum (spectral index) are almost identical whether the observer is looking parallel or perpendicular to the mean magnetic field. We performed a least-squares linear fitting and found the differences of the slopes in every case are well within the fitting uncertainties.

We conclude that VCA is applicable to anisotropic turbulence studies.

## 6 Summary and conclusions

In this paper we have tested the VCA technique numerically. We used compressible MHD simulations, and an ensemble of Gaussian fields to emulate spectroscopic observations (PPV data cubes). We computed the power spectra in velocity channels to test the analytical predictions in LP00.

To obtain the PPV data we modified the output of the simulations to force power-law indices to obtain density fields with spectral indices of $`2.5`$ and $`11/3`$, and velocity fields with spectral indices of $`3.3`$ to $`4.0`$). The spectral indices derived from our emulated observations are within a $`15\%`$ error of the analytical predictions. We confirmed trends predicted in LP00, like the steepening of the emissivity index as we increase the width of the velocity channels, or that the velocity fluctuations dominate the emissivity for thin velocity channels. This last result is in agreement with a previous study by Pichardo et al. (2000), and is important because warns observers against interpret any structure in PPV as density enhancements or ‘clouds’ (see also discussion in LP00).

We showed that the deviations from power-law behaviour at the largest spatial scales (small $`k`$ numbers) are due to lack of good statistics at those scales. This was comparing the results obtained with PPV produced with a constant density and the modified velocity field of spectral index $`11/3`$, to those using an ensemble of $`60`$ Gaussian cubes with the same index.

An important source of error addressed is a flattening of the power spectrum at large wave-numbers. This is due to the limited number of emitters along the LOS. We illustrated this point comparing the results obtained with the ensemble of $`60`$ cubes with $`216`$ emitters along the LOS, and an ensemble (two-deimensional) with $`1536`$ emitters. We also try including an intrinsic width to the emitters along the LOS to smooth the sharp edges arising from the limited number of emitters. This procedure reduced the flattening of the spectrum at large $`k`$ numbers. But, also increases the effective width of the channels, reducing our ability to make use of the scales on wich we reduced such flattening of the power spectrum. This problem was found to be the most restrictive to our work. However, is exclusive of numerical simulations, and is not expected to be important in real observations (like those presented in Stanimirović & Lazarian 2001), where the number of emitting elements is essentially infinite.

We introduced linear shear to our data, and showed that its effect only becomes important when the magnitude of the shear is comparable to the local turbulent shear. We showed, that the effect of large scale shear is marginal at small scales. And that avoiding the largest scales to computed the power-spectrum allowed to recover the turbulence statistics. This shows that for ordinary observational situations we can safely neglect shear in the VCA.

We analysed the cross-correlation for 2D maps of integrated intensity and emissivity at the centroids of the lines, and showed that they are anisotropic (elongated), with symmetry axis defined in the direction of the mean magnetic field. This constitutes a new technique to study the direction of the magnetic field. We hope that the anisotropies will reveal the magnetic field within dark clouds, where grain alignment and therefore polarimetry fails. From the point of view of VCA, we illustrated that the effect of the anisotropy of the turbulence is marginal.

## Acknowledgments

We thank the referee for valuable suggestions that improved this work. AE acknowledges financial support from CONACyT (Mexico). Research by AL and JC is supported by NSF grant AST 01-25544.