Issue 
A&A
Volume 501, Number 3, July III 2009



Page(s)  999  1011  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/200811553  
Published online  27 May 2009 
A differentially rotating disc in a highmass protostellar system
M. R. Pestalozzi^{1}  M. Elitzur^{2}  J. E. Conway^{3}
1  Dept. of Physics, University of Gothenburg, 412 96, Göteborg, Sweden
2 
Dept. of Physics and Astronomy, Univ. of Kentucky, Lexington, KY 405060055, USA
3 
Onsala Space Observatory, 439 92 Onsala, Sweden
Received 19 December 2008 / Accepted 23 April 2009
Abstract
Context. A strong signature of a circumstellar disc around a highmass protostar has been inferred from high resolution methanol maser observations in NGC 7538IRS1 N. This interpretation has however been challenged, with a bipolar outflow proposed as an alternative explanation.
Aims. We compare the two proposed scenarios for best consistency with the observations.
Methods. Using a newly developed formalism, we model the optical depth of the maser emission at each observed point in the map and LOS velocity for the two scenarios.
Results. We find that if the emission is symmetric around a central peak in both space and LOS velocity, then it has to arise from an edgeon disc with sufficiently fast differential rotation. Disc models successfully fit independent measurement points in positionvelocity space with 4 free parameters to an overall accuracy of 34%. Solutions for Keplerian rotation require a central mass of at least 4 .
Close to bestfitting models are obtained if Keplerian motion is assumed around a central mass equaling 30 ,
as inferred from other observations. In contrast, we find that classical bipolar outflow models cannot fit the data, although it could be applicable in other sources.
Conclusions. Our results strongly favour the differentially rotating disc hypothesis to describe the main feature of the 12.2 (and 6.7) GHz methanol maser emission in NGC 7538 IRS1 N. Furthermore, for Keplerian rotation around a 30
protostar, we predict the position and velocity at which tangentially amplified masers should be detected in high dynamic range observations. Also, our model predicts the amplitude of the proper motion of some of the maser features in our data. Confirmation of a large central mass would strongly support the idea that even the highest mass stars (>20 )
form via accretion discs, similarly to lowmass stars. Finally we note that our new formalism can readily be used to distinguish between discs and outflows for thermal emitting line sources as well as masers.
Key words: stars: formation  radio lines: ISM  masers  stars: circumstellar matter  stars: abundances
1 Introduction
1.1 Highmass star formation through accretion discs
While it is generally accepted that lowmass stars (M < 8 ) form via accretion discs, the situation is less clear for more massive stars (see e.g. Zinnecker & Yorke 2007, for a review). Highmass protostars produce large radiation pressures which are able to reverse accretion flows and prevent the growth in mass of the central object. The simplest models suggest that stars more massive than 8 cannot form by standard accretion. To explain the observation of more massive stars, alternative mechanisms have to be invoked, such as competitive accretion, mergers of lower mass stars (see e.g. Bonnell et al. 2004) or accretion through a ``trapped'' hypercompact H IIregion (see e.g. Keto 2003, and references therein). It is therefore of great astrophysical interest to determine observationally whether highmass protostars have accretion discs or not.
In recent years a number of claims of compact discs (<1000 AU in radius) surrounding highmass stars have been made (see the review of Cesaroni et al. 2007). These are based both on radio and millimetre continuum emission (see e.g. G192.163.82, Shepherd et al. 2001; GL490 Schreyer et al. 2006; and OrionI Reid et al. 2007) as well as millimetre spectral line emission (IRAS 20216+4104 Cesaroni et al. 1997,2005; and Cep A, Patel et al. 2005). Maser emission observations include the equatorial disc outflow in Orion KL (Greenhill et al. 2004,1998) and OH masers in IRAS 20216+4104 (Edris et al. 2005). Despite these results, no disc has yet been observed around a protostar with a mass larger than 20 (Cesaroni et al. 2007). Detailed numerical simulations invoking nonspherical accretion via dense discs are now able to model the formation of stars up to this mass, but it is still an open question as to whether disc accretion is a valid formation mode for more massive stars (Zinnecker & Yorke 2007).
1.2 Methanol masers, discs and outflows
Methanol maser emission has been recognised as one of the best tracers of highmass star forming regions and has been extensively searched for in the last two decades across the Milky Way. The most recent of these searches, after having covered 60% of the galactic plane, has yielded some 800 sources (see e.g. Green et al. 2009). Of these, 520 objects were already reported in the literature and show a variety of characteristics (Pestalozzi et al. 2005). All the known methanol masers are associated with highmass star forming regions, most of them in very early stages of evolution, prior to the creation of an Ultra Compact (UC) H IIregion (see e.g. Walsh et al. 1998; Ellingsen 1996). A considerable number of these methanol masers are rich in spectral features that very often align both in space and in lineofsight (LOS) velocity (Norris et al. 1998). These lines were interpreted as coming from edgeon rotating discs since in such systems, the conditions for building up long, velocity coherent amplification paths are naturally met.
Although the disc interpretation was the first model suggested for lines of methanol masers, other models have subsequently been proposed. It has been argued for instance that they arise in jets or outflows. In W3(OH) for instance, Moscadelli et al. (2002) model the methanol maser line as tracing the surface of a bipolar cone spiraling with a constant velocity around the cone axis. Another model proposes that the linear 6.7 GHz methanol masers are produced by planar shocks propagating nearly perpendicular to the line of sight in a rotating cloud (Dodson et al. 2004). Attempts have been made to distinguish between disc and outflow models by comparing the directions of the lines of masers and the H_{2} outflow axis in sources where both are detected (De Buizer 2003). Among the sources that produce a clear signature, in 6 the two axes are parallel, indicating an outflow origin for the masers, and in 2 they are perpendicular, in agreement with the disc interpretation. Six additional sources are classified as ``likely'', though not definitively, parallel. While these results strengthen the outflow hypothesis, they indicate a disc origin in some cases.
1.3 Previous observations and interpretation of NGC 7538 IRS1 N methanol maser emission
Methanol maser emission in the highmass star forming region NGC 7538 at 6.7 and 12.2 GHz was discovered by Menten (1991) and Batrla et al. (1987) respectively. Both masers were found to be coincident with the IR source NGC 7538 IRS1 N. Subsequently, a number of high spatial resolution observations have been performed using the European VLBI Network (EVN^{}) and the Very Long Baseline Array (VLBA^{}). From these observations, very accurate maps of the spatial and dynamical structure of the main spectral maser feature have been made. Recent MERLIN observations at lower spatial resolution have revealed that two other weak maser spectral features are associated with objects 1.5 arcmin to the south of the main feature (close to IR sources IRS 9 and NGC 7538S). These results suggest that several massive stars are forming in the NGC 7538region (Pestalozzi et al. 2006).
The main maser emission feature in NGC 7538 IRS 1 (feature ``A'', see Pestalozzi et al. 2006 for the nomenclature) is 2 km s^{1} wide and has a peak flux of 350 Jy at 6.7 and 120 Jy at 12.2 GHz. It is seen projected on a UC H IIregion with a brightness temperature of 10 00015 000 K (Gaume et al. 1995; Campbell 1984). The source powering the UC H IIregion appears to be an O7 star (Akabane & Kuno 2005). Minier et al. (1998) recognised, in the very linear shape, both in space and LOS velocity, of this spectral feature at both 6.7 and 12.2 GHz, the first hint for a rotating disc seen edgeon. Assuming the maser spots to lie on the outer radius of a thin disc around a central star of 10 , the authors suggested the radius at which the masers occurred to be several hundred AU. The same data was modelled in detail by Pestalozzi et al. (2004) assuming for the first time masing methanol over a range of radii (3501000 AU), distributed in a disc in Keplerian rotation around a 30 star^{}.
Challenging the above disc model, midIR observations have been interpreted as showing dust emission at different temperatures, tracing cavities excavated by an outflow (De Buizer & Minier 2005). Since these cavities are oriented roughly parallel to the line of methanol masers, the maser emission has been interpreted as arising within the cavities of a collimated jet/outflow rather than a disc. A controversy then exists: the methanol maser emission toward NGC 7538 IRS 1 seems to be equally validly interpreted as arising from a rotating disc seen edgeon or from a bipolar outflow, depending on the data considered for the analysis.
In even more recent work, the outflow hypothesis for the methanol maser has been weakened. Kraus et al. (2006) find that the large scale CO outflow originating in NGC 7538 IRS 1 is probably precessing. Detailed precession models can be made consistent with a disc origin for methanol masers and the large scale outflows with cavities traced by CO and IR emission.
In this paper, we further develop the formalism used in Pestalozzi et al. (2004) with the aim to test the competing disc/outflow explanations of the methanol masers in this object. We show that the bipolar outflow geometry is not able to reproduce the emission shown in the data. The only way to obtain a symmetric emission around a maximum both in space and LOS velocity is to model the maser as emitted by a differentially rotating disc seen edgeon.
2 Methanol maser data in NGC 7538 IRS1 N
The data modelled in this paper are shown in Figs. 1 and 2, where the maser optical depth is presented. This is obtained from the flux density maps using Eq. (1) below. The 6.7 GHz data were taken in February 2001 using the EVN, the 12.2 GHz in March 2005 using the VLBA. The resolution in the data is indicated by the size of the beam in the lower right corner of the figures and is about 5 and 2 mas for the 6.7 and 12.2 GHz data respectively. The velocity resolution is 0.088 and 0.048 km s^{1} for the 6.7 and 12.2 GHz data respectively. Both the velocity integrated maps (Fig. 1) and the positionvelocity diagrams (Fig. 2) contain a continuous smooth central feature extending over some 40 mas in RA. At the resolution of our VLBI observations, this translates into 8 to 20 independent measurements across the central maser feature at 6.7 and 12.2 GHz, respectively. Such a large amount of data allows our models to be highly constrained. At displacements greater than 20 mas, some weak disconnected regions of emission are visible. These outliers were already seen in the early 6.7 GHz interferometry data but were detected at 12.2 GHz only in the most recent observations.
Figure 1: Map of the velocity integrated optical depth of the maser emission in NGC 7538, obtained from Eq. (1). From top to bottom: 6.7 GHz data, 12.2 GHz data, 12.2 GHz data convolved at the resolution of the 6.7 GHz data. The peak value of is 18 and 16 for the 6.7 and 12.2 GHz respectively, in the outliers it reaches 11 and 14. The beam size is indicated in the lower right corner of each panel. Contours are at 70, 80, 90, 95, 99% of the peak in all panels. Note that the lowest contour level is set by the dynamic range limit of our observations. 

