# Turbulent Bipolar Outflows in Young Stellar Objects: Multifractal Universality Classes and Generalized Scaling

## 1 Introduction

The presence of bipolar outflows in the vicinity of young stellar objects (YSOs) is a generic feature of early stellar evolution. The outflows consist mostly of cold molecular gas, ionized material, and dust grains ejected with velocity in the range 100-300 km/s. A second characteristic of the circumstellar environment is a Keplerian accretion disk surrounding the protostar; the existence of these optically thick disks has been suggested by various observations, such as the high degree of linear polarization shown to be the result of multiple scattering by material in the disk (Bastien & Ménard 1990), and the absence of a red-shifted component in forbidden-line emission (Appenzeller et al. 1984; Edwards et al. 1987).

The observation that the thrust of the outflowing material cannot be radiatively driven (Bally & Lada 1983; Cabrit 1989) suggests that the flows are hydromagnetic in nature. Moreover, experimental evidence of correlations between the physical properties of the disk and those of the ejected material hint that the accretion disk is a necessary ingredient in the driving mechanism. For example, continuum observations have indicated that YSOs with outflows of larger spatial extent correspond to stronger millimeter fluxes, which in turn can be related to the presence of a more massive disk (Cabrit & André 1991; Bontemps et al. 1996). Furthermore, there are correlations between the forbidden-line emission profiles which are generated in the outflows, and near-infrared excess believed to be related to the accretion flow evolution (Gomez de Castro & Pudritz 1992; Edwards et al. 1987).

Current models for the origin of the outflows involve magnetohydrodynamics (MHD) winds driven by the accretion disk, and can be divided into two categories. The first class of models postulates that the winds are driven by a magnetic field threading through the disk (Pelletier & Pudritz 1992; Ferreira & Pelletier 1993a,b, 1995): particles rotating in the vicinity of the disk are subjected to a centrifugal force whose component perpendicular to the magnetic field lines may be neglected (assuming the magnetic field to be sufficiently strong); this results in a (centrifugal) acceleration along the field lines, thus generating the outflows. A second popular interpretation is the so-called X-celerator model (Shu et al. 1988, 1994a,b); in the region of the stellar equatorial plane where the radial acceleration is approximately null (due to counterbalance of the gravitational and centrifugal forces), the strong magnetic field of the YSO takes over the accretion mechanism, thus generating outflows. For the sake of briefness, we shall refer to these two categories of models as the PP and Shu models, respectively.

It is clear that further experimental data and analysis methods are needed to clarify the situation. In this spirit, we believe that a study of the scaling properties of the outflows could provide complementary information on the dynamics involved. The search for scaling in the outflows is motivated by the observation that over all scales, the equations of MHD for a non-dissipative medium present no characteristic length (e.g., Carbone 1993); in addition, since the law of gravity is a scale invariant (power law) form, the gravitational field of the disk-protostar system will not introduce a characteristic scale, and hence a priori it is possible to obtain scale invariance over all scales ranging from the macroscopic size of the lobes down to the much smaller dissipation scale (i.e., the inner scale of the inertial range). A fundamental observation in the analysis of three-dimensional MHD turbulence is the existence of three physical quantities, namely the total energy (i.e. the sum of the magnetic and the kinetic energies), the cross-helicity, and the magnetic helicity <sup>1</sup><sup>1</sup>1The cross-helicity and magnetic helicity densities are defined by $`\stackrel{}{v}\stackrel{}{B}`$ and $`\stackrel{}{A}\stackrel{}{B}`$, respectively, where $`\stackrel{}{v}`$ is the velocity field, $`\stackrel{}{B}`$ is the magnetic field, and $`\stackrel{}{A}`$ is the vector potential., which are conserved by the non-linear terms in the ideal equations (i.e., without forcing terms, nor dissipation). Considering the disk-protostar system as injecting the corresponding fluxes at large scales, they subsequently propagate (or cascade) towards smaller scales until the scale where dissipation becomes important is reached. Furthermore, such propagation mechanisms are more efficient between scales of similar magnitude (usual MHD turbulence is correspondingly “local” in Fourier space; e.g., Biskamp 1993; Carbone et al. 1996), and over the inertial range, the (statistical) laws governing the propagation of flux are scale invariant. The combination of conserved fluxes and local interactions (in scale) is the basis for phenomenological cascade models of turbulence. In the case of YSOs, since MHD turbulence implies the existence of three conserved fluxes in the inertial range, the dynamics of the outflows can approximatively be described by three non-linearly coupled cascade processes (see Schertzer & Lovejoy 1995; Schertzer et al. 1997b), and, in the approximation that dust grains constitute a passive scalar quantity (i.e., they are advected by the velocity field without disturbing it), there will be an additional coupled cascade corresponding to the conserved flux of passive scalar variance (Obukhov 1949; Corrsin 1951).

While the idea of cascades in hydrodynamic turbulence was introduced by Richardson (1922), the first explicit cascade models were not developed until the 1960’s (Novikov & Stewart 1964; Yaglom 1966; Mandelbrot 1974), and have become since then the basic tools for studying turbulent intermittency. Developments in the following two decades led to the current understanding that cascade models generically produce multifractals (Schertzer & Lovejoy 1985, 1987b), hence establishing their relevance in analyzing and modeling scale invariant multifractal fields. Indeed, multifractals have already been applied to various astrophysical problems such as the large-scale structure of the universe (e.g., Wiedenmann et al. 1990; Coleman & Pietronero 1992; Borgani et al. 1993; Martinez & Coles 1994; Garrido et al. 1996; Sylos Labini & Montuori 1998; Sylos Labini et al. 1998; Lovejoy, Garrido & Schertzer 1999), Ly$`\alpha `$ clouds (Carbone & Savaglio 1996), the cosmic microwave background radiation (Pompilio et al. 1995), photospheric magnetic fields (Cadavid et al. 1994), photometric data of NGC 4151 (Longo et al. 1996), and the solar wind (Carbone 1993, Politano & Pouquet 1995).

