1. Introduction
Pre-protostellar cores (PPCs) are starless cores that sample the earliest stages of low-mass star formation. Due to their lack of a meaningful internal heat source, they hold the promise of providing a unique picture of the physics of the initial conditions leading to star formation. These regions have been significantly studied in molecular lines, such as CO, CS, N2H+, NH3, HCO+, and their isotopes [1] -[6] . These studies show relatively quiescent cores that are dynamically young and chemically old as evidenced by depletion of gas molecules onto grain surfaces at the low temperatures in these cold cores.
Studies of the structures of these regions have focused on far-infrared and sub-millimeter continuum observations of emission from dust grains. In the case of L1498, a meaningful history of such studies has taken place. Langer and Willacy [7] used semi-analytic modeling to show that a spherical source with a broken power law fit the source well-namely n(r) ~ r−m, where m is the power law index. They found that the data could be fit with a relatively flat density inner core, a steeper middle region having m ~ 1.8, and a sharp outer edge with m ~ 4. Tafalla et al. [8] [9] and Shirley et al. [10] used radiative transfer modeling to study the same source. Tafalla et al. found that a density profile of the form n(r) ~ 1/[1 + (r/r0)m] fits the data, with m ~ 3.5. Shirley et al. utilized a Bonner-Ebert sphere, but also found that the external radiation field needed to be somewhat decreased (s = 0.5 - 1, where s is the ratio of the local radiation field to the standard [11] radiation field strength) and attenuated (τISRF ~ 1).
Unfortunately, the analyses above either do not take in account the non-spherical nature of the source [8] -[10] , or have not completed detailed, self-consistent radiative transfer modeling [7] . In particular, it is interesting to note that the spherical models do not require the sharper edge inferred by the semi-analytical modeling. This is perhaps not a surprise as models based upon a spherical assumption will necessarily smooth-out gradients in the observed radiation field. Previous work in multi-dimensional radiative transfer modeling for a different source, L1544 [12] , showed that some differences can arise even in gross properties such as mass, when considering the non-spherical nature of the source.
In this paper we report a study comparing existing data with detailed, three-dimensional, continuum radiative transfer models to constrain the structure and surroundings of L1498. In Section 2, we describe L1498 and previously inferred properties. In Section 3, we briefly describe the model we apply. We compare this model to the observations in Section 4. Finally, we conclude in Section 5.
2. L1498
L1498 is a nearby molecular core, located in the Taurus molecular cloud [13] at a distance of D = 140 pc. Ward-Thompson et al. [14] classified it as a PPC due to the lack of an IRAS 100 μm source in combination with its detection in the submillimeter. The source has been well observed ([8] [10] [14] [15] ), with fluxes reproduced in Table1
The fluxes (SED) described in Table 1 can be used to help constrain the source structure. However, the constraint is not unique (see [16] -[18] ). This is due to the fact that the SED primarily samples the amount (mass) of dust at each temperature, modified by the amount of absorption by intervening dust. As a result, there are multiple possible internal structures that could yield the same integrated SED.
Table 1. Observed fluxes toward L1498.
*fluxes used in fitting the SED.
On the other hand, the spatial intensity distribution in the form of spatially resolved maps provides significantly more data with which to constrain models. L1498 was observed using SCUBA by Shirley et al. [10] . They used a chop of 120" perpendicular to the major axis of the core. While the 450 μm maps are difficult to fit due to a low S/N at the higher contour levels, the 850 μm maps are much cleaner. The reduced map contours were fit with ellipsoidal profiles. The resulting sizes of the contours along each of the axes are given in Table2
3. Model
We utilize a detailed, self-consistent, three-dimensional, dust continuum radiative transfer model. The model uses a monte-carlo method combined with approximate lambda iteration to ensure true convergence even at high optical depths. The model has been tested against existing 1D [19] and 2D [20] codes with good success.
Based upon the input parameters discussed below, we solve for the dust temperature and radiation field at each point in the cloud. The emergent radiation is then convolved with the appropriate telescope beam [21] .
For L1498, the chopping process might be important. Even a 120" chop may still include self-subtraction of the source from itself. As a result, we directly simulate the chop process by simulating observations at the source position, as well as the chop positions, and then subtracting the simulated chop observations from the “on-source” data. While this is a somewhat tedious process, the direct simulation of the observing process gives the greatest fidelity between the observations and simulations, especially as regards the relative fall-off of intensity near the source edge. This is an important point that needs to be considered when using spatial map data of extended sources.
3.1. Input Data
We adopt a piecewise triaxial ellipsoid density structure similar to utilized in [12] . In particular, in each section a density structure of the form
is used. The sections are joined by choosing n0 in each section so that the densities are the same at the matching point. In this case, the x axis is defined along the long axis of the cloud, the y-axis along the short axis, with the z-axis perpendicular to both. As noted earlier, previous work has suggested that a broken power law can fit the data, with low values of m in the interior, and potentially larger values of m in the exterior.
Following our previous work on L1544, we reduce our density structure to a prolate spheroid, with symmetry in the yand z-directions. Mathematically, this is done by taking b = c, and defining the axis ration q = a/b. Based upon the observations, we would expect 1 < q < 3.
The interstellar radiation field is the dominant heat source. We adopt the ISRF of [11] (hereafter MMP). The strength of the ISRF is parameterized by sISRF = (adopted ISRF)/(MMP ISRF). We expect sISRF to be in the range 0.1 - 5 or so. We also consider the effect of a luminous internal heat source by considering 10−6 < L*/Lsun < 1. The ISRF is allowed to be attenuated by the surrounding natal cloud. This attenuation is parameterized by τISRF and is expected to be in the range 0 < τISRF < 3 or so.
The asymmetry in the observed intensity map may be due to an asymmetric external radiation field, as opposed to being due to the density structure itself. In order to allow for this possibility, we adopt a plane-parallel external radiation field. This plane-parallel radiation field is taken to have the same spectral distribution as the
Table 2 . Contour size at 850 mm along each azis for L1498.
Uncertainties given in square brackets, so that 36.7 [2.5] = 36.7 ± 2.5.
ISRF, but to come from only a single direction. Its strength and attenuation are specified by sPPRF and τPPRF in a manner similar to the ISRF above. The angle from which the PPRF is incident is also specified by the parameter θPPRF.
We adopt the dust opacities from [22] , column 5—commonly called the OH5 dust properties. They have been successful in modeling both high mass and low-mass forming regions, and should be applicable to the L1498 PPC.
Finally, the total grain density is parameterized by the optical depth in the visible part of the spectrum, τv. Based upon previous work, and the opaque nature of the cloud, this is expected to be in the range 3 < τv < 50.
3.2. Grid of Models
The models above are run for various sets of values for the input parameters, and compared with observations via a χ2 (see discussion in Section 4). The overall best fit is determined through a global minimization using Powell’s method [23] . We find that different initial conditions yield the same best fit for multiple different initial conditions.
Once a best fit is identified, grids of models are run in the vicinity of the best fit in order to best understand the parameter space topology. More importantly, these grids allow us to place uncertainties/constraints on the model parameters. Between the minimization scheme and the local grids, over 5000 models were run as part of this approach.
3.3. Comparison with Observations
In order to compare the simulated and actual observations, a reduced χ2 statistic was utilized. For the 850 μm isophote contours,
where i corresponds to the isophote contour number. In this expression, r is the radial location of the isophote and σ is the uncertainty in that location, from Table2
We can likewise determine the quality of the SED fit in an analogous manner
where i is the number of the data point, S is the flux (either observed or modeled), and σ is the uncertainty in the flux, from Table1
The overall quality of fit is then taken to be the average of these two reduced chi-squares (assuming that they are independent), as
This choice for overall fit has the advantage that neither one of the individual chi-squares generally dominate the other. Furthermore, the individual chi-squares tend to be sensitive to different parts of the structure—the SED chi-square is sensitive to total amount of dust and the strength of the external radiation field, while the 850 μm chi-square is more sensitive to the details of the density profile. As a result, a best fit for this total chisquared gives the best simultaneous fit to all data, and will be the chi-squared value that we discuss below.
4. Results
4.1. Best Fit
The best fit results for the 850 μm isophotes and the flux are shown in Figure 1 and Figure 2. Figure 1 shows the simulated contours (solid curved lines) and the best fit ellipses on each axis as points with 3σ error bars. All contours are given as a fraction of the maximum intensity in the image. Note that the data are well fit by the simulations along each axis. Also note the East-West and North-South asymmetries in the data; in particular, the
Figure 1. Comparison of model (curved lines) and observed (points with 3-sigma error bars) isophotes at 850 microns toward L1498.
Figure 2. Comparison of simulated (solid line) and observed (symbols with errorbars) fluxes toward L1498. While the simulated data are calculated on a finer wavelength grid than shown, only the points corresponding to the observational data are used in the plot so that the most direct comparison can be made.
observed contours are wider spread to the West and the South than they are to the East and the North. These asymmetries are accounted for in the simulations. In Figure 2, the solid line represents the best fit for the flux. The symbols with error bars are the observations from Table1 The best fit simulation clearly fits the observed flux data well. The properties of the best fit model with uncertainties are reported in Table3
4.2. Central Luminosity
Chi-squared contours in the external radiation field strength (sISRF)—luminosity (L*) plane near the best fit are
Table 3 . Properties for best fit model for L1498.
shown in Figure 3. Here the strength of the external radiation field is specified by sISRF = (strength of external radiation field)/(standard ISRF), so that a value of 1 means an external radiation field of the same intensity as the standard [11] . Likewise, the central luminosity is specified in solar luminosities. It has been noted previously that the temperature of the central source is unimportant in modeling the infrared and submillimeter emission from dust clouds [18] [24] , and this is true in our modeling also.
As can be seen in the figure, the central luminosity of the source is less than 0.001 solar luminosities. This is consistent with earlier findings [14] where a lack of a mid-infrared excess suggested that a central source had not yet formed.
4.3. External Radiation Fields
As noted above, the central luminosity of this source is quite low—less than 0.001 Lsun. As a result, most if not all of the heating of the dust is contributed by an interstellar radiation field. Following [10] , we allow for an ISRF of variable strength and attenuation. The strength of the ISRF (relative to the standard) is given by sISRF as noted above. The attenuation of the ISRF is given by an optical depth τISRF. These are adjusted as part of the fitting procedure described in Section 3.2 above. As before, a grid of models around the best fit is constructed so that the topology of the region and quality of the fit can be understood.
The results of the grid in the τISRF – sISRF plane are shown in Figure 4. As can be seen, there is some degeneracy between the two values. This is to be expected, since increasing the attenuation will commensurately require an increase in the intrinsic strength of the radiation field. However, since the attenuation is not the same at all wavelengths, the degeneracy is not complete.
For a power-law ellipsoidal density distribution discussed previously, a uniform external radiation field would lead to symmetric isophotes. This is not what is seen in the observational data. As a result, the density distribution is either asymmetric or the external radiation field contains an asymmetric component. Since young clouds are expected to be quiescent and relatively undisturbed, and since this source is part of a larger star-forming region, we view the possibility of an asymmetric external radiation field as much more likely. As a result, we adopt an additional external plane-parallel radiation field (PPRF). Similar to the ISRF above, we characterize the PPRF by its strength relative to the ISRF (sPPRF), an attenuation (τPPRF), and also an angle of incidence θPPRF.
In Figure 5 and Figure 6 we present the chi-squared contours for the cuts through the sPPRF, τPPRF, and θPPRF parameter space. In Figure 5 the chi-squared contours through the τPPRF − sPPRF space. As with the ISRF in Figure 4, there is a degeneracy between the attenuation and the strength of the external field. What can be seen is that while the intrinsic strength of the PPRF can vary significantly depending upon attenuation, a PPRF is
Figure 3. Chi-squared contours as a function of the strength of the external radiation field and the central luminosity.
Figure 4. Chi-squared contours as a function of the attenuation and strength of the interstellar radiation field.
Figure 5. Chi-squared contours of quality of fit as a function of the attenuation and strength of the external plane-parallel radiation field.
Figure 6. Chi-squared contours of quality of fit as a function of the angle and strength of the external plane-parallel radiation field.
necessary. The unattenuated strength of the PPRF is about equal to the normal ISRF (1.0 ± 0.5). A best fit attenuation of 1.25 ± 0.75 is found as well.
In Figure 6, the chi-squared contours of the PPRF are shown as the angle of the PPRF and its strength are varied. As can be seen, the PPRF does not impinge along an axis, but rather impinges at θππρφ = (21/16)π ± π/8 radians—this is from the North-East at nearly 45 degrees from the North-South axis.
4.4. Density Distribution
Perhaps the most significant clue as to the environment, origin, and evolution of star forming cores is their density distribution. Most density distributions are parameterized as n(r) = n0 (r0/r)m. In the case of the singular isothermal sphere, m = 2. Shu characterized an infalling sphere as having m tending to 1.5 inside the infall radius. In the case of a self-gravitating cloud at early times and small radii, m tends to zero. In each of these cases, the cloud is assumed to be spherical and either supported by thermal pressure, or collapsing as self-gravity overcomes thermal pressure.
Another possibility is a cloud bounded on the outside and constrained by external pressure—a situation also known as pressure confinement. This has been discussed by numerous authors [25] -[29] . Pressure confinement mayhelp explain observational data [30] , including obervations of 30 Doradus [31] , the Pipe Nebula [32] [33] , and high-mass infrared dark clumps [34] (Pitann et al. 2013).
As a result, our fitting procedure allowed for a broken power law, with a density power law index m1 for r < rcut, and m2 for r > rcut, where rcut is the cutoff between the two power laws. Our best fit finds m1 = 0 ± 0.25, m2 < −10, and rcut = 0.073 ± 0.005pc. A plot of the chi-squared contours in the m2 – rcut plane are shown in Figure 7.
This density profile is extremely steep. Previous authors have found that a Bonner-Ebert sphere fit the profile [10] , that a modified power law with external power law m = 3.5 fit the profile [9] , and a broken power law with an exterior having m > 4 [7] . While the density profile found in this model is steeper than that found in previous work, it is not necessarily surprising. First, as seen above, previous work has suggested a strong outer density gradient—stronger than a canonical m = 1.5 or 2 profile. Second, the previous works all assumed spherical symmetry for the density distribution while we have undertaken a true three-dimensional representation of the density distribution. An ellipsoid, by its elongated nature, will tend to have more mass at larger radii than an equivalent sphere. As a result, the ellipsoidal power law will be steeper at large radii than the sphere in order to maintain the same mass at large radii.
The density profile itself is shown in Figure 8. The left hand axis is the number density of dust grains. Taking the grains to be of size a = 0.1 μm, this density corresponds to an H2 density of n(H2) ~ 5 × 104 cm−3. This density is consistent with the densities predicted to be necessary to be near collapse [35] .
Figure 7. Contours of chi-squared for the total fit as a function of the outer density power law and the cutoff radius in the region near the best fit.
Figure 8. Dust density distribution for the best fit model. Notice the flat interior and extremely fast fall-off in the exterior.
The density distribution here is potentially quite important. This fall-off is precipitous. It is as if there is a sharp “wall” between the nearly uniform density interior and the exterior of the cloud. This is exactly what one would expect in a pressure confinement scenario, where the cloud was being confined by the exterior medium. The surrounding gas is expected to have a number density on the order of 1 - 5 × 103 cm−3, and could thus provide that confinement. Indeed, it has been suggested [36] that the outer envelope may be still growing by accretion from the surrounding cloud. This would be consistent with the pressure bound core that has not yet achieved sufficient mass to collapse under its own gravity.
5. Conclusions
We have constructed models for the non-spherical preprotostellar core L1498 utilizing a fully three-dimensional continuum radiative transfer model. After convolving the output with the actual SCUBA beam profile and including the beam chopping directly in the simulated observations, we are able to compare our results with existing observations, including both SEDs and sky maps. We find that:
1) The observational data (SEDs and sky maps) can be well fit by a three-dimensional prolate spheroid density distribution. This fit includes the effects of beam chopping in the observations. [Section 4.1]
2) The central source has a luminosity < 10 - 3 Lsun, consistent with the lack of an mid-infrared excess. [Section 4.2]
3) The main source of heating is an external radiation field. The first is a uniform field of strength sISRF = 0.5 ± 0.25 of the standard ISRF, accompanied by an attenuation τISRF = 1 ± 0.25, consistent with the fact that the source is embedded in a surrounding lower density cloud. The second radiation field is a plane-parallel field of strength sPPRF = 1.0 ± 0.5, and attenuation τPPRF = 1.25 ± 0.75, impinging at an angle of θPPRF = (21/16)π ± π/8. The second field is consistent with the notable asymmetry in the sky map isophote contours. [Section 4.3]
4) The density distribution in the source is consistent with a constant density central region (density power law m1 = 0.0 ± 0.25) out to a cutoff radius rcut = 0.073 ± 0.008 pc. Outside of this, the density drops off precipitously, with a density power law m2 > 10. This density profile is interpreted as being consistent with a pressure-confined cloud which has not yet achieved sufficient mass to collapse under its own gravity. [Section 4.4]
Acknowledgements
We would like to thank Yancy Shirley for providing copies of his data for L1498 and for useful discussions. We would also like to thank Brant Bowers for his help in preparing Figure 7.