Open with DEXTER 
Figure 2: LOS velocity versus spatial displacement (diagram) of of the main maser emission feature in NGC 7538. The top panel shows 6.7 GHz data, the bottom panel the 12.2 GHz data. The peak value of is 18 and 16 for the 6.7 and 12.2 GHz respectively, in the outlier features at 40 mas and 80 mas and it reaches 11 and 14. The beam size is indicated in the lower right corner of each panel. Contours are at 70, 80, 90, 95, 99% of the peak in both panels. As for Fig. 1, the lowest contour level is set by the dynamic range limit of our observations. 

Open with DEXTER 
From phasereference astrometry, it is seen that the 6.7 and 12.2 GHz maser emissions are perfectly cospatial within 2 mas (see Pestalozzi et al. 2004). The cospatiality of the two emissions is also supported by the strong similarity between the maps in the first and third panels of Fig. 1. The former shows the 6.7 GHz emission and the latter the 12.2 GHz data convolved to the resolution of the 6.7 GHz map. This similarity is remarkable given that that the two maser lines are emitted by two different isomers of the methanol molecule, and hence are expected to trace different regions.
The superposition of the maser emission on bright centimetre continuum emission, together with the cospatiality of maser emission at two frequencies, leads to the assumption that the maser is the result of the amplification of a background continuum. In that case, the brightness at displacement (x,y) from the centre of each frequency plane with LOS velocity
is expressed as^{}:
where I_{B} is the background continuum and the (negative) maser optical depth. This expression describes foreground amplification by any maser, irrespective of its degree of saturation (e.g., Elitzur 1992). The optical depth is therefore defined at every point in the cube and, knowing the brightness temperature of the background, it is a measured quantity, . Considering K and the maser peak brightness per unit area I(0,0,0) at the centre we have that and 18 for the 12.2 and 6.7 GHz respectively. In the outliers we have and 11. This drop in optical depth between centre peak and outliers (30%) corresponds to a dynamic range in brightness of 100:1. Note that the shape of the contours in Figs. 1 and 2 does not depend on the value of the background radiation I_{B}. The background brightness is only needed to obtain the value of the maser optical depth at the peak.
There are two further important characteristics of these data. The first is that the bulk of the maser emission ( mas) is symmetric about a maximum of emission both in space (Fig. 1) and in LOS velocity (Fig. 2). The second is that in the diagrams, the emission is seen to bend away from the overall gradient at a displacement close to 1015 mas. These two facts together will be shown to be naturally explained by emission from a differentially rotating disc seen edgeon, whilst they are virtually impossible to reproduce with an outflow geometry. In the following discussion, we will concentrate on the central part of the 12.2 GHz data ( mas), as these show the highest spatial resolution.
3 The formalism: fundamental expression
Assuming that the observed maser emission is generated from amplified continuum radiation, then the observed radiation is completely characterised by the maser amplification at each position (see Eq. (1)). Introduce Cartesian axes with z toward the observer and denote by
the (negative of the) linecenter maser absorption coefficient at an arbitrary position in the source and by
the local line profile. The maser optical depth is then
,
where
is the LOS component of the local bulk velocity at (x,y,z). Introduce
,
the linecenter optical depth at some conveniently chosen fiducial point
in the plane of the sky. Then the maser optical depth can be written
,
where
and where contains all the effects of temperature, density and inversion gradients within the masing cloud. The dimensionless function T, which is normalised such that , contains all the spatial and velocity information in the maser data cube. Because the data in Fig. 1 show a remarkably high aspect ratio, we assume that the masers arise from a flat structure. We therefore rotated the cubes by and integrated them along the height axis to obtain displacementLOS velocity diagrams (diagram, see Fig. 2). Hence, in what follows we concentrate on modelling the function . In the cases of interest here, the spatial amplification profile is primarily a function of , the radial distance from the center. The size scale of the system can be specified by some characteristic length r_{1}, the radial coordinate at some convenient point in the source; that is, the profile is a mathematical function of the scaled dimensionless variable . Taking for a Gaussian profile with a width finally yields
This general expression describes the optical depth of any maser, whether saturated or not. Saturation implies that the population inversion is affected by the propagating maser radiation, becoming one of the factors that shape the radial profile . Here this is a moot issue because we limit ourselves to a parametric description of the profile without attempting to uncover its physical foundation. Saturation also causes broadening of the frequency profile of the maser absorption coefficient. This typically requires 15 across the entire maser, and the profile broadening is then proportional to length in excess of this threshold (Elitzur 1992). With an optical depth of 1618, saturation can be expected to have some effect here, but only at the very outer segments of the source where a decrease in temperature could have an opposite, offsetting effect. The width is therefore taken as constant, an approximation that is not expected to affect our fundamental conclusions.
While describes the spatial distribution of the maser amplification, contains the details of the dynamics of the system. Maser emission requires velocity coherence along the LOS  in order to participate in the maser action at a certain frequency v, the difference in LOS velocity of two points along a given direction cannot exceed the thermal width . That is, once the geometry and dynamics of the system are defined, the variation of along any LOS determines the coherence length that controls the maser optical depth for that x. Figure 2 shows that the maser emission peaks at the center of the diagram. The challenge is to identify distributions and that reproduce this property.
In the next section we consider the application of the above formalism to edgeon discs. In Sect. 5 we fit the disc model to the NGC 7538 data and derive the best fitting disc parameters. In Sect. 6 we apply the formalism for the case of outflow geometry, and show why this model cannot fit the data shown in Figs. 1 and 2.
4 Edgeon disc models
In an edgeon rotating disc, the LOS velocity
and rotational velocity
of a point at distance
from the center and azimuthal angle
(see Fig. 3) are related according to:
where is the local angular velocity. Significantly, the LOS velocity is determined by the angular velocity, not the rotational velocity itself. Denote by the rotational velocity at r_{1} and by the corresponding angular velocity. A fast rotating large disc and slowly rotating compact one will have widely different values of and r_{1} yet they will produce the same LOS velocity field if they have the same and radial variation of .
In this section we show from some simple considerations how the morphology of the data in Figs. 1 and 2 constrains the disc properties. In particular, we demonstrate how the simultaneous appearance of a maximum of emission in the centre of both the map and the diagram can naturally be reproduced only by sufficiently fast differential rotation.
Figure 3: Geometry for the disc model, with the rotation velocity. The observer is in the plane. 