A basic difficulty in multifractal analysis and modeling is that at a purely general level, multifractals implicitly involve an infinite number of parameters (e.g., the codimension function), and hence would be unmanageable if no further simplification could be made. Fortunately, there exist multifractal universality classes, which are stable attractors of multiplicative cascades (Schertzer & Lovejoy 1987b, 1989a,b, 1991, 1997a). Universality is of practical importance, since it reduces the number of parameters required for the description of the scaling function of multifractal fields to only three.

Although the cascades and corresponding multifractals usually discussed in the literature involve isotropic scaling, physical systems generally exhibit preferred spatial directions (for example, almost all scaling geophysical systems are strongly stratified due to gravity). The need to handle scaling freed from the constraints of isotropy led to the development of the formalism of Generalized Scale Invariance (hereafter GSI, Schertzer & Lovejoy 1985, 1987a,b, 1989a,b, 1991; Lovejoy & Schertzer 1985) which is the most general framework for describing anisotropic scaling. In the case of bipolar nebulae, the application of GSI to the study of their scaling properties is motivated by the observation that the outflowing material is not isotropically ejected from the disk-protostar system, and the anticipation that the direction and strength of the anisotropy may vary with scale.

In this paper, we present a multifractal analysis of near-infrared light scattered from dust grains in bipolar outflows, using images of V380 Orionis, V645 Cygni (GL 2789), LkH$`\alpha `$ 101/NGC 1579, LkH$`\alpha `$ 233, PV Cephei (for a discussion of their physical parameters, see Bastien & Ménard 1990), V633 Cassiopeiae (e.g., Asselin et al. 1996), and GGD 18 (Gyulbudaghian et al. 1978). An immediate issue in interpreting scattered light field is the relationship of radiative transfer to the density field of the lobes; indeed, this basic remote sensing problem of radiative transfer through multifractal clouds constitutes an important application of multifractals (e.g., Lovejoy et al. 1995; Lovejoy & Schertzer 1995). Since the radiative transfer equation has no characteristic scale, if the spatial distribution of scatterers is scale invariant so will the related radiation field. It is therefore reasonable to infer that if the scattered light field displays scale invariance over a wide range of scales, then the underlying matter distribution will also possess scale invariant statistics. We shall come back to previous attempts of quantifying the correlations between the respective multifractal parameters of the scattered light and matter fields.

Information concerning the acquisition of CCD images is provided in section 2. In section 3, the basic concepts of universality are presented, while section 4 summarizes the results of the universal parameters for the seven objects in the ensemble. An overview of the GSI formalism is the subject of section 5, while the measured anisotropy parameters of GGD 18 are presented in section 6. Finally, we discuss in section 7 the potential implications of similar analyses on current and future models of star formation.

## 2 Observations

Images of V380 Orionis, V645 Cyg (GL 2789), LkH$`\alpha `$ 101/NGC 1579, LkH$`\alpha `$ 233, PV Cep, and V633 Cas were obtained at the f/15 focus of the 1.6m Ritchey-Chrétien telescope at the Observatoire Astronomique du Mont Mégantic (hereafter OMM). The first five objects were observed on 1997 October 23 in the I bandpass using a 2048x2048 pixel<sup>2</sup>, 16-bit Loral CCD camera with a scale of 0.13” per pixel; with this optical configuration and at the time of observation, the seeing was of the order of 1.5”. On the other hand, V633 Cas was observed on 1989 September 25 in the I bandpass using a 512x320 pixel<sup>2</sup> RCA chip with 0.48” per pixel. The seeing for these images was of the order of 1.3”. Images of GGD 18 were obtained at the f/8 focus of the 3.6m Canada-France-Hawai Telescope (hereafter CFHT) with the RCA2 1024x640 pixel<sup>2</sup> CCD for a scale of 0.108” per pixel, during the period 1987 December 23-28; the seeing for these images was 0.5”. All images were processed to correct for cosmetic defects of the CCDs, for reading noise, and for cosmic rays. The resolution of images of GGD 18 was increased by deconvoluting the images with the method of maximal entropy (Gravel 1990), resulting in an effective seeing of 0.39”. Acquired images and exposures for all the afore-mentioned objects are summarized in Table 1.

As YSOs generically occupied only a fraction of the region of observation, it was necessary to cut the images into sections for the purpose of our analysis. Being physically distinct from the dust shells, background stars as well as the protostar were always excluded from these sub-images since their presence would bias the analysis of scattered light. Sections were also chosen to be far enough from the protostar to avoid problems of large intensity gradients associated with the protostar, such that the cloud–radiation physics (and hence statistics) for a given sub-image could be considered approximately invariant under translations. Finally, regions where the intensity was not significantly greater than the background noise were also excluded. A contour plot of a sample CCD image acquired at the OMM, namely PV Cep, is shown in Figure 1, where the box delimits the section used for the analysis.

Figure 2 presents a contour plot of the region of GGD 18 analyzed, divided into sub-regions that will be used for the subsequent multifractal analysis (sections 4 and 6).