Open with DEXTER 
4.1 The norotation case
In a nonrotating disc, . Therefore the maser molecules are velocitycoherent along the whole source for every LOS, and the maser optical depth obeys , where is the linecenter optical depth at position x. The positional variation of the amplification is controlled solely by the length of the amplifying column because the velocity profile is the same at each position. The diagram of the amplification contours in this case is shown in the top panel of Fig. 4. The contours show two distinct peaks, displaced symmetrically from the center at the two inner tangents, while the center is a local minimum. In any nonrotating disc with an inner radius, the central LOS has a lower opacity than the LOS that is tangent to the inner radius, and this holds no matter what the distribution is. This happens because, just like the LOS to the centre, the LOS tangent to the inner radius samples all , but each interval of radius corresponds to a longer path along the LOS than for the path toward the centre. Nonrotating discs will therefore never produce the central peak of the diagram evident in the data shown in Fig. 2.
4.2 Solidbody rotation
In solidbody rotation, the angular velocity does not depend on radius, therefore is constant along each LOS (Eq. (4)). Similar to the norotation case, the material remains fully velocity coherent along each LOS, only the velocity profile is now centered on instead of v = 0 so that the optical depth obeys . As a result, the structure of amplification contours in the diagram remains the same, only rotated by the angle tan (and slightly stretched to maintain the peak positions at the two tangents to the inner radius), as shown in the two bottom panels of Fig. 4. The center, x = 0, remains a local minimum, in analogy with the nonrotating disc. Solidbody rotation, too, will never produce the central peak observed in the diagram.
Figure 4: diagram of maser optical depth (Eq. (2)) for edgeon discs with constant amplification ( ) between inner radius r_{1} and outer radius 3r_{1}. Contour values are normalised to the value at the diagram center and are spaced by 0.1, with the first contour at 0.6. Black contours indicate values greater than the central value. Top: a nonrotating disc. The emission maxima are along the tangents to the inner edges , the longest paths through the disc. Mid and Bottom: same configuration as top panel only the disc rotates as a solid body with the angular velocity , as indicated. Notice that all the maps are qualitatively identical except for the tilt, i.e. norotation and solid body rotation are two aspects of the same case. 

Open with DEXTER 
4.3 Differential rotation
The data shown in Fig. 2 cannot be reproduced by either a nonrotating disc or one rotating as a solidbody. Both scenarios fail because the disc geometry implies that the LOS through the disc center is a local minimum of pathlength. The introduction of differential rotation changes the situation fundamentally. Although the amplification path is geometrically short, it maintains its full velocity coherence across the entire disc, irrespective of the rotation law, because the motion is perpendicular to the LOS and hence . In contrast, along any other path, the rotation velocity has a finite LOS component, with the consequence that segments of the path can move out of velocity coherence when varies with radius. In this situation the optical depth at any LOS will no longer depend exclusively on the geometrical pathlength across the disc but will be limited by velocity coherence. This situation is represented in the left panel of Fig. 5, where coloured areas illustrate how the length of coherence paths changes with displacement in Keplerian differential rotation.
Differential rotation alone will not produce a maximum in the centre of the diagram if it is not sufficiently different from the two previously discussed cases. In general, differential rotation that is either too slow or too weakly dependent on will produce a twinpeak diagram, similar to the norotation or solidbody cases. For a given differential rotation law, a transition from a twopeaked to a single, central peaked diagram is expected only when the rotation velocity increases above a certain threshold. This is illustrated in Fig. 6, which shows the evolution of the diagram with increasing rotation velocity for a Keplerian disc with constant . In this particular case, the threshold rotation velocity in the inner radius lies roughly between 4 and 6 .
Figure 5: Contours of constant LOS velocity () and the velocity fields that produce them. The contours are spaced by 0.1 of the highest value and are also colourcoded to indicate the Doppler velocity. The velocity vector field is indicated by arrows. Left: Keplerian rotating disc with masing material confined within inner and outer radii with ratio 1:3. Right: a bipolar flow that produces the same diagram as the disc geometry (see Eq. (16) below), which has outflowing gas on the right side of the central LOS and inflowing gas on the left side, as indicated by the arrows. The radial velocity was multiplied by a factor of 2.5 for clarity. The masing material has been confined within a cone of opening angle 20. 

Open with DEXTER 
Figure 6: Same as Fig. 4, only the disc is in Keplerian rotation, with angular velocity on the inner radius. Notice that a maximum at the centre is obtained when the rotation is fast enough ( bottom panel). 

Open with DEXTER 
4.4 The p, diagram ``spine''
In each of the displayed diagrams, the locus of points with strongest amplification at all LOS stands out as a clearly visible feature. We term this feature the spine of the diagram. In the case of solidbody rotation, it is the straight line , reverting to for the nonrotating disc. Differential rotation introduces a curvature in the spine, evident in the varying orientation of highestamplification contours in Fig. 6. A similar curvature is clearly visible also in the data in Fig. 2, shown by the change in velocity gradient at peak amplification between the inner and outer parts of the central maser feature. Spine analysis provides useful insight into the effect of the rotation on maser amplification. In particular, it provides an answer to the question: How fast must the differential rotation become for the central emission to turn from a local minimum to a local maximum?
The condition
determines the relation between v and x on the spine, defining the spine curve
in the diagram. Maser amplification along this curve is the spine amplification
.
For a general spine analysis
we express the rotation velocity as a power law, so that
Here is the rotation velocity at r_{1} and is the angular velocity there. Keplerian rotation corresponds to while solidbody rotation is . The spine condition then yields
The integral on the right differs from only by the additional term in the integrand (cf. Eq. (3)); when = 1 the two integrals are identical, and we recover the solidbody result ; the spine in this case is a linear structure tilted by angle . Consider now the centre of the diagram. At that point we always have T(0,0) = 1, while the integral on the right of Eq. (6) is , where we introduced the notation for the moments of the profile. Because angular velocity along the spine is , we have
where is the slope of the spine at the origin of the diagram. This important result relates the measured quantity to the free parameters of the model. It is used extensively in our data analysis, as explained in Sect. 5.
Equation (6) defines the spine only in an implicit form because
enters on both sides. Straightforward series expansion in powers of
produces in 2nd order
where
This shows that the amplification along the spine will increase (decrease) with x when b is positive (negative), producing a local minimum (maximum) at the center of the diagram. A nonrotating disc has , therefore and the origin is always a local minimum in this case. Similarly, solid body rotation has so that again , identical to the norotation case. A local maximum for at the origin requires b < 0, yielding
This is the condition for a local peak at the origin of the diagram. It can be fulfilled only when both , i.e., differential rotation, and exceeds by some dependent prescribed factor, i.e., sufficiently fast rotation.
In addition to its potential to produce a central peak in the diagram, another fundamental property of differential rotation is the spine curvature. Similar to Eq. (8), a 2nd order expansion for the angular velocity along the spine yields
where
and where is given in Eq. (7). Solidbody rotation ( ) yields c = 0 as it should (the angular velocity is constant on the spine), but in the case of differential rotation . Then the angular velocity on the spine deviates markedly from its value at the origin, and the spine curvature becomes noticeable, when the quadratic term becomes significant and the velocity coherence prescription begins to play a role in the determination of the amplification path length. The presence of a clear curvature in our data provides further strong evidence for the differentially rotating disc hypothesis.
5 Fit to the data
In this section, we describe the results of the multidimensional fitting of disc models to our data. We first identify the natural scales in our problem for space and velocity with r_{1} and , respectively. In addition to these scales, the problem contains the two unknown functions and . Depending on their parametrisation, these functions will add to the total number of free parameters to be fitted.
Two quantities can be derived directly from the data without detailed analysis of the full diagram. The first one is the velocity scale = 0.432 km s^{1}, obtained from Gaussian fitting to the data at . The other quantity is the spine slope at the origin of the diagram, = 0.064 km s^{1} mas^{1}, where D (=2.7 kpc; see Blitz et al. 1982; Moscadelli et al. 2009) is the distance to the NGC 7538 region. Through Eq. (7), provides a constraint on the free parameters we use in the detailed modeling.
5.1 The function
From Eqs. (8) and (10) it is evident that both the shape of the spine and the value of the normalised optical depth T on the spine depend on the dynamics (represented by the value of )
and on low order moments of the function
.
It follows that the exact form of
will be unimportant in fitting the data; only its overall properties such as width and kurtosis (lobsidedness) likely will matter. This fact, together with the need
to minimise the number of fitted parameters, suggest that a twoparameter function for
is a
practical choice. Pestalozzi et al. (2004) used a powerlaw function between inner and outer radii, obtaining good results. Although this function has only two parameters (ratio of radii and powerlaw exponent) it also has unphysical sharp edges at the inner and outer radii, which we wish to avoid in this work. After some experimentation, we settled on a modified lognormal function of the form
where A is determined from the normalisation . To avoid blowup when , must be positive and the argument of the exponential must remain negative, therefore p must be an even integer. Further, without loss of generality, we choose for r_{1} the location of the maximum of , that is, the maximum of is at . This yields the following relation between p, and :
This implies that the argument of the exponential term in Eq. (11) becomes . The profile is then completely defined by the two parameters p and . Figure 7 shows three examples for different combinations of p and . While the function is continuous, there usually exist steeply rising ``shoulders'' of the profile. Large values of p produce profiles with very steep shoulders, while small values give smoothly decreasing functions. At a given p, the profile is very peaked around , similar to a thin ring, when (small ; see Eq. (12)) and is very extended when .
Figure 7: Radial variation of the linecenter absorption coefficient (Eq. (11)) for values of p and of three of the bestfit solutions in Table 1 with (solid), 0 (dashed) and 0.3 (dasheddotted); the legend lists also the corresponding width in each case (Eq. (12)). All profiles are normalised to unit area and have their maximum at = 1 (r = r_{1}), though r_{1} is different in each case (see Table 1). 

Open with DEXTER 
5.2 The function and dynamical mass
From Eq. (5), the function for the LOS velocity
adds two free parameters to the problem,
and ,
the angular velocity at r_{1} and the index of the radial variation of
,
respectively. The latter is set by the centripetal acceleration produced by the gravitational force of the mass distribution. Spherical distributions give a reasonable estimate of this dynamical mass even for the disc geometry because material inside and outside any r has opposite effects in the two geometries that tend to cancel each other. The interior mass would produce a smaller centripetal acceleration when arranged as a sphere because some of the gravitational force goes into components perpendicular to the rotation plane. On the other hand, the mass outside r would decrease the centripetal force by its outward pull when arranged in a disc but will have no effect when distributed in a spherical shell. The net effect is that at a given radius, a spherical mass produces rotational velocity similar to that of flat disc models having the same total mass (see e.g. Fig. 1 in Toomre 1963). Considering a spherical mass distribution and using Eqs. (5) and (7), the dynamical mass inside radius is
In Keplerian rotation, M does not depend on , i.e., the cause of the dynamics can be considered a central point mass. Formally, the mass diverges when for any but this divergence is meaningless in practice because the data probes only finite radii (although the mathematical expression for has no formal bounds). Differentiating the above expression for shows that the mass density is proportional to . Since the density must be positive, must be larger than ; Keplerian rotation provides the lower limit on .
Table 1: Bestfit models of edgeon differentially rotating discs. The parameters = 0.432 km s^{1} and = 0.064 km s^{1} mas^{1} are determined directly from the data and are common to all models. The distance to the source D = 2.7 kpc.
Figure 8: Data ( left) and modeling ( middle) of the bulk of the 12.2 GHz maser emission (within 20 mas of the feature centre). Contours are at 70, 80, 90, 95 and 99% of the peak optical depth (). Red dots show the amplification spine (see Sect. 4.4). Right: map of the model relative error at every position in the diagram, corresponding to the individual terms in the sum in Eq. (14). Note the color bar on the right for the error magnitude in percent. The error averaged over the entire map is 3.07%. The model corresponds to an edgeon Keplerian disc around a 31 star (see Table 1 for details). 

Open with DEXTER 
5.3 Free parameters and fitting procedure
The masing absorption coefficient profile is described by the two free parameters p and , and the LOS velocity by the two additional ones and . Directly from the data, the spine tilt at the origin of the diagram is = 0.064 km s^{1} mas^{1}, providing a constraint that reduces by one the number of free parameters: given , p and , the angular velocity is set by Eq. (7). Two additional quantities are introduced by the fundamental expression (Eq. (3)) for the amplification map: the independent velocity scale and spatial scale r_{1}. The first is determined directly from the data, = 0.432 km s^{1}, while the second adds another free parameter, bringing their total number to four: , r_{1}, p and . Here we carry out for the first time a fully unbiased search over these four parameters. This represents a significant improvement over Pestalozzi et al. (2004), where Keplerian rotation was assumed, a central mass of 30 was taken from other observations and the best fit was found by optimising the two parameters defining , which was chosen to be a single powerlaw between sharp cutoffs.
For a given set of the four free parameters, we determine the quality of the model by comparing the optical depth determined from Eq. (3) with the data at every point in the central part of the diagram (within 20 mas of the origin) where a measurement above the lowest detected opacity contour exists. Since spectral resolution is about 2 pixels in the velocity direction (corresponding to 0.1 km s^{1}) and the FWHM beam size is 5 pixels or about 2 mas, there are 100 data measurements within the lowest contour.
As a quality estimator of each model we use the relative difference between model and data averaged over the detected part of the source:
where N is a normalisation parameter equalling the number of pixels in the fitted diagram; hence X is a measure of the mean fractional error of the model opacity compared to that observed in the data. We chose this procedure after some experimentation. Although the measurement errors generally increase toward the edge of the central feature, all points were given equal weight in this averaging procedure because of the difficulty in estimating the noise in data determined from the logarithm of the intensity map. For this reason we also did not attempt a formal reduced analysis. The choice of the absolute value of the difference between model and data is in general preferable to the square of the difference, which tends to emphasise noisy points (see Press et al. 2007).
While the quality estimator X includes only the central feature ( 20 mas), larger displacements from the origin play an equally important role in constraining disc models. As discussed in Sect. 4.3, and displayed in Fig. 6, differential rotation that is not sufficiently fast will produce tangential amplification features at displacements 5 times larger than the central feature size, i.e., 100 mas. Our data show that such features do not exist to within the sensitivity of the observations, whose dynamic range is 100:1. This implies that maser optical depth in the tangential features must be less than 70% of , its value at the origin of the diagram. Acceptable models must meet this criterion, in addition to producing a suitably small X.
5.4 Results
Table 1 presents a series of best fits. Starting from the different values of listed in the first column, we found the combination of r_{1}, p and that minimised X (Eq. (14)) and produced acceptable models. The quality of all tabulated fits is high. Figure 8 shows the detailed comparison of model and data for the tabulated Keplerian model with 31 . This model produces a formal error of X = 3.07%, and the figure's right panel shows that it describes the data perfectly at the center of the diagram and to better than 10% for the large majority of calculated points within 20 mas of the origin. As is evident from the table, all rotation curves in the range are capable of producing a highquality fit with an average error of only 34%. The corresponding amplification profiles vary from broad distributions at large to peaked shapes as decreases toward Keplerian rotation. This is evident from the few examples plotted in Fig. 7 as well as the tabulated values of the radii and where decreases to 5% of the peak value. Moreover, the few entries listed for Keplerian rotation show that an unlimited number of models are producing essentially the same quality fits in this case. These different entries have nearly the same and , differing from each other only in the size scale r_{1}, which can be increased without bound. The same behaviour has been found for models with +0.4. For those models, Table 1 contains only the solutions with minimal r_{1}.
Figure 9: Full diagram of the best fitting Keplerian model with a 31 central protostar. Contours are plotted at 1, 3, 10, 30, 50, 60, 70, 80, 90, 95 and 99% of the peak optical depth = 16. The red contours outline the region visible in our data (>70% of ); this is the region shown in Fig. 8. Notice the small fraction of the full extent of the masing disc covered by the available data. The green dots indicate the positions of the outliers. 