The image of GGD 18 was the only one in the ensemble with sufficient intensity and spatial resolutions to allow a GSI analysis (see section 6), and it is thus worthwhile to describe some of its known features. This object is located at 30.25” NW of the binary GL 961 (Cohen 1973), and was first identified as a possible candidate for a HH object (Gyulbudaghian et al. 1978). A few years later, Lenzen et al. (1984) discovered an IR source embedded in the nebulous region of GGD 18 (it corresponds to the maximum of intensity of Figure 2, identified by a small cross), but did not speculate on its nature. Polarization maps from deconvoluted images revealed (Gravel 1990) a centrosymmetric pattern centered on the GGD 18 source showing that it was a YSO, along with a region of aligned vectors near the source indicating the presence of an optically thick disk (the annular structure near the source and along the NE and SW directions is believed to correspond to a fraction of the disk that is visible); the inclination of the latter with respect to the line of sight was estimated at $`30^o`$ (the northern lobe points toward the observer), and the symmetry axis of the nebula points at $`30^o`$ NW (Gravel 1990). Finally, spiral-like polarization patterns were discovered in the vicinity of the source, indicating a complex distribution of matter that may be the result of rotation in the material near the source, combined with the effects of possible local inhomogeneities in the disk. An unresolved issue is the extent at which the dynamics of the outflows are influenced by the proximity of the binary GL 961; evidence that the CO bipolar jet of GL 961 reaches and influences the outflow dynamics of GGD 18 was provided by Gravel (1990), but a quantitative study of such interactions is still needed. To summarize his argument, we notice that both lobes are more extended in the western direction, and the apparent contraction on the eastern side is believed to result from interactions with the jet; the latter appears as a faint island (approximately 12” in length) located in (1,1), and points in the NW direction. It is more visible in (1,2) and enters both (1,3) and (2,3), where it turns around towards the SW direction; such deviations could be the result of interactions with nebulous structures to the north of GGD 18. The local maximum in section (0,0) corresponds to the western component of the binary GL 961; since GL 961 is an intense source of unscattered light, the region (0,0) will be excluded to allow an unbiased analysis of scattered light in the outflows. For similar reasons, section (2,0) will also be excluded from the analysis.

## 3 Review of Multifractal Processes

### 3.1 Universal Multifractals

It is well known that the Navier-Stokes (NS) equations are invariant under the rescaling $`xx\lambda ^1`$, $`vv\lambda ^H`$, and $`tt\lambda ^{H1}`$. Assuming that $`ϵ`$, the energy flux to smaller scales, is a scale invariant quantity, it is found that $`H=1/3`$, and dimensional analysis leads to the famous scaling law $`E(k)k^{5/3}`$ for the energy density in momentum space (Komolgorov 1941; hereafter K41). While the equations of MHD satisfy similar scaling relations, there is no consensus on what the analogue of K41 should be (the corresponding dimensional analysis no longer gives a unique dimensional combination, hence unique exponent). In terms of Elsasser variables, eddies fall into two classes depending on their direction of propagation along the magnetic field lines: interactions between eddies belonging to different classes are less likely, thus weakening the energy transfer (Iroshnikov 1963; Kraichman 1965; hereafter IK). Consequently, the characteristic interaction time $`\tau _{eddy}`$ (i.e. the eddy turnover time) is increased to $`(\tau _{eddy}/\tau _A)^a`$, where $`\tau _A`$ is the characteristic time for Alfvén waves, and $`a`$ is some positive constant (Politano & Pouquet 1995). The scaling relation becomes:

$$E(k)k^{\left(1+\frac{2}{a+3}\right)},$$
(1)   
with $`a=1`$ corresponding to the IK theory.

It is found experimentally that in hydrodynamical media K41 is generally not respected for individual realizations – even in estimates of the ensemble average, the exponent differs from 5/3; the discrepancy can be attributed to intermittency (fluctuations in $`ϵ`$ due to small-scale non-linear structures). As discussed in the introduction, multiplicative cascades model the propagation of conserved fluxes in turbulent intermittent media, and in general present the expected characteristics of a fully developped turbulent field. The general outcome of such cascades is a multifractal field which is scale invariant over the inertial range, and whose flux density at resolution $`\lambda `$, denoted by $`ϵ_\lambda `$, is described by (Schertzer & Lovejoy 1987b):

$$ϵ_\lambda ^q=\lambda ^{K(q)},$$
(2)   
where the brackets indicate an average over many realizations, and $`K(q)`$ is the moment scaling function. The arbitrariness of $`H`$ in the rescaling of the NS or MHD equations allows the possibility of scaling of different moments of the intensity spectrum, and this feature is what equation 2 describes.

While in general $`K(q)`$ need only be convex, cascades possess stable, attractive universality classes (Schertzer & Lovejoy 1987b, 1997) whose description requires only three parameters, namely $`\alpha `$, $`C_1`$ and $`H`$, with the corresponding moment scaling function determined by the universality relation:

$$K(q)qH=\frac{C_1}{\alpha 1}(q^\alpha q).$$
(3)   
The significance of each of the three universality parameters on the multifractal field can be described as follows:

* $`C_1`$ corresponds to the codimension of the mean field, and thus distinguishes between a field whose mean is dominated by a few localized intense peaks (large $`C_1`$), and one with a mean dominated by a larger proportion of its surface (small $`C_1`$ — for non-fractal such as white noise, $`C_1=0`$);
* $`H`$ is a measure of the degree of (scale by scale) non-conservation of the field, or qualitatively a measure of its smoothness (with large values of $`H`$ corresponding to smoother fields, see eqs. 4 and 9). For example, in usual hydrodynamic turbulence the energy flux to smaller scales is conserved ($`H`$=0) whereas the velocity shears have the Komolgorov value $`H=1/3`$;
* $`\alpha `$ is the degree of multifractality, i.e., a measure of the deviation from the monofractal case. As $`\alpha `$ is the Lévy index of the multifractal generator, we have the restriction $`0\alpha 2`$, with $`\alpha =0`$ and $`\alpha =2`$ corresponding to monofractal ($`\beta `$ model) and log-normal models, respectively <sup>2</sup><sup>2</sup>2Note that the frequently used expressions “log-Lévy” and “log-normal” are rather misleading because of the divergence of high order statistical moments; the statistics will only be approximately “log-normal” and “log-Lévy” up to the given critical order of divergence of the moments..

### 3.2 Analysis of Physical Fields

A preliminary verification of the existence of scale invariance is that the spectral energy density satisfies a general (isotropic) scaling law:

$$E(k)k^\beta ,$$
(4)   
where $`\beta `$ is the scaling exponent, or spectral slope.