Open with DEXTER 
The reason for the degeneracy among so many highly successful models is that the data sample only a small portion of the full diagram generated by the disc. As an example, Fig. 9 shows the full diagram for the 31 Keplerian model, drawn to a level of 1% of the peak value. While the map extends to mas, the data cover only the inner mas (red contours in the figure). This has a profound impact on the sensitivity of the model fitting. In fact, as is evident from the fundamental expression Eq. (3), the dependence on the scale r_{1} comes only from the integration lower limit and the curvature term in the denominator, and in both of them it enters as . While the data extend only to 20 mas, r_{1} is at least 110 mas for Keplerian rotation (and 220 mas for the model displayed in Fig. 9), so and is even smaller still at all measured points. Therefore, the entire data cube could be adequately described by a series expansion limited to the first few powers of , similar to the spine analysis presented in Sect. 4.3. The expansion coefficients involve only moments of the amplification profile , which can be easily adjusted with a suitable choice of its parameters; this is evident from the expressions for the spine (Eqs. (8) and (10)), the data most prominent feature. Furthermore, because is so small, all acceptable models are practically indistinguishable from the limit , which can be formally obtained by taking while holding constant. This formal limit implies , i.e., an infinite mass (see Eq. (13)). As can be seen from Table 1, every Keplerian entry with 12 and above has the same profile (same p and ) and the same , differing only in the scale size r_{1}; that is, the infinite mass limit is already reached at 12 in the case of Keplerian rotation. Similar behavior is found for the other tabulated values with . Although r_{1} decreases as increases, the profile becomes much more extended and the bulk of the integration originates from radii much larger than the observed displacements, yielding a similar behavior to the Keplerian case, which is concentrated around .
As this discussion shows, analysis of the available data is quite insensitive to the free parameter r_{1}, which is equivalent to (= ). The only meaningful constraint on the latter is a lower limit on to ensure that the origin of the diagram is a local maximum (Eq. (9)) and that the tangential emission is sufficiently suppressed. Removal of the degeneracy among the bestfit models requires interferometry with a higher dynamic range, capable of detecting the tangential features. Current data do show emission at displacements between 30 mas and 80 mas (see Figs. 1, 2 and 9), and we term these features ``outliers''. We do not think these outliers are the tangential features predicted by the models because their displacements are much smaller than the disc outer radius and they are not located quite on the extended spine, although close to it. Instead, the outliers are most likely regions of chance enhanced LOS coherence close to the spine, as suggested by their location and by the fact that they occur on only one side of the disc. Because of the maser exponential amplification, an enhancement of only 10% in the local value of would suffice to produce these outliers. Also, small deviations from axial symmetry (outliers only on one side) could be easily generated by, e.g., spiral density waves or a small warp that brings the East side of the maser disc slightly closer to edgeon than the West side, producing better LOS alignment of maser molecules (and hence outliers only on one side of the disc). Such a warp might be consistent with the disc precession in NGC 7538 IRS 1, inferred by Kraus et al. (2006).
5.5 Dynamical considerations
Edgeon discs in differential rotation fully capture the structure of maser amplification in the limited region of the diagram covered by our data. Most disc parameters are largely irrelevant as long as the rotation is sufficiently fast. This degeneracy makes it impossible to determine the disc properties purely from a bestfit analysis of the maser data. We must invoke additional considerations to narrow down the range of acceptable disc parameters.
For each successful model, Table 1 lists the dynamical masses and , calculated from Eq. (13), contained within and , the respective radii where decreases to 5% of its peak value on each side of the peak at r_{1} (). These radii effectively mark the inner and outer boundaries of the disc maser region. As noted above, most sets of parameters allow r_{1}, and the corresponding , to increase indefinitely; both and are left intact, though. As is evident from Eq. (13), both and then increase without bound but their ratio remains constant. The tabulated show that all models with are essentially devoid of any mass interior to the maser region. These models correspond to selfgravitating discs without any central object, either a star or just a central bulge, and thus can be discarded as unlikely to arise in realistic situations. Models with do contain a sizeable central object and thus are more likely to correspond to possible configurations. However, in each case the disc contains a significant fraction of the full mass, an inherently unstable situation: as shown by Adams et al. (1989), such discs are unstable to growth of eccentric distortions arising from small deviations between the positions of the star and the system centre of mass. Such unstable systems are unlikely to harbour the remarkably smooth, regular structure observed in the maser central feature. The problem is avoided only in the Keplerian models, where the mass in the disc is negligible in comparison with that in the central object.
Keplerian models are the only tabulated ones to offer stable physical systems, yielding a lower bound of 4.1 on the central mass; all larger masses produce equally successful fits while smaller masses produce tangential features that conflict with the data. Therefore, we can conclude with reasonable confidence that the maser central feature arises from an edgeon Keplerian disc around a central star heavier than 4 . Kraus et al. (2006) conclude that the most plausible explanation for the jet precession around NGC 7538 IRS1 is a circumbinary disc around a binary pair separated by 7 mas (19 AU), which is much smaller than the minimal value of 110 mas for r_{1}. The maser disc fully encompasses the binary pair, whose total mass therefore must exceed 4 . This lower limit on the central mass is consistent with the 30 inferred from the O7 spectral type of IRS1 (Akabane & Kuno 2005) but the degeneracy of our model fits prevents us from conclusively identifying this star with the center of the maser disc. However, the apparent alignment of the disc with the ``waistline'' of the associated ultracompact H II region strongly suggests that this is indeed the case. The alternative would require a chance coincidence of this star with a foreground disc around a lower mass star, which seems less likely.
6 Bipolar outflow
In the previous two sections, we considered the application of the general opacity formalism of Eq. (3) to the specific case of an edgeon rotating disc. However, as noted in Sect. 1.3, it has been proposed that the linear maser feature in NGC 7538 IRS 1 is instead due to an outflow. This alternative can readily be modelled using the same general expression by modifying the velocity field. Consider a flow with radial velocity
.
The LOS velocity of a point at distance
from the center and azimuthal angle
is now
(see Fig. 10). Other than a different relation between velocity and its LOS component (cf. Eq. (4)), the outflow is handled in exactly the same way to the edgeon disc. Is it possible devise a radial configuration that produces the observed central maximum in the diagram? In this section we first demonstrate that, when measuring only Doppler velocity components, every disc model can be formally considered as if arising from radial motion, but this formal equivalence is unphysical as it requires alternating infall and outflow along adjacent radial rays (Sect. 6.1). We then show that physical outflow models cannot reproduce the observed central peak (Sect. 6.2).
6.1 Transferring disc solutions to radial flow models
Comparison of Eqs. (4) and (15) shows that, given an edgeon disc rotating with velocity
,
we can always define an equivalent radial flow with the velocity
and the two will have the exact same LOS velocity field. As long as we measure only the LOS component of the velocity, it would be impossible to distinguish between the two configurations. The right panel of Fig. 5 shows the radial flow velocity field equivalent to the Keplerian disc shown in the figure's left panel. Since these two very different morphologies will produce the same LOS pattern, no direct observations can differentiate between them without measuring proper motions. However, although they will produce identical maps and diagrams, the physical foundations of the two configurations differ significantly. The Keplerian field arises naturally from rotation around a central star. In contrast, the equivalent radial motion does not arise from any plausible physical scenario, requiring zero velocity along the central LOS, outflow to one side of this direction and inflow to the other. This peculiar morphology is essential to reproduce the central peak and ordered gradients in the observed diagram. Therefore, although radial flow configurations can be mathematically devised to fit the data as well as edgeon rotating discs, they represent highly contrived velocity fields that are not produced by any physical model.
Figure 10: Geometry of a physical outflow model whose axis is inclined at an angle i to the LOS. Distance along the outflow is measured by . The radius marks an injection point for material into the outflow (see text). Beyond this point the outflow has a position dependent bulk flow velocity (dashed arrows) and mass flux conservation is assumed. In general the outflow has a crosssectional radius , which for the conical flow illustrated is proportional to . 

Open with DEXTER 
6.2 A physical outflow scenario
Rather than starting from the disc solution as in Sect. 6.1, we can ab initio consider possible collimated outflows that could fit our data. Consider a narrow central symmetric outflow having crosssectional radius w, density n and velocity , all only depending on the distance from the centre (see Fig. 10). If the inclination of the outflow axis to the LOS i is significantly greater than the outflow opening angle, then every LOS samples only one value of the LOS component of , . From the diagram of such an outflow, it is then possible to directly derive the outflow velocity function . Within this framework, our data in Fig. 2 show constant acceleration for mas and then a reduction in acceleration at larger displacements, indicated by the ``bend'' of the overall gradient. This gives in general that . Mass flux conservation within the outflow implies and the total gas column density across the jet varies as . Assuming that the maser optical depth is proportional to this column density, even without having any constraints on the function w we can say that will decline at least as . This dependence becomes even stronger if we further assume that the jet cross section w is increasing with (e.g. conical outflow). The relationship between and position along the outflow cannot obviously be extended to , and in any case, any outflow must have a reservoir from which outflowing material is drawn. We can simulate such a reservoir with an injection point , indicated by the thick lines in Fig. 10. If the reservoir is an accretion disc, would correspond to half of the disc thickness. At , the mass flow rate is conserved, inside this is no longer the case due to complex dynamics. This fact introduces a discontinuity in the outflow geometry that is not present when modelling an edgeon disc. This discontinuity is inherent in the nature of the bipolar outflow geometry and could therefore be detected in our observations, provided observations are made at a high enough spatial resolution. Our observations suggest that if such a discontinuity is present in our data, it has to be on a scale that is smaller than half of our resolving beam width, i.e. 1 mas.
From our data, we can express a strict constraint on the orientation of the outflow axis. This comes from making the (reasonable) assumption that the outflow velocity at the injection radius is approximately equal to the internal velocity dispersion . The inclination angle is then obtained from the ratio of the LOS velocity at a certain and . As noted above, the injection point lies at less than a beam FWHM projected distance from the centre (i.e. 1 mas) and therefore has a LOS velocity of 0.064 km s^{1} (see estimate of the overall gradient in the data). Using km s^{1} we have that , i.e. the outflow is almost in the plane of the sky.
Figure 11 shows an example diagram of a bipolar outflow illustrating the effect of a varying inclination i to the LOS. The dependence makes decline very quickly with displacement, in strong conflict with our data, which show a smooth structure over about 10 beams. We denote by the ratio in maser optical depth between one beam out, i.e. , and ten beams out, i.e. mas. In a constant crosssection outflow, from the ratio of the corresponding LOS velocities. An increasing crosssection would make the ratio larger still  a conical outflow with its apex on the central star would give . In contrast, the observations show that since only decreases from 15.5 at one beam FWHM from the peak to 11 at ten beams away from the centre. It is impossible to reconcile such extended maser emission with the rapid density decline in outflows.
Figure 11: diagram of maser optical depth for the biconical outflow sketched in Fig. 10. Each panel shows a different inclination i of the outflow axis to the LOS. Contours are drawn at 70, 80, 90, 95 and 99% of the peak optical depth. The outflow has a constantly accelerating velocity, starting from at the injection radius . Maser opacity beyond is proportional to the total gas column density. The opacity and dynamics inside the injection radius are illdefined: in this example the opacity there is taken as zero; quiescent maser emission from material within would add in each panel a horizontal strip about the displacement axis, similar to the top panel of Fig. 4. 

Open with DEXTER 
An important assumption in our model is that maser opacity scales with total gas column density. While the factors that affect maser opacity can be very complex, we argue that the change in total gas column density in an outflow is so rapid (at least ) that this must strongly affect the maser opacity. The only alternative to mitigate this effect is a finetuning of the maser opacity per total gas column density to almost exactly compensate for the dependence of the total gas column density. We consider this highly unlikely.
A final possibility to consider is that a bipolar outflow could produce maser emission from interaction with ambient clumps of the ISM. In this case, the fall off in gas column density within the outflow itself is not directly relevant, the important quantity is the gas column density of the ambient ISM. Such models will produce a string of maser features, symmetrically displaced around a central minimum in the diagram. Although numerous sources do display such structures and thus could be explained with bipolar outflows, these models cannot explain the smooth central peak observed in NGC 7538 IRS1 N.
7 Discussion and conclusions
7.1 Discussion
The dominant methanol maser feature in NGC 7538 IRS 1 N is remarkable amongst observed maser structures for its smoothness and symmetry. At 12 GHz (see Fig. 1), the overall linear structure extends over 50 beams. The central feature extends over 20 beams and shows a clear ``S'' shaped symmetry in its diagram. Such a structure must come from some coherent geometry which controls the gas flow; both discs and outflows have been proposed to provide this controlling geometry. In this paper, we have applied a general formalism to simulate these two cases and we find that discs can readily fit the data. In contrast, outflows predict gaps in the emissivity at the centre of the diagram and rapid falloff of the maser opacity with projected distance, all contrary to our observations. Only extremely contrived outflow models are able to fit the data and therefore we rule these out.
The fundamental reason why discs are superior to outflows to fit the observations is the different functional forms of the opacity versus projected distance which comes from their different symmetries. At the front of a disc, the column density per velocity coherence length scales as 1 + kx^{2} (see Sect. 4.4 and Eq. (8)). For a sufficiently rapid, differentially rotating disc, k can be negative, in which case the centre becomes a local maximum in opacity, as observed. The absolute value of k, however, remains small enough compared to the range of so that the opacity can, as observed, change only by a few tens of percent for changes of a factor of 10 in . In contrast, for an outflow the opacity per unit velocity width scales as down to some minimum where n is at least 1 (see Sect. 6.2). This means that the change in opacity with projected distance is very fast and can only be reconciled with the data if one proposes a very contrived behavior for the ratio of maser opacity to total column density. An additional fundamental property of outflows is that there must be discontinuity in its properties at the minimum ; a discontinuity which we do not see.
We believe that our NGC 7538 data provide perhaps the best purely dynamical evidence to date for a rotating disc in a highmass star forming region, from any spectral line observation, thermal or maser. Amongst maser sources, the detection of a single smooth continuous symmetrical structure spanning 20 beams is extremely unusual as other maser sources in starforming regions generally consist of a series of disconnected spots. Furthermore, this number of beams across the source is much larger than for thermal emission observations: these usually span only a few beams, only providing average velocity gradients. In contrast, in our diagram data we find unambiguous evidence for a differentially rotating disc seen edgeon, supported by a number of independent facts. These include the symmetry and decrease in intensity from a central peak as well as the curvature of the diagram spine. Our modelling fits a 100 point data matrix with only four free parameters to an overall accuracy of better than 34%. This is the highest ratio of observables to fitted parameters that we are aware of for any dynamical study claiming discs around a highmass protostar.
Based on our present maser data, our modelling gives estimated central masses for Keplerian discs in the range of 4 to infinity (see Sect. 5.4). As noted in Sect. 5.5, other observations of outflows in radio continuum and millimeter lines and of the total bolometric luminosity suggest a mass >8 with the best estimate being 30 . It seems clear that our differentially rotating gas disc surrounds a highmass star but it is not yet clear what the mass of the central object is. Confirmation of a mass >20 would be highly significant because it would be beyond the mass threshold of the highest mass stars so far produced via accretion discs in hydrodynamic simulations (see e.g. Cesaroni et al. 2007; Zinnecker & Yorke 2007).
It is difficult to accurately constrain the central mass with our present maser data because we only detect emission from a small part of the disc (see Sect. 5.4). However, prospects are good for making mass estimates with future deeper observations because for different central masses we can predict the full diagram and compare it with these new observations. Figure 9 shows the full diagram distribution for Keplerian rotation around a central mass of 31 ; interestingly, this already suggests some agreement with the data that were not part of the original fitting. In the figure, the red contours correspond to the extent of the fitted data, the green dots indicate the position of the outliers at displacements around 40 and 80 mas (see Fig. 1). Even though these outliers were not considered in making the fit, the spine of the model passes close to at least one of them. This model also predicts local maxima in the maser emission at the edge of the disc (socalled tangential maser emission) at displacements 220 mas, at levels just below our present observational limits. Deeper observations can search for this emission. Such a detection would be a very powerful confirmation of the validity of our model and would more tightly constrain the central mass estimate. Another approach to constraining the mass is to look for the proper motions of distinct components relative to the central peak. In fact, in the framework of an edgeon rotating disc, if the outliers are indeed independent inhomogeneities orbiting the central star in a Keplerian fashion, then they should show a transverse motion either toward or away from the central LOS, marked by the brightest emission. In particular, if the outlier closest to the main feature is at a radius of some 600 AU from a 30 central protostar, we expect for Keplerian dynamics a proper motion of some 0.20.3 mas yr^{1}, which should be detectable within the multiple epoch data that are presently being analysed (Jerkstrand et al., in prep.).
Finally, further investigations could cast light on the real shape of the profile. As stated in Sect. 3, the shape of the function contains details of inversion, temperature and density gradients, which we do not aim to discuss in this paper. Hints for tackling the exact modelling of the maser pumping come from the fact that 6.7 and 12.2 GHz maser maps are strikingly similar (see Fig. 1). This is in fact slightly surprising as the two masers are produced by two different molecules that in turn are expected to be excited differently and hence trace different regions. High resolution (millimetre interferometric) observations of thermally emitting transitions of methanol sharing levels with maser transitions will allow us to directly estimate the level of amplification of both masers and possibly the level of inversion. These projects are ongoing and will be reported on in a later publication.
7.2 Summary and conclusion
In this paper, we have carried out a detailed analysis of the main methanol maser feature in NGC 7538 IRS 1 N. We have extended the formalism and edgeon disc model first presented by Pestalozzi et al. (2004). In addition, we have considered alternative outflow models. Our main conclusions are:
 The main feature of the 12.2 GHz (and 6.7 GHz) methanol maser emission toward NGC 7538 IRS 1 N is well modelled by a sufficiently fast, differentially rotating disc seen edgeon.
 In contrast, outflow models are unable to explain the central peak seen in the diagram and the gradual falloff in opacity versus projected distance. They also do not naturally explain the