After the existence of a scaling regime is established, one can proceed to compute the scaling function $`K(q)`$ and test for universal multifractal behavior. An efficient technique for that purpose is the Double Trace Moment (DTM) method (Lavallée 1991; Lavallée et al. 1991, 1993). A new function $`K(q,\eta )`$ is first defined similarly to $`K(q)`$ in equation 2, but for a field $`ϵ_\mathrm{\Lambda }`$ at its maximal resolution $`\mathrm{\Lambda }`$, raised to the $`\eta `$ power (i.e., $`ϵ_\mathrm{\Lambda }ϵ_\mathrm{\Lambda }^\eta `$), and renormalized by spatial averaging. Hence, writing $`(ϵ_\mathrm{\Lambda }^\eta )_\lambda ^q`$ to indicate the $`\lambda `$-resolution $`q^{\mathrm{th}}`$ moment of $`ϵ_\mathrm{\Lambda }^\eta `$, we obtain the following generalization of equation 2:

$$(ϵ_\mathrm{\Lambda }^\eta )_\lambda ^q=\lambda ^{K(q,\eta )}.$$
(5)   
While it can be shown (Lavallée 1991) that $`K(q,\eta )`$ and $`K(q)`$ are related by

$$K(q,\eta )=K(q\eta )qK(\eta ),$$
(6)   
the advantage of the DTM technique for testing and characterizing universality compared with other methods (e.g., Schmitt et al. 1995) is realized when $`K(q)`$ is universal (eq. 3), in which case the $`\eta `$-dependence factorizes:

$$K(q,\eta )=\eta ^\alpha K(q).$$
(7)   
The DTM method allows the computation of $`K(q,\eta )`$, and assuming universality, the Lévy index $`\alpha `$ can be deduced from equation 7. The remaining parameters $`C_1`$ and $`H`$ are also determined from a knowledge of $`K(q,\eta )`$. Explicitly, one finds:

$$C_1=(\alpha 1)\frac{K(q,1)}{q^\alpha q},$$
(8)   
and

$$H=\frac{\beta 1}{2}C_1\frac{2^{\alpha 1}1}{\alpha 1}.$$
(9)   
We conclude this section with a few comments concerning the range of validity of the above equations when analyzing images of physical fields. In practice, the finite size of the sample implies that sufficiently high order moments are dominated by the largest value assumed by the field, and therefore underestimate the true ensemble moments. Beyond some threshold $`q_s`$, which for a single realization is given by

$$q_s=\left(\frac{d}{C_1}\right)^{1/\alpha },$$
(10)   
equation 3 is no longer expected to hold true, and $`K(q)`$ becomes linear. Here, $`d`$ is the dimension of the space over which the analysis is made ($`d=2`$ for the images discussed here). One therefore encounters a “multifractal phase transition” (Schertzer et al. 1992).

## 4 Double Trace Moments Results

As discussed in section 3, scale invariance is a necessary condition for a physical field to be multifractal, and its existence should be verified before computing the universal parameters. A first analysis is to consider the isotropic power spectrum since the spectral energy density $`E(k)`$ of an isotropic MHD turbulent medium is expected to obey a scaling law with exponent $`\beta `$ (see eq. 4). Examples of the power spectra obtained are shown in Figure 3a) and b) for a sub-region of GGD 18 and LkH$`\alpha `$ 101/NGC 1579, respectively, where $`E(k)`$ has arbitrary units.

In each case, we note a break in scaling occurring at small scales (i.e., large wavenumber), followed by a regime of constant or increasing energy density. We shall argue that the scale at which scaling breaks corresponds to the resolution at which structures of the intensity field become dominated by noise. This is easily seen in the case of GGD 18, where the break occurs at a resolution of approximately 0.3”, while the seeing at the time of data acquisition was estimated at 0.4”. On the other hand, the high frequency behavior of $`E(k)`$ for the spectrum of LkH$`\alpha `$ 101 is linear ($`\beta 1`$) in wavenumber, as expected for Gaussian white noise in 2D space <sup>3</sup><sup>3</sup>3In the case of LkH$`\alpha `$ 101, the fact that the break is not a manifestation of the seeing results from the limited exposures compared to those of GGD 18.. These two sample power spectra illustrate a feature common to all power spectra in the ensemble, namely isotropic scaling down to the scale where noise dominates the statistics (the scaling anisotropy – see sections 5 and 6 – is removed by the angular integration). Least-squares fits over the respective linear regions of Figure 3a) and b) yield $`\beta `$=2.4$`\pm `$0.3 and 2.0$`\pm `$0.2, respectively, where the uncertainties are estimated from the fitting procedure.

With the existence of scaling established, a DTM analysis was performed on each sub-image in the ensemble using four values of $`q`$, namely $`q`$=0.5, 0.6, 0.75, and 0.9. As an example of the DTM results, Figure 4 presents $`K(\eta ,q)`$ with $`q=0.5`$ for a sub-image of LkH$`\alpha `$ 233, from which we note a power-law dependence over the range 0.06$`\eta `$1, in agreement with equation 7.

The slope of the linear region yields $`\alpha =1.93`$, while the intercept $`K(q,1)`$ gives $`C_1=0.023`$ (eq. 8). Finally, there is a departure from the power-law dependence in the range $`4\eta 10`$, beyond which $`K(\eta ,q)`$ becomes independent of $`\eta `$. As explained in section 3, equation 7 breaks down beyond max$`(q\eta ,\eta )q_s`$ due to the finite size of analyzed samples. To confirm that this is indeed the cause for the multifractal phase transition, the substitution of the above values of $`\alpha `$ and $`C_1`$ in equation 10 gives $`q_s10`$; on the other hand, $`K(\eta ,q)`$ becomes horizontal at $`\eta 10`$, which implies that max$`(q\eta ,\eta )=10`$ (recall that $`q=0.5`$), as expected.

Variations of the universality parameters over different sub-regions of a given nebula were in general observed to be within the known uncertainties, and this was explicitly verified in the case of GGD 18. Consequently, it was reasonable to average the universal parameters over all sub-regions for each object, and the results are summarized in Table 2.

Note that the accuracy on the parameters of individual objects can be roughly estimated from the accuracy of the DTM method (measured from numerical simulations in Lavallée et al. 1991) to be $`\pm 0.1`$ for $`\alpha `$, $`\pm 0.05`$ for $`C_1`$, and $`\pm 0.2`$ for $`H`$. An immediate observation from Table 2 is that $`\alpha `$ is very close to 2 for every object in the ensemble, a value which corresponds to the highest degree of multifractality. A second observation is the uniformity of the parameters over the ensemble, allowing one to compute ensemble averages listed in Table 2, where the uncertainties quoted for the ensemble averages correspond to one standard deviation from the mean. While these parameters describe the statistics of the field of scattered light, it is important to consider their relationship with those of the underlying field of matter; simulations (Naud et al. 1996) and simple theoretical arguments (Schertzer et al. 1997) have suggested that only $`H`$ is significantly affected by the scattering process, with the radiative value observed to be larger (corresponding to a smoother texture).

Recall from section 3 that equation 7 follows from equation 6 provided that the scaling function $`K(q)`$ is universal. To test the validity of this assumption, one can compute $`K(q)`$ directly and compare the result with the prediction of the universality relation (eq. 3). Such a comparison is illustrated in Figure 5 for a sub-region of GGD 18.

According to Figure 5, we note that the two curves agree up to $`q5`$, beyond which $`K(q)`$ becomes linear; from equation 10 we find $`q_s6`$, indicating that the linear dependence is another manifestation of the finiteness of the sample. Similar agreement was observed for the other images in the ensemble.

## 5 Generalized Scale Invariance

Our discussion has so far been constrained to the case of isotropic, or self-similar, scale invariance, for which structures at different scales are related by isotropic magnifications. However, more general anisotropic scaling “zooms” are possible and may arise as a result of forces inducing stratification and differential rotation in the dynamics, for instance. Rather than imposing, a priori, an isotropic notion of scale, the latter may be determined by the nonlinear dynamics; in this section, we attempt to empirically characterize this anisotropic scaling.

The need for a more general framework led to the Generalized Scale Invariance (GSI) formalism, in which three ingredients are necessary for the description of scaling: (i) a unit ball $`B_1`$ consisting of all vectors of unit length, which can implicitly be defined by:

$$B_1\{\stackrel{}{x}:g_1(\stackrel{}{x})1\},$$
(11)   
where $`g_1`$ determines the notion of unit length; (ii) a scale changing operator $`T_\lambda `$ mapping a vector between two scales of ratio $`\lambda `$; once a unit ball is specified, all other scales can be identified by repeated applications of $`T_\lambda `$, thus generating a family of balls, $`B_\lambda `$; (iii) a definition of a measure of the $`B_\lambda `$, such as the volume or a power of the volume.

It follows from the definition of $`T_\lambda `$ that these operators form a one-parameter multiplicative group, with corresponding generator $`G`$ given by:

$$T_\lambda =\lambda ^G.$$
(12)   
While equation 12 allows $`G`$ to be non-linear (Schertzer & Lovejoy 1985), in order to simplify the computations we shall make a linear approximation by assuming that $`G`$ and $`T_\lambda `$ are real matrices with constant coefficients, which is equivalent to assuming that the anisotropy of the scales is translationally invariant within the image analyzed. Since GSI analyses are conventionally performed in Fourier space, the operator of interest is the Fourier analogue of $`G`$, which shall also be denoted $`G`$ by an abuse of notation (in the case of linear GSI, these matrices are transpose of each other). Expanding the matrix generator in terms of pseudo-quaternions, we write:

$$G=1+f\sigma _xie\sigma _y+c\sigma _z=\left(\begin{array}{cc}1+c& fe\\ f+e& 1c\end{array}\right),$$
(13)   
where $`1`$ denotes the unit matrix, $`\sigma _x`$, $`\sigma _y`$, and $`\sigma _z`$ are the $`SU(2)`$ generators in the Pauli representation, and $`c,e,f`$ are real (note that since the coefficient of the identity matrix would correspond to an isotropic magnification, it has been set equal to unity without loss of generality). This choice of basis is particularly convenient for GSI analysis as it decouples rotation and stratification. It follows from equations 12 and 13 that, within the approximation of linear GSI, scale transformations are completely determined by the specification of $`c`$, $`e`$, and $`f`$.

It is instructive at this point to discuss a few examples of $`G`$, and corresponding family of balls obtained by repeated applications of $`T_\lambda `$ on the unit ball. For the purpose of this discussion, we assume that $`G`$ is a 2x2 matrix, corresponding to the analysis of two-dimensional data with translationally invariant statistics. The simplest example is the case where $`G`$ is the identity matrix, corresponding to $`T_\lambda =\lambda ^1`$, or self-similar scaling. In addition, if the unit ball is chosen to be the unit circle, the family of balls generated by $`T_\lambda `$ are circles as well. More generally, the scaling may be different in two or more preferred directions. For simplicity, suppose there are two such directions assumed to coincide with the x and y axes, in which case the scaling is said to be self-affine: the off-diagonal elements of $`G`$ remain null, but the diagonal entries may be different from unity. If the unit ball is taken to be the unit circle, the corresponding balls are ellipses with principal axes pointing in the x and y directions. Furthermore, as illustrated in Figure 6a), one finds that the balls are horizontally elongated for large wavenumbers ($`k>1`$) and vertically elongated for small ones ($`k<1`$).

An example of physical systems presenting approximately such scaling are vertical cross-sections of the atmosphere where self-affinity is caused by the (stratifying) gravitational field of the Earth (Pflug et al. 1991, 1993). Finally, a matrix generator with non-zero off-diagonal elements indicates differential rotation in the scaling (see Figure 6b)); in the atmosphere for instance, the observed differential rotation can be generated by the Coriolis force.

The numerical calculation of the anisotropy parameters $`c`$, $`e`$, and $`f`$ is performed using the Scale Invariant Generator (henceforth, SIG) technique (Lewis et al. 1999). Let us first define $`P(\stackrel{}{k})`$ to be the modulus squared of the Fourier transform of the field at wave vector $`\stackrel{}{k}`$. It follows from this definition that an angular integration of $`P(\stackrel{}{k})`$ yields $`E(k)`$, as defined in section 3. As a generalization of equation 4, we have the scaling relation:

$$P(\stackrel{}{k})\stackrel{}{k}^s,$$
(14)   
where the norm is with respect to $`G`$, the brackets indicate an average over all realizations, and $`s`$ is a generalized scaling exponent given in (isotropic) 2-D by $`s=\beta +1`$. It follows that $`P(\stackrel{}{k})`$ is in fact a function of $`G`$, $`B_1`$, $`s`$, and $`\stackrel{}{k}`$, that is, $`P(\stackrel{}{k})P(c,e,f,B_1,s,\stackrel{}{k})`$. If the analyzed sample consists of $`N`$ data points denoted by $`\stackrel{}{k}_i`$, $`i=1,\mathrm{},N`$, we can define an error function $`E_r`$ by:

$$E_r^2(G)\frac{1}{N}\underset{i,j}{\overset{N}{}}\left[\mathrm{ln}P\left(\lambda _i^G\stackrel{}{k}_j\right)+s\mathrm{ln}\lambda _i\mathrm{ln}P(\stackrel{}{k}_j)\right]^2,$$
(15)   
where the sum is over the data points ($`\stackrel{}{k}_j`$) and scale ratio ($`\lambda _i`$). The anisotropic parameters are then estimated by minimizing this error function.

## 6 Scale Invariant Generator Results

The evaluation of the anisotropy exponents with the SIG technique requires data over a wide range of scales (sufficient resolution), and good signal-to-noise ratio, such that the (approximate) isotropic scaling is valid over many scales. As mentioned in section 2, the images obtained at the OMM did not show sufficiently good statistics for the type of analysis described in this section.

Figure 7 displays the results of the SIG analysis performed on the sub-images of GGD 18 (defined in fig. 2), where the quoted uncertainties on the parameters $`c`$, $`e`$, and $`f`$ were estimated from the numerical optimization of the error function $`E_r`$ defined in equation 15.

Since the error surface near its minimum was typically wider in the $`e`$ direction, this exponent is known with less accuracy than the other two. Note that boxes (2,0) and (0,0) were left empty in Figure 7 as they contain the source of GGD 18 and the western component of the binary GL 961. An immediate observation is the substantial variations of the parameters from one sub-region to the next; while the choice of linear GSI simplifies the numerical calculations, the fact that the generator $`G`$ is not constant over the region analyzed suggests that the framework of non-linear GSI may be more appropriate in this case. We shall nevertheless assume that linear GSI is a good approximation over each subimage.

The parameters $`e`$ and $`f`$ which determine the off-diagonal elements of $`G`$ (eq. 13), are non-zero for most sub-regions, hence providing evidence for the existence of differential rotation in the outflow dynamics. They are comparatively larger in sections (1,0), (2,1) and (3,1), which cover the portions of the southern and northern lobes that are the closest to the protostar–disk system. As discussed in section 2, boxes (1,2) and (2,3) contain parts of a jet of matter presumably emitted by GL 961 at some point in the past. The large values of the off-diagonal elements of $`G`$ in these regions suggest that there is some rotation induced in the interaction region of this jet with GGD 18 material. However, it is not clear at this moment where the rotation exactly occurs (e.g. in the jet, the perturbed material, or both) since the linear SIG technique doesn’t resolve variations of the exponents within subimages.

It has been noted by Gravel (1990) that the eastern edge of the northern lobe of GGD 18 appears to be pushed westward presumably by the jet of material of GL 961. It is interesting to note that the region (1,1) where the jet and the northern lobe presumably comes the closest, is the only one with a positive value of $`e`$.

Let us divide the image into two subregions (see Figure 7) with region I covering what is morphologically identifiable as the lobes of the nebula, and region II covering the rest of the observed field (with the exception of the emission features). The assignation of subimages to either region was performed using the following criteria: subregions belonging to region I (i) were close to the $`30^o`$ symmetry axis of the nebula, (ii) the corresponding light intensity was more important, and (iii) their intensity contours had conical shapes (indicative of ejection). With this classification, there appears to be a systematic difference in the mean value of $`c`$ between regions I and II: quantitatively, the difference in mean values, $`c_Ic_{II}`$ is found to lie 4.0 standard deviations away from the case $`c_I=c_{II}`$. Although a precise identification of the stratifying factors is beyond the scope of our phenomenological framework, the observed stratification appears primarily dynamical in nature since $`c`$ is not observed to decrease monotonically with distance from the protostar (as would be expected for gravitationally generated stratification), and since values of $`c`$ present clear differences within and outside the lobes of the nebula.

## 7 Conclusion

The statistics of scattered light in a small ensemble of YSOs have been shown to obey different scaling relations depending on the moment $`q`$ (“multiscaling”), with the relation between scaling laws and statistical moments given by $`K(q)`$. It was also shown that these $`K(q)`$ functions fit reasonably well into a multifractal universality class, and that the corresponding universal parameters are fairly uniform over the ensemble (with $`\sigma _\alpha /\alpha 1\%`$, $`\sigma _{C_1}/C_147\%`$, and $`\sigma _H/H28\%`$). Although universality may be thought of as describing the attractors of multiplicative cascades, and consequently is not sufficiently sensitive on the initial conditions of the cascade to allow the observation and characterization of fine details in the dynamics of a given multifractal field, the reasonable uniformity of the results presented suggests that the seven objects in the ensemble have similar (presumably turbulent) dynamics, as expected on theoretical grounds.

Finally, although obvious differences (e.g. the existence of a magnetic field) exist between the dynamics of the radiative fields of YSOs and that of terrestrial water clouds, our results suggest that the corresponding statistics are similar. Indeed, the spectral slope and the universality parameters of atmospheric clouds ($`\beta 2.2,\alpha =1.79,H=0.63`$ and $`C_1=0.061`$ – Sachs et al. 1999) are close to that of bipolar nebulae (see table 2). Since empirically cloud liquid water behaves statistically approximately like a passive scalar (Lovejoy & Schertzer 1995), its statistical similarity with YSOs would make sense if dust grains in bipolar outflows constituted a passive scalar as well. Furthermore, the physics of radiative transfer is also similar since in each case it is dominated by scattering rather than absorption/emission processes.

All the objects in the ensemble presented reasonable isotropic scaling, that was systematically broken near the scale where noise becomes dominant in the measured signal. Only GGD 18 presented sufficient resolution to allow a GSI analysis of its statistics. Most of the sub-regions analyzed had matrix generators with non-zero off-diagonal elements, revealing the existence of differential rotation. The origin of the latter and its influence on mass ejection mechanisms should be accounted for in models of star formation; an obvious source of rotation in YSOs is the rotation of the disk-protostar system, and some of the properties of the coupling of this rotation to the ejected material might be studied by techniques similar to those used in this work. All sub-images presented a non-zero value of $`c`$, indicating a possible stratifying force in the outflow mechanism. It was argued in section 6 that such stratification would not be primarily gravitational in nature, but instead could result from dynamical pressure gradients related to physical forces, such as the centrifugal acceleration in the PP models, or the magnetic force in the Shu models. Finally, an important issue in the dynamics of GGD 18 is its interactions, if any, with the neighbouring YSO GL 961. Showing that such interactions are indeed involved would confirm that GGD 18 is located at approximately the same distance from the Earth as GL 961, namely 1.6 kpc. We found that the generator of (1,1) was the only one to have off-diagonal elements of opposite sign; as this box is believed to contain the region of closest approach between the northern lobe of GGD 18 and the jet from GL 961, the peculiarity of its anisotropy parameters could be a sign of interactions between the two YSOs.

It should be kept in mind that linear GSI is probably not accurate enough to probe fine details in the dynamics, and is in fact increasingly understood as measuring local multifractal textures (Pecknold et al. 1996, 1997). While developments in non-linear GSI or models involving non-scalar cascades, along with an increased ensemble of sufficient resolution, are probably necessary to obtain statistically robust statements concerning the outflow dynamics, it is hoped that our analysis has provided a foretaste of the vast possibilities of multiscaling analyses.

###### Acknowledgements.

P.B. thanks the director of the Canada-France-Hawai telescope for a generous time allotment. S.L. was partly supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada.  REFERENCES  Appenzeller, I., Jankovics, I., & Ostreicher, R. 1984, A&A, 141, 108   
Asselin, L., Ménard, F., Bastien, P., Monin, J.-L., Rouan, D. 1996, ApJ, 472, 349A   
Bachiller, R. 1996, ARA&A, 34, 111   
Bally, J., & Lada, C. J. 1983, ApJ, 265, 824   
Bastien, P., & Ménard, F. 1990, ApJ, 364, 232   
Beckwith, S. V. W., Sargent, A. I., Koresko, C. D., & Weintraub, D. A. 1989, ApJ, 343, 393   
Biskamp, D. 1993, Nonlinear Magnetohydrodynamics (Cambridge: Cambridge University Press)   
Blitz, L. & Thaddeus, P. 1980, ApJ, 241, 676   
Bontemps, S., André, P., Terebey, S., & Cabrit, S. 1996, A&A, 311, 858   
Borgani, S., Murante, G., Provenzale, A., Valdarnini, R. 1993, Phys. Rev. E, 47, 3879   
Cabrit, S. 1989, in ESO Workshop on Low Mass Star Formation and Pre-Main Sequence Evolution, ed. B. Reipurth (Garching: ESO), p. 119   
Cabrit, S., & André, P. 1991, ApJ, 379, L25   
Cadavid, A. C., Lawrence, J. K., Ruzmaikin, A. A., Kayleng-Knight, A. 1994, ApJ, 429, 391   
Carbone, V., 1993, Phys. Rev. Lett., 71, 1546   
— & Savaglio, S. 1996, MNRAS, 282, 868   
—, Veltri, P., Bruno, R. 1996, Nonlin. Proc. Geophys., 3, 247   
Cohen, M. 1973, ApJ, 185, L75   
Coleman, P. H., & Pietronero, L. 1992, Phys. Rep., 213, 311   
Corrsin, S. 1951, J. Applied. Phys., 22, 469   
Edwards, S., Cabrit, S., Strom, S. E., Heyer, I., Strom, K. M., & Anderson, E. 1987, ApJ, 321, 473   
Ferreira, J. & Pelletier, G. 1993a, A&A, 276, 625   
— & — 1993b, A&A, 276, 637   
— & — 1995, A&A, 295, 807   
Garrido, P., Lovejoy, S., Schertzer, D. 1996, Physics A, 225, 294   
Gomez de Castro, & A. I., Pudritz, R. E. 1992, ApJ, 395, 501   
Gravel, P. 1990, M.Sc. thesis, Université de Montréal, Montréal, Québec, Canada   
Gyulbudaghian, A. L., Glushkov, Y. I., Demisyuk, E. K. 1978, ApJ, 224, L137   
Halsey, T. C., Jensen, M. H., Kadanoff. L. P., Procaccia, I., Schraiman, B. 1986, Phys. Rev. A, 33, 1141   
Iroshnikov, P. 1963, Astron. Zh., 40, 742   
Kraichnan, R. H. 1965, Phys. Fluids, 8, 1385   
Lada, C. J. 1985, ARA&A, 23, 267   
—, Gauthier III, T. N. 1982, ApJ, 261, 161   
Lavallée, D. 1991, Ph.D. thesis, McGill University, Montréal, Canada   
—, Schertzer, D., Lovejoy, S. 1991, in Scaling, Fractals and Non-Linear Variability in Geophysics, eds. D. Schertzer and S. Lovejoy (Klumer), p. 99   
—, Lovejoy, S., Schertzer, D., Ladoy, P. 1993, in Fractals in Geography, eds. L. De Cola, N. Lam (Prentice Hall), p. 158   
Lenzen, R., Hodapp, K. W., Reddman, T. 1984, A&A, 137, 365   
Lewis, G., Lovejoy, S., Schertzer, D., Pecknold, S. 1999, Computers in Geophys (in press)   
Longo, G., Vio, R., Paura, P., Provenzale, A., Rifatto, A. 1996, A&A, 312, 424   
Lovejoy, S., & Schertzer, D. 1985, Wat. Resour. Res., 21, 1233   
—, — 1990, Physics in Canada, 46, 4, 46   
—, Watson, B., Schertzer, D., Brosamlen, G. 1995, in Particle Transport in Stochastic Media, ed. L. Briggs (Portland: American Nuclear Society), p. 750   
—, & — 1995, in Fractals in Geoscience and Remote Sensing, ed. G. Wilkinson (Luxembourg: Office for Official Publications of the European Communities), p. 102   
—, Schertzer, D., Tessier, Y., 1998, Int. Journal. of Remote Sensing (submitted)   
—, Garrido, P., Schertzer, D. 1998, Physica A (submitted).   
Mandelbrot, B. B. 1974, J. Fluid Mech., 62, 331   
Martinez, V. J., Coles, P. 1994, ApJ, 437, 550   
Naud, C., Schertzer, D., Lovejoy, S. 1996, in Stochastic Models in Geosystems, eds. S. A. Molchansov and W. A. Woyczynski (Spinger-Verlag), p. 239   
Novikov, E. A., & Stewart, R. 1964, Izv. Akad. Nauk. SSSR. Ser. Geofiz., 3, 408   
Obukhov, A. 1949, Izv. Akad. Nauk. SSSR. Ser. Geogr. I Geofiz, 13, 55   
Parisi, G., & Frisch, U. 1985, in Turbulence and Predictability in Geophysical Fluid Dynamics and Climate Dynamics, eds. M. Ghil, R. Benzi, and G. Parisi (North-Holland), p. 72   
Pecknold, S., Lovejoy, S., Schertzer, D. 1996, in Stochastic Models in Geosystems, eds. S. A. Molchansov and W. A. Woyczynski (Spinger-Verlag), p. 269, 85   
Pecknold, S., Lovejoy, S., Schertzer, D., Hooge, C. 1997, in Scale in Geophysical Information Systems, eds. D. Quattrochi and M. F. Goodchild (Florida: CRC Press), p. 361   
Pelletier, G., & Pudritz, R. E. 1992, ApJ, 394, 117   
Pflug, K., Lovejoy, S., & Schertzer, D. 1991, in Nonlinear Dynamics of Structures, eds. R. Z. Sagdeev, U. Frisch, A. S. Moiseev, and A. Erokhim (World Scientific), p. 72   
—, —, & — 1993, J.Atmos.Sci, 50, 538   
Politano, H., & Pouquet, A. 1995, Phys. Rev. E, 52, 636   
Pompilio, M. P., Bouchet, F. R., Murante, G., Provenzale, A. 1995, ApJ, 449, 1   
Richardson, L. F. 1922, Weather prediction by numerical process (Cambridge University Press)   
Sachs, Lovejoy, S., Schertzer, D. 1999 (in preparation for Fractals)   
Schertzer, D., & Lovejoy, S. 1983, in Forth Symposium on Turbulent Shear Flows, Karlshule, Germany   
—, & — 1985, Phys. Chem. Hydrodyn., 6, 623   
—, & — 1987a, Ann. Sci. Math. du Québec, 11, 139   
—, & — 1987b, J.Geophys.Res.D., 92, 9693   
—, & — 1989a, in Fractals: Physical Origin and Consequences, ed. L.Pietronero (Plenum), p. 49   
—, & — 1989b, Pageoph, 130, 57   
—, & — 1991, in Scaling, Fractals and Non-Linear Variability in Geophysics, eds. D. Schertzer and S. Lovejoy (Kluwer), p. 41   
—, & — 1995, in Space/time Variability and Interdependence for various hydrological processes, ed. R. A. Feddes (Cambridge: Cambridge University Press), p. 153   
—, —, Schmitt, F., Chigirinskaya, Y., Marsan, D. 1997, Fractals 5, 427   
—, & — 1997a, J. Applied Meteorology, 36, 1296   
—, & — 1997b, ARM Proceedings   
Schertzer, D., Schmitt, F., Naud, C., Marsan, D., Chigirinskaya, Y., Margeurite, C., Lovejoy, S., 1997, 7$`\mathrm{th}`$ Atmos. Rad. Meas. Sci. Team Meeting, 327-335   
Schmitt, F., Lavallée, D., Schertzer, D., Lovejoy, S. 1992, Phys. Rev. Lett., 68, 305   
Schmitt, F., Lovejoy, S., Schertzer, D. 1995, Geophys. Res. Lett., 22, 1689   
Shu, F. H., Lizano, S., Ruden, S. P., Najita, J. 1988, ApJ, 328, L19   
—, Najita, J., Ostricker, E., Wilkin, F., Ruden, S., Lizano, S. 1994a, ApJ, 429, 781   
—, —, Ruden, S., Lizano, S., 1994b, ApJ, 429, 797   
Sylos Labini, F., Montuori, M. 1998, A&A, 331, 809   
Sylos Labini, F., Montuori, M., Pietronero, L., 1998, Phys. Rep., 293, 726   
Turner, D. G. 1976, ApJ, 276, 65   
Wiedenmann, G., Atmanspacher, H., Scheingraber, H. 1990, Can. Journ. Phys., 69, 9, 827   
Yaglom, A. M. 1966, Sov. Phys. Dokl., 2, 26