curvature seen in the p,vspine data. Only by resorting to extremely finetuned models can outflow models duplicate the observed features in the data.
 For Keplerian rotating discs, we find good fits for central masses ranging from 4
to infinity. Models with central masses of 30 ,
as required by the bolometric luminosity, provide excellent fits.
 Because of the high ratio of observables to unknowns in our modelling, we argue that the NGC 7538 IRS 1 N maser emission provides one of the best pieces of evidence to date for a compact accretion disc around a highmass protostar. Our results further bolster the case for massive stars forming via the same accretion disc mechanism at work in lowmass stars.
 For different central masses, our disc model makes clear predictions for the opacity in the full diagram and proper motions. These predictions can be tested against future observations to more accurately constrain the central mass.
An important final point is that the application of our formalism is not restricted to masers but can also be applied to optically thin thermal spectral line emission. For such emission, the brightness is directly proportional to the optical depth ( ), while for masers . This paper has investigated model distributions from edgeon rotating discs, which were compared to the natural logarithm of the brightness of our observed maser emission. For optically thin thermal edgeon discs^{}, the same disc model can instead be compared directly to the intensity. One subtle point should be noticed however: the average over a varying exponential function in p,vspace is not the same as taking the average of opacity and then calculating the exponential of this average. This issue has to be taken into account if we compare models and data in maser systems that we do not resolve in space or velocity, or when we consider calculating integrated spectra of thermal and maser sources. In the latter case, for the same disc geometry, kinematics and , the total integrated spectra arising from maser emission would not in general equal the exponential of the spectrum from optically thin thermal emission. With the above practical caveats, the results from this paper can be applied to (nearly) edgeon thermal discs. These results include the equations for opacity in the diagram (Eq. (3)) and opacity spine (Eq. (6)), the inner gradient of this spine (Eq. (7)) and its intensity variation (Eq. (8)) as well as its curvature (Eq. (10)). The condition for having a central peak in opacity (Eq. (9)) also applies.
Acknowledgements
This work is the result of a fruitful collaboration involving people from several institutes in Europe, USA and Chile aiming at studying the highmass star forming region NGC 7538 (see also http://physics.gu.se/~micpes). The fitted 12 GHz data were reduced in collaboration with A. Jerkstrand. Michele Pestalozzi thanks Gothenburg University for support through a postdoctoral position. Moshe Elitzur acknowledges NSF support through grants AST0507421 and AST0807417. John Conway acknowledges partial support by a Swedish Research Council grant.
References
 Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959 [NASA ADS] [CrossRef] (In the text)
 Akabane, K., & Kuno, N. 2005, A&A, 431, 183 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Batrla, W., Matthews, H. E., Menten, K. M., & Walmsley, C. M. 1987, Nature, 326, 49 [NASA ADS] [CrossRef] (In the text)
 Blitz, L., Fich, M., & Stark, A. A. 1982, ApJS, 49, 183 [NASA ADS] [CrossRef]
 Bonnell, I. A., Stephen, G. V., & Matthew, R. B. 2004, MNRAS, 349, 735 [NASA ADS] [CrossRef] (In the text)
 Campbell, B. 1984, ApJ, 282, L27 [NASA ADS] [CrossRef]
 Cesaroni, R., Felli, M., Testi, L., Walmsley, C. M., & Olmi, L. 1997, A&A, 325, 725 [NASA ADS]
 Cesaroni, R., Galli, D., Lodato, G., Walmsley, C. M., & Zhang, Q. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 197 (In the text)
 Cesaroni, R. nad Neri, R., Olmi, L., Testi, L., Walmsley, C. M., & Hofner, P. 2005, A&A
 De Buizer, J. M. 2003, MNRAS, 341, 277 [NASA ADS] [CrossRef] (In the text)
 De Buizer, J. M., & Minier, V. 2005, ApJ, 628, L151 [NASA ADS] [CrossRef] (In the text)
 Dodson, R., Ojha, R., & Ellingsen, S. P. 2004, MNRAS, 351, 779 [NASA ADS] [CrossRef] (In the text)
 Edris, K. A., Fuller, G. A., Cohen, R. J., & Etoka, S. 2005, A&A, 434, 213 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Elitzur, M. 1992, Astronomical Masers (Kluwer Academic Publishers) (In the text)
 Ellingsen, S. P. 1996, Ph.D. Thesis, University of Tasmania, Hobart
 Gaume, R. A., Goss, W. M., Dickel, H. R., Wilson, T. L., & Johnston, K. J. 1995, ApJ, 438, 776 [NASA ADS] [CrossRef]
 Green, J. A., Caswell, J. L., Fuller, G. A., et al. 2009, MNRAS, 392, 783 [NASA ADS] [CrossRef] (In the text)
 Greenhill, L. J., Gwinn, C. R., Schwartz, C., Moran, J. M., & Diamond, P. J. 1998, Nature, 396, 650 [NASA ADS] [CrossRef]
 Greenhill, L. J., Reid, M. J., Chandler, C. J., Diamond, P. J., & Elitzur, M. 2004, in Star Formation at High Angular Resolution, ed. M. Burton, R. Jayawardhana, & T. Bourke, IAU Symp., 221, 155
 Keto, E. 2003, ApJ, 599, 1196 [NASA ADS] [CrossRef] (In the text)
 Kraus, S., Balega, Y., Elitzur, M., et al. 2006, A&A, 455, 521 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Menten, K. M. 1991, ApJ, 380, L75 [NASA ADS] [CrossRef] (In the text)
 Minier, V., Booth, R. S., & Conway, J. E. 1998, A&A, 336, L5 [NASA ADS] (In the text)
 Moscadelli, L., Menten, K. M., Walmsley, C. M., & Reid, M. J. 2002, ApJ, 564, 813 [NASA ADS] [CrossRef] (In the text)
 Moscadelli, L., Reid, M. J., Menten, K. M., et al. 2009, ApJ, 406
 Norris, R. P., Byleveld, S. E., Diamond, P. J., & co authors 1998, ApJ, 508, 275 [NASA ADS] [CrossRef] (In the text)
 Patel, N. A., Curiel, S., Sridharan, T. K., et al. 2005, Nature, 437, 109 [NASA ADS] [CrossRef] (In the text)
 Pestalozzi, M. R., Elitzur, M., Conway, J. E., & Booth, R. S. 2004, ApJ, 603, L113, erratum: 2004, ApJ, 606, L173 [NASA ADS] [CrossRef] (In the text)
 Pestalozzi, M. R., Minier, V., & Booth, R. S. 2005, A&A, 432, 737 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Pestalozzi, M. R., Minier, V., Motte, F., & Conway, J. E. 2006, A&A, 448, L57 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes, Third Edition (Cambridge University Press) (In the text)
 Reid, M., Menten, K. M., Greenhill, L. J., & Chandler, C. J. 2007, ApJ, 664, 950 [NASA ADS] [CrossRef] (In the text)
 Schreyer, K., Semenov, D., Henning, T., & Forbrich, J. 2006, ApJ, 637, L129 [NASA ADS] [CrossRef] (In the text)
 Shepherd, D. S., Claussen, M. J., & Kurtz, S. E. 2001, Science, 292, 1513 [NASA ADS] [CrossRef] (In the text)
 Toomre, A. 1963, ApJ, 138, 385 [NASA ADS] [CrossRef] (In the text)
 Walsh, A. J., Burton, M. G., Hylard, A. R., & Robinson, G. 1998, MNRAS, 301, 640 [NASA ADS] [CrossRef]
 Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481 [NASA ADS] [CrossRef] (In the text)
Footnotes
 ... (EVN^{}
 The EVN is a joint facility of European, Chinese, South African and other radio astronomy institutes.
 ... (VLBA^{}
 The VLBA is operated from the National Radio Astronomy Observatory's Array Operation Centre in Socorro, NM.
 ... star^{}
 Note that in Pestalozzi et al. (2004) the value of 30 was adopted for consistency with observations at other wavelengths and was not obtained from the model.
 ... expressed as^{}
 For the model fit to the radio interferometry data we adopt the radio convention to convert frequencies into LOS velocities: , where f is frequency.
 ... discs^{}
 For nonedgeon optically thin thermal discs, the formalism of this paper can be applied to 2D positionvelocity data created from the full data cube by integrating emission along the minor axis of the projected disc.
All Tables
Table 1: Bestfit models of edgeon differentially rotating discs. The parameters = 0.432 km s^{1} and = 0.064 km s^{1} mas^{1} are determined directly from the data and are common to all models. The distance to the source D = 2.7 kpc.
All Figures
Figure 1: Map of the velocity integrated optical depth of the maser emission in NGC 7538, obtained from Eq. (1). From top to bottom: 6.7 GHz data, 12.2 GHz data, 12.2 GHz data convolved at the resolution of the 6.7 GHz data. The peak value of is 18 and 16 for the 6.7 and 12.2 GHz respectively, in the outliers it reaches 11 and 14. The beam size is indicated in the lower right corner of each panel. Contours are at 70, 80, 90, 95, 99% of the peak in all panels. Note that the lowest contour level is set by the dynamic range limit of our observations. 

Open with DEXTER  
In the text 
Figure 2: LOS velocity versus spatial displacement (diagram) of of the main maser emission feature in NGC 7538. The top panel shows 6.7 GHz data, the bottom panel the 12.2 GHz data. The peak value of is 18 and 16 for the 6.7 and 12.2 GHz respectively, in the outlier features at 40 mas and 80 mas and it reaches 11 and 14. The beam size is indicated in the lower right corner of each panel. Contours are at 70, 80, 90, 95, 99% of the peak in both panels. As for Fig. 1, the lowest contour level is set by the dynamic range limit of our observations. 

Open with DEXTER  
In the text 
Figure 3: Geometry for the disc model, with the rotation velocity. The observer is in the plane. 

Open with DEXTER  
In the text 
Figure 4: diagram of maser optical depth (Eq. (2)) for edgeon discs with constant amplification ( ) between inner radius r_{1} and outer radius 3r_{1}. Contour values are normalised to the value at the diagram center and are spaced by 0.1, with the first contour at 0.6. Black contours indicate values greater than the central value. Top: a nonrotating disc. The emission maxima are along the tangents to the inner edges , the longest paths through the disc. Mid and Bottom: same configuration as top panel only the disc rotates as a solid body with the angular velocity , as indicated. Notice that all the maps are qualitatively identical except for the tilt, i.e. norotation and solid body rotation are two aspects of the same case. 

Open with DEXTER  
In the text 
Figure 5: Contours of constant LOS velocity () and the velocity fields that produce them. The contours are spaced by 0.1 of the highest value and are also colourcoded to indicate the Doppler velocity. The velocity vector field is indicated by arrows. Left: Keplerian rotating disc with masing material confined within inner and outer radii with ratio 1:3. Right: a bipolar flow that produces the same diagram as the disc geometry (see Eq. (16) below), which has outflowing gas on the right side of the central LOS and inflowing gas on the left side, as indicated by the arrows. The radial velocity was multiplied by a factor of 2.5 for clarity. The masing material has been confined within a cone of opening angle 20. 

Open with DEXTER  
In the text 
Figure 6: Same as Fig. 4, only the disc is in Keplerian rotation, with angular velocity on the inner radius. Notice that a maximum at the centre is obtained when the rotation is fast enough ( bottom panel). 

Open with DEXTER  
In the text 
Figure 7: Radial variation of the linecenter absorption coefficient (Eq. (11)) for values of p and of three of the bestfit solutions in Table 1 with (solid), 0 (dashed) and 0.3 (dasheddotted); the legend lists also the corresponding width in each case (Eq. (12)). All profiles are normalised to unit area and have their maximum at = 1 (r = r_{1}), though r_{1} is different in each case (see Table 1). 

Open with DEXTER  
In the text 
Figure 8: Data ( left) and modeling ( middle) of the bulk of the 12.2 GHz maser emission (within 20 mas of the feature centre). Contours are at 70, 80, 90, 95 and 99% of the peak optical depth (). Red dots show the amplification spine (see Sect. 4.4). Right: map of the model relative error at every position in the diagram, corresponding to the individual terms in the sum in Eq. (14). Note the color bar on the right for the error magnitude in percent. The error averaged over the entire map is 3.07%. The model corresponds to an edgeon Keplerian disc around a 31 star (see Table 1 for details). 

Open with DEXTER  
In the text 
Figure 9: Full diagram of the best fitting Keplerian model with a 31 central protostar. Contours are plotted at 1, 3, 10, 30, 50, 60, 70, 80, 90, 95 and 99% of the peak optical depth = 16. The red contours outline the region visible in our data (>70% of ); this is the region shown in Fig. 8. Notice the small fraction of the full extent of the masing disc covered by the available data. The green dots indicate the positions of the outliers. 

Open with DEXTER  
In the text 
Figure 10: Geometry of a physical outflow model whose axis is inclined at an angle i to the LOS. Distance along the outflow is measured by . The radius marks an injection point for material into the outflow (see text). Beyond this point the outflow has a position dependent bulk flow velocity (dashed arrows) and mass flux conservation is assumed. In general the outflow has a crosssectional radius , which for the conical flow illustrated is proportional to . 

Open with DEXTER  
In the text 
Figure 11: diagram of maser optical depth for the biconical outflow sketched in Fig. 10. Each panel shows a different inclination i of the outflow axis to the LOS. Contours are drawn at 70, 80, 90, 95 and 99% of the peak optical depth. The outflow has a constantly accelerating velocity, starting from at the injection radius . Maser opacity beyond is proportional to the total gas column density. The opacity and dynamics inside the injection radius are illdefined: in this example the opacity there is taken as zero; quiescent maser emission from material within would add in each panel a horizontal strip about the displacement axis, similar to the top panel of Fig. 4. 

Open with DEXTER  
In the text 
Copyright ESO 2009
